The seagrass meadow studied in the present study is situated at the estuary of the Yellow River in Dongying, Shandong Province, China (37° 51.070’ N, 119° 06.184’ W; Figure S1 ). The region features a temperate monsoon continental climate with pronounced seasonal variations. The meadow extends for 25-30 km, with a distribution width of 200-500 m from the shore to the sea, and encompasses an area of more than 1000 ha. It is the largest seagrass bed discovered in China to date (Yue et al., 2021). Field observation found that the eelgrass continuously distributed on the seaward side of the habitat of cordgrass (Sporobolus alterniflorus) and partially mixed with cordgrass. We collected seagrass-vegetated samples solely from an area covered by seagrass, without any cordgrass mixing, at a distance of approximately 200 m from the shore during the lowest tide period. Additionally, the unvegetated samples were collected from an unvegetated area located ~100 m further towards the seaward side, where seagrass growth was absent throughout the whole year of 2018.

Sampling was conducted during January, April, July, and September of 2018. Surface-layer (top 0-5cm) sediment samples were collected from a seagrass-covered area and an adjacent barren area, using a custom-made corer (10 cm in length and 1.5 cm in diameter). Triplicate samples, spaced roughly 5 meters apart, were pooled to form a single sample for both vegetated (V) and unvegetated (U) regions, resulting in a total of 8 samples. The samples were immediately transported in an ice box to the laboratory within an hour. Subsamples were taken from the well-mixed sediments and stored at −80°C for DNA extraction, while the remaining sediments were stored at −20°C for subsequent physicochemical determination.

The parameters of the overlying water, including temperature, salinity, dissolved oxygen (DO) concentrations, and pH, were measured using a YSI 556 (YSI, USA) ( Table S1 ). The homogenized sediments were subjected to freeze-drying, grinding, and sieving in preparation for the determination of sediment parameters. The grain sizes of sediments were measured using a Mastersizer 2000F laser diffraction particle size analyzer (Malvern, England). The concentrations of nitrate (NO₃ ⁻-N), nitrite (NO₂ ⁻-N), and ammonium (NH₄ ⁺-N) were analyzed using a Nutrient AutoAnalyzer (Seal, Germany). The sulfate (SO₄ ²⁻) concentrations were determined by an ion chromatography Dionex ICS 3000 (USA). The contents of total organic nitrogen (TON) and total organic carbon (TOC) were measured with a Vario Micro Cube elemental analyzer (Elementar, Germany). The concentrations of metals were determined using an ELAN DRC II plasma mass spectrometer (PerkinElmer, Hong Kong) ( Table S1 ).

Genomic DNA of triplicates was extracted from 0.5 to 1.0 g of sediment using a FastDNA® Spin Kit for Soil (MP Biomedical, USA) following the manufacturer’s protocol. The quantity of DNA was determined by a NanoDrop ND-3000 spectrophotometer (Thermo Scientific, NC, USA). The sequencing library was constructed with NEBNext Ultra II DNA PCR-free library prep kit (New England Biolabs, Ipswich, MA, USA). Shotgun sequencing was performed on the Illumina Hiseq 4000 platform with 2×150-bp paired-end reads at Majorbio Bio-Pharm technology company (Shanghai, China).

The raw metagenomic reads were preprocessed using KneadData v0.7.2 ( Table S2 ), which integrates Trimmomatic (Bolger et al., 2014) and Bowtie2 (Langmead and Salzberg, 2012) to trim low-quality sequences (parameters: ‘SLIDINGWINDOW:4:20 MINLEN:50’) and remove contaminant reads originating from the seagrass Zostera marina reference genome (GCA_001185155.1 assembly) (parameters: ‘-very-sensitive-dovetail’). Each sample’s metagenomic reads were assembled individually into contigs using MEGAHIT v1.1.3 with default parameters (Li et al., 2016). To yield more genic and ecological (biomass) information of microorganisms, all clean reads were combined and co-assembled with MEGAHIT. Assembly qualities were evaluated using QUAST v5.0.2 (Gurevich et al., 2013) and are summarized in Table S3 .

Protein-coding open reading frames (ORFs) were predicted from all contigs, which were at least ≥ 500 bp in length, using Prodigal v2.6.3 (Hyatt et al., 2010) with parameters: ‘-p meta -m’. The extracted ORFs, longer than 100 bp, were further clustered at a criterion of 95% identity (Sunagawa et al., 2015; Wang et al., 2022) to generate a non-redundant (NR) gene catalog (5,004,158 ORFs) using CD-HIT-EST v4.8.1 (Li and Godzik, 2006) with parameters ‘-T 86 -c 0.95 -M 0 -n 10 -G 0 -aS 0.9 -g 1 -d 0’. These NR ORFs were functionally annotated with KOfamscan v1.3.0 (Aramaki et al., 2020) and taxonomically annotated with DIAMOND-BLASTp v0.8.36.98 (parameters: ‘-evalue 1e-5’) against NCBI nr database (25 October 2019). The coverage and abundance (RPKM, reads per kilobase per million sample reads) were calculated by mapping high-quality reads to the NR gene set using BBMap v38.84 with default settings (Bushnell, 2014).

The functions of genes involved in nitrogen- and sulfur-cycling were identified based on the assignment of KEGG Orthology (KO) (see Table S4 ). For each predicted gene, CD-HIT was utilized to further eliminate redundant sequences with a similarity of 90%. The diversity of each gene was calculated based on the number of non-redundant sequences. Nitrogen-cycling (hao, napB, narG/nxrA, narI/narV, nirS, norB, nosZ, nrfA) and sulfur-cycling genes (aprB, dsrA, fccA, phsC, sat, soeC, soxY, sqr) were selected according to following criteria: (i) sequence length ≥ 100, (ii) genes involved assimilatory reduction and mineralization processes were excluded, and (iii) the top 8 genes were with a minimum sequence count of 173 for sulfur genes and 117 for nitrogen genes ( Table S4 ). Taxonomic annotation of these functional genes was carried out against NCBI nr database using BLASTp v.2.9.0 with an E-value cutoff of 1e-5. The full taxonomic lineages were obtained by querying the “taxids” against the NCBI taxonomy using “lineage” function in TaxonKit v0.10.1 (Shen and Ren, 2021). The composition of different phyla or classes for each functional gene was calculated as the sum of the abundance of sequences affiliated with that phylum or class.

The SSU rRNAs obtained from assembled contigs were screened using Metaxa2 (Bengtsson-Palme et al., 2015). These SSU rRNAs were further subjected to a denoising workflow that involved filtering, dereplication, chimera removal, and merging of paired-end reads. This workflow was carried out using the DADA2 pipeline (Callahan et al., 2017) in R version 3.6.3 (R Core Team, 2018), as previously described by (Liu et al., 2021). Taxonomic assignment was performed against the SILVA database (v138) and the community composition was visualized using the ggplot2 package in R software (Wickham, 2011).

Contigs with length ≥ 1000 bp were retained and binned into metagenome-assembled genomes (MAGs) utilizing the function modules of MetaWRAP v1.2.1 pipeline (Uritskiy et al., 2018) and three binning tools: MetaBAT2 v2.12.1 (Kang et al., 2019), MaxBin2 v2.2.6 (Wu et al., 2016), and CONCOCT v1.1.0 (Alneberg et al., 2014). The resulting MAGs were then refined using the ‘Bin_refinement’ module in MetaWRAP, and their completeness and contamination were evaluated by running CheckM v1.0.12 (Parks et al., 2015). The sets of MAGs from all samples were then combined and de-replicated using dRep v.2.6.2 (Olm et al., 2017) with default parameters, yielding 33 high-quality MAGs (completeness ≥ 90%, and contamination ≤ 5%) and 50 medium-quality MAGs (completeness ≥ 60%, and contamination ≤ 10%).

The final MAGs were re-annotated more thoroughly, with taxonomic classification performed using GTDB-Tk v1.2.0 (Chaumeil et al., 2020). Genes prediction was carried out using Prodigal, and protein functions were annotated through a combination of InterProScan v5.34 (Jones et al., 2014) and KOfamscan. Identification of the types and quantity of tRNA sequences was performed using tRNAscan-SE v2.0.5 (Meier-Kolthoff et al., 2013). The presence of subunit rRNA genes was screened using Metaxa2 v2.2.3 (Bengtsson-Palme et al., 2015). The coverage of each MAG was calculated as the average coverage of all contigs belonging to that MAG and normalized based on sequencing depth. A summary of the taxonomic classification, genome size, completeness, contamination, strain heterogeneity, N50, GC content, and coverage of recovered MAGs is presented in Table S7 .

The phylogeny of the final MAGs was inferred using the UBCG v3.0 pipeline (Na et al., 2018), which was based on 92 housekeeping genes. A concatenated maximum-likelihood phylogenetic tree was constructed using IQ-TREE (Nguyen et al., 2015) with 1000 bootstraps and was visualized using Interactive Tree of Life (iTOL, v5.0) platform (Letunic and Bork, 2021).

The abundance of each gene was calculated as the sum of the coverage of all ORFs assigned to that gene. Student’s t-test was performed to test the difference in relative abundance of major taxonomic groups of functional genes between seagrass-vegetated and unvegetated samples. To explore the relationships between diversity and abundance of functional genes and environmental parameters, linear regression and Spearman’s rank correlation were conducted using SPSS v.13.0 (IBM, USA).

A total of ~418 million clean metagenomic reads was retained after quality control, with an average of 52.3 ± 7.1 million reads per sample ( Table S2 ). Following assembly, 831,728 contigs with length > 1 kb were obtained, ranging from 34,537 to 149,379. The N50, representing the length at which 50% of the assembled contigs were equal to or longer than, ranged between 699 and 895 bp, with an average G+C content of 54.5%. Co-assembly resulted in a total of 1,423,701 contigs (>1 kb), with an N50 of 1,072 bp and an average G+C content of 55.6% ( Table S3 ).

The taxonomy of prokaryotic community in 8 samples was analyzed by screening SSU rRNAs from the assembled contigs. Our analysis identified the presence of 26 prokaryotic phyla, as well as the prominent subclasses within Proteobacteria ( Figure S2 ). A large proportion of reads (approximately 58.1% on average) was assigned to Proteobacteria, with variations ranging from 52.8% to 64.5% in individual samples. Within Proteobacteria, γ- (14.9%-34.5%), δ- (9.4%-30.4%), and α-proteobacteria (6.6%-19.2%) were the most abundant. While ε-, and β-proteobacteria were present in relatively lower abundance, with average proportions of less than 1%. The second most abundant bacterial phylum was Bacteroidetes, accounting for 15.3% of reads, followed by Firmicutes and Actinobacteria, which accounted for 6.0% each. The rest (e.g., Planctomycetes, and Verrucomicrobia) appeared to be minor. Compared to the bacteria, archaea were less abundant, accounted for 3.5% of the total sequences, with the majority belonging to Euryarchaeota phylum ( Figure S2 ).

After clustering protein sequences at a 90% similarity level, 4,612 ORFs were annotated as sulfur-cycling-related genes ( Table S4 ). The gene fccA, which encodes sulfide dehydrogenase, dominated the pool of sulfur metabolism genes, particularly in the summer sediments that were seagrass-vegetated ( Figure 1A ). Additionally, the abundances of dissimilatory thiosulfate reduction gene phsC (average RPKM, 89.8 vs. 66.5) and assimilation gene cysC (30.1 vs. 14.9) were notably higher in vegetated sediments compared to unvegetated sediments (P < 0.05). However, genes involved in thiosulfate oxidation (soxXYZ, 73.1 vs. 110.7), sulfate reduction (aprA, 53.4 vs. 77.2), and tetrathionate reduction (ttrB, 53.8 vs. 72.6) were more abundant in unvegetated areas, indicating a potentially different mechanism of sulfur metabolism between the two sediment types. Regarding seasonal comparisons, fccA, and sseA peaked during summer in both sediment types, while sat and aprB genes were generally higher in the autumn ( Figure 1A ).

For the nitrogen cycles, 2,343 ORFs were obtained ( Table S4 ). The highest abundance of nitrogen-cycling-related genes was hao, encoding hydroxylamine dehydrogenase, with RPKM ranging from 61.8 to 189 ( Figure 1B ), followed by narI/narV encoding nitrate reductases. The proportions of the nitrogen fixation gene (nifK, 16.5 vs. 1.5) and hydroxylamine oxidation gene (hao, 167.6 vs. 102.1) were significantly higher in vegetated sediments compared to unvegetated areas. However, the relative abundances of nitrite reduction gene (nirS, 62.8 vs. 110.7) and nitrous-oxide reduction gene (nosZ, 39.1 vs. 67.5) were notably higher in unvegetated sediments. In seagrass-covered area, nitrogen metabolism genes (hao, nifH, narI/narV, and hzo) were more abundant in spring and winter compared summer and autumn. In contrast, the barren area generally had higher gene abundance for hao, narI/narV, and gdhA during summer and autumn ( Figure 1B ).

The taxonomic affiliation and niche distribution of microbial taxa involved in sulfur and nitrogen cycles were analyzed, as depicted in Figures 2 and S3 . In the seagrass meadow, sulfur cycling was predominantly governed by proteobacterial lineages, constituting an average relative abundance of 80.7% ( Figure 2A ). Within the Proteobacteria phylum ( Figure S3 ), δ- and γ-proteobacteria played key roles in sulfate reduction, including aprB (85%-95%), dsrA (83%-99%), phsC (19%-96%), and sat (58%-87%). The δ-proteobacteria were mainly represented by Desulfosarcina genus, while γ-proteobacteria were dominated by genera Ca. Thiodiazotropha and Sedimenticola ( Table S5 ). Conversely, α- and γ-proteobacteria were predominantly associated with sulfide oxidation (fccA, soxY, soeC, and sqr), accounting for over 80% of the abundance. Predominant members of α-proteobacteria included Hyphomonas (fccA and soxY), Roseobacter (soeC), and Pseudoruegeria (sqr) ( Table S5 ). Apart from Proteobacteria, other functional groups were also identified in sulfur cycle. Acidobacteria and Chloroflexi were important sources of genes related to sulfur reduction (aprB, dsrA, phsC, and sat), with an average abundance of 8.1%. Bacteroidetes contributed to both sulfide-oxidation (sqr, 20.3%) and thiosulfate-reduction (phsC, 5.2%) processes. Additionally, Planctomycetes and Gemmatimonadetes participated in sulfite oxidation (soe) and sulfide oxidation (fccA), respectively.

Regarding nitrogen cycling, Proteobacteria was also the dominant participant of nitrogen-cycling (70.6% on average) affiliated microbial taxa in the Z. japonica system ( Figure 2B ). The gene encoding hydroxylamine oxidoreductase (hao) was mainly assigned to δ- (62.8%, mainly Desulfosarcina) and γ -proteobacteria (16.2%, mainly Ca. Thiodiazotropha), as well as Planctomycetes (8.3%). The α-, δ-, and γ-proteobacteria collectively contributed significantly to the abundance of napB, norB, and narI/narV, with proportions ranging from 52.3% to 94.5%. Furthermore, napB could also be categorized to phylum Acidobacteria (7.9%) and Bacteroidetes (10%), while norB could partially be assigned to Spirochaetes (19.4%). Contig sequences aligning to gene narG/nxrA were mostly affiliated with α- (53.9%, mainly Aestuariivita), δ-proteobacteria (22.2%, mainly Desulfatitalea and Desulfotignum), and Actinobacteria (11.7%). nrfA sequences belonging to Proteobacteria were exclusively from δ-proteobacteria (mainly Desulfuromonas and Malonomonas), while Bacteroidetes, Planctomycetes, and Chloroflexi were also its carriers. Gene nirS was affiliated with γ- (32.5%, mainly Thiohalomonas), α-proteobacteria (12.1%, mainly Roseovarius), and, Gemmatimonadetes (6.3%). The main contributors of nosZ were α- (22.3%, mainly Roseovarius), γ-proteobacteria (18.5%, mainly Sulfurivermis), Bacteroidetes (34.7%), Chloroflexi (8.6%), and Armatimonadetes (2.1%; Table S6 ).

We utilized linear regression analysis ( Figures 3 , 4 ) and Spearman’s rank correlation analysis ( Figures S4 , 5 ) to identify key geochemical parameters influencing the diversities and abundance of sulfur- and nitrogen-cycling genes. The diversities of sulfur-cycling genes, particularly gene sat (R = 0.882, P = 0.004), increased with warming ( Figure 3A ). Spearman’s rank correlation further highlighted that temperature significantly promoted sat (ρ = 0.946, P < 0.001) and sqr (ρ = 0.776, P < 0.05) diversities, while salinity and specific metal elements (Mn, As, and Cd) substantially reduced the diversities of sulfate-reducing genes dsrA and aprB (ρ < -0.731, P < 0.05; Figure S4A ).

In contrast, no significant linear relationship was observed between sulfur-genes abundance and temperature ( Figure 3B ). Spearman’s correlation revealed a positive correlation between dsrA abundance and DO (ρ = 0.714, P < 0.05) but a negative correlation with NO₃ ^- content (ρ = -0.738, P < 0.05). Additionally, gene aprB abundance was negatively correlated with Cd (ρ = -0.714, P < 0.05; Figure S4B ).

As for nitrogen-cycling genes, particularly gene narG/nxrA (R = 0.778, P = 0.023), their diversities increased with rising temperature ( Figure 4A ). Spearman’s coefficients indicated a significant positive correlation between temperature with narG/nxrA diversity (ρ = 0.795, P < 0.05), a positive correlation between NH₄ ⁺ contents and norB diversity (ρ = 0.714, P < 0.05), and a negative relationship between salinity and narI/narV diversity (ρ = -0.719, P < 0.05; Figure S5A ). Although no significant linear relationship was found between nitrogen genes abundance and temperature ( Figure 4B ), a significantly positive correlation was observed between NH₄ ⁺ and norB abundance (ρ = 0.762, P < 0.05; Figure S5B ).

We employed a genome binning strategy to comprehensively explore the genomic potential for sulfur and nitrogen cycling, resulting in 83 medium- to- high quality metagenome-assembled genomes (MAGs) after removing duplicates. In accordance with recent MAGs standards (Bowers et al., 2017), all MAGs in this study exhibited >80% completeness with < 10% contamination. Detailed statistics for each MAGs are provided in the Supplementary Table S7 . In summary, these 83 MAGs exhibited genome size ranging from 1.9 to 9.4 Mbp, GC content spanning 34.3 to 41.4%, and contained 2,102 to 8,478 ORFs. Using the GTDB-Tk tool, we classified the 83 MAGs in to 12 different phyla ( Figure 5 ; Table S7 ), with the majority classified into Proteobacteria (30 MAGs, predominantly Halieaceae, Rhodobacteraceae, and Sedimenticolaceae), Bacteroidota (21, mainly Flavobacteriaceae), and Desulfobacterota (16, mainly Desulfocapsaceae). Smaller numbers of MAGs were assigned to various other phyla, including Actinobacteriota (4, Acidimicrobiia), Desulfuromonadota (2, Desulfuromonadales), Acidobacteriota (2, Pyrinomonadaceae), Chloroflexota (2, Anaerolineaceae) and Gemmatimonadota (2, Gemmatimonadetes), Nitrospirota (1, Thermodesulfovibrionales), Planctomycetota (1, Phycisphaerae), Cyanobacteria (1, Cyanobiaceae), and BMS3Abin14 (1).

Among the 83 MAGs, each featured at least one sulfur-cycling gene ( Figure 5 ; Table S8 ). Notably, a greater proportion of MAGs contained genes for sulfide oxidation (64 MAGs) and dissimilatory sulfite reduction (60 MAGs) compared to other sulfur-cycling pathways ( Figure 6 ). Various MAGs encompassed genes encoding sulfide-oxidizing enzymes, specifically sulfide:quinone reductase (sqr, 53 MAGs) and fccAB (38 MAGs; Figure 5 ). Specific MAGs assigned to α- (Bin1) and γ-proteobacteria (Bin22, 24-25, 27-29, 66, 72, 75-78, and 82) exhibited the dual capability of encoding both sulfide dehydrogenase (fccA) and sulfide:quinone oxidoreductase (sqr) ( Figure 5 ; Table S8 ), indicating a pronounced potential for sulfide oxidation. Moreover, three MAGs within γ- (Bin27) and α-proteobacteria (Bin37 and 38) were able to encode a complete Sox-enzyme system ( Figure 5 ; Table S8 ). Additionally, 15 MAGs assigned to Desulfobacterota (Bin14-15, 19, 20, 40, 51, 56-60), Nitrospiraceae bacterium UBA6898 (Bin62), and γ-proteobacteria (Bin66, 76, and 82), were identified as hosting the full canonical pathway for sulfate reduction to sulfide (sat+qmoABC+aprAB+dsrAB) ( Figure 5 ; Table S8 ). The gene encoding 3-mercaptopyruvate sulfurtransferase (TST), involved in the mineralization of organic sulfur compounds, exhibited widespread prevalence across these MAGs (58 MAGs) ( Figure 5 ; Table S8 ).

For the nitrogen cycle, the gene encoding glutamate dehydrogenase (E1.4.1.4), which is involved in the nitrogen-mineralization, was the most widely distributed and detected in 42 MAGs ( Figure 5 ; Table S8 ), followed by narI/narV. The complete pathways of assimilatory nitrate reduction, nitrogen fixation, DNRA, and denitrification were identified in the 83 MAGs ( Figure 5 ), and the number of MAGs containing nitrate reduction genes was more abundant (33 MAGs) than other nitrogen-cycling steps ( Figure 6 ). Along the denitrification pathway, MAGs associated with Desulfobacterota (NO₃ ^- → NO₂ ^- and NO → N₂O), γ-proteobacteria (NO₂ ^- → NO), and Bacteroidota (N₂O → N₂) were dominant in different steps, respectively ( Figure 5 ; Table S8 ). Additionally, five MAGs from Desulfobacterota (Bin48 and 58), Gemmatimonadota (Bin63), and γ-proteobacteria (Bin72 and 76) contained both dissimilatory nitrate- and nitrite-reduction gene clusters (narGHI/napAB+nirBD/nrfAH), which have the potential to reduce nitrate to ammonium. Desulfobacterota was found to be the main nitrogen-fixing bacteria ( Figure 5 ; Table S8 ).

In all metagenomes, we identified gene sets responsible for dissimilatory sulfate reduction and sulfide oxidation, highlighting the ecological significance of sulfate respiration and sulfide removal in the seagrass meadow ecosystem. Notably, the gene fccA was found to be the most abundant among sulfur-cycling gene pool in the seagrass sediments. The prevalence may be attributed to the provision of aerobic conditions created by seagrass, which are favorable for SOB. In contrast, sqr, the homolog of fccA, appeared relatively scarce. The gene fcc encodes flavocytochrome c, typically utilized for sulfide oxidation under low sulfide concentrations and aerobic conditions, while sqr encodes enzyme sulfide:quinone oxidoreductase, primarily employed for oxidizing sulfide under high sulfide concentrations (Cúcio et al., 2018). The observed high abundance of sulfide-oxidation genes potentially signifies a rapid recycling mechanism for reduced sulfides in sediments colonized by seagrass. This aligns with previous research by Holmer et al. (2003), which highlighted an accelerated turnover of reduced sulfides in these seagrass-covered sediments.

Our study observed that seagrass-vegetated sediments specifically promoted the abundance of thiosulfate reductase gene phsC, which is responsible for sulfide production ( Table S5 ). This result is in line with previous reports of enrichment of sulfate-reduction bacteria in seagrass-vegetated sediments (Holmer et al., 2005; Sun et al., 2015; Cúcio et al., 2018). These findings indicate that seagrass vegetation may stimulate the microbial capacity to reduce thiosulfate, with thiosulfate serving as a major electron acceptor in sedimentary microbial communities involved in sulfur cycle. In general, SRB utilizes thiosulfate as an electron acceptor to completely oxidize various fatty acids and alcohols to CO₂ (Brysch et al., 1987) or incompletely to acetate (Brenner et al., 2005). Thus, these thiosulfate-utilizing bacteria aid in the mineralization and recycling of organic matter in productive seagrass-vegetated sediments, and also provide the seagrass with an alcohol detoxification process as a compensation for the sulfide they produce (Cúcio et al., 2016). In contrast, the abundances of genes related to sulfur oxidation (soxXYZ), sulfate reduction (aprA), and tetrathionate reduction (ttrB) were significantly higher in the samples from the unvegetated areas and are involved in the formation of sulfur cycle intermediates (SCIs), including sulfite, thiosulfate, tetrathionate, elemental sulfur, polysulfides. SCIs play a significant role in ecological and biogeochemical perspective in marine environments, serving as an energy source for microbial communities (Vigneron et al., 2021) and creating shortcuts in tropical sulfur metabolism (Jørgensen et al., 2019). Hence, the enrichment of genes involved in SCIs production may provide additional energy for microbial guilds in the oligotrophic unvegetated environments.

Taxonomic annotation of sulfur-related genes revealed that Ca. Thiodiazotropha and Sedimenticola, members of γ-proteobacteria, accounted for a large portion of the genes involved in dissimilatory sulfate reduction (aprB, dsrA, phsC, and sat; Table S5 ). Ca. Thiodiazotropha, prominent member of γ-proteobacteria in the Z. japonica system, is typically known as a symbiotic bacteria of lucinids, which are frequently found in seagrass sediments (Osvatic et al., 2021). Notably, a comprehensive global survey of rhizosphere microbial diversity across 14 seagrass species revealed that Ca. Thiodiazotropha is abundant in seagrass environments worldwide, regardless of the presence of lucinids. This suggests that Ca. Thiodiazotropha offers a mechanism for seagrasses to alleviate sulfide stress (Martin et al., 2020). Additionally, members of δ-proteobacteria, particularly Desulfosarcina significantly contribute to dissimilatory sulfate reduction ( Table S5 ), consistent with our metagenome annotation ( Table S8 ). The genus Desulfosarcina predominates among δ-proteobacteria in the rhizosphere of various seagrass species, including Cymodocea nodosa, Zostera marina, Zostera noltii, Halophila ovalis, Halodule wrightii, Syringodium isoetifolium, Thalassia hemprichii, and Thalassia testudinum (Cúcio et al., 2016; Wang et al., 2020; Aires et al., 2021; Ling et al., 2021b). Previous study have demonstrated that members of genera Desulfosarcina are capable of utilizing a diverse range of electron donors, such as ethanol, benzoate, lactate, formate, H₂, propionate, lactate, and 1-propanol, for sulfate reduction (Kuever et al., 2015; Cúcio et al., 2018). Consequently, the genus Desulfosarcina was conjectured to contribute significantly for nutrient availability and remineralization in seagrass sediments (Holmer et al., 2001). In seagrass-vegetated sediments, high concentrations of various carbon substrates supplied from root exudates and organic matter decomposition provide a competitive advantage for the metabolically versatile SRB (Smith et al., 2004).

Sulfide oxidation processes (fccA, soxY, soeC, and sqr) are primarily carried out by genus Hyphomonas, Roseobacter, and Pseudoruegeria ( Table S5 ). These genera share a common ecological strategy, where Hyphomonas utilizes a cylindrical prostheca as a holdfast to attach to the surface of seagrass leaves (Halophila stipulacea) (Weidner et al., 2000). Similarly, certain Roseobacter and Pseudoruegeria strains are reported to form biofilms on surfaces of marine algae, invertebrates, and particles (Ding et al., 2023). Above-mentioned these gene-level findings aligns well with genome binning analysis, emphasizing the pivotal role of gproteobacteria in sulfur metabolism in Z. japonica meadows, participating in both sulfate reduction and sulfide oxidation.

In this study, we observed an increase in diversity of sulfur-cycling genes, specifically the gene sat, with higher temperature ( Figure 3A ). Previous studies have demonstrated that increased temperature can affect diversity and intensify sulfate reduction reactions of SRB (Knoblauch and Jørgensen, 1999; Robador et al., 2016). For instance, Ma et al. (2021) found that the diversity of dsrB gene was highest observed at mid-temperature in hot springs. Given the established association between seagrass loss with elevated temperatures and sulfides accumulation (Koch et al., 2007; García et al., 2012; García et al., 2013), further research on sulfur biogeochemical cycles is imperative to understand the impact of global warming on seagrass ecosystems.

The gene hao, responsible for hydroxylamine oxidation to nitrite, was abundant in vegetated sediment compared to unvegetated sediment ( Figure 1B ). Hydroxylamine has been identified as a key nitrogen metabolite for microbial interactions within microbial communities and engineered systems (Soler-Jofra et al., 2021). Adding external hydroxylamine to biofilm system has been shown to disperse microcolonies into scattered cells (Kindaichi et al., 2004; Harper et al., 2009). Therefore, the prevalence of hao gene suggests a significant role in shaping the microbial communities and the biofilm dynamics in seagrass ecosystems. Additionally, we observed a higher level of the nitrogenase gene (nifK) in the vegetated sediments ( Figure 1B ), consistent with the findings of Sun et al. (2015). Nitrogen fixation by rhizosphere microbiomes is considered a major contributor to seagrasses’ nitrogen requirements (Welsh, 2000). The presence of diazotrophic assemblages has been characterized on seagrass surfaces through sequencing nitrogenase gene (nifH) (Sun et al., 2015) and 16S rRNA amplicons (Jensen et al., 2007). In contrast, the low abundance of nitrogen-fixing encoding genes (nifK) and the high abundance of denitrification genes (nirS and nosZ) in the unvegetated sediments imply that nitrogen was not a limiting nutrient in the unvegetated environments.

In agreement with our expectations, the majority of nitrogen-cycling genes and MAGs derived from metagenomic data exhibited affiliation with δ- and γ-Proteobacteria ( Figure S3B and Figure 5 ). This suggests that Proteobacteria play a crucial role in the nitrogen cycling process, with various lineages of this phylum contributing differently to the abundance of different nitrogen-cycling genes. For example, genera that are potentially associated with nitrogen cycling were mainly affiliated with Desulfosarcina, Desulfuromonas, and Malonomonas ( Table S6 ), have the capacity to use sulfate as a terminal electron acceptor (Mills et al., 2012).

In concordance with sulfur-cycling genes, we observed that the diversity of nitrogen cycle genes increased with increasing temperature ( Figure 4A ), especially for nitrification gene nxrA, which is responsible for nitrate (NO₃ ^–N) and nitrite-nitrogen (NO₂ ^–N) conversion. This corresponds with previous findings that temperature is positively correlated with nitrification rates in estuaries or coastal systems (Dai et al., 2008; Damashek et al., 2016). Elevated nitrification may result from increased ammonium regeneration due to enhanced microbial respiration at high temperatures (Caffrey, 2004), or it may be due to higher growth rates of nitrifying bacteria at higher temperatures. Additionally, denitrification genes narI/narV, encoding nitrate reductase responsible for the transformation of nitrate to nitrite (Sohaskey and Wayne, 2003), exhibited reduced diversity in response to salinity ( Figure S5A ). This observation is consistent with the widely accepted concept that salinity acts as a stressor on denitrification process by inhibiting microorganism activities (Miao et al., 2015; Zhai et al., 2020).

Approximately 68% of the draft genomes were found to contain genes responsible for partial nitrogen and sulfur metabolisms, indicating tight couplings between sulfur and nitrogen cycles. The most abundant MAGs classified as Proteobacteria suggest that this taxonomic group may represent the most abundant active population. Among the members of Proteobacteria, γ- and α-proteobacteria were genetically equipped to perform denitrification and oxidation of reduced inorganic sulfur compounds. Different taxa were equipped with genes for different metabolic steps, implying a cooperative effort among microbial species. The dissimilatory nitrate reduction (to ammonium or denitrification) coupled with sulfur oxidation is a significant energy source for inorganic carbon fixation (Nakagawa and Takai, 2008; Li et al., 2018) and utilizes nitrate and nitrite as electron acceptors. For instance, the frequently detected SOB Sulfurovum sp. in seagrass systems (Sun et al., 2020; Zhang et al., 2020; Martin et al., 2022; Zhang et al., 2022), have been shown to use nitrate as an electron acceptor (Li et al., 2018).

It has been reported that fix nitrogen is widespread among SRB, which were previously considered the dominant diazotrophs in seagrass-colonized sediments (Welsh, 2000; Smith et al., 2004). Studies have shown that the nitrogen fixation fulfilled by the sulfate-reducing community could supply 60-95% of the nitrogen requirements of seagrasses (Welsh, 2000; Smith et al., 2004). Similarly, the dominant members of the nifDHK host guild in the Z. japonica system were found to be sulfate-reducing Desulfobacterota, including the families Desulfococcaceae, Desulfocapsaceae, and Desulfosarcinaceae ( Figure 5 ; Table S6 ). This finding provides molecular evidence for the tight coupling of nitrogen fixation and sulfate reduction in seagrass meadows (Welsh, 2000; Smith et al., 2004; Ling et al., 2021a). Additionally, we also identified the potential of the SRB genomes from Desulfobacterota and γ-proteobacteria participate in the reductive sulfur cycle processes by oxidizing hydroxylamine or nitrite ( Figure 5 ). This suggests that they may facilitate the connection between sulfate reduction and nitrification processes. Alternatively, the presence of diverse functional profiles in these microorganisms allows for metabolic switch between available compounds in seagrass ecosystems (Cúcio et al., 2018). Our results suggest a tight interdependence between sulfur- and nitrogen-metabolite exchange, with cooperation among various microbial species being crucial for the completion of benthic nitrogen and sulfur cycles in seagrass meadows.

This study sheds light on the unknown seasonal variations in sulfur- and nitrogen-cycling genes between seagrass-vegetated and unvegetated sediments, thus supporting the hypothesis that seagrass colonization could affect microbial-mediated biogeochemical cycles. We found that genes involved in sulfide oxidation (fccA) and hydroxylamine oxidation (hao) dominated the gene pools of sulfur- and nitrogen cycles. Microorganisms residing in vegetated sediments exhibited higher proficiency in sulfite reduction (via phsC) and nitrogen fixation (nifK), whereas those in unvegetated sediments displayed enhanced activity in sulfur oxidation (soxXYZ) and denitrification (nosZ and nirS). We provided a finer taxonomic resolution, revealing that microbial participants in different processes are highly diverse and their relationship with environmental drivers such as temperature and salinity. The dominance of proteobacterial lineages with versatile functions in both sulfur and nitrogen cycles, indicating tight coupling between these processes in seagrass ecosystems. This was reinforced by annotation of 83 MAGs recovered from our metagenomic dataset. Taken together, our work provides valuable insights into the complex and interconnected biogeochemical processes within seagrass ecosystems, emphasizing the importance of microbial communities in mediating biogeochemical cycling in seagrass Z. japonica habitats.