Sustainability of the rice-crayfish co-culture aquaculture model: microbiome profiles based on multi-kingdom analyses

While the rice-crayfish culture (RCFP) model, an important aquaculture model in Asia, is generally considered a sustainable model, its sustainability in terms of microbial community profiles has not been evaluated. In this study, multi-kingdom analyses of microbiome profiles (i.e., bacteria, archaea, viruses, and eukaryotes) were performed using environmental (i.e., water and sediment) and animal gut (i.e., crayfish and crab gut) microbial samples from the RCFP and other aquaculture models, including the crab-crayfish co-culture, crayfish culture, and crab culture models, to evaluate the sustainability of the RCFP systematically. Results showed that RCFP samples are enriched with a distinct set of microbes, including Shewanella, Ferroplasma, Leishmania, and Siphoviridae, when compared with other aquaculture models. Additionally, most microbes in the RCFP samples, especially microbes from different kingdoms, were densely and positively connected, which indicates their robustness against environmental stress. Whereas microbes in different aquaculture models demonstrated moderate levels of horizontal gene transfer (HGT) across kingdoms, the RCFP showed relatively lower frequencies of HGT events, especially those involving antibiotic resistance genes. Finally, environmental factors, including pH, oxidation–reduction potential, temperature, and total nitrogen, contributed profoundly to shaping the microbial communities in these aquaculture models. Interestingly, compared with other models, the microbial communities of the RCFP model were less influenced by these environmental factors, which suggests that microbes in the latter have stronger ability to resist environmental stress. The findings collectively reflect the unique multi-kingdom microbial patterns of the RCFP model and suggest that this model is a sustainable model from the perspective of microbiome profiles.

on experience and lacks systematic evaluation protocols. Because microbial communities play a basic material and energy cycle driving role in various ecosystems [8][9][10], evaluating the sustainability of different aquaculture models in terms of their microbial profiles is necessary. The sustainability of an aquaculture model at the microbial community level could be represented and measured by the active interactions of microbes in the community, the frequency of gene transfer, especially of antibiotic resistance genes (ARGs), and ability of the communities to withstand environmental stress.
Crayfish farming is an important form of aquaculture in Asia, because crayfish is an excellent source of protein and essential amino acids [11]. Efforts to improve the sustainability of crayfish farming have led to the development of the rice-crayfish co-culture (RCFP) model. This model is a dominant model of crayfish farming and contributes ~ 90% of the total crayfish production [12] by taking advantage of the synergistic effects of co-cultured species [13]. It also improves rice yields and produces extra-economic profits [14]. In the RCFP model, alternating water and lower inputs of pesticides and chemical residues provide a green production environment (e.g., higher water quality, soil fertility, and dissolved oxygen contents), thereby reducing the risk of disease [13][14][15][16]. However, how the microbial profile of this model differs from that of other aquaculture models, as well as its robustness against the surrounding environment, is as yet unknown.
The microbial communities in an aquaculture co-culture model usually consist of bacteria, archaea, viruses, and eukaryotes. Further, their interactions include predation (e.g., some protists feed on bacteria), pathogenicity (e.g., microbes could interact with pathogens), and parasitism (e.g., some viruses live by parasitizing bacteria) [17]. These multi-kingdom species and their interactions jointly maintain the stability of the microbial community. However, the current knowledge on aquaculture microbial communities is mainly based on bacterial, fungal, or viral communities or the combination of two kingdoms [6,[18][19][20]. Studies profiling aquaculture microbial communities from the perspective of multiple kingdoms, including bacteria, archaea, viruses, and eukaryotes, are scarce. These knowledge gaps limit our in-depth understanding of microbial profiles in aquaculture models at the multi-kingdom level.
Aquaculture microbial communities are strongly influenced by environmental factors, including human activities, which often manifests as an increase in antimicrobial resistance [21]. When ARGs are transmitted into humanassociated pathogens, they may pose a great environmental risk [22]. Many ARGs and other potentially harmful functional genes are spread through the environment by horizontal gene transfer (HGT) [23][24][25]. However, the current knowledge on HGT events in aquaculture models, especially those involving multi-kingdom microbes and ARGs, is limited. Moreover, the combined influence of multiple environmental factors on microbial communities in different aquaculture models remains unclear.
In the present study, we analyzed the microbial communities of various aquaculture models to achieve a systematic assessment of the sustainability of the RCFP and other aquaculture models (Fig. 1). We focused on aquaculture model-specific microbial community patterns (Fig. 1C), multi-kingdom interactions (Fig. 1D), HGT events (Fig. 1E), especially those involving ARGs, and environmental factors driving microbial community divergence (Fig. 1F). To this end, we collected environmental (i.e., water and sediment) and animal gut (i.e., crayfish and crab gut) microbial samples from the crabcrayfish culture (CCFP), crab culture (CP), crayfish culture (CFP), and RCFP models (Fig. 1A, B). Specifically, we aim to answer four questions: (1) What are the distinct microbial community patterns of different aquaculture models at the multi-kingdom level? (2) How do bacteria, archaea, viruses, and eukaryotes interact with each other in different culture models? (3) How do HGT events, especially those involving ARGs, differ among different aquaculture models? (4) How do environmental factors influence the microbial communities in different aquaculture models? We found that the unique microbial profiles of the RCFP model maintain its stability when faced with environmental pressure, as reflected by its water, sediment, and crayfish gut microbial communities. The results confirm that the RCFP model is a sustainable model from the perspective of microbiome profiles and provide new insights into sustainable aquaculture.

Sample collection
Samples from water, sediment, crayfish, and crab were collected from Honghu farms (29.92° N, 113.49° E), in Hubei Province, China. A total of 20 water and 20 sediment samples were collected from four different types of aquaculture models, including the CCFP, CFP, CP, and RCFP models (Fig. 1A, B). Five parallel samples were collected for each culture model (Figs. 1A, B, 2A and Additional file 1). The water and sediment samples were divided into two parts. One part was used to measure the effects of environmental factors according to our previous study [26]. The other part was utilized for metagenomic sequencing. Ten crab samples were collected from the CP (n = 5) and CCFP (n = 5) models, and 17 crayfish samples were collected from the CCFP (n = 5), CFP (n = 6), and RCFP (n = 6) models.
The intestinal contents of the animals were aseptically extracted using the conventional anatomical method and placed in a sterile centrifuge tube (5 mL). In total, 67 samples were prepared and stored at -80 ℃ for metagenomic sequencing.

Sample processing and sequencing
The total DNA of water filter membranes and sediment samples was extracted using the FastDNA ™ SPIN Kit (MP, USA) following the manufacturer's instructions. DNA extraction from crab gut and crayfish gut samples The workflow for systematically assessing the sustainability of the RCFP and other aquaculture models from the perspective of aquaculture model-specific microbial community patterns, multi-kingdom interactions, horizontal gene transfer (HGT) events, and environmental factors driving microbial community divergence. A Water, sediment, crayfish gut, and crab gut samples were collected from four representative aquaculture models: crab-crayfish culture model (CCFP), crab culture model (CP), crayfish culture model (CFP), and rice-crayfish co-culture model (RCFP). B The inhabiting habitats for the microbial community. C The aquaculture model-specific microbial community patterns of water, sediment, crayfish gut, and crab gut habitats. D Comparison of multi-kingdom interactions (bacteria, archaea, viruses, and eukaryotes) between RCFP and other aquaculture models in water, sediment, and crayfish gut habitats, respectively. E HGT events detecting across kingdoms between RCFP and other aquaculture models in water, sediment, and crayfish gut habitats, respectively. F The response of the multi-kingdom microbial community to environmental factors Fig. 2 Microbial community composition and diversity. A Overview of the sampling habitats and aquaculture models. B Viral community composition. C Bacterial community composition at the genus level. Unweighted paired-group method with arithmetic means (UPGMA)-based hierarchical clustering (Bray-Curtis distance) was used to analyze the microbial community structure. Each column of Fig. 1B is vertically aligned to the columns of Fig. 1C. D Bacterial, archaeal, viral, and eukaryotic diversities across the crab gut, crayfish gut, sediment, and water habitats. E Comparison of microbial diversity across habitats via PCoA using Jaccard coefficients as the distance measurement. The 90% confidence intervals of each group are also shown in the background. Samples and tree branches colored in green, orange, blue, and yellow represent water, sediment, crayfish gut, and crab gut samples, respectively. CP: crab culture model; CFP: crayfish culture model; CCFP: crab-crayfish mixed culture model; RCFP rice-crayfish co-culture model; "uc'': unclassified. "*":p < 0.1; "**":p < 0.05; "***":p < 0.01; "****":p < 0.001; ns: not significant   (~ 1 g) was performed using the QIAamp DNA Stool Mini Kit (Qiagen, Germany). The extracted DNA was treated using the NEBNext Ultra DNA Library Prep Kit for Illumina (NEB, USA) for whole-genome amplification and library preparation. Finally, paired-end sequencing was performed on an Illumina HiSeq X Ten platform. Except for 9 samples that could not be sequenced, 58 metagenomic samples were sequenced for metagenomic analysis (Additional file 1: Table S1). A total of 2.52 billion pair-end reads of raw sequencing data (Additional file 1: Table S2) were first evaluated using FastQC (version 0.11.6) [27]. Low-quality reads and adaptors were then trimmed by Trimmomatic (version 0.38) [28] to eliminate reads less than 100 bp in length, adapters, leading or trailing bases with Phred base quality (BQ) scores of < 20, and strings of every five bases with an average BQ score of < 25. After quality control, we obtained 2.41 billion high-quality pair-end reads (Additional file 1: Table S2).

Microbial classification and diversity
To obtain the taxonomical composition of different models, we annotated high-quality reads (94.65% of the raw reads) from each sample by using MetaPhlAn2 (version 2.6.0) [29]. For each kingdom annotation, we filtered other kingdoms using "-ignore_bacteria -ignore_ archaea -ignore_viruses -ignore_eukaryotes" items. Alpha diversity based on the Shannon index was determined using the diversity() function in R "vegan" package (version 2.6-2). Principle coordinate analysis (PCoA) based on Jaccard coefficients as a distance measurement was used to cluster samples according to their group with 90% confidence intervals.

Detection of indicator microbes for different culture models
The indicator value of microbes in each culture model was calculated to characterize aquaculture model-specific microbial compositions by using the indval() function in R "labdsv" package (version 2.0-1) [30]; here, the frequency and relative abundance of each microbe were considered. In our study, only the microbes of an aquaculture model with the highest indicator values at p < 0.05 were considered indicator microbes for this model.

Measurement and analysis of environmental factors
Environmental factors influencing the water and sediment samples were determined following our previously published study [26]. These environmental factors included physicochemical factors, such as temperature, total nitrogen (TN), total phosphorus (TP), sediment volume-weight (Sed_VM), and moisture content (MC), as well as antibiotic factors, such as roxithromycin (ROX), erythromycin (ERY), and erythromycin derivative 1 (ERY1). In total, we detected 18 physicochemical data and 5 antibiotic factors for water samples, and 9 physicochemical factors and 4 antibiotic factors for sediment samples (Additional file 1: Figs. S1, S2).

Analysis of multi-kingdom co-occurrence networks
Microbes with a relative abundance of ≥ 2.00% and coverage of > 20% samples were selected for multi-kingdom interaction analysis to investigate microbial correlations among bacteria, archaea, viruses, and eukaryotes in different culture models. We used Spearman's correlation analysis to calculate correlations among the four kingdoms. Only Spearman correlations of ≥ 0.65 or ≤ − 0.65 with p < 0.05 were considered as strong correlations and visualized in Cytoscape (version 3.8.1) [31].

Analysis of horizontal gene transfer events
We investigated variations in ARGs and HGT events across different habitats and aquaculture models. We mainly used two methods to detect genes transferred among different kingdoms in water, sediment, and crayfish gut. In the first method, ARGs were detected using DeepARG (version 1.0.2) [32] and then mapped to microbial contigs using BLASTN (version 2.7.1+) with "-task megablast -evalue 1e-10. " These contigs were annotated for taxonomy assignment using Kraken2 (version 2.0.8beta) [33]. In the second method, the genes involving HGT events were detected using MetaCHIP (version 1.10.4) [34] with "-r pcofgs" to predict reference-independent HGT events at the community level. Here, highquality bins were built using MetaWRAP (version 1.2.2) [35]. A total of 496 bins (0.90 Gb) with an average length of 1,806,454 per bin were assembled from the 58 contigs (4.80 Gb; an average length of a contig: 2306). Then, all bins obtained from the metagenomic data, their taxonomic classifications, and their group information were inputted into MetaCHIP.

Analysis of the effect of environmental factors on microbial profiles
We used detrended correspondence analysis (DCA) to judge the major axis length. Because the lengths of the first four major axes were less than three, redundancy analysis (RDA) was used to reveal the effects of the physicochemical and antibiotic factors of water and sediment habitats on microbial community ordination using the R package "vegan" (version 2.6-2). All environmental factors were tested in RDA analysis using the envfit() function with 999 permutations in R "vegan" package (version 2.6-2), while those with p < 0.05 were considered as significant factors influencing the corresponding microbial community.

Statistical analysis
Differences among the water, sediment, crayfish gut, and crab gut habitats, as well as the RCFP and other culture models (i.e., the non-RCFP group), were determined using the Wilcoxon test in R "stats" package (version 4.0.5). Variations in habitat among different aquaculture models were tested by analysis of variance with the leastsignificant difference test. The p value was adjusted using the Benjamini and Hochberg method.

Differences in environmental and animal gut microbial community patterns at the multi-kingdom level
Microbial community samples from different habitats, including water, sediment, crayfish gut, and crab gut, were collected from the CCFP, CFP, CP, and RCFP models ( Fig. 1A, B, and 2A). After demultiplexing, adaptor trimming, and quality control, 2.41 billion high-quality pair-end reads were used for downstream analysis (Additional file 1: Table S1 and Fig. S3). The microbial community composition varied across the different habitats. Through multi-kingdom analysis, we detected an average abundance of 61.84% bacteria, 4.62% archaea, 30.13% viruses, and 2.94% eukaryotes in each sample. The proportion of clean reads that could be assigned to different taxonomic levels in each habitat was also provided in Additional file 1: Table S3. We also found that environmental (sediment and water) microbial communities were dominated by bacteria, whereas animal gut microbial communities of crayfish and crab were dominated by both bacteria and viruses ( Fig. 2 and Additional file 1: Figs. S4-S6). Water and sediment habitats harbored more archaea than those of crab gut and crayfish gut habitats, whereas animal gut habitats contained more eukaryotes than those of water and sediment habitats (Additional file 1: Figs. S4-S6). All detected microbial taxa in each habitat across these four kingdoms were also provided in Additional file 1: Fig. S7 and Table S4. Unweighted paired-group method with arithmetic means (UPGMA) analysis illustrated that the water, sediment, crayfish gut, and crab gut samples were respectively clustered together ( Fig. 2 and Additional file 1: Fig. S6), thereby indicating that the microbial community composition largely depends on the habitat. In terms of microbial alpha diversity, water and sediments showed higher bacterial and viral diversities compared with animal guts. However, the diversities of archaeal and eukaryotic organisms did not differ among these four habitats ( Fig. 2D and Additional file 1: Fig. S8). Beta diversity also revealed that the animal gut microbial communities differed from those obtained from the environment (Fig. 2E).

Multi-kingdom co-occurrence networks in different aquaculture models
To compare multi-kingdom co-occurrence networks in the RCFP models with those in other aquaculture models, we categorized all samples into two groups, namely, RCFP and non-RCFP; here, the RCFP group includes samples from the RCFP, while the non-RCFP group includes samples from all other aquaculture models.
The RCFP network was dominated by dense and positive multi-kingdom interactions in both water and sediment habitats. The RCFP model demonstrated higher network complexity in water habitats than the non-RCFP models. Specifically, in water habitat, the RCFP model showed 66.67% more nodes (40 vs. 24), 436.36% more edges (177 vs. 33), 83.3% greater network density (0.11 vs. 0.06), and 221.82% greater node connectivity (node degree: 8.85 vs. 2.75; all microbes were positively correlated) compared with the non-RCFP group. Details of the network properties obtained in water habitat were provided in Fig. 3A and Additional file 1: Table S9. We also observed that Polynucleobacter and Lactobacillus were positively correlated with most microbes, including bacteria, archaea, viruses, and eukaryotes, in the RCFP model, while the weak correlations were found in these microbes in the non-RCFP group in water habitat ( Fig. 3A and Additional file 1: Fig. S9A). In sediment, although the number of microbes (46 nodes) used for RCFP network construction was ~ 6.52% fewer than that used for the non-RCFP group (49 nodes), the resulting network of the RCFP model was denser (0.134 vs. 0.04), and closer (average node degree: 11.00 vs. 2.00; the number of edges: 231 vs. 26) compared with that of the non-RCFP group (Fig. 3B, Additional file 1: Fig. S9B, and Additional file 1: Table S9). Moreover, compared with the RCFP group, we found weaker interactions between Corynebacterium and other microbes in the non-RCFP group. Taken, together, these results reflect the dense and close multi-kingdom interactions in the RCFP model in sediment habitats, consistent with the interactions found in water habitats.
The RCFP group also showed dense and positive networks across multiple kingdoms in the crayfish gut habitat. The network density of the RCFP group (network density: 0.24, node degree: 21.78) was approximately 2.5 times greater than that of the non-RCFP group (network density: 0.07, node degree: 6.25) (Fig. 3C, Additional file 1: Fig. S9C, and Additional file 1: Table S9). Compared with those in the RCFP network, the correlations among Lactobacillus, Streptococcus, and other microbes decreased in the non-RCFP group (Fig. 3C and Additional file 1: Fig. S9C). We found that Citrobacter was positively correlated with other kingdoms in the RCFP models, but these interactions were reduced or negative in the non-RCFP group. Besides, in the non-RCFP network, we determined that Mycoplasma was densely correlated with bacteria, archaea, viruses, and eukaryotes, which was not observed in the RCFP network ( Fig. 3C and Additional file 1: Fig. S9C). Moreover, the network analysis has elucidated the specificities of the viral host in different aquaculture models (Additional file 1: Fig. S10). Collectively, these results indicated that the close and dense multi-kingdom interactions in the RCFP model may lead to the improved robustness of its microbial communities against environmental stress.

Insights into HGTs within microbial communities ARG transfer among multiple kingdoms
The RCFP microbial communities demonstrated few ARG-associated HGT events in water and sediment habitats. We categorized samples into the RCFP and non-RCFP groups. HGT analysis showed that many ARGs are transferred among different microbes in water (5689 ARGs involved in HGT events, accounting for 14.00% of the total ARGs identified in water samples) and sediment (5091 ARGs involved in HGT events, accounting for 12.54% of the total ARGs identified in these samples) habitats ( Fig. 4 and Additional file 1: Figs. S11-S12). In water habitats, 60 genes belonging to 8 ARG types were transferred across bacteria, archaea, and viruses, accounting for 34 microbes and 73 HGT events (Fig. 4A and Additional file 1: Fig. S11). In sediment habitats, four ARGs participated in 4 HGT events between bacteria and archaea (Additional file 1: Fig. S12). Among the aquaculture models studied, the RCFP model contributed the lower frequencies of HGT events both detected in water (1.4 HGT events per sample in the RCFP group; 4.7 HGT events per sample in the non-RCFP group) and sediment (zero HGT events per sample in the RCFP group; 0.29 HGT events per sample for the non-RCFP group) habitats ( Fig. 4A and Additional file 1: Fig. S12).
RCFP microbial communities also contributed few ARG-involved HGT events in the crayfish gut habitat. Among 782 detected ARG-associated genes, three were annotated as MLS|macB and involved in two HGT events across kingdoms (Fig. 4B). Specifically, MLS|macB transferred between Methanosphaera and Brevibacillus in the RCFP group, as well as between Methanosphaera and Spiroplasma in both RCFP and non-RCFP groups. Together, these results implied that HGT events might play an important role in balancing the function of microbial communities and that RCFP microbial communities were less affected by ARG-associated HGT events than non-RCFP microbial communities.

Functional genes participating in HGT events
The RCFP microbial community contributed few HGT events in both water and sediment habitats. Using MetaCHIP analysis, we detected 1.2 and 1.3 HGT events per sample in the RCFP and non-RCFP groups, respectively, in water habitats (Additional file 1: Fig. S13). In the RCFP group, Polynucleobacter contributed the largest number of HGT events, followed by uc_Rhodocyclaceae (Additional file 1: Fig. S13A). We also found that Polynucleobacter contributed to fewer HGT events (2 HGT events in all samples) in the non-RCFP group (Additional file 1: Fig. S13B). Moreover, we found that Aeromonas, a virulent and antibiotic-resistant pathogen [36,37] . 3 Multi-kingdom co-occurrence networks of samples in the RCFP and non-RCFP groups. Networks in A water, B sediment, and C crayfish gut habitats are presented, respectively. Microbes with a relative abundance of ≥ 2% and coverage of > 20% samples were used to construct the multi-kingdom co-occurrence networks. Only Spearman correlations of ≥ 0.65 or ≤ -0.65 with p < 0.05 were considered strong correlations and visualized in the network. The nodes in green, blue, purple, and orange represent bacteria, archaea, viruses, and eukaryotes, respectively. Red edges indicate positive correlations, whereas blue edges reflect negative correlations events in the non-RCFP group (Additional file 1: Fig.  S13B).
In the crayfish gut habitat, the RCFP group contributed fewer HGT events compared with the non-RCFP group. Using MetaCHIP analysis, we detected only 4 HGT events across Shewanella, Enterococcus, and other unclassified microbes in all samples of non-RCFP model at the genus level (Additional file 1: Fig. S14). We also found that Shewanella and Enterococcus were among the top 10 bacterial hosts of ARGs (Additional file 1: Fig. S11D), and thus, could potentially accelerate ARG spreading. Taken together, the limited number of HGT events, especially ARG-associated HGT events, of RCFP microbial communities in both environmental and animal gut habitats suggested that the RCFP model was an environment-friendly aquaculture model.

Effect of environmental factors on microbial communities
The RCFP was characterized by low nitrogen levels and antibiotic concentrations, especially in water habitats. In this study, the effects of environmental factors, including physicochemical parameters and antibiotics, on the microbial community structures were assessed. In water habitats, the RCFP model had higher salinity and oxidation-reduction potential (ORP) compared with other aquaculture models; by contrast, the concentrations of nitrate-nitrogen (NO 3 − -N), TN, turbidity, and erythromycin derivative 2 (EYR2) of the RCFP model were lower than those of other aquaculture models (p < 0.05) (Additional file 1: Fig. S1). These findings suggested low levels of nitrogen and antibiotic concentrations in the RCFP model. In sediment habitats, the pH of the RCFP model was higher compared with other aquaculture models (p < 0.05), but no significant differences in other environmental factors were found across the aquaculture models (Additional file 1: Fig. S2).
The effects of environmental factors on microbial communities across different aquaculture models significantly differed. RDA revealed that, in water habitats, temperature (r 2 = 0.510, p = 0.005) was a primary environmental factor affecting non-RCFP bacterial communities; other environmental factors, including pH and nitrogen level, also significantly influenced non-RCFP bacterial communities (Fig. 5A). By comparison, ORP (r 2 = 0.546, p = 0.004) was the primary factor influencing the RCFP bacterial community (Fig. 5A). No other environmental factors significantly influenced (p < 0.05) the RCFP bacterial community. These environmental factors did not significantly affect the microbial community in sediment habitats across aquaculture models. Collectively, the results indicated that the microbial community of the Fig. 4 ARG-associated HGT events among kingdoms. ARGs transferring in A water and B crayfish gut habitats are shown. Each band in the inner or outer circle represents a microbe at the genus level, the name of which is colored according to the kingdom type: bacteria (black), archaea (red), and viruses (blue). Bands among bacterial, archaeal, and viral genomes mean the ARGs are involved in HGT events across kingdoms. Different band colors represent different microbes bearing genes involved in HGT events RCFP was more stable under environmental stresses than those of other aquaculture models.
The microbial community of the RCFP model in the crayfish gut was only slightly affected by the studied environmental factors. In the crayfish gut habitat, chlorophyll (r 2 = 0.268, p = 0.065) was the main factor influencing non-RCFP bacterial community structures (Fig. 5B). By contrast, we did not detect any environmental factor significantly affecting the RCFP microbial community (Fig. 5B,  D). These results indicated that the microbial community in the RCFP model was more stable under environmental stresses in the crayfish gut than the microbial communities of other aquaculture models.

Differences in microbial community patterns across aquaculture models
To the best of our knowledge, this study is the first to assess the microbial community of freshwater aquaculture wetlands from the perspective of multi-kingdoms. In terms of relative abundance, our research showed that bacteria are the most abundant microbes, followed by C Effect of water environmental factors on gut viral community compositions in crayfish and crab habitats, respectively. D Effect of sediment environmental factors on gut viral community compositions in crayfish and crab habitats, respectively. All environmental factors were tested in RDA analysis using the envfit() function with 999 permutations, and those with p < 0.1 were shown in the figure, while those with p < 0.05 were considered as significant factors influencing the corresponding microbial community. RCFP: rice-crayfish co-culture model; non-RCFP: all other aquaculture models excluding RCFP; ORP: oxidation-reduction potential; TN: total nitrogen; NH4 + -N: ammonium nitrogen concentration; NO 3 − -N: nitrate-nitrogen concentration; Sed_VM: sediment volume-weight; MC: moisture content; ERY1: erythromycin derivative 1; ERY2: erythromycin derivative 2. Arrows indicate the contribution of environmental factors to the microbial community. Purple symbols indicate microbes in the RCFP group, and red symbols represent microbes in the non-RCFP group viruses, archaea, and eukaryotes. This result is consistent with previous studies on the microbial structure of the ocean environment at the multi-kingdom level [10,38].
Specifically, our results showed that the RCFP model has a distinct set of microbes that play important ecological roles. In a water habitat, Desulfotomaculum is enriched in the RCFP model compared with other aquaculture models. Desulfotomaculum arcticum and Desulfotomaculum geothermicum are capable of degrading acetone and sulfate [39]. Polydnaviridae, a virus enriched in the non-RFCP model of water habitats, could encode virulence genes and infect the cells of the caterpillar host [40], which is harmful to crayfish farming. In sediment habitats, Ferroplasma, an indicator of RCFP, could remove inorganic sulfur compounds by mediating extracellular electron transfer [41]. Verticillium has been detected as an indicator of CFP and reported as a pathogen [42]. In the crayfish gut habitat, Methanopyrus, an obligate chemolithoautotrophic methanogen [43,44], is enriched in the RCFP model and plays an important role in the remineralization and recycling of organic matter [45]. Millerozyma has been reported as a pathogen, and its enrichment in non-RCFP environments is detrimental to crayfish farming [46]. Such results suggest that compared with the RCFP model, other aquaculture models are enriched in many adverse microbes, which increases the opportunistic infection risk for crayfish [42,46]. Moreover, in the RCFP model, many indicator microbes, including bacteria, archaea, viruses, and eukaryotes, play important roles in maintaining the balance of the aquaculture ecology [47,48].

Multi-kingdom interactions in the RCFP model
In this study, the multi-kingdom correlation network of the RCFP system was fairly complex, with a large number of nodes, edges, and degrees and high network density, which means it has a fairly stable microbial community structure, regardless of the habitat examined. Previous researchers showed that complex microbial network interactions reflect a stable microbial community structure [49,50]. Because microbes are the basic driving force of environmental materials and energy flow [8][9][10], stable microbial structures often have better ability to resist external interferences, including human interference [51,52]. Therefore, RCFP ecosystems can be considered to have good ecological sustainability. Such high ecological sustainability may be related to low exogenous pollution inputs (e.g., antibiotics and bait) and strong pollution self-purification ability (e.g., the purification role of rice) [12,14,52]. These results reveal the potential advantages of the RCFP model from the perspective of microbes.

ARG exchange across kingdoms in the RCFP model
Given the extensive use and abuse of antibiotics in recent years, ARGs have gradually become a research hot spot [52][53][54]. In fact, ARGs are the inherent genes of microbes, and these genes have been found in permafrost microbes [53,54]. In addition to resisting the stress of antibiotics, resistance genes often participate in information exchange among microbes and other important functions [55]. We found that ARG-associated HGT events play an important role in multi-kingdom microbial interactions in aquaculture models, i.e., multikingdom interactions existing in both environmental and animal habitats, and many of these interactions are mediated by ARGs. For example, in crayfish gut habitats, the ARG gene MLS|macB is transferred between Methanosphaera and Brevibacillus in the RCFP model. Methanosphaera is a methanogen [56], while Brevibacillus could degrade vinyl acetate in methanogenic reactors [57], which is beneficial for organic matter recycling [45]. Another example is Polynucleobacter in the RCFP model, which has been reported as the host of most ARGs [52,58]. The functional genes of Polynucleobacter exert a positive selection process across different environments [59], which may help other microbes adapt to changing environmental conditions through HGT to acquire new functions [55]. Additionally, Rhodocyclaceae, which was also identified in the RCFP model, can degrade toluene and aromatic compounds under hypoxic [60]. The genes of these microbes participating in HGT events could help microbes rapidly adapt to changing environments through the acquisition of new functions [25,61]. However, ARG-associated HGT events varied across different aquaculture models. For instance, Citrobacter was positively correlated with other kingdoms in the RCFP model, while interactions associated with Citrobacter were much less or negative in non-RCFP models. Citrobacter is an opportunistic pathogen in aquatic animals; in non-RCFP models, larger inputs of antibiotics and pesticides may induce the pathogenicity of Citrobacter, potentially causing high mortality in aquatic animals [62][63][64]. Collectively, these phenomena suggest that multi-kingdom interactions, especially those mediated by ARGs, shape the microbial composition across different aquaculture models.
Microbial antibiotic resistance and ARG pollution are considered among the most important global public threats in the twenty-first century [65]. In this study, we also discovered that non-RCFP aquaculture, particularly in the water habitat, was predicted to have a higher frequency of putative HGT events than that of the RCFP model. To seek high production and economics from aquaculture using fewer limited resources, a greater amount of pesticides were input into the non-RCFP model for treatment and prophylaxis [12,66], which may induce the ARG spearing and contamination, thereby threatening the human and the surrounding environment [66,67]. These antibiotics are directly entered into the water, while the water also acts as inhabiting filtering, thereby accumulating more ARGs [13,14,16,68]. Our previous study also discovered more ARGs in the non-RCFP model in sediment compared with water habitats [52], while in this study, we found a higher frequency of HGT events in water habitat than that of sediment in the non-RCFP models, especially for the ARG-mediated HGT events, which might be attributed to a higher fluidity of water compared with sediment. Therefore, we believe that the lower frequency of ARGs present in the RCFP compared to non-RCFP, especially those transferred across microbes, is a strong indication of the sustainability of this model for several possible reasons. First, in the RCFP model, the presence of crayfish is linked to fewer antibiotic inputs [12,14], and alternating water and high sediment quality could reduce the accumulation and spreading of ARGs [14]. Second, the water and sediment in the model could play as inhabiting filtering, such as ARG filtering, thereby reducing the occurrence of antibiotic resistance selection in the microbial community of crayfish [13,14,16,68]. Finally, the high salinity in the RCFP model could reduce the relative abundance of ARGs [69]. Taken together, these phenomena suggest that the RCFP model is an environment-friendly and ecologically balanced aquaculture model.

Effect of environmental factors on microbial communities
Microbial communities are affected by a variety of environmental factors [26,70]. Previous studies mainly focused on the impact of different environmental factors on bacterial communities [71,72], but an in-depth understanding of the impact of environmental factors on microbial communities at the multi-kingdom level in aquaculture models is still lacking. We found that, at the multi-kingdom level, the microbial communities of non-RCFP models are mainly affected by temperature, pH, and nitrogen level. This result is in line with the results of earlier multi-kingdom studies on the global ocean microbiome [10], which also showed that temperature and pH strongly influence microbial community structures and functions. In addition, previous studies reported that salinity, rather than temperature, could explain a significant portion of the distribution patterns of microbial communities [70]. In the present study, our results indicated that salinity (r 2 = 0.28, p = 0.07) is indeed important but temperature (r 2 = 0.510, p = 0.005) could explain a larger portion of the distribution patterns of the water microbial communities in non-RCFP aquaculture models. This finding may be attributed to higher temperatures generally being able to stimulate microbial growth and metabolism [73], which could increase the stability of the microbial community. Thus, temperature plays an important role in shaping the microbial composition structure in non-RCFP models. Our results indicated that RCFP microbial communities are less affected by environmental factors than non-RCFP microbial communities. In fact, we determined that the bacterial community of the RFCP model in water habitats was significantly influenced by ORP (r 2 = 0.546, p = 0.004) only, mainly because ORP could adjust microbial metabolism via modulating the intracellular redox balance [74]. This result is consistent with a previous study that found that ORP exerts an important effect on the succession of both bacterial and fungal communities in the sludge composting process [74].
We further found that compared with the water microbial community, the crayfish gut microbial community is less affected by environmental factors, especially in the RCFP model. These results indicate that the microenvironment of animal guts is stable and fairly resistant to the effects of the external environment. Crayfish, for example, is able to resist exogenous environmental stress, which is related to the environmental adaptability of crayfish [16,52]. In the RCFP model, the external environment appears to exert minimal effects on the microbial community in crayfish gut, because the rice-crayfish co-culture environment is relatively clean [12,14,52,68] and can indirectly influence the crayfish gut microbial community [16].

Conclusions
This study quantified the sustainability of the RCFP model at the multi-kingdom level on the basis of microbial communities. The sustainability of different aquaculture models has previously been examined but rarely quantified at this level. The multi-kingdom microbial profiles of different culture models were investigated from multiple aspects, including the complexity of the microbial community, network interactions, and the HGTs of functional genes, especially ARGs. Our results clearly illustrated that microbial communities from the RCFP model have unique indicator microbes, such as Shewanella, Ferroplasma, Leishmania, and Siphoviridae. Moreover, the RCFP microbes were densely and positively connected in the network, which suggests that these microbial communities are resilient to environmental stress. Our results further illustrated that HGT events of functional genes, especially ARGs, among bacteria, archaea, and viruses have lower frequencies in the RCFP model than that in other aquaculture models. Finally, environmental factors, such as pH, ORP, temperature, and TN, could substantially shape the microbial communities in the aquaculture environments, although the microbial communities from the RCFP model, especially in the crayfish gut, are less affected by these factors.
The results collectively indicate that the RCFP possesses specific patterns, including distinct microbial community structures, densely and positively connected microbial networks, lower frequency of HGT events, and robustness against environmental factors. The findings quantitatively confirm the sustainability of the RCFP culture model. Quantification of the sustainability of the RCFP on the basis of microbial profiles could provide a deeper understanding of the links between microbial communities in different aquaculture models and the environmental factors influencing these communities. This work represents one of the first studies on the microbial community of freshwater aquaculture wetlands from the multi-kingdom perspective and provides new insights into sustainable aquaculture.