Metagenomic study on functional characteristics of soil carbon, nitrogen, and phosphorus cycling under rice-crayfish farming

Sustainable agricultural systems are crucial for balancing productivity and environmental conservation^[1,2]. As a prominent integrated farming model, rice-crayfish (RC) co-cultivation has been widely promoted in China due to its ecological and economic benefits^[3]. Crayfish activities such as bioturbation and excretion introduce additional organic matter into paddy soils, altering the biogeochemical environment^[4,5]. Soil microorganisms are core drivers of element cycling and energy flow, regulating organic matter mineralization and nutrient transformation to maintain ecosystem stability^[6,7].

Emerging evidence suggests that RC systems could enhance soil quality and microbial activity; for instance, the increased soil organic matter characteristic of RC systems provides more active carbon sources, supporting greater microbial diversity^[4,8]. Specific shifts in microbial functional groups have also been reported, such as a lower abundance of methanogens coupled with a higher abundance of methanotrophs compared to rice monoculture (RM), and notable alterations in the community structure of denitrifying microbes harboring the nirK gene, alongside increased abundances of the nosZ and AOA/AOB amoA genes^[9]. Despite these valuable insights, significant knowledge gaps persist, constraining a mechanistic understanding of nutrient cycling in RC systems.

Currently, existing research has largely focused on individual elemental cycles, particularly C or N, in a fragmented manner^[10,11]. A comprehensive, simultaneous investigation into the interconnected microbial-driven mechanisms of C, N, and P cycling is notably lacking. This singular focus prevents a holistic view of how these nutrient cycles are coupled and co-regulated within the integrated system^[12]. Secondly, there is a pronounced temporal and spatial bias in current research^[13,14]. Most studies are confined to the rice-growing season, overlooking the critical non-rice period, which is essential for crayfish cultivation^[15,16]. During this phase, substantial feed input provides rich nitrogen sources, and the distinct, perpetually flooded conditions create a biogeochemical environment vastly different from that of drained fields in RM, likely fostering unique microbial processes^[17]. Furthermore, the trench (RCT) area, which represents a unique habitat that serves as the main arena for crayfish activity, a repository for residual feed, and a zone of persistent deep flooding, has been largely ignored^[18]. The adaptive responses and functional roles of microbial communities, especially those involved in nitrification, denitrification, and phosphorus transformation, in this distinctive niche remain virtually unexplored^[19]. Consequently, the microbial driving mechanisms driving transformation processes, particularly for phosphorus transformations, are still poorly understood, limiting the ability to predict and manage the overall efficiency and environmental impact of RC systems.

This study aimed to: (1) systematically characterize the functional genes and co-occurrence patterns of C-, N-, and P-cycling microorganisms in RC systems; (2) explore potential functional couplings among these cycles; and (3) identify key environmental factors regulating elemental cycling. It was hypothesized that: (1) RC farming would significantly enhance the abundances of genes involved in methane oxidation, N degradation, denitrification, anammox, organic P mineralization, and P transport compared to RM; (2) soil nutrients and their availability would be the primary environmental factors regulating these cycling processes.

The experiment was established in 2014 at the rice-crayfish farming experimental base of Hunan Agricultural University, located in Wangcheng District, Changsha City, China (Supplementary Fig. S1). The site has a subtropical monsoon climate with an annual average temperature of 17.5 °C Hunan Agricultural University, located in Wangcheng District, Changsha City, China would significantly enhance the abundances of genes involved in methane oxida experiment (2013), the soil was classified as a clay loam (29% sand, 21% silt, 50% clay) with the following chemical properties: soil organic C = 22.7 g kg⁻¹; ammonium N = 11 mg kg⁻¹; available potassium (AK) = 35.3 mg kg⁻¹; pH = 6.15.

Experimental design and soil sampling

This study employed a completely randomized design with three replications, focusing on RM and RC systems. Within the RC system, two spatially and functionally distinct zones were defined: the central paddy field (RCP) and the surrounding forage ditch (RCT). Consequently, three primary treatments were established for comparison: RM (rice monoculture), RCP (paddy field within the RC system), and RCT (forage ditch within the RC system).

In the RC systems, each 660 m² (20 m × 33 m) experimental field was surrounded by a circular water ditch. This ditch measured 1.5 m wide and 1.2 m deep, covering approximately 66 m². Its primary function was to provide habitat and refuge for crayfish. Adjacent plots were separated by soil barriers covered with plastic sheeting (40 cm × 40 cm) to prevent the exchange of moisture and nutrients between plots. To prevent crayfish escape, all RC experimental fields were surrounded by a 1.5 m embankment and an additional 0.3 m-high nylon net.

The rice cultivar used was Nongxiang 42. In early June, at the beginning of the rice growing season, the paddy field areas in all RC and RM plots were leveled, compacted, and received a basal fertilizer application. Fertilizer management during the rice growing season was identical for both treatments. All fertilizers were applied exclusively in this season at the following rates: urea (N 46.0%) supplying 100.0 kg N hm⁻², calcium superphosphate (P₂O₅ 12.0%) supplying 30.0 kg P₂O₅ hm⁻², and potassium chloride (K₂O 60.0%) supplying 50.0 kg K₂O hm⁻². Four days after fertilization and flooding, rice seedlings were transplanted at a density of 160,000 hills hm⁻². One week prior to rice transplanting, juvenile crayfish (average weight 7.0 ± 1.0 g) were stocked into the RCT at a density of 1.2 × 10⁵ individuals hm⁻². Besides, the planting field and trench were initially connected to allow water exchange. After the basal fertilizer application and field flooding, the connection was closed until the drainage period before rice harvest. Throughout the growth period, crayfish could freely move across the low barrier between the trench and the paddy field for foraging. Crayfish were harvested one week prior to rice harvest.

Soil samples were collected at a depth of 0−20 cm after rice harvest. Six subsamples were collected per plot to cover heterogeneity, thoroughly mixed to form one composite sample, and stored at −80 °C for subsequent analysis. Methods for determining soil physicochemical characteristics are detailed in Supplementary Text 1.

DNA extraction and sequencing

Total genomic DNA was isolated from 0.5 g of fresh soil using the EZNA® Soil DNA Kit (Omega Bio-Tek, Norcross, GA, USA). DNA concentration and quality were evaluated using a Thermo Fisher Scientific spectrophotometer and 1% agarose gel electrophoresis. Qualified DNA samples were fragmented using a Covaris M220 focused-ultrasonicator. Sequencing libraries were constructed with the NEXTFLEX Rapid DNA-Seq Kit and subjected to paired-end sequencing on the Illumina NovaSeq 6,000 platform (Illumina, San Diego, CA, USA) by Shanghai Majorbio Bio-pharm Technology Co., Ltd (Shanghai, China). All raw sequencing data were deposited in the NCBI Short Read Archive (SRA) for public access.

Raw metagenomic sequences were quality-controlled using Fastp, retaining high-quality paired-end and single-end reads^[20,21]. Sequences were assembled at varying depths using MEGAHIT (https://github.com/voutcn/megahit)^[22]. Using MEGAHIT v1.1.2, the optimally assembled sequences were selected for ORF prediction^[23]. For subsequent translation and analysis, predicted genes were retained only if they met the length threshold of ≥ 100 bp^[23]. After aligning high-quality reads to genes in SOAPaligner v2.21, gene abundance was quantified by normalizing the mapping read count for each gene according to its length^[24]. Based on the NCBI-NR and KEGG databases, Diamond (v0.8.35) was employed to classify and annotate gene sequences for functional assignment, thereby generating their KEGG functional profiles^[25]. Based on prior research, core genes involved in C, N, and P cycling (present in all samples) were identified, and the abundance of each microbial group was defined as the sum of abundances of all genes assigned to that taxonomic unit^[26,27].

Statistical analysis

Statistical analyses were performed using R software (v4.3.2). Redundancy analysis (RDA) was applied to assess the impact of environmental variables on microbial communities, with a variance inflation factor (VIF) threshold of 10 to mitigate multicollinearity. The Kruskal-Wallis test identified significant differences in microbial functional gene abundances among treatments. The Mantel test evaluated correlations between environmental factors and microbial functional genes. Co-occurrence network analysis was conducted on all samples (n = 9) using a stringent Spearman correlation criterion (|r| > 0.6, p < 0.05), with visualization via Gephi software (https://gephi.org).

An overview of the metagenomic sequencing and assembly for nine soil samples across the RM and RC systems is provided in Supplementary Table S1. After quality control, an average of 48.3 million clean reads and 211,000 contigs per sample were obtained, identifying 8,386 KEGG orthologs (KOs) in total. Microbial taxonomic diversity in the RCP area exceeded that of RM (Supplementary Fig. S2a), though functional diversity did not differ significantly (Supplementary Fig. S2b). Moreover, RCP outstripped RCT in both taxonomic (max: 5.37) and functional (max: 7.63) diversity. A significant positive correlation existed between these two diversity types (p < 0.05) (Supplementary Fig. S2c).

The rice-crayfish system significantly altered the microbial functional profiles for C, N, and P cycling (Fig. 1), with 49, 36, and 41 related genes predicted. PCoA confirmed a significant shift (p < 0.05) in these genes compared to monoculture, indicating that co-culture restructured the microbial community and its C, N, and P metabolism.

Figure 1. 

Microbial functional characteristics of (a) soil C, (b) N, and (c) P cycles identified by principal coordinates analysis (PCoA) based on Bray-Curtis distances. RM: rice monoculture; RCP: paddy field of rice-crayfish farming system; RCT: the peripheral trench of rice-crayfish farming system.

The effect of RC farming on soil C cycling

In this study, the functional microbes are characterized as associated with soil carbon cycling under the rice-crayfish system. The related functional genes are summarized in Supplementary Table S2 and categorized into carbon fixation, carbon degradation, and methane metabolism, respectively. Among these, the rTCA cycle accounted for 39.05%–41.23% of the relative functional gene abundance. Carbon degradation was classified into labile and recalcitrant carbon fractions, encompassing starch, hemicellulose, pectin, cellulose, chitin, and lignin (p < 0.05; Supplementary Fig. S3e). Genes associated with starch catabolism exhibited the highest relative abundance, ranging from 36.36% to 41.42%.

The carbon cycling process in the RC system (Fig. 2) was evaluated. Wilcoxon rank-sum tests revealed significant alterations in relevant gene abundances. In methane metabolism, the abundances of methane-oxidizing genes pmoA, pmoB, and pmoC were significantly enriched in the RCP and RCT areas, but the abundances of methanogenic genes were unaffected (p < 0.05). Conversely, abundances of korA and korB (rTCA cycle) genes were increased, while rbcL and cbbR (Calvin cycle) genes were decreased in the rice-crayfish farming system. Regarding carbon degradation, chi (chitin), malQ, and malZ (starch) genes were significantly increased (p < 0.05). Notably, the RCT area had a lower abundance of the chi gene than the RCP area (p < 0.05).

Figure 2. 

Diagram depicting the different carbon cycling processes based on metagenomic sequencing. Red indicates that gene abundances were significantly and consistently increased in the former treatment; blue indicates that gene abundances were significantly and consistently decreased in the former treatment. RM: rice monoculture; RCP: paddy field of rice-crayfish farming system; RCT: the peripheral trench of rice-crayfish farming system.

The relative abundance of the phylum-level microbial community in the rice-crayfish system was shown in Supplementary Fig. S4. Taxonomic annotation of C-cycle genes revealed Proteobacteria as the dominant phylum across all treatments, with a mean relative abundance of 40.07%, followed by Acidobacteria (17.37%), and Chloroflexi (11.84%) (Supplementary Fig. S4a). Significant shifts in the abundances of major bacterial phyla were observed in the RCP area compared to the monoculture. These shifts were characterized by an enrichment of Verrucomicrobia and Actinobacteria alongside a depletion of Proteobacteria and Acidobacteria (p < 0.05). Furthermore, the RCT area exhibited a significant increase in Bacteroidota and Actinobacteria, alongside a decrease in Proteobacteria and Chloroflexi (p < 0.05) (Supplementary Fig. S4b).

The effect of RC farming on soil N cycling

Functional gene analysis revealed that nitrogen degradation (55.88%−59.01%) and denitrification (25.35%−27.03%) were the two most dominant processes in the soil N cycle of the rice-crayfish system (Supplementary Fig. S5,Supplementary Table S3).

Differential nitrogen cycling processes in the RC system were shown in Fig. 3. Wilcoxon tests confirmed significant alterations in N-cycling genes. Specifically, compared with rice monoculture, in the RCP area, increased abundances were recorded for glnA (5.15%), gltB (6.57%), gltD (13.62%), and arcC (11.37%) (p < 0.05). Significant elevations were also observed for denitrification genes (narG, narH, narI, nirK, nirS, norB, norC) and nitrification genes (amoA, amoB, amoC) (p < 0.05). In the RCT area, increases were measured for gltB (6.75%) and gltD (5.39%), along with narG, nirK, and nirS (11.59%, 5.32%). The napB gene was significantly increased by 6.77%, whereas the amoA gene was significantly reduced by 7.27% (p < 0.05).

Figure 3. 

Diagram depicting the different nitrogen cycling processes based on metagenomic sequencing. Red indicates that the gene abundances were significantly and consistently increased in the former treatment, blue indicates that the gene abundances were significantly and consistently decreased in the former treatment. RM: rice monoculture; RCP: paddy field of rice-crayfish farming system; RCT: the peripheral trench of rice-crayfish farming system.

Taxonomic annotation of N-cycling genes revealed that Proteobacteria, Acidobacteria, Actinobacteria, Chloroflexi, and Bacteroidota accounted for 81.66%−88.02% of the total community (Supplementary Fig. S6a). Compared to RM, RC systems showed increased mean relative abundances of Acidobacteria and Bacteroidota, with reduced Proteobacteria (p < 0.05) (Supplementary Fig. S6b).

The effect of RC farming on soil P cycling

Functional microbes associated with soil P cycling were characterized (Supplementary Table S4). As shown in Supplementary Fig. S7, genes pertaining to the P transport system were most abundant (54.82%–56.72%). Those contributing to inorganic P solubilization were the next highest, at 21.02% of the P (p < 0.05). The rice-crayfish system elevated the abundance of inorganic P solubilization genes by 3.34%–6.22% while depressing that of P transport system genes by 2.70%–3.35% (p < 0.05), respectively.

Phosphorus cycling processes in the RC farming system were shown in Fig. 4. Significant P-cycling gene alterations were induced by rice-crayfish farming. In organic P mineralization, abundances of appA, phnA, phnM, and phnX were increased in RCP and RCT, while phoN and phnG were decreased (p < 0.05); RCT showed stronger effects than RCP. For inorganic P solubilization, the abundance of ppa was reduced with RCP > RCT. In P regulation, the abundance of phoB was decreased (p < 0.05). In P transport, abundances of phnC and ugpC were increased, while pstB, pstS, and phnD were decreased; RCP showed higher pstB, pstS, and phnD but lower phnC and ugpC than RCT (p < 0.05).

Figure 4. 

Difference of soil P cycle gene in different rice farming systems. RM: rice monoculture; RCP: paddy field of rice-crayfish farming system; RCT: the peripheral trench of rice-crayfish farming system.

The P-cycling genes were predominantly distributed among Proteobacteria (24.09%−30.14%), Acidobacteria (18.78%−22.93%), Chloroflexi (15.29%−21.96%), Candidatus_Rokubacteria (8.56%−9.84%), and Actinobacteria (5.38%−6.89%) (Supplementary Fig. S8a). Significant enrichments of Acidobacteria and Actinobacteria were observed in RCP, while Acidobacteria, Chloroflexi, and Actinobacteria were significantly enhanced in RCT compared to monoculture (p < 0.05) (Supplementary Fig. S8b).

The relationship between microbial functional genes and the environment

Mantel tests revealed significant relationships between environmental factors and C/N/P cycling components in the rice-crayfish system (Fig. 5). Soil TP, DOC, and NH₄⁺-N significantly influenced functional genes, while soil SW, CEC, SOC, TN, TP, DOC, NH₄⁺-N, MBC, and MBN significantly affected the relevant microbial communities key to biogeochemical cycling processes involving C, N, and P.

Figure 5. 

Environmental drivers of (a) soil microbial functional genes, and (b) taxa in different rice farming systems. RM: rice monoculture; RCP: paddy field of rice-crayfish farming system; RCT: the peripheral trench of rice-crayfish farming system.

Soil environmental factors differentially influenced specific functional genes involved in C, N, and P cycling (Supplementary Fig. S9). For carbon cycling, significant negative correlations were observed between soil BD and vanA and accA genes, while soil pH showed positive correlations with pel, xylH, FUCA, and MAN2C1. Soil SW was positively correlated with vanB, coxM, korA, icd, lacA, abfA, and malQ. Soil SOC and CEC exhibited positive correlations with cbbR, fhs, and yvak. A significant positive correlation was observed between soil DOC and the genes mttB, chi, and malZ, whereas a negative correlation was found between DOC and pmoC, cbhA, and bglB (Supplementary Fig. S9a). In nitrogen cycling, significant negative correlations were identified between soil BD and pmoC, as well as between soil NAG and hao. Soil SW was positively correlated with narH and nxrB. Soil CEC and NH₄⁺-N were negatively correlated with norB. Soil TN and MBN were negatively correlated with norC but positively correlated with ureA (Supplementary Fig. S9b). Soil SW and NAG were negatively correlated with gltB. For phosphorus cycling, several soil phosphorus forms (TP, Olsen-P, MBP, AP) were inversely linked to phnM but positively linked to pstB. Separately, soil BD was correlated negatively with phoU and ugpA, and soil SW showed a similar negative relationship with phnX and phnA (Supplementary Fig. S9c).

In C cycling, negative correlations of Proteobacteria with pH, Olsen-P, NH₄⁺-N, TP, TN, MBC, MBN, MBP, and AP; negative correlations of Acidobacteria with CEC and SOC; positive correlations of Bacteroidota and Firmicutes with BD; negative correlations of Candidatus_Rokubacteria with pH and SW; positive correlation of Gemmatimonadetes with SOC (Supplementary Fig. S10a). Within nitrogen cycling, Proteobacteria exhibit a negative correlation with pH, while Candidatus_Rokubacteria display negative correlations with both pH and MBC (Supplementary Fig. S10b). Within phosphorus cycling, Proteobacteria exhibit a positive correlation with BD. Acidobacteria are positively correlated with both BD and DOC. Chloroflexi show positive correlations with SW and NAG. In contrast, Candidatus_Rokubacteria display negative correlations with SW, TP, TN, MBN, MBP, BG, NAG, and AP. Bacteroidota and Gemmatimonadetes are positively correlated with multiple factors, including SW, CEC, SOC, TN, TP, NH₄⁺-N, MBC, MBN, MBP, BG, NAG, and AP (Supplementary Fig. S10c).

As illustrated in Fig. 6, this study examined the interactions between genes participating in the C, N, and P cycling processes of the rice-crayfish co-culture system. The top 20 genes with the highest degree centrality scores included five C-cycling genes (korA, pccA, coxS, mttB, coxM), eleven N-cycling genes (nirB, nifK, nifD, norB, nifH, nirD, narG, nrfA, gltD, napA, GDH2, narH), and four P-cycling genes (phoR, pstS, pstC, phnE). According to Spearman correlation analysis, significant correlations were observed among C-, N-, and P-related genes. Among the high-degree centrality genes, korA (rTCA cycle) showed significant negative correlations with nifK and norB (denitrification genes) as well as with nifD, nifK, and nifH (nitrogen fixation genes). mttB (methanogenesis gene) was significantly negatively correlated with GDH2 (nitrogen degradation gene). Additionally, korA (rTCA cycle gene) was significantly positively correlated with phoR (P regulation gene) but negatively correlated with pstC and phnE (P transport genes). gltD and gltB showed significant negative correlations with pstC and phnE, whereas nifK, norB, narG, and narH (denitrification genes) showed significant negative correlations with phoR and pstS.

Figure 6. 

Network analysis between genes related to (a) carbon, nitrogen, and phosphorus cycling genes, and (b) degree centrality scores for different rice farming systems. RM: rice monoculture; RCP: paddy field of rice-crayfish farming system; RCT: the peripheral trench of rice-crayfish farming system.

This study explored the RC farming system, aiming to deeply investigate its sustainability from a deep perspective on biogeochemical cycles. Simultaneously, it revealed the key microbial driving mechanisms driving the cycling processes in this integrated system. Soil microbiota serves as a vital indicator of ecosystem health, where rich microbial diversity helps maintain ecosystem stability and promotes efficient nutrient cycling^[28,29]. The significant increase in soil microbial diversity within the RC system suggests that the observed structural and functional shifts in the microbiota are a direct result of the agricultural system conversion, rather than minor environmental fluctuations (Fig. 1; Supplementaty Fig. S2). The study found that rice-crayfish farming involves diverse nutrient inputs, including straw incorporation, feed supply, crayfish excretion, and fertilizer application, which significantly enhance soil organic matter content and subsequently provide energy for microorganisms, promoting microbial metabolism and reproduction^[30]. This predation pressure may also result in a simplified community composition of soil invertebrates and protists in the integrated system compared to monoculture^[6,27]. This top-down effect can propagate through the food chain, reducing predation pressure on microorganisms and, consequently, enhancing the diversity of fungi and archaea at lower trophic levels in the rice-crayfish system^[31]. The significantly higher microbial taxonomic and functional diversity in RCP soils compared to RCT areas (Figs 2−4). It is mainly attributed to long-term deep flooding, which enhances soil-reducing conditions, decreases aerobic microbial abundance, and ultimately reduces microbial taxonomic diversity in RCT areas.

The rice-crayfish farming altered soil C degradation and fixation processes

Soil microorganisms actively participate actively in organic matter decomposition and transformation via diverse metabolic pathways, playing a pivotal role in the carbon cycle of paddy ecosystems^[28,29], and serving as a key component of the system's carbon budget that modulates soil carbon storage and turnover^[17]. Notably, the relative abundance of chitin-degradation genes was significantly elevated in the rice-crayfish system relative to the rice monoculture (Fig. 2; Supplementary Fig. S2). This increase is primarily attributed to the frequent molting of crayfish over ten times during growth, as their chitin-rich shells stimulate the proliferation of chitin-degrading microorganisms, thereby elevating the abundance of related functional genes^[30−32]. The distinct chitin (chi) degradation gene abundance between RCP and RCT areas further reflects the regulatory effects of crayfish activity and dissolved oxygen levels (Fig. 2). In contrast, the lower relative abundance of starch degradation genes (malQ and malZ) in the rice-crayfish system (Fig. 2) may stem from uneven soil organic matter distribution caused by crayfish feeding and excretion, which alters the ecological niches of starch-degrading microorganisms^[33,34]. Regarding methane oxidation, all rice cultivation systems showed higher relative abundance of aerobic methane oxidation genes than anaerobic ones (pomA, pomB, pomC) (Fig. 2); however, continuous deep flooding in the rice-crayfish system enhanced anaerobic methane oxidation gene abundance, especially in RCP areas^[35,36]. In the rice-crayfish system, the reductive acetyl-CoA pathway and CO oxidation were suppressed in the rice-crayfish system (Fig. 2). This may be attributed to heterotrophic microorganisms rapidly utilizing available carbon from crayfish excrement and rice straw, thereby reducing heterotrophic carbon fixation^[37,38]. Conversely, these organic inputs significantly increased functional gene abundance associated with the rTCA pathway, as the heterotrophic rice-crayfish system provides sufficient carbon and energy sources to promote autotrophic carbon metabolism^[39,40].

The rice-crayfish farming altered soil N and P cycling processes

Across different rice cultivation systems, nitrogen degradation accounts for the largest proportion of the soil nitrogen cycle (Supplementary Fig. S4), indicating the potential of paddy soils to convert both fertilizer- and inherent-soil nitrogen into NH₄⁺-N for crop uptake^[41−43]. The rice-crayfish system exerted a divergent effect on the nitrogen-cycle gene pool, fostering the abundance of genes for degradation, denitrification, and DNRA at the expense of those responsible for nitrification (Fig. 3; Supplementary Fig. S4). This suggests that the integrated system accelerates nitrogen turnover in paddy fields^[43,44]. The elevated abundances of glnA, gltB, gltD, and arcC in the rice-crayfish system promote soil nitrogen degradation and increase NH₄⁺-N production (Fig. 3). Under high NH₄⁺-N conditions, this facilitates amino acid formation, which microorganisms can subsequently assimilate^[45,46]. Continuous deep flooding in the rice-crayfish system significantly reduced the relative abundance of pomA, pomB, and pomC genes, thereby inhibiting the conversion of NH₄⁺-N to NO₃⁻-N and promoting the accumulation of soil NH₄⁺-N^[47]. The enhanced denitrification in the rice-crayfish system also contributed to increased N₂O emissions^[48]. Studies have found that nitrogen loss through soil denitrification is generally comparable to nitrogen emissions as N₂O^[49]. Furthermore, the relative gene abundances of napA, napB, norB, and nirK were significantly increased in the rice-crayfish system, primarily due to abundant exogenous nitrogen inputs and enhanced reducing conditions^[6,50]. The DNRA process converts nitrogen from its highest oxidation state (NO₃⁻-N) to its lowest reduced state (NH₄⁺-N) while simultaneously harvesting energy from soil organic carbon. These key mechanisms effectively improve the retention of available nitrogen in soil, thereby promoting its absorption and utilization by crops^[51].

Different farming systems influenced soil P availability (Supplementary Table S5), with soil microorganisms playing a major role in regulating phosphorus cycling^[52]. Compared to rice monoculture, the rice-crayfish system appeared to enhance the genetic potential for mineralizing organic phosphorus, indicating a greater microbial capacity to access this nutrient pool (Supplementary Fig. S6). Elevated soil C and N levels in the integrated system could stimulate microbial investment into the genetic machinery required for organic P mineralization (Supplementary Figs S7 & S8). The mineralization of organically bound phosphorus is often enhanced under conditions of elevated carbon and nitrogen availability, as supported by multiple studies^[53,54]. Another study indicated that additional carbon inputs to soil can also increase available P content^[55,56]. Furthermore, elevated soil C and N can create a favorable niche for bacteria involved in P transformations, thereby selecting for a more abundant community (Supplementary Fig. S9). Genes critical for organic phosphorus mineralization, particularly phoD, phnA, and phnP, were markedly more abundant in the rice-crayfish system than in monoculture soils (Supplementary Fig. S7). The phoD, phnA, and phnP genes are widely recognized as key drivers of the organophosphorus mineralization process^[36,57], while phoN, phnA, and phnG influence phosphatase activity, thereby affecting soil available phosphorus content^[58]. The rice-crayfish system exhibited the highest abundance of inorganic P solubilization genes, indicating a greater advantage in inorganic phosphorus dissolution^[59,60]. Existing research has shown that inorganic phosphorus solubilization is an essential pathway for increasing bioavailable phosphorus in soil^[61]. The functional divergence in inorganic P solubilization was attributable to changes in gcd and gppa gene abundance^[62], which were significantly elevated in the rice-crayfish system (Supplementary Fig. S7). These genes are considered hallmark genes for this specific microbial function^[63]. The machinery for phosphorus assimilation centers on the phoB two-component system, alongside the high-affinity PstSCAB system and complemented by the low-affinity inorganic phosphate transporter (pit)^[64]. Although the abundance of pstB was significantly elevated in rice-crayfish farming soils, its corresponding phoR did not exhibit a corresponding change (Fig. 4). This imbalance may have led to impaired signal transduction between pathways^[53], thereby affecting phosphorus regulation efficiency. This phenomenon may stem from differences in their sensitivity to soil available phosphorus concentrations^[27,65].

Responsive patterns of C, N, and P cycle-related microbial taxa under rice-crayfish farming

The functional potential of microbial communities is highly correlated with microbial taxa (Figs 5 & 6). Here, bacteria dominated the soil microbiota in the rice-crayfish system, while other microorganisms (archaea and fungi) were non-dominant, consistent with the findings of Wu et al.^[1]. According to the oligotrophic-copiotrophic theory, the observed increase in functional microbe abundance in this study may be attributed to enhanced copiotrophic bacterial populations, such as phoD-carrying Actinobacteria^[66]. Furthermore, the predicted functional gene classification highlighted the significance of Proteobacteria and Actinobacteria in C, N, and P cycling during rice-crayfish co-culture^[67]. Members of Proteobacteria are recognized for their high nutrient-use potential and serve as major drivers of functional changes in fertile soils^[68]. Some studies indicate that the r- and K-selection growth strategy classification may not apply to Actinobacteria, making their response to labile carbon changes difficult to predict^[69]. Previous studies have demonstrated Actinobacteria enrichment in nutrient-rich soils^[70]. Acidobacteria are typically considered oligotrophic bacteria^[16], yet, in this study, Acidobacteria carrying C-, N-, and P-cycling functional genes showed significantly higher abundance in the nutrient-richer rice-crayfish soils compared to monoculture (Fig. 5). This is primarily because the rice-crayfish system could maintain neutral soil pH conditions, which favor Acidobacteria growth^[63]. Studies have identified Acidobacteria as thriving in neutral pH environments^[62].

Soil TN, TP, and DOC were identified as primary factors influencing C, N, and P cycling processes in rice-crayfish system, including both functional genes and microbial taxa (Fig. 5). This is primarily attributed to the input of diverse nutrients during rice-crayfish co-culture—including straw, fertilizers, feed, and crayfish excreta—which provide abundant exogenous carbon, nitrogen, and phosphorus, significantly enhancing soil TN, TP, and DOC concentrations^[71]. Soil carbon, nitrogen, and phosphorus are pivotal regulators of ecosystem productivity, functioning through their influence on microbial dynamics, organic matter turnover, and the cycling of essential nutrients^[72]. Functioning as a central node, soil TN modulates carbon sequestration and decomposition processes in parallel with N and P cycling^[73]. In the context of rice-crayfish co-culture, P availability emerges as a key modulator of this integrated nutrient network^[74]. Based on these findings, it was hypothesized that soil TN, TP, and DOC likely constitute key environmental factors that microbially drive C, N, and P cycling in the rice-crayfish farming system.

The genetic machinery mediating the transformations of C, N, and P demonstrates tightly coupled interactions, characterized by both synergistic and antagonistic relationships (Fig. 6). Nodes with higher degree centrality scores were identified as dominant elements maintaining network connectivity^[49,75]. Despite their lower numerical representation compared to nitrogen cycle genes, genes involved in carbon and phosphorus cycles held a higher degree of centrality within the co-occurrence network. This indicates they exerted a disproportionately strong influence, underscoring their critical role in regulating nitrogen transformations^[49]. In the rice-crayfish system, NADH encoded by GDH not only participates in soil N mineralization but may also promote methane oxidation^[76]. Additionally, reducing nitrogen oxides through denitrification has been reported to enhance carbon fixation under anaerobic conditions, while this study detected negative correlations between denitrification-related N-cycle genes and soil carbon sequestration-related C-cycle genes in the rice-crayfish system^[77]. Microorganisms exhibiting high denitrification and carbon fixation capabilities through the rTCA cycle may indicate synergistic interactions between these processes^[76]. The coordinated expression of genes linked to carbon and nitrogen cycles fine-tunes nutrient and metabolite levels. This co-regulation serves to maintain metabolic homeostasis, ensuring that the assimilation of nitrogen is synchronized with the availability of replenished carbon^[70]. Coupling between the C and P cycles occurs principally through inositol phosphate metabolism, which produces acetyl-CoA (catalyzed by mmsA) and glyceraldehyde 3-phosphate (catalyzed by TPI). These products subsequently enter the core pathways of central carbon metabolism^[50,78]. A critical step in the C-P cycle linkage is the hydrolysis of inositol phosphates by phytase (phy, appA), which generates bioavailable phosphate for plant and microbial use, thus acting as a key nexus between these elemental cycles^[79].

In the rice-crayfish farming system, the coupled cycling of C, N, and P is governed by the functional genetic potential and the taxonomic composition of the soil microbiota. Specifically, a marked increase was observed in the abundance of genes governing methane oxidation, nitrogen degradation, denitrification, anammox, organic phosphorus mineralization, and phosphorus transport within the integrated system. Across the integrated field, the system acted to enhance carbon retention, which in turn substantially supported the concurrent cycling and transformation of nitrogen and phosphorus. The findings indicate that the rice-crayfish farming system could be a promising, environmentally friendly strategy, capable of enhancing carbon sequestration and increasing the availability of NH₄⁺-N for plants by stimulating aerobic DNRA and denitrification processes. The results of this study pinpoint soil TN, TP, and DOC as the principal environmental determinants behind the observed biogeochemical transformations of C, N, and P. A tight functional coupling was observed among the genes involved in these cycles, and their interrelationships were altered by variations in nutrient input. Consequently, these findings offer a theoretical foundation for improved management of fertilization and the microbial regulation of soil nutrient cycles.