Novel insights into the rhizosphere and seawater microbiome of Zostera marina in diverse mariculture zones

Background Seagrasses offer various ecosystem services and possess high levels of primary productivity. However, the development of mariculture has affected the homeostasis of seagrass meadow ecosystems. Plant-microbiome associations are essential for seagrasses health, but little is known about the role of environmental microbiomes and how they affect seagrass in a mariculture environment. In this study, we investigated the influence of mariculture on the rhizosphere and seawater microbiome surrounding Zostera marina and focused on the bacterial, eukaryotic, and fungal components in the composition, diversity, metabolism, and responses to mariculture-related environmental factors. Results Significant differences in the composition, richness, diversity, and internal relations of the bacterial community between the seawater and rhizosphere sediment surrounding Z. marina were observed, while differences in the eukaryotic and fungal communities were less significant. More complex bacterial and fungal co-occurrence networks were found in the seawater and rhizosphere sediment of the Saccharina japonica (SJ) and sea cucumber (SC) culture zones. The seawater in the SJ zone had higher levels of dissimilatory and assimilatory nitrate reduction, denitrification, and nitrogen fixation processes than the other three zones. The assimilatory sulfate reduction enzymes were higher in the rhizosphere sediments of the SJ zone than in the other three zones. Tetracycline, sulfonamide, and diaminopyrimidine resistance genes were enriched in the mariculture SJ and SC zones. Conclusions Our findings might contribute to a better understanding of the effects of mariculture on the seagrass and the meadow ecosystems and thus revealing their potential operating mechanisms. These insights may serve to raise awareness of the effects of human activities on natural ecosystems, regulation of antibiotic usage, and environmental restoration. Video Abstract Supplementary Information The online version contains supplementary material available at 10.1186/s40168-024-01759-3.


4]
. Recent studies have revealed that plant-microbe associations, which facilitate several metabolic exchanges and biogeochemical transformations thus altering the homeostasis of their physiological and biological functions in seagrass meadow ecosystems, are essential for seagrass health [4][5][6][7].It has been reported that the presence of microbial community is crucial for developing seagrasses, including seed germination, phytohormone production, and the formation of defense mechanisms against pathogens [8,9].
The rhizosphere has a rich diversity of microorganisms, which are beneficial to plants through the actions of pathogenic infection suppressions and help in acquiring nutrients [10].Seagrass rhizosphere sediments and surrounding seawaters are hotspots of nitrogen and sulfur cycles, processes that are governed by specific microorganisms, particularly bacteria [11,12].These processes include nitrogen fixation and nitrification [13], ammonium production [14], sulfide oxidation, and sulfate reduction [4,13].Research has demonstrated that fungi are crucial for the fitness and health of terrestrial plants and act as pathogens or endophytes, but the role of fungi in marine habitats has been poorly studied [15].It is thought fungi potentially play important roles in the cycle of organic matter and the dynamics of food webs in marine ecosystems, such as seagrass meadows [16,17].Recently, studies have investigated the diversity of the seagrass-associated fungal community and emphasized the importance of understanding the variables that influence the dynamics of these communities [18,19].Seagrass meadows commonly co-exist with algae, especially in tropical oceans [20], in which affect the stability and primary productivity of the seagrass community, and ultimately affect their carbon sink [21,22].Previous studies have shown that the red Smithora algae persists with the seagrass Z. marina and influences the seagrass-associated assemblages [23].Seagrass-associated microorganisms, including bacteria, fungi, and eukaryotic algae, have increasingly been recognized as key components of seagrass ecology [8,24].
Maricultures, constructed on coastal lines, are crucial for supplying marine seafood for humans, raising the incomes of coastal communities, and promoting the development of a marine economy.However, it will inevitably impact the original ecological homeostasis [25].China has the largest mariculture area and production in Asia, which accounts for ~71% of the global mariculture production [26].Over the past few decades, mariculture has grown at an annual rate of ~6% in China [27].The rapid development of mariculture has unfortunately resulted in environmental issues, such as marine pollution, eutrophication, antibiotic residues, and food safety concerns [28].The residues of antibiotics enrich the antibiotic resistance genes (ARGs) in the water columns and sediments, change the microbial communities and biogeochemical cycles, and eventually threaten the lives of marine organisms and human health [29,30].Some bacterial taxa, such as Vibrio and Mycobacterium, are strongly linked to the ARGs and contribute to the ARG enrichment [31,32].
The microbiome in seawater and rhizosphere sediment of seagrass may change when exposed to different interfering factors.Sea cucumbers are associated with seagrass meadows in the natural marine environment and form complex interrelationships with the seagrass [33].Seagrasses could protect sea cucumbers from predators, while sea cucumbers, in turn, may act as bioturbators or ecosystem engineers to promote the transformation of nutrients in the rhizospheres of seagrasses [34].However, sea cucumbers mariculture may alter the microbiomes and metabolic processes of natural seagrass meadows [33].Active seagrass restoration is an important tool for facilitating seagrass meadow recovery, as such restoration produces fewer physical stressors and provides more positive intraspecific interactions that may enhance seagrass growth [35].Positive interactions, such as mutualistic effects, are common in natural seagrass meadows [36], but those in the seagrass restoration areas may be different.Some studies reported that kelp cultivation did not impact Z. marina biomass [37], but the environmental microbiome in the seagrass meadow ecosystems will change with any type of interference (such as nutrient input) from kelp mariculture.Currently, studies are primarily focused on the mariculture economy or antibiotic residues in the environment [27,29].However, the influence of different mariculture on seagrass environmental microbiome remains unknown.
The seagrass-associated microbiomes in seawaters and sediments have been proven to be different [11,38].Our research interests in the microbiome of seawaters and rhizosphere sediments surrounding Z. marina were motivated by the effects of mariculture on seagrass and its microenvironment.Although mariculture affects the metabolism of the seagrass, little is known about the roles of microbial communities in this process.To fill this knowledge gap, we collected seawater and rhizosphere sediment samples from natural (N) sea, Saccharina japonica (SJ) culture, sea cucumber (SC) culture, and seagrass restoration (SR) zones, which were surrounding the Z. marina, and performed high-throughput sequencing of 16S and 18S rRNA and internal transcribed spacer (ITS) genes in addition to performing metagenomic sequencing.Our study had several aims: (1) investigate the differences in bacterial, eukaryotic, and fungal communities in the seawaters and rhizosphere sediments of the Z. marina in different mariculture zones; (2) analyze the metabolic pathways and functional genes using metagenomic sequencing; and (3) elucidate the potential influencing mechanisms of mariculture on the environmental microbiome surrounding Z. marina.

Sampling and environmental parameters
Surrounding seawater and rhizosphere sediment samples of Z. marina were collected from natural and mariculture zones (N, SJ, SC and SR) in Yantai and Weihai of Shandong Province, China, in June 2021 (Fig. S1).The seawater samples from NS, SJ, SC and SR zones were named NSW, SJSW, SCSW and SRSW, respectively, and rhizosphere sediment samples of these four zones were named NRS, SJRS, SCRS and SRRS, respectively.At each sampling location, four replicate samples (~50 m apart each) were collected.One liter of seawater was filtered through 0.22-μm polycarbonate membranes (Millipore Corporation, Billerica, MA, USA) for the sequencing of bacterial, eukaryotic, and fungal communities.Ten liters of seawater was filtered through 0.22-μm polycarbonate membranes (Millipore Corporation, Billerica, MA, USA) for metagenomic sequencing.The Z. marina seagrass was shaken to remove the loose soil and scraped using a plastic spoon to collect the corresponding rhizosphere sediment samples.Membranes and rhizosphere sediment samples were frozen in liquid nitrogen immediately, and then stored at -80 °C in the laboratory.
Nutrients in seawater samples, including nitrate (NO 3 ), were measured using a nutrient auto-analyzer (AA3, Seal Analytical Ltd, UK) after filtering with 0.45-μm cellulose acetate membranes [39].The mean (MD) and mid-value (D50) diameters of rhizosphere sediment were measured with a laser particle analyzer (Malvern Instruments, UK) after removing the large particles using a 2-mm sieve.The rhizosphere sediment samples were dried at 60 °C, sieved through 149-μm sieve and ground using a porcelain mortar for the detection of the following environmental parameters.For nutrients detection, 25 mL KCl (2 M) was added to 2.5 g rhizosphere sediment sample and shaken for 1 h.After filtering with 0.45-μm cellulose acetate membranes, NO ), chloride (Cl − ), and bromide (Br − ), ultrapure water was added to the rhizosphere sediment at a ratio of 5:1 (w/w) and shaken for 1 h.The Cl − , Br − , and SO 4 2− concentrations were detected using Ion Chromatograph (Dionex ICS-3000, USA) after filtered with 0.45-μm cellulose acetate membranes.For the detection of total carbon (TC), total organic carbon (TOC), and total organic nitrogen (TON), 3 mL hydrochloric acid (HCl, 1 M) was added to a brown bottle containing 2 g rhizosphere sediment, shaking for 1 h and discarding supernatant.Adding 1 M HCl into bottle until no bubbles generated, discarding the supernatant after shaking for 1 h.Three milliliters of ultrapure water was added into the brown bottle, mixed and dried in a 60 °C oven.The TC, TOC, and TON were detected using vario MICRO cube elemental analyzer (Elementar, Germany) [40].

DNA extraction
Total DNA of seawater samples, including samples for 16S rRNA, 18S rRNA, ITS sequencing, and metagenomic sequencing, was extracted using E.Z.N.A. ® Soil DNA Kit (Omega Bio-Tek, Georgia, USA) following the manufacturer's protocol [41].The DNA of rhizosphere sediment samples was extracted using the FastDNA TM Spin Kit for Soil (MP Biomedicals, California, USA) [42].All the DNA was dissolved in 30 μL DNase-free water to increase DNA concentration and sent to Majorbio Bio-Pharm Technology Co. Ltd. (Shanghai, China) for sequencing.

Bacterial, eukaryotic, and fungal community structure analyses
The bacterial 16S rRNA gene was amplified using 338F and 806R primers (Table S1) [43].The PCR reaction mixture (20 μl) contained 2× Pro Taq (10 μl), forward and reverse primers (4 μM), and template DNA (200 ng).PCR cycling conditions were as follows: 1 cycle of 3 min at 95 °C, 30 cycles of 30 s for denaturation at 95 °C, 30 s for annealing at 53 °C and 45 s for elongation at 72 °C, 1 cycle of 10 min at 72 °C.The amplified PCR products were purified using an AxyPrep DNA Gel Extraction Kit (Axygen Biosciences, CA, USA) and quantified using QuantiFluor ™ -ST (Promega, USA).Purified amplicons were paired-end sequenced (2 × 300 bp) on the Illumina MiSeq platform (Illumina, San Diego, USA) by Majorbio Bio-Pharm Technology Co. Ltd. (Shanghai, China).The eukaryotic 18S rRNA gene was amplified using TAReu-kFWD1F and TAReukREV3R primers (Table S1) [44], and the PCR condition was the same with 16S rRNA gene except for the 55 °C annealing temperature and 35 cycles.The fungal ITS gene was amplified using ITS1F and ITS2R primers (Table S1) [45], and the PCR conditions were the same with 18S rRNA gene amplification.
Raw data were filtered using QIIME2 (version 2022.2) [46] and then trimmed with FASTP to remove lowquality reads (< 100 bp) [47].Paired-end sequences were joined with at least a 10-bp overlap and < 5% mistakes using FLASH [48].Then, the amplicon sequence variants (ASVs) of 16S rRNA, 18S rRNA, and ITS were clustered using QIIME2 with 100% similarity.Barcodes, chimeras, primers, and low-quality reads were removed using the DADA2 plug-in within QIIME2 software [49].The random rarefaction of each sample was performed to obtain an equal sequencing depth (i.e., to obtain the same minimum number of sample sequences).The taxonomic classification of the 16S rRNA gene sequence for all ASVs was analyzed against the Silva 138 16S rRNA database using a classify-sklearn classifier with a confidence threshold of 70%.
Alpha diversity indices, including Good's coverage, Chao 1, and Shannon indices, were calculated using the MOTHUR to estimate the richness (Chao 1) and diversity (Shannon) of the bacterial, eukaryotic, and fungal communities [41].The linear discriminant analysis effect size (LEfSe) was performed to assess the taxa enrichment in these four zones with linear discriminant analysis (LDA) values > 4, and taxa that were significantly higher in all four zones were set as the enriched taxa (biomarkers) [50].The nonmetric multidimensional scaling (NMDS) analysis was used to evaluate ASVs differences.Canonical correspondence analysis (CCA) was used to analyze correlations between bacterial communities and environmental factors.Both the NMDS and CCA analyses were performed on R software (v3.3.1) using "vegan" package.The variation partition analysis (VPA) was performed on R software (v3.3.1) using "vegan" package to investigate the relative contribution of different environmental parameter groups to bacterial, eukaryotic, and fungal communities.The correlations between different environmental parameters as well as the bacterial, eukaryotic, and fungal taxa were calculated using Mantel test on R software (v3.3.1) using "vegan" and "ggcor" packages [51].The co-occurrence network was analyzed on R software (v3.3.1) using "Hmisc" package and plotted using Gephi software (v0.9.2).
The initial paired-end reads were co-assembled using IDBA-UD [55].Subsequently, these reads were mapped to the final assembled contigs with Burrows-Wheeler Alignment tool [56].Contigs were grouped into metagenome assembled genomes (MAGs) based on empirical probabilistic distances of genome abundance and tetranucleotide frequency using program MetaBAT2 [57].Completeness and contamination of MAGs were assessed using CheckM (v1.2.2) [58], and MAGs with a completeness above 80% and contamination lower than 10% were considered for further analysis.The assignment of MAGs taxonomy and the subsequently phylogenomic tree construction were conducted using GTDB-Tk tool (v2.3.2) based on GTDB (https:// gtdb.ecoge nomic.org) database [59].Genome annotation of the MAGs was carried out using the RAST server (https:// rast.nmpdr.org).
Representative sequences of non-redundant genes were annotated against the NR database (v20200604) in NCBI to obtained the taxonomy assignment using Diamond (v0.8.35) with an e-value cutoff of 10 −5 (blastp) [60].The Kyoto Encyclopedia of Genes and Genomes (KEGG) annotation was conducted against the KEGG database (version 94.2) using the Diamond software (v0.8.35) with an e-value cutoff of 10 −5 [61].ARGs annotation was conducted using Diamond (0.8.35) against the Comprehensive Antibiotic Resistance Database (CARD) with an e-value cutoff of 10 −5 [62].The relative contributions of each taxon to KEGG and CARD functions were calculated as described previously [63].The LEfSe was performed to assess the enriched functions in these four zones with LDA values > 2. Functions that were significantly higher in all of these four zones were set as enriched functions (biomarkers).

Statistical analyses
Unless otherwise stated, the difference in environmental parameters, the bacterial, eukaryotic, and fungal community (composition and α diversity), and the metagenomic analysis (KEGG and CARD analyses) between these four zones were tested using the Wilcoxon rank-sum test.The correlations between different environmental parameters and taxa of bacteria, eukaryota, and fungi were analyzed using the Spearman correlation test.All statistical analyses were performed using SPSS version 24.0 (SPSS Inc., Chicago, IL, USA).P < 0.05 and P < 0.01 were considered statistically significant.

Environmental parameters
The environmental parameters of the seawater and rhizosphere sediment samples surrounding Z. marina are shown in Tables S3 and S4.Generally, the seawater samples of these four zones had a higher level of SiO , and MD) samples.

Bacterial, eukaryotic, and fungal community compositions
The sequencing of 16S and 18S rRNA and ITS genes was performed on seawater and rhizosphere sediment samples from four diverse zones surrounding Z. marina.For the 16S rRNA sequencing, we obtained 966,249 sequences after "denoising" the system using DADA2, ranging from 19,831 to 46,237.The average length of these samples was 417 bp.A total of 32,182 ASVs were assigned after random rarefaction.Good's coverage values of all samples ranged from 98.39 to 100.00%, indicating a good representativeness of the bacterial communities in these samples (Table S5).Generally, a higher ASV numbers were found in the rhizosphere sediment samples versus the seawater samples (P = 0.021 for all the four zones), indicating a higher bacterial community richness in the rhizosphere sediments.The alpha indices also showed that the richness (Chao 1 index) and diversity (Shannon index) of the bacterial communities in the rhizosphere sediment of these four zones were higher than in seawater samples (P = 0.030 for all the four zones, Fig. S2A and D).Besides, the richness (Chao 1 index) and diversity (Shannon index) of the bacterial communities in the SJSW samples were higher than in the NSW and SRSW samples (P = 0.030 for both, Fig. S2A and D).
For the 18S rRNA sequencing, we obtained 1,971,641 sequences after applying a "denoising" protocol using DADA2 that ranges from 32,051 to 84,069.The average length of these samples was 378 bp.After random rarefaction, a total of 19,123 ASVs were assigned.The Good's coverage values of all samples ranged from 99.50 to 100.00%, indicating a good representativeness of the eukaryotic communities (Table S5).The ASV numbers in the SRRS samples were higher than in the SRSW samples (P = 0.021), but no significant differences in the ASV numbers between the seawater and rhizosphere sediment samples in the NS, SJ, and SC zones (P = 0.248, 0.149 and 0.773, respectively) were detected.The alpha indices also indicated that only the eukaryotic community in the SRRS samples had a higher richness (Chao 1) and diversity (Shannon) than in the SRSW samples (Fig. S2B and  E).Furthermore, the eukaryotic community of SRSW samples had the minimum richness and diversity among the seawater samples of all these four zones (Fig. S2B and  E), but no significant differences were found between the rhizosphere sediment samples (P > 0.05).
For ITS gene sequencing, we obtained 1,531,903 sequences after applying a "denoise" protocol using DADA2, ranging from 26,120 to 71,292.The average length of these samples was 242 bp.After random rarefaction, a total of 11,642 ASVs were assigned.Good's coverage values of all samples ranged from 99.86 to 100.00%, indicating good representativeness of the fungal communities (Table S5).The ASV numbers in NRS and SRRS samples were higher than those in corresponding seawater samples (P = 0.021 for both), but no such difference was found in the SJ and SC zones (P = 0.248 and 0.083, respectively).The richness (Chao 1) of the fungal community was also higher in the NRS and SRRS samples than in seawater samples, but not in the SJ and SC zones (Fig. S2C).As for the diversity (Shannon), only the SRRS samples were lower than the SJRS and SCRS samples (P = 0.030 for both).
To seek the diversity and potential functions of uncultured novel bacteria, we constructed 109 high-quality MAGs with a contaminant threshold of 10% and a completeness threshold of 80% (Table S6).These MAGs included 24 Alphaproteobacteria (mainly the Rhodobacteraceae), 12 Gammaproteobacteria, 36 Bacteroidia (all the Flavobacteriaceae), 21 Verrucomicrobiae (potential probiotics), and some other bacteria such as Desulfobacteria.Interestingly, ~87% MAGs were obtained from the seawater samples, which was ~6.7-fold higher than the rhizosphere sediment samples (Fig. S5), although these samples had similar clean data (16.26 to 18.26 Gb).The results of the phylogenomic tree analysis suggested that the MAGs from NSW samples were mostly Alphaproteobacteria and Bacteroidia, while MAGs from SRSW were mostly Alphaproteobacteria. MAGs in the SJSW and SCSW had higher diversity than the former two sample groups, and most of them belonged to Verrucomicrobiae and Bacteroidia (Fig. S5).

Interactions and diversity of bacterial, eukaryotic, and fungal communities
The co-occurrence network analysis was performed using the top 100 ASVs to study the interactions of bacterial, eukaryotic, and fungal communities in these four zones (Fig. 3, S6, and S7).The nodes in the co-occurrence network analysis were similar and ranged from 92 to 97, 90 to 93, and 94 to 97 in bacterial, eukaryotic, and fungal communities, respectively (Table 1).The edges in the bacterial community in the SJSW and SRSW samples were higher when compared with the NSW samples (1.9-and 1.7-fold, respectively), and the proportions of positive and negative correlations were equal (Fig. 3A, B, D).The SCRS samples also had a higher edge number than the NRS samples (1.5-fold), which was mainly due to positive correlations (Fig. 3G).As for the eukaryotic community, no significant difference in the number of edges between these four zones was found, but the proportions of positive correlations in SJSW (89.69%) and SRSW (72.46%) samples were higher than in NSW (52.51%) samples (Fig. S6B and D).When compared with the eukaryotic community in the NRS samples, only the SCRS samples had a higher edge number, and the proportions of positive correlations in the SJRS, SCRS, and SRRS samples (65.25%, 69.13%, and 75.72%, respectively) were lower than in the NRS samples (88.15%,Fig. S6E-H).The edge numbers in the fungal communities of the SJSW and SRSW samples were higher than in the NSW samples, and the SJSW, SCSW, and SRSW samples showed a higher proportion of positive correlations (88.86%, 83.24%, and 97.35%, respectively) than the NSW samples (67.83%,Fig. S7A-D).The edge numbers of fungal communities in the rhizosphere sediment were similar except for the SCRS samples (60.00%), which was lower than the other three zones (Fig. S7E-H).
The NMDS and CCA analyses were performed to investigate the dissimilarity of bacterial, eukaryotic, and fungal communities (Figs.S8 and S9).The results showed significant differences in the bacterial, eukaryotic, and fungal communities between seawater and rhizosphere sediment samples, particularly in the bacterial community (red and blue spots, Fig. S8).The Anosim analysis also confirmed significant differences in these three communities between the seawater and rhizosphere sediment samples (R = 0.999, 0.962 and 0.925 for bacterial, eukaryotic, and fungal communities, respectively, P = 0.001 for all of them).The bacterial, eukaryotic, and fungal communities in the NSW, SJSW, and SCSW were similar, but those in the SRSW samples differed from the other three zones.However, no such significant trend was observed in the rhizosphere sediment samples for the bacterial, eukaryotic, and fungal communities.

Correlations between environmental parameters and communities of bacteria, eukaryotes, and fungi
When combined with the environmental parameters, the CCA revealed that the correlations between environmental factors and communities of bacteria, eukaryotes, and fungi were similar both in seawater and rhizosphere sediment samples (Fig. S9).The NO 3

−
, NO 2 − , and SiO 3 2− concentrations were the main environmental factors that affected the bacterial, eukaryotic, and fungal communities in the SJSW samples, and NH 4 + and PO 4 3− concentrations affected the bacterial, eukaryotic, and fungal communities in the NSW and SRSW samples, respectively (Fig. S9A-C).As for the rhizosphere sediment samples, NO 3 − , NO 2 − , NH 4 + , Br − , and MD were the important factors for the SJRS samples, and TOC concentration may have affected the bacterial, eukaryotic, and fungal communities in the SCRS samples (Fig. S9D-F).
To evaluate the influence of different environmental parameter groups on those communities, N (including NO 3 .PO 4 3− contributed more to the eukaryotic community of seawater samples than to the bacterial and fungal communities (Fig. S10A-C).As for the rhizosphere sediment samples, the C group was the main contributor (29.25%) for the bacterial community, followed by the D group a with 16.73% contribution (Fig. S10D).All of the C, N, and D groups contributed less to the eukaryotic communities in the rhizosphere sediment samples (8.81%, 2.55%, and 3.08%, respectively, Fig. S10E).The N group was the main contributor to fungi and accounted for 14.48% of the total community (Fig. S10F).
The Mantel test was performed to analyze the correlations between different environmental factors and elucidate their interactions with distinct microorganism groups of the bacterial, eukaryotic, and fungal communities in the seawater and rhizosphere sediment samples (Fig. 4).Significantly positive correlations were observed in the seawater samples between NO 3 − , NO 2 − , and SiO 3 2− (Fig. 4A-C, P = 0.001 for all, Pearson's correlation test).For the rhizosphere sediment samples, significantly positive correlations were observed between particle diameter (D50 and MD) and some other environmental parameters, including NH 2− and Br − (P = 0.001), D50 and MD (P = 0.001), TON and TC (P = 0.009), TON and TOC (P = 0.027), and TC and TOC (P = 0.001) in the rhizosphere sediment samples (Fig. 4D-F, Pearson's correlation test).A significantly negative correlation was only found in the rhizosphere sediment samples between D50 and TC (P = 0.026, Pearson correlation test).
In seawater and rhizosphere sediment samples, the top 50 ASVs of the eukaryotic community were grouped into Chlorophyta (10 ASVs), Arthropoda (7 ASVs), Dinophyceae (Alveolata clade, 5 ASVs), and others (28 ASVs) clusters at the phylum level.The number of significant correlations between environmental parameters and eukaryotic clusters was less than those in the bacterial community both in the seawater and rhizosphere sediment samples (Fig. 4B, E).Significant correlations were observed between Chlorophyta and the concentrations of PO 4 3− (P = 0.027) and NH 4 + (P = 0.035) in addition to between Arthropoda and PO 4 3− (P = 0.011) concentrations in seawater samples (Fig. 4B).The Chlorophyta in the rhizosphere sediment samples significantly correlated with NO 2 − (P = 0.007), D50 (P = 0.002), and MD (P = 0.001), while Dinophyceae only significantly correlated with TON (P = 0.010, Fig. 4E).No significant correlations were found between Dinophyceae and environmental factors in seawater samples.In addition, no significant correlations were found between Arthropoda and environmental factors in the rhizosphere sediment samples.
In seawater and rhizosphere sediment samples, the top 50 ASVs of the fungal community were grouped into Ascomycota (13 ASVs), Basidiomycota (8 ASVs), Chytridiomycota (three ASVs), and others (26 ASVs) clusters at the phylum level.The number of significant correlations between environmental parameters and fungal clusters was less than in the bacterial community both in the seawater and rhizosphere sediment samples, which was similar to the eukaryotic community (Fig. 4C, F).Ascomycota and Basidiomycota in seawater samples only correlated significantly with PO 4 3− concentration (P = 0.026 and 0.049), and Chytridiomycota in seawater samples correlated with NH 4 + (P = 0.001), NO 2 − (P = 0.010), and SiO 3 2− (P = 0.046) concentrations (Fig. 4C).As for the rhizosphere sediment samples, both Ascomycota and Chytridiomycota were significantly correlated with TON (P = 0.001 and 0.001), TC (P = 0.004 and 0.001), and TOC (P = 0.010 and 0.002) concentrations.Besides, Ascomycota also correlated significantly with NH 4 + (P = 0.010), and Chytridiomycota also significantly correlated with MD (P = 0.049).Basidiomycota in rhizosphere sediment samples correlated significantly with particle diameters (D50 and MD, P = 0.040 and 0.013) as shown in Fig. 4F.

Metagenome and metabolism analysis
The corresponding metagenomic sequencing (32 samples) was performed on seawater and rhizosphere sediment samples.The results of KEGG functions showed that the global and overview maps of metabolism (such as metabolic pathways, biosynthesis of secondary metabolites, biosynthesis of amino acids), nucleotide metabolism (both purine and pyrimidine metabolism), and genetic information processing (such as aminoacyl-tRNA biosynthesis, ribosome) were higher in seawater while the global and overview maps of metabolism (such as microbial metabolism in diverse environments, carbon metabolism), signal transduction (two-component system), cellular community of cellular processes (quorum sensing), carbohydrate metabolism (such as glycolysis/ gluconeogenesis, citrate cycle, pyruvate metabolism), and energy metabolism (such as oxidative phosphorylation, methane metabolism) were higher in rhizosphere sediments (Fig. S11C, P = 0.001).
The contribution of various species to the KEGG functions in seawater and rhizosphere sediment samples was evaluated to investigate the different roles of bacteria, eukaryota, archaea, and viruses in terms of the KEGG functions in seawater and rhizosphere sediment samples (Fig. 5).The main contributor to biosynthesis of secondary metabolites, microbial metabolism in diverse environments, biosynthesis of amino acids, and carbon metabolism was bacteria, followed by archaea and eukaryota.However, the contribution of viruses to those functions was lower especially in the rhizosphere sediment samples (Fig. 5).As for the functions of ABC transporters, two-component system, quorum sensing, and pyruvate metabolism, bacteria were still the main contributors.Interestingly, viruses contribute more to purine and pyrimidine metabolism, homologous recombination, mismatch repair, DNA replication, RNA degradation, and one carbon pool by folate than the bacterial, eukaryotic, and fungal communities (Fig. S12).Archaea were found to play essential roles in the rhizosphere sediment samples with higher contributions than in seawater samples for most of these KEGG functions (Fig. S12).However, the contributions of archaea to purine and pyrimidine metabolism, mismatch repair, DNA replication, one carbon pool by folate, and nucleotide excision repair were opposite and were higher in the seawater samples than in rhizosphere sediment samples (Fig. 5  and S12).The contribution of archaea to carbohydrate metabolism in the SRSW samples was lower than those in the other seawater samples and rhizosphere sediment samples, such as pyruvate, butanoate, propanoate, citrate cycle, pentose phosphate pathway, glyoxylate, and dicarboxylate metabolism.

Nitrogen and sulfur metabolism
Considering that the nitrogen metabolism pathway was significantly enriched in the SJRS samples, we analyzed the dissimilatory reduction, assimilatory nitrate reduction, denitrification, nitrogen fixation, nitrification, and anammox pathways in the seawater and rhizosphere sediment samples of those four zones (Fig. 7).Generally, the levels of most nitrogen metabolic proteins, especially those related to the dissimilatory nitrate reduction, denitrification, and nitrogen fixation pathways, were higher in the rhizosphere sediment samples than in seawater samples (Fig. S13).When compared with the NSW, SCSW, and SRSW samples, the SJSW samples had a higher proportion of dissimilatory nitrate reduction and assimilatory nitrate reduction metabolism.The abundances of the peptidoglycan-associated protein NapA (0.87% ± 0.15%), nitrate reductase NirBD (5.35% ± 0.98%), and the cytochrome c nitrite reductase NrfA (0.23% ± 0.07%) proteins, which function in the dissimilatory reduction of nitrate to ammonia, were significantly higher in SJSW samples than in the other three seawater samples (Fig. 7A).The abundances of NarB (0.53% ± 0.14%) and NirA (0.38% ± 0.05%) were significantly higher in the SJSW samples than in the other three seawater samples (not significant between SJSW and SRSW for NirA), indicating higher assimilatory reduction of nitrate to ammonia (Fig. 7A).The denitrification pathway from nitrite to nitric oxide (NirK, 0.75% ± 0.17%), nitrous oxide (NorBC, 0.61% ± 0.15%), and nitrogen (NosZ, 0.49% ± 0.14%) were also higher in the SJSW samples.The nitrogen fixation (NirDKH and AnfG, 0.23% ± 0.06%) in the SJSW samples was only significantly higher than in the NSW samples.The first step of nitrification (AmoABC) Fig. 7 The nitrogen metabolism in seawater (A) and rhizosphere sediment (B) samples from the N, SJ, SC, and SR zones.Wilcoxon rank-sum test is used to analyze the difference in nitrogen metabolic genes in the seawater and rhizosphere sediment samples between the N, SJ, SC, and SR zones.P < 0.05 is marked with "*" in the SJSW (0.03% ± 0.02%), SCSW (0.02% ± 0.01%), and SRSW (0.02% ± 0.03%) was higher than in the NSW samples (0.01% ± 0.01%), but no significant difference was observed in the subsequent processes.No significant enrichment of the anammox pathway was observed in the seawater samples.For the rhizosphere sediment samples, no consistent trend was found in any of the six nitrogen metabolism pathways (Fig. 7B).The dissimilatory reduction of nitrate to nitrite was lower in the SJRS (NarGHI, 4.25% ± 0.25%) and SCRS (NapA, 3.66% ± 0.41%) samples, whereas the assimilatory reduction of nitrate to nitrite was higher in the SJRS (NarB, 0.30% ± 0.12%) and SCRS (NR, 0.05% ± 0.03%) samples.Interestingly, the nitrogen fixation pathway in the SJRS was lower than the other three rhizosphere sediment samples, which was opposite to the seawater samples (Fig. 7B).

Discussion
The interactions between seagrass and microbial communities are essential for seagrass adaptation and the homeostasis of seagrass meadows through metabolic exchanges and biogeochemical transformations [4,5].Microorganisms in seawater and rhizosphere sediment benefit the seagrass regarding nitrogen and sulfur cycling [4,64], which may be affected by environmental factors.We investigated the impacts of mariculture on the environmental microbiome surrounding Z. marina in seawater and rhizosphere sediments, focusing on community composition, metabolic pathways, and functional genes.

Rhizosphere sediment has higher levels of bacterial diversity, internal relations, and nitrogen metabolic pathways than seawater
The most striking result was the significant difference in community diversity and richness between seawater and rhizosphere sediment samples in all four zones.Especially for the bacterial community, the rhizosphere sediment samples possessed higher diversity (Shannon index) and richness (OTU number and Chao 1 index) as shown in Fig. S2.These results agree with Shang et al., who found a higher Shannon index for the bacterial community in sediment rather than in seawater but not for the fungal community [65].This finding may have been because the rhizosphere sediments contain more nutrients than the corresponding seawaters and could provide adequate substrates for seagrass growth [66].The NMDS results showed that the bacterial, eukaryotic, and fungal communities in seawaters were separated from the rhizosphere sediments, revealing the distinct difference between these two biotopes (Fig. S8).Previous studies addressing the microbial community also revealed the significant difference between seawaters and sediments in the coastal and ocean environments [67,68].The Mantel test results revealed that the top three abundant families/phyla of bacterial, eukaryotic, and fungal communities were significantly correlated with diameters in the rhizosphere sediment samples (Fig. 4), which may be due to the effect of nutrient access, transport, and respiration in this habitat [69,70].Zai et al. found that root diameter differentially affected rhizospheric bacterial communities and implied a close relationship between the bacterial microbiome, root functions, and soil properties, and this finding confirmed our results that diameters of the rhizosphere were correlated with bacterial, eukaryotic, and fungal communities [71].PO 4 3− concentration may be the key factor influencing these top three abundant families/phyla, which has been proven in previous studies [72][73][74][75][76][77].Besides, in seawater samples, the NH 4 + concentration was another important factor that affected the top three abundant families of bacterial communities, and the majority of these taxa were involved in nitrogen metabolism [78,79].
No significant differences in the phylum of Proteobacteria in addition to the eukaryotic taxa were found in our study, a result that was different from Shang et al. who found a significantly higher level of Proteobacteria in rhizosphere sediments [65].These differences may result from the various environmental conditions in which our two seagrass meadows (SJ and SC zones) were influenced by mariculture.The region-specific characteristics of eukaryotic and fungal communities in both seawater and rhizosphere sediment may reveal their susceptibility to environmental impacts in seagrass meadows.For example, Trebouxiophyceae in a previous study was related to clay and silt [80], but no such enrichment of this taxon was found in the zones with higher clay and silt level (SRRS samples, 11.47% and 65.46%, respectively).Chlorophyta (mainly the Trebouxiophyceae and Ulvophyceae) were found to be related to PO 4  3− and NH 4 + concentrations (Fig. 4B, E).This finding may suggest their demands for the phosphate and ammonium nutrients in these natural and mariculture environments, which had been reported to influence the health of Trebouxiophyceae and Ulvophyceae, such as the fatty acid proportions [75,81], the photosynthesis rates [82], and the morphology of algae [83,84].SRRSW samples may be the main contributors to these correlations considering the highest phosphate concentration.Significant correlations were observed between Arthropoda (mainly Hexanauplia, the crustaceans) and PO 4 3− concentrations.Calcium phosphate was crucial for the formation of mineralized chitinous exoskeleton [85], and its lowest concentration in the SJSW may suggest the limitation for the Hexanauplia growth in the SJSW samples (Fig. 1B).
The significantly enhanced KEGG pathways in seawaters were mainly the global and overview maps of metabolism, but the rhizosphere sediments possessed diverse enhanced pathways such as carbon metabolism, two-component system, quorum sensing, and oxidative phosphorylation (Fig. S11C).Bacteria were found to possess the highest contributions to most of the KEGG metabolic pathways, especially the secondary metabolites, microbial metabolism in diverse environments, biosynthesis of amino acids, and carbon metabolism, while the contributions of viruses were the lowest especially in rhizosphere sediments.It is reported that seagrassrelated bacteria could increase the nitrogen availability for uptake by seagrass, mineralizing amino acids, thus enhancing the growth and productivity of seagrass meadows [86].These higher bacterial contributors may suggest their essential roles in maintaining healthy seagrass growth in our study.Seagrass rhizomes and roots could release significant amounts of dissolved organic carbon and stimulate microbial activity, including sulfate reduction in the rhizosphere in turn [24].However, viruses played important roles in nucleotide metabolism, replication and repair, which was even higher than bacteria (Fig. S12).Some studies have revealed that viruses could cause alterations in the nucleotide metabolism pathways and DNA repair of hosts [87][88][89], which confirmed the genetic information related role of viral components in the environment surrounding seagrass.These findings suggest that bacteria may be the vital component in the microbiome of the seagrass surrounding environment, and viruses may affect genetic information processing.
The nitrogen metabolic proteins, especially those in the dissimilatory nitrate reduction, denitrification, and nitrogen fixation pathways, were higher in rhizosphere sediment samples than in seawater samples (Fig. S13), indicating that the rhizosphere was a hotspot of nitrogen metabolism.These findings are consistent with previous knowledge that higher fixation levels [90,91], denitrification [92,93], and dissimilatory nitrate reduction [94] were found in the rhizosphere or bulk sediment of seagrass meadows than in the overlying seawater.Denitrification could reduce NO 3 − and NO 2 − to N 2 thus removing reactive nitrogen from a biological system, and dissimilatory nitrate reduction reduces NO 3 − to NH 4 + and contributes to the nitrogen retention of seagrass [95].
The competition between denitrification and dissimilatory nitrate reduction could regulate nitrogen elimination/retention dynamics in the seagrass holobiont in our study.This finding was confirmed by Aoki et al. who revealed the importance of denitrification and dissimilatory nitrate reduction to ammonium formation in the seagrass root zone using a push-pull incubation method [94].Assimilatory sulfate reduction and DMSP degradation were higher in seawaters, while the dissimilatory sulfate reduction and oxidation were higher in the rhizosphere sediments (Fig. S14).It is known that an increase in sulfite concentration could harm seagrass health and result in dieback [96].The variation in the sulfate reduction in the seawater and rhizosphere sediment of the Z. marina seagrass would impact its homeostasis.

S. Japonica culture enhances bacterial internal relations, nitrogen metabolism, and sulfate reduction levels
The richness and diversity of the bacterial community in SJSW samples were higher than in the NSW and SRSW samples, indicating that the culture of S. japonica may cause an increase in the complexity of the bacterial community in the seawater environment (Figs.S2A and 2D).This finding was confirmed by the bacterial co-occurrence network that the SJSW samples had more edges and positive correlations (Fig. 3).Studies have shown that kelp cultures could enhance coastal biogeochemical cycles by maintaining bacterioplankton richness [97], and dissolved organic carbon released by S. japonica may be the key matter that enhances the stability of the bacterial community [37,98].Psychromonas was enriched in the SJSW samples, which had previously been found in the epiphytic bacteria of S. japonica [99].It is an important primary degrader of both proteins and lipids and could produce alginate lyases [100,101].Approximately 50% of brown seaweed cell walls and intracellular material are alginate; thus, the enriched Psychromonas in the SJSW samples may indicate the degradation demands of S. Japonica.Members of Flavobacteriaceae and Rhodobacteraceae were enriched in SJRS samples, which were mainly detected in the diseased S. japonica [102].This enrichment may be related to the degradation of disease kelp and subsequently decomposition in the SJ zone.Desulfosarcinaceae (Desulfobacteraceae) and Desulfobulbaceae are known as the dominant sulfate-reducing bacteria [4], and these enriched families revealed a potentially higher level of sulfate-reducing process in the SJRS samples.Besides, Liu et al. found an increasing sulfate level in kelp blanching water [103], and the high SO 4 2− (9.11 ± 4.49 mmol/kg) concentrations in SJRS samples may be due to dead kelp degradation.Thus, these sulfate-reducing bacteria cause an increase and then help reduce SO 4 2− .Nitrogen metabolism was significantly enriched in the SJSW samples based on the LEfSe results (Fig. 6A).After an in-depth analysis, we found higher dissimilatory and assimilatory nitrate reduction processes in the SJSW samples when compared with the other three kinds of seawaters surrounding seagrass.Nitrogen fixation and the first step of nitrification were also higher in the SJSW samples when compared with the NSW samples (Fig. 7A).These findings revealed that the nitrogen metabolism in the seawater surrounding seagrass in the SJ zone was higher than in the other three zones.The cultivation of S. japonica needs a sufficient nitrate supply which limited their productivity [104,105], and sufficient nitrate enhances the nitrogen metabolism in SJ zone.Wang et al. found that nitrate reductase activity and nitrate absorptivity were found to be fairly high in the S. japonica sporophyte [106], which support the higher nitrogen metabolism in the SJ zone.No consistently significant differences were found for the nitrogen metabolism between SJRS and the other three kinds of rhizosphere sediment samples, indicating that the impact of S. japonica culture on seagrass rhizosphere was less than that of the seawater.Sandy marine sediments could provide a significant source of nitrogen to the water column in which it was difficult to preserve nutrients [107]; thus, the seawater may have higher nitrogen metabolism rates than in the sediment environment considering the highest sand proportions in SJRS samples (71.87-95.82%).Furthermore, the highest NO 2 − and NO 3 − concentrations were observed in the SJSW and SJRS samples, which also confirmed the presence of active nitrogen metabolism in the SJ zone (Tables S3 and S4).The CCA results suggested that the bacterial, eukaryotic, and fungal communities in the SJSW and SJRS samples were both influenced by NO 3 − and NO 2 − concentrations, especially the bacterial and fungal communities in seawater.This finding further confirmed the nitrogen metabolism results above (Figs.S9 and S10).Many of the metabolic pathways for S. japonica substrates, including the metabolism of amino acids, lipids, and xenobiotics, were enriched in the SJRS samples, and these enrichments may due to the demand of the S. japonica residues degradation process [108].The abundance of dissimilatory nitrate reduction proteins was lower in the SCSW and SCRS samples compared to the NSW and NRS samples.On the contrary, the assimilatory nitrate reduction proteins were higher in the SCSW and SCRS samples than in the NSW and NRS samples.Zheng et al. found nitrate and ammonium uptake, which helps remove nitrogen from the aquaculture systems, in the sediment of the sea cucumber culture pond [109].
Nitrogen in the sea cucumber culture is enriched due to animal feces, extra feed, and fertilizers [110].Thus, the excess nitrate in the sea cucumber culture may be removed from the seagrass meadows via assimilatory nitrate reduction in the SC zone.
The assimilatory nitrate reduction from PAPS to sulfide was generally higher in the SJSW, SJRS, and SCRS samples, which indicates that the S. japonica and sea cucumber culture may have brought sulfide pressure to the seagrass.Furthermore, Sulfurovum could oxidize sulfide and thiosulfate to sulfate [111,112], and the enrichment of this genus in the SCRS samples may have helped mitigate the sulfide pressure of seagrass in the SC zone.Choi et al. revealed that sulfur metabolism was a key driver of refractory dissolved organic carbon production during the biodegradation of S. japonica [113], and the higher level of sulfate reduction and in the SJ zone may have been due to the requirements of S. japonica degradation, which may bring sulfite pressure for seagrass in the meantime.The level of sulfite pressure may increase in the S. japonica and sea cucumber culture (SJ and SC) zones as revealed by an increase in SO 4 2− concentration, assimilatory sulfate reduction, and sulfate-reducing related bacteria.Thus, artificial breeding scales should be under control and supervision to preserve seagrass meadows from degradation risk.
Considering the enriched microorganisms in all these natural and maricultural zones, it was found that most of them were related to sulfate and nitrogen metabolism, which were shown to assist seagrass in nitrite removal and nitrogen fixation [114].Members of Sulfitobacter (enriched in NSW) are specializing in sulfite oxidation, and these marine bacteria are abundant in the seagrass tissues and could protect algae such as Emiliania huxleyi from pathogen infection [115].Sulfurovum sp.(enriched in SCRS) could use oxygen and nitrate as electron acceptors and gain energy by oxidizing reduced sulfur compounds through the sulfur-oxidizing (Sox) pathway [116]; thus, this enrichment may be a response to sea cucumber culture for maintaining the seagrass health in our study.MAGs of Desulfobacterales and Desulfuromonadales, mainly the sulfur-reducing bacteria, were obtained from our metagenome data, and these bacteria were commonly associated with anoxic, sulfide-rich seagrass rhizosphere [117].The annotation of these MAGs also revealed the sulfur-reducing ability in these bacteria (DsrC, DsrD), and maybe potential seagrass probiotics (Table S6).These associations between seagrass and sulfide-oxidizing bacteria could alleviate toxic sulfate accumulation, produce higher biomass and form more complex rhizome structures [118].Besides, some other enriched bacteria, such as Psychromonas (enriched in SJSW) [119] and Rhodobacteraceae (enriched in SJRS) [120], are capable of nitrite reduction, nitrate assimilation and/or nitrogen fixation, which may aid Z. marina in maintaining health under kelp culture pressure.MAGs of Rhodobacteraceae were also obtained from the metagenomic data, and these bacteria had potential nitrogen fixation ability with NifU annotated.Further, MAGs in the family Akkermansiaceae contained genes for dissimilatory nitrate reduction, which was consistent with Weigel's work [121], and these bacteria might potentially promote the growth of Z. marina considering its nitrate reduction ability.Although there was no study revealed the effects of Lentimonas and Polaribacter to the growth of seagrass, genes involving in fucoidan degradation [122] (alpha-l-fucosidase) were observed in MAGs of these two genera.Considering that fucoidan is mainly derived from kelp, bacteria of Lentimonas and Polaribacter could provide nutrients for seagrass growth through degrading dead kelp, being seagrass probiotics.
Different amounts of antibiotics have been used in the mariculture and are released into the marine and terrestrial environments [123].Although the predominant ARGs in these four zones were similar, the ARGs against bicyclomycin, pleuromutilin, tetracycline, sulfonamide, and diaminopyrimidine were significantly enriched in the mariculture zones (SJ and SC zones, Fig. S16).Bicyclomycin, sulfonamides, and tetracyclines have been widely used in mariculture to treat bacterial and protozoan infections [124,125].ARG enrichment indicates the impact of antibiotic addition to the seagrass surrounding environment, which may enhance their ARG levels compared to the natural sea.Effective supervision measures should be taken to reduce the addition of antibiotics in mariculture to protect seagrass meadow ecosystems.

Conclusions
This study analyzed the impacts of mariculture on the environmental microbiome surrounding the seagrass Z. marina, focusing on the bacterial, eukaryotic, and fungal components in the composition, diversity, metabolism, and responses to the mariculture environment.We found significant differences in the composition, richness, diversity, and internal relations of the bacterial community between the seawater and rhizosphere sediments surrounding Z. marina, while the eukaryotic and fungal communities were less significant.More complex bacterial and fungal co-occurrence networks were found in the seawater and rhizosphere sediments of the SJ and SC zones.The SJSW samples had higher levels of dissimilatory and assimilatory nitrate reduction, denitrification, and nitrogen fixation processes than the other three zones.Proteins related to assimilatory sulfate reduction were higher in the rhizosphere sediments of the SJ zone.Tetracycline, sulfonamide, and diaminopyrimidine resistance genes were enriched in the mariculture zones.Our study showed the impact of mariculture on the microbiome of seagrass Z. marina and could enhance the understanding of the mariculture impacts on seagrass meadow ecosystems.These novel insights may help evoke the emphasis on the human productive activity impacts on natural ecosystems.

Fig. 1
Fig. 1 The bacterial (A), eukaryotic (B), and fungal (C) community composition in seawater and rhizosphere sediment samples from the N, SJ, SC, and SR zones at class level.The top 20 classes are shown

Fig. 3
Fig. 3 Co-occurrence network diagram of bacterial communities at ASV level in seawater (A, B, C, and D) and rhizosphere sediment (E, F, G, and H) samples from the N, SJ, SC, and SR zones.ASVs of the top 100 abundance in each sample are selected.The colors of spots indicate different bacterial phyla, and lines between spots indicate correlations between ASVs.Red lines, positive correlations.Blue lines, negative correlations

3 − , NO 2 − , NH 4 +
samples and C (including TC and TOC), N (including NO , and TON), and D (including MD and D50) groups for rhizosphere sediment samples were chosen to conduct the VPA analysis as shown in Fig.S10.The N group contributed the highest variation (43.67%) for the total bacterial community in the seawater samples, which was the same for the eukaryotic community (32.62%).However, the contribution from the N group to the fungal community in seawater samples was lower at 8.20% in addition to the contributions from PO 4 3− and SiO 3 2−

4 +Fig. 4
Fig. 4 The Mantel test of bacterial, eukaryotic, and fungal communities in seawater (A, B, and C) and rhizosphere sediment (D, E, and F) samples from the N, SJ, SC, and SR zones.ASVs of the top 50 abundance in each sample are selected.Families (bacteria) or phyla (eukaryota and fungi) of top three abundance are merged respectively as Spec01, Spec02, and Spec03.Line colors and thickness represent the P and R values of the Mantel test

Table 1
The properties of the co-occurrence network for bacterial, eukaryotic, and fungal communities in seawater and rhizosphere sediment samples of natural and mariculture zones