Chemical pollution drives taxonomic and functional
shifts in marine sediment microbiome, influencing
benthic metazoans

Lin-Lin Xu, Shelby E Mcilroy, Yueqiong Ni, Isis Guibert, Jiarui Chen, Ulisses

Rocha, David M Baker, Gianni Panagiotou

► To cite this version:

Lin-Lin Xu, Shelby E Mcilroy, Yueqiong Ni, Isis Guibert, Jiarui Chen, et al.. Chemical pollution drives taxonomic and functional shifts in marine sediment microbiome, influencing benthic metazoans. ISME Communications, 2025, 5 (1), 10.1093/ismeco/ycae141. [https://doi.org/10.1093/ismeco/ycae141](https://doi.org/10.1093/ismeco/ycae141). [hal-05449739](https://hal.science/hal-05449739v1)

HAL Id: hal-05449739

https://hal.science/hal-05449739v1

Submitted on 9 Jan 2026

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Chemical pollution drives taxonomic and functional shifts in marine sediment microbiome, influencing benthic metazoans

Lin-Lin Xu<fnref n="1" />, Shelby E. McIlroy<fnref n="2" />,<fnref n="3" />,<fnref n="4" />, Yueqiong Ni<fnref n="1" />,<fnref n="5" />,<fnref n="6" />, Isis Guibert<fnref n="2" />,<fnref n="3" />, Jiarui Chen<fnref n="7" />, Ulisses Rocha<fnref n="8" />, David M. Baker<fnref n="2" />,<fnref n="3" />,*, Gianni Panagiotou<fnref n="1" />,<fnref n="5" />,<fnref n="9" />,†

## Abstract

Microbial communities in marine sediments contribute significantly to the overall health and resiliency of marine ecosystems. However, increased human disturbance undermines biodiversity and, hence, natural functionality provided by marine sediments. Here, through a deep shotgun metagenomics sequencing of the sediment microbiome and COI metabarcoding of benthic metazoans, we demonstrate that >50% of the microorganisms' and metazoan's taxonomic variation can be explained by specific chemical pollution indices. Interestingly, there was a significant correlation between the similarity in microbiome communities' taxonomical and functional attributes and the similarity of benthic metazoans community composition. Furthermore, mediation analysis was conducted to evaluate the microbiome-mediated indirect effect, suggesting that microbial species and functions accounted for 36% and 26%, respectively, of the total effect of pollution on the benthic metazoans. Our study introduces a multi-level perspective for future studies in urbanized coastal areas to explore marine ecosystems, revealing the impact of pollution stress on microbiome communities and their critical biogeochemical functions, which in turn may influence macrofaunal composition.

**Keywords:** shotgun metagenomics, pollution, marine sediments, DNA metabarcoding, mediation analysis

## Introduction

Coastal marine habitats, including a large variety of habitats such as estuaries, mangroves, and coral reefs, are the world's most heavily used and commercially important ecosystems [1]. These ecosystems provide numerous benefits, including food, recreation, protection from storms and erosion, and economic opportunities. For example, marine fisheries provide >260 million livelihoods, generating sizeable revenues for many countries, including US$ 80 billion dollars in export revenues for developing countries [2]. Despite the important services that coastal marine ecosystems bring to humanity, they face various anthropogenic threats [1]. Particularly, heavy metal pollution in coastal marine environments from domestic and industrial effluents has

drawn attention for decades [3] and continues to be a focal point, especially in developing countries [4, 5]. Additionally, polycyclic aromatic hydrocarbons (PAHs), which exist as hazardous complex mixtures in nature and accumulate in sediment environments, pose potential risks to benthic organisms [6] and human health [7].

Sedimentary habitats span the majority of the ocean floor and make up more than two-thirds of the Earth's surface, support diverse benthic communities. Microbial communities in marine sediments play an important role in biogeochemical cycling [8–10], and provide chemical cues that are essential for the recruitment of marine invertebrates [11]. Urban development has increased drastically in the last few decades, bringing pollution

<footnote n="1">¹Microbiome Dynamics, Leibniz Institute for Natural Product Research and Infection Biology—Hans Knöll Institute, Adolf-Reichwein-Strasse 23, Jena, Thuringia, 07745, Germany</footnote>

<footnote n="2">²The Swire Institute of Marine Science, The University of Hong Kong, Cape D'Aguilar Road, Shek O, Hong Kong SAR, P.R. China</footnote>

<footnote n="3">³School of Biological Sciences, The University of Hong Kong, Pok Fu Lam, Hong Kong SAR, P.R. China</footnote>

<footnote n="4">⁴Simon F. S. Li Marine Science Laboratory, School of Life Sciences, The Chinese University of Hong Kong, Shatin, Hong Kong SAR, P.R. China</footnote>

<footnote n="5">⁵Cluster of Excellence Balance of the Microverse, Friedrich Schiller University Jena, Fürstengraben 1, Jena, Thuringia, 07743, Germany</footnote>

<footnote n="6">⁶Shanghai Key Laboratory of Diabetes Mellitus, Department of Endocrinology and Metabolism, Shanghai Diabetes Institute, Shanghai Clinical Center for Diabetes, Shanghai Sixth People's Hospital Affiliated to Shanghai Jiao Tong University School of Medicine, 227 South Chongqing Road, Shanghai, 200025, P.R. China</footnote>

<footnote n="7">⁷Guangdong Provincial Key Laboratory of Food, Nutrition and Health, and Department of Nutrition, School of Public Health, Sun Yat-Sen University, Zhongshan Rd 74, Guangdong, Guangzhou, 510080, P.R. China</footnote>

<footnote n="8">⁸Department of Applied Microbial Ecology, Helmholtz Centre for Environmental Research—UFZ GmbH, Permoserstrasse 15, Leipzig, Saxony, 04318, Germany</footnote>

<footnote n="9">⁹Faculty of Biological Sciences, Friedrich Schiller University, Fürstengraben 1, Jena, Thuringia, 07743, Germany</footnote>

<footnote>*Corresponding authors: David M. Baker, School of Biological Sciences, The University of Hong Kong, Kadoorie Biological Sciences Building 6S-14, Pokfulam Road, Lung Fu Shan, Hong Kong, P.R. China. Email: dmbaker@hku.hk and Gianni Panagiotou, Microbiome Dynamics, Leibniz Institute for Natural Product Research and Infection Biology Hans Knöll Institute, Adolf-Reichwein-Strasse 23, Jena, Thuringia, 07745, Germany. Email: Gianni.Panagiotou@leibniz-hki.de</footnote>
and impacting the benthic marine ecosystem, including marine sediment microbiome and benthic metazoans. Pollutants have been demonstrated to affect both composition and function of sediment microbial communities, including nitrogen cycling [12] and antibiotic resistance genes (ARGs) [13]. Furthermore, nematodes, benthic metazoans commonly found in marine sediments, exhibit bioaccumulation of heavy metals, which can lead to biomagnification and increase the transfer of contaminants along food webs to larger organisms [14]. Evidence indicates that human-introduced antibiotics alter environmental microbiomes, shifting zooplankton microbiomes and affecting their growth [15]. Additionally, microbial biomass likely serves as a crucial food source for hadal meiofauna, driving increased abundance and variations in their densities within sediments [16]. However, while changes resulting from pollution stressors on marine sedimental microbial diversity and function are well established, their impact on macrofauna has remained largely unexplored. Due to our limited understanding of the essential microbial functions in the marine ecosystem, our capacity to benefit from such interactions in conservation practices is severely constrained.

Marine animals live in close relationships with microorganisms that modulate the environmental conditions animals are exposed to, impacting the resilience and adaptation of marine holobionts [17]. The extent to which the marine sediment microbiome taxa and function may affect macrofauna, however, has not been thoroughly investigated. Here, we examined the role of the marine microbiome in the coastal maritime zones across various pollution sources by integrating sediment quality indices, marine sediment microbiomes, and benthic metazoans. To assess benthic metazoan communities, we used Autonomous Reef Monitoring Structures ARMS, well-established, standardized, passive samplers which, when combined with metabarcoding approaches, are ideal to characterize the largest portions of marine biodiversity [18, 19]. Deep metagenomics sequencing of the marine microbiomes allowed the construction of high-quality assemblies for functional characterization of the microbiota and revealed the shifts triggered by pollution stressors. By combining a variety of methodological approaches, we demonstrated a strong association between pollution-driven microbiome shifts and benthic metazoans. The identified microbiome-benthic metazoan signatures not only further characterized the coastal marine ecosystem under pollution stress, but also underscored the necessity for more holistic approaches in conservation practices. The comprehensive understanding of marine ecosystems from multiple perspectives will shed light on protecting marine biodiversity, ecosystem services, and the sustainable use of marine resources in the context of a rapidly developing metropolis's coastal shorelines.

## Materials and methods
### Sampling sites and sample collection

Seven locations [San Shek Wan (SSW), Peng Chau (PC), Center Island (CI), Sham Wan (SW), Cape D'aguilar (CDA), Bluff Island (BI), Tung Ping Chau (TPC)] in Hong Kong's major coastal marine were studied. ARMS were deployed near the shore at ~4 meters depth in Jan 2018 (n=3 per site) and July 2018 (n=3 per site), and retrieved 12 months later in Jan 2019 and July 2019, respectively. The three ARMS were placed ~5 meters apart at each site. Just prior to the retrieval of each ARMS, 50 g of sediments from the top 10 cm of seafloor, ~10 cm from ARMS was scooped into a sterile 50 mL falcon tube. Sediments were

snap frozen in liquid nitrogen once on board and stored at -80°C upon return to the laboratory for subsequent metagenomic analyses.

## Pollution data and benthic metazoans data collection

Geochemical indicators across in year 2019 were obtained from the Hong Kong Environmental Protection Department's periodic sediment monitoring database (https://cd.epic.epd.gov.hk/EPICRIVER/marine/?lang=en). Monitoring stations were selected based on their proximity to the ARMS deployments and only the nearest date to sediment collection was considered. The average distance between the pollution monitoring site and the sampling locations across all seven sites is 4.4 kilometers (Table S1). The details of ARMS processing, metazoans data collection and processing are available in McIlroy et al [20]. Briefly, each ARMS is a ~22 x 22 x 25 cm structure consisting of 9 stacked polyvinyl chloride (PVC) plates with alternating open and partially closed spaces between. Following field deployment, during which small motile and sessile organisms naturally recruit and settle onto ARMS, the ARMS are removed and taken apart. All motile organisms between 100 and 500 um are collected by sieving, and the bulk sessile communities are scrapped from the plate and homogenized. Those samples were preserved in 90% ETOH, extracted using Qiagen Power Soil kits, and COI sequences are amplified and sequenced using standard metabarcoding techniques [21]. The forward and reverse sequence reads were processed with Cutadapt [22] and DADA2 pipeline [23] for quality control. Sequences were further aligned to Moorea BIOCODE library with MASCE [24]. Only sequences that had zero stop codons, zero frameshifts, zero insertions and no more than three deletions were retained and the remaining amplicon sequencing variants (ASVs) were clustered using VSEARCH [25] into operational taxonomic units (OTUs) based on 97% sequence similarity. OTUs were assigned to taxonomic groups using a custom-built database including the CO-ARBitrator database [26], sequences from GEOME [27] and 735 barcodes generated from large motile macrofauna of the same ARMS using VSEARCH within QIIME2 [28]. A 95% sequence similarity threshold was set for species-level taxonomic assignments where feasible. Unassigned OTUs underwent a second round of identification with an 80% threshold, followed by BLAST searches against the entire National Center for Biotechnology Information (NCBI) NT database (as of May 2021; word size=7; max e-value=5e-13). Taxonomy is assigned based on the lowest common ancestor (LCA) of the top 100 hits.

## Deoxyribonucleic acid extraction and metagenomic sequencing

Deoxyribonucleic acid (DNA) from the sediment samples was extracted using the Mag-Bind Stool DNA 96 Kit following the manufacturer's instructions. DNA concentrations were verified using the Qubit Fluorometer and DNA quality was assessed via agarose gel electrophoresis (using Takara λ-Hind III digest and the Tiangen D2000 marker) prior to library construction. DNA libraries were constructed using the BGI in-house library preparation kit. Metagenomic shotgun sequencing was then performed on an Illumina HiSeq Xten (150 bp PE) platform.

## Taxonomic profiling

Human contamination reads were removed after mapping raw sequencing paired-end reads to human genome sequences with bwa mem using default parameters. After that, we used
pipelineForQC.pl [29] to perform quality control and remove adapter regions, low-quality reads, and duplicate reads. DIA-MOND v2.0.14.152 BLASTX [30] was used to map these reads to the NCBI nonredundant (nr) database (as of February 2024) using an e-value <1e-10 cutoff. The LCA algorithm was used to estimate the taxonomic composition of the samples using the daa-meganizer from MEGAN6 [31] with default parameters (blast hits with a bit score <50, and hits outside the top 10% of the highest bit score were excluded). Reads with bacteria or archaea information were used for further analyses. Prokaryotic community profiles were constructed at each taxonomic level for further statistical analysis.

## Microbial community composition

R package vegan [32] was used to compute alpha diversity indices describing each sample's microbial community composition. The Shannon and Simpson indices were used to assess alpha diversity based on the relative abundance. Bray-Curtis distances calculated by R package vegan based on the relative abundance of each taxon at different levels were used to estimate community dissimilarities. Adonis from the vegan package in R was used as a permutational multivariate analysis of variance (PERMANOVA) test to determine the significance of a variable in determining distance variation (the number of permutations were set as 999). Microbiome vectors were overlaid onto the PcoA using the envfit function within the vegan package and the false discovery rate (FDR) is adjusted by Benjamini-Hochberg's (BH) multiple-testing correction method ($P < .05$, $FDR < 0.1$). Distance-based redundancy analyses (dbRDA) were performed with R package vegan, using Bray-Curtis distances. dbRDA analyses were first performed using the capscale function from vegan with 35 pollution indices collected from sediment. A nonredundant contribution to variation of these variables was calculated using forward stepwise variable selection via the ordiR2step function from vegan. Univariate analyses with both chemical pollution indices and six environmental variables were calculated and model p values were corrected using BH multiple-testing correction to select variables with BH-adjusted P values <.05. Hierarchical clustering of pollution indices, selected through stepwise dbRDA, was performed using the ComplexHeatmap package (v2.6.2) in R. The clustering was based on the Ward.D method with Euclidean distance.

## Metagenome data assembly and functional annotation

Clean reads passed quality control were used for assembly by metaSPAdes from metaWRAP v1.2.3 [33]. The resulting assemblies were subjected to eukaryotic contig assignment using EukRep v0.6.7 [34]. Eukaryotic contigs were filtered out from the original assemblies, resulting in individual "prokaryotic-enriched" assembly files for each sample. Metagenemark v3.38 [35] was used to predict genes from each assembled genome. Predicted genes from each assemblies were aligned to the database of Clusters of Orthologous Genes (COG) [36] and the KEGG Orthology (KO) database [37] built by KOBAS v3.0 [38] using DIAMOND [30] v0.9.19.120 BLASTP (-e 1e-10, best hits reserved). The identified COGs were further annotated into different pathways and categories based on predefined collections in the COG database and were quantified by copies per million (CPM). The identified KOs were quantified by CPM and a curated heavy metal resistance KO list [39] based on BacMet database [40] were used to retrieve metal resistance KOs.

## Data analyses and statistics

### Differential analyses

Differential abundance analyses of prokaryotic- and benthic metazoans community profiles were performed using R package MaAsLin2 [41] for prokaryotic taxa profile (compound Poisson linear model: CPLM model), prokaryotic function profile (linear model: LM model; log transform), and benthic metazoans' taxa profile (linear model: LM model) to identify pollution-associated microbes (season as fixed effect) and season-associated microbes (sampling sites as random effect).

### Correlation analyses

Inter-dependencies among the chemical pollution indices and environmental variables were evaluated using Pearson correlation with R package Hmisc. Correlations between pollution-associated microbes and chemical pollution variables were calculated using partial Spearman correlation analysis, adjusted for season, with the R package ppcor. Correlations between heavy metals and heavy metal resistance KOs markers were assessed using Spearman correlation with the R package ggpubr.

### Richness comparison

Pearson correlation between community richness was calculated by R package ggpubr. Mantel test with function mantel from R package vegan, using Spearman's correlation coefficient with 1000 times permutations was used to analyse the associations of dissimilarity between microbiome community and benthic metazoans community. Bray-Curtis dissimilarity matrices for the microbiome community and Jaccard dissimilarity matrix for the benthic metazoans community were computed to perform this test.

### Mediation analyses

Mediation associations on the community level were evaluated by R package method for multivariate omnibus distance mediation analysis (MODIMA) [42] with Euclidean distance for chemical pollution indices, Bray-Curtis distance for the microbiome community and Jaccard distance for the benthic metazoans. R package LDM-med [43] is used to validate the mediation associations on the community level. Individual-level mediation analysis was carried out with the mediate function from the R package mediation (version 4.5.0) [44], using pollution-associated features identified with R package MaAsLin2 ($P \le .05$, $FDR \le 0.2$), including microbial species, COG pathways, and benthic metazoans OTUs. The BH procedure was applied to control FDR.

## Results

### Specific pollutants explain microbiome community variation in marine sediments

Seven sites across Hong Kong water zone (Fig. 1A) chosen for this study were sampled in 2 seasons (2019: summer and winter, n = 37). Thirty-five sediment related pollution indices and six environmental variables were gathered from the monitoring stations close to each sampling site's depth and geographical location (Table S1). To characterize the sediment microbiome community, deep metagenomics sequencing (to facilitate assembly in the downstream analysis of the microbiomes) of the sediment samples generated on average 70 million reads (standard deviation: 10 million reads) of ~150 nucleotide length per sample after quality control (Table S2). To reveal the composition of the
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Sedimental pollution shapes microbiome community. (A) Map of seven sampling sites in Hong Kong: SSW, PC, SW, CDA, CI, BI, and TPC. Map of the geographic distribution of the samples were constructed using the R package ggmap. (B) Taxonomic profile of the prokaryote community at the genus level. (C) dbRDA plot based on Bray-Curtis distance. The red dashed line represents the cut-off for significant nonredundant contribution to the multivariate model. (D) Heatmap of microbiome associated pollution variables in seven sites with hierarchical clustering. Averaged and scaled value of pollution variables are used as heatmap colors reflecting standard deviations from the mean. Clustering is based on Euclidean distances. (E) Beta diversity of prokaryote community at the species level using Bray-Curtis distance.</figcaption></figure>
marine sediment prokaryotic communities, we extracted bacte-
rial and archael reads among samples at different taxonomic
levels from kingdom to species (Table S3). The relative abun-
dance of bacteria and archaea in the sediments was 97.5% and
2.5%, respectively. At the phylum level, the most abundant bac-
teria and archaea were *Proteobacteria* (59.0 ± 5.9%), *Planctomycetes*
(9.9 ± 2.1%), *Chloroflexi* (5.2 ± 4.4%), *Bacteroidetes* (5.1 ± 2.1%), and
*Acidobacteria* (3.6 ± 0.7%). The most abundant genus and species
were *Candidatus Entotheonella* (9.8 ± 7.0%) and *Gammaproteobacteria*
*bacterium* (13.3 ± 2.4%; Fig. 1B; Fig. S1).

We explored potential associations between the species-level
microbiota composition and pollution indices. More than 50%
of the variance in the species-level microbiota composition
was significantly explained by all covariates together (dbRDA,
R² = 76.6%, P = .001). In univariate analysis, 11 of these covariates
showed a significant correlation to microbiota composition
(dbRDA, FDR < 0.05). These significant covariates were classified
as heavy metal, nutrient, PAHs, and inorganic pollutants
(Table S4). Considering the potential multicollinearity between
covariates (Fig. S1), we further applied stepwise dbRDA to identify
the group of pollution indices that could explain variation in
the microbiome community without redundancy. Out of all
covariates, arsenic, total sulfide, ammonia nitrogen, pyrene, total
cyanide, lead, copper, and boron were found to explain both
significant and nonredundant fractions of compositional micro-
biome variation (Fig. 1C, stepwise dbRDA, R² = 70.4%, P = .001).
Environmental variables were assessed; however, none showed
significant contributions in the univariate dbRDA analysis, nor
were they selected in the stepwise dbRDA model (Table S4). These
microbiome community related pollution indices selected by
stepwise dbRDA were further used for the hierarchical clustering
of the seven sites revealing complex and diverse pollution sources
for each site (Fig. 1D). Sampling sites (—TPC, BI, and CDA—) were
clustered together showing relatively low pollution level, whereas
sampling sites (—SSW, PC, and CI, and SW—) showing the opposite.
SSW, PC, and CI are significantly impacted by human activities [45,
46], leading to high nutrient loads and pollution.

To evaluate the differences in community diversity between
high- and low-pollution sites, as clustered in Fig. 1D, we calculated
the alpha- (Shannon and Simpson diversity) and beta-diversity.
Beta diversity (Bray–Curtis distance) highlighted significantly dif-
ferent microbiome communities among the high and low pollu-
tion sites at the species and genus (Fig. 1E; Fig. S2) level (PER-
MANOVA test, species: R² = 0.15, P = .003, genus: R² = 0.18, P = .001).
However, no significant difference between seasons (summer and
winter) were observed in beta-diversity (PERMANOVA test: P > .05)
or alpha-diversity (Wilcoxon signed-rank test: P > .05; Fig. S2).

Specific pollutants associated with abundance
changes of bacteria taxa in marine sediments

Since our sampling included both the summer and winter periods, we built linear mixed-effects models with sampling sites as a random factor to find significant changes in the abundance of bacteria and archaea associated with season (Fig. S3). In total, we found 3 genera and 14 species significantly affected by the season using CPLM model implemented in MaAsLin2 (P ≤ .05, FDR ≤ 0.2), in which *Ilumatobacter* and *Synechococcus* are more abundant in summer, whereas *Candidatus Promineofilum* is enriched in winter. Therefore, when we subsequently explored the microbiome relative abundance changes at different taxonomic levels between high- and low- pollution sites, we adjusted for the season using CPLM model implemented in MaAsLin2 (P ≤ .05, FDR ≤ 0.2). We found 19 phyla, 95 genera and 93 species were

significantly affected by the pollution (Fig. 2A and B; Fig. S3;
Table S5). Six of the nineteen phyla have the most species
(number > 3) significantly affected by pollution, including both
pollution-sensitive phyla (*Proteobacteria*, *Planctomycetes*, and
*Candidatus Tectomicrobia*) and pollution-tolerant phyla (*Chloroflexi*,
*Bacteroidetes*, and *Nitrospinae*). Consistent with previous studies
[13] the phylum *Planctomycetes* showed a significant decreasing
abundance from low to high-polluted sites (Fig. 2A), which our
data suggest is season-independent. In contrast, the phylum
Nitrospinae, including *Nitrospina* and *Candidatus Nitronauta*, showed
an opposite trend (Fig. 2A). Although *Nitrospinae* is important to
the nitrogen cycle in marine ecosystems, it can also contribute
to the eutrophication [47]. In addition, we found a decreased
abundance of genus *Ilumatobacter* and its species *Ilumatobacter
nonamiensis* and *Ilumatobacter coccineus* associated with the
elevated pollution level (Table S5), which may be related to the
sensitivity of those microorganisms to metal contamination [48].
*Proteobacteria*, known for driving carbon cycling through various
processes, such as primary production (sulfur and methane
oxidation [49]), heterotrophy [50], and photosynthesis [51], was
found in our study to include both pollution sensitive and
pollution tolerant species (Fig. 2B).

To delineate the associations between pollution concentrations
and microbial abundances in marine sediments we applied partial
Spearman correlation analysis (P < .05, FDR < 0.2, adjusted for the
season; Fig. 2C). From the significantly affected by the pollution
taxa, our analysis showed that *Planctomycetes* (including genus—
*Pirellula*, *Mariniblastus*, and * Maioricimonas*—and species—*Pirellula
staleyi*, *Mariniblastus fucicola*, and *Lignipirellula cremea*) was signifi-
cantly negatively correlated with heavy metals and PAHs. In addi-
tion, *Nitrospinae* (including genus *Nitrospina* and species *Nitrospina
gracilis*) was significantly positively correlated with PAHs. In addi-
tion, we observed *Ilumatobacter* and more specifically *Ilumatobacter
coccineus* to be significantly negatively correlated with heavy metal
(arsenic and iron) and PAHs pollution indices (Fig. S4).

A functional adaptation of the marine
microbiome to pollution-associated stress

Following quality control, ~2.6 billion filtered reads were used to assemble contigs for each individual sample, resulting in an average of 2.7 billion bp in assembled metagenome (average N50: 261 bp) per sample (Table S6). After removing eukaryotic contigs (EukRep), we predicted genes on constructed prokaryote-enriched metagenome assemblies for functional comparisons. A total of 4646 COG terms, ranging from 4227 to 4415 in each sampling site, were identified and categorized into 63 COG pathways. Shannon diversity index was used to calculate the alpha diversity of COGs. Significant differences between pollution sites were found in alpha diversity based on COG terms (Wilcoxon test, P < .05, Fig. 3A).

Beta diversity calculated by Euclidean distance of COG terms and COG pathways revealed significantly different functional profiles between high and low pollution sites (PERMANOVA test, P < .05, Fig. 3B). We did not observe any significant seasonal differences in alpha or beta diversity in COG terms and COG pathways (Wilcoxon signed-rank test: P > .05, PERMANOVA test: P > .05). The overall difference in COG terms across pollution sites appears to be significantly driven by different microbes (envfit from R package vegan, P < .05, FDR < .05, Fig. 3B), the abundance of which varies significantly depending on the pollution. Phylum *Planctomycetes* and *Proteobacteria* were significantly associated with the functional signature of the microbiome in the low pollution
A

<figure><img src="image_3.png" /></figure>

B

<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

C

<figure><img src="image_8.png" /><figcaption>**Figure 2.** Microbial abundance changes related to pollution. (A) Significantly differentially abundant microbial phyla between low- and high-pollution sites identified by R package *MaAsLin2* (*P* ≤ .05, FDR ≤ 0.2). (B) Significantly differentially abundant microbial species between low- and high-pollution sites identified by R package *MaAsLin2* (*P* ≤ .05, FDR ≤ 0.2). (C) Partial Spearman correlations (adjusted for the season) between pollution-associated microbes and environmental pollution variables. Correlations with raw *P*-value <.05 were labeled with a circle. Correlations with raw *P*-value <.05 and BH *P*-value <.2 were labeled with a triangle.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 3.** Microbial functional adaption to pollution stress. (A) Comparison of Shannon diversity based on COG terms between low- and high-pollution sites. (B) PCoA plot based on COG terms and COG pathways using Euclidean distance. Confidence ellipses: 95% confidence level. Individual lines with an arrow show vector matrices of the significant pollution-related microbial taxa. (C) Significantly changed COG categories between low- and high-pollution sites identified by R package *MaAsLin2* (*P* ≤ .05, FDR ≤ 0.2). B, chromatin structure and dynamics; D, cell cycle control, cell division, chromosome partitioning; F, nucleotide transport and metabolism; I, lipid transport and metabolism; J, translation, ribosomal structure, and biogenesis; K, transcription; L, replication, recombination, and repair; N, cell motility; O, posttranslational modification, protein turnover, chaperones; Q, secondary metabolites biosynthesis, transport and catabolism; R, general function prediction only; S, function unknown; V, defense mechanisms; Z, cytoskeleton.</figcaption></figure>

assemblage, whereas phylum *Nitrospinae* and *Chloroflexi* were associated with high pollution assemblages.

To identify the functional categories of the prokaryotic communities mostly affected by the pollution, we examined the COG categories, COG terms and pathways that showed overall significant differences by R package *MaAsLin2* (*P* ≤ .05, FDR ≤ 0.2). In total 14 COG categories (Fig. 3C) showed significant differences, suggesting a substantial influence of pollution on functional profiles of the microbiome in marine sediments. Significantly more genes related to “Replication, recombination and repair (L)” were found in the high pollution sites as previously shown in heavy-metal contaminated sediments [52]. We also observed significantly more genes related to “Translation, ribosomal structure and biogenesis (J)” in high pollution sites. The relative abundance of 10 COG pathways and 751 COG terms was significantly reduced as the pollution increased, while 29 COG pathways and 1126 COG terms had an opposite trend (Fig. S5 and Table S7). Several metabolic pathways, including fatty acid biosynthesis, biotin biosynthesis, asparagine biosynthesis, and serine biosynthesis were found to decrease at high-polluted sites. Most of the increased pathways are related to COG categories J (translation, ribosomal structure and biogenesis), E (amino acid transport and metabolism), and F (nucleotide transport and metabolism). In highly polluted areas, we discovered an increase in genes related to FOF1-type ATP synthase, which is important for the synthesis of ATP [53]. To assess whether heavy metals specifically influence microbiome function, we examined a set of marker KOs for metal resistance mechanisms [39, 40]. We found that several metal

resistance KO markers were significantly positively associated with their corresponding heavy metals, such as arsenic, cadmium, copper, lead, and zinc (Fig. S5).

## Sessile benthic metazoan structure affected by pollution

Next, we compared microbial and pollution data to COI metabarcoding results of all <2 mm motile organisms and homogenized sample of all sessile scrapings from ARMS deployed for 12 months at the seven sites. After quality filtering and keeping only taxa that are known benthic infaunal metazoans (Table S8) [54–62], in total 571 OTUs were included in the downstream analysis (Table S9). Beta diversity of the benthic metazoans assessed using the Jaccard index revealed significantly different communities across the pollution sites at the species level (PERMANOVA, *P* = .001). To understand how benthic metazoans’ biodiversity is related to environmental pollution stress, we constructed a dbRDA with sediment pollution indices. The assemblages formed by high and low-pollution sites was clearly separated. All covariates together explained 60.0% of the benthic metazoan variance at the species level (dbRDA, *R*² = 60.0%, *P* = .001), which is strikingly similar to the microbiome variance explained by pollution indices. Twelve covariates together explained both significant and nonredundant fractions of variation (Fig. 4A, stepwise dbRDA, *R*² = 58%, *P* = .001). Thirty-one pollution indices significantly explained the variance of the benthic metazoans in univariate analysis (dbRDA, FDR < 0.05), with heavy metals (copper, arsenic, zinc, and lead), PAH pollutants [Benzo(b)fluoranthene and Benzo(k)fluoranthene]
<figure><img src="image_10.png" /><figcaption>**Figure 4.** Benthic metazoans affected by pollution. (A) dbRDA plot based on Jaccard distance. Confidence ellipses: 95% confidence level. The dashed line represents the cut-off for significant nonredundant contribution to the multivariate model. (B) Significantly differentially abundant benthic metazoans between low- and high-pollution sites identified by R package MaAsLin2 (P ≤ .05, FDR ≤ 0.2).</figcaption></figure>

and total sulfide explained the most of variance. Out of all nonredundant covariates, heavy metal pollution indices such as copper, lead, arsenic, and zinc were associated with high-polluted assemblages. Total sulfide, copper, lead, arsenic, pyrene, and total

cyanide were found to account for the nonredundant variance in both the microbiome and benthic metazoans.

We identified significantly enriched benthic metazoans for each season using MaAsLin2 (P ≤ .05, FDR ≤ 0.2) based on linear
mixed-effects models (Fig. S6). More benthic metazoan OTUs were enriched in the summer (35) compared to winter (19), in which OTUs from phyla *Bacillariophyta*, *Gastrotricha*, and *Nematoda* are only enriched in summer. Similar analysis adjusting for season revealed pollution-associated benthic metazoans by R package MaAsLin2 ($P \le .05$, $FDR \le 0.2$), including forty-seven pollution-tolerant and forty-four pollution-sensitive OTUs (Fig. 4B). *Bacillariophyta* exhibited the most pollution-sensitive OTUs. Most pollution-associated OTUs belong to the phyla *Annelida*, *Bacillariophyta*, and *Arthropoda*. Specifically, pollution-tolerant OTUs included *Boccardiella hamata* and *Lumbrineris perkinsi* from *Annelida*, as well as *Alpheus bisincisus* and *Alpheus edwardsii* from *Arthropoda*. In contrast, pollution-sensitive OTUs comprised *Parougia albomaculata* from *Annelida*, *Synalpheus streptodactylus* from *Arthropoda*, and *Mitra luctuosa* from *Mollusca*.

## Microbiome mediates the pollution effect on sessile benthic metazoans

By integrating the data from the microbiome and benthic metazoans, we found that the richness of the microbiome taxonomic and functional profiles was significantly positively correlated with the richness (at the OTU level) of benthic metazoans (Fig. 5A, species: $R=0.54$, $P=.0006$; function: $R=0.35$, $P=.032$). Furthermore, we found statistically positive correlations between the similarity in benthic metazoans and the dissimilarity in bacterial community taxonomic composition and functional composition (Fig. 5B, Mantel test, $P < .01$). To further investigate whether the microbiome community acts as a mediator in the influence of pollution indices on benthic metazoans, we conducted mediation analyses at both the community and individual levels. Mediation analysis is a statistical method that elucidates how an independent variable impacts a dependent variable through an intermediary, known as a mediator. This technique is essential not only for assessing the presence of an effect but also for understanding its underlying mechanism. Notably, mediation analysis with MODIMA showed that microbial communities can mediate the pollution impact on benthic metazoans (Fig. 5C). The panel of microbial species accounted for 36% ($P < .01$) of the total effect of pollution quality on the benthic metazoans, whereas the panel of microbial function accounted for 26% ($P < .01$). We further validated the significant mediation effect (global mediation: $P < .01$) of microbial species and functions with another mediation analysis method LDM-med [43] (Fig. S7).

Subsequently, we performed individual-level mediation analysis based on pollution-associated microbiomes, microbial functions, and benthic metazoans to evaluate the specific microbial species and functions that can mediate the impact of pollution on the benthic metazoans ($FDR < 0.05$). In total, 142 mediation linkages via microbial species and 20 mediation linkages via microbial functions revealed the effects of pollution on pollution-sensitive benthic metazoans (Tables S10 and S11). Among these, microbes from the phylum *Planctomycetes* were the primary mediators (81 out of 142 mediation linkages). Notably, PAHs exert significant effects on pollution-sensitive phyla such as *Annelida* and *Bacillariophyta* primarily through microbial mediation by *Planctomycetes*, which accounted for 60 linkages. Specifically, our mediation analysis suggested that Benzo(k)fluoranthene and Benzo(b)fluoranthene contribute to a decline of *Dodecaceria_sp.* AMW23009 and *Bacillariophyta* OTUsq858, by decreasing *P. staleyi* and *Lignipirellula cremea* levels, respectively (Fig. 5D). Conversely, Nitrospinae, the phylum containing pollution-tolerant microbes, mediated the effects (19 linkages) of pollution on sensitive benthic metazoans including *Annelida* and *Nematoda* (Fig. 5D).

From a mechanistic perspective, we identified several microbial functions to mediate the impact of pollution on pollution-sensitive *Annelida*, mainly including Na + -translocating NADH dehydrogenase, FoF1-type ATP synthase, pyrimidine salvage, and molybdopterin biosynthesis (Table S11). Furthermore, we investigated the potential effects of pollution on microbial species mediated by benthic metazoans (Table S12). In total, we identified 102 mediation linkages, 79 of which were mediated by benthic species from the phylum *Annelida*. The mediation effect of microbes was primarily driven by PAH pollution, accounting for 81 out of 142 linkages, while the mediation effect of benthic metazoans was mainly associated with heavy metal pollution, representing 60 out of 102 linkages.

## Discussion

Recent studies have shown that coastal marine microbiomes are impacted by pollutants resulting from fast urbanization and human activities, such as metal contamination, plastic leachates, and organic pollutants [63–65]. Besides the multiple roles attributed to the environmental microbiome, it serves also as a source for the microbiomes associated with metazoans. For instance, the carapace and gut of fiddler crabs are colonized by different pools of microbial colonists acquired from the environment [66], whereas, seawater and marine sediments were identified as microbial sources for the colonization of fish mucus [67]. However, studies examining holistically and in high taxonomic resolution the micro-organisms and benthic metazoans are required to determine how microbial processes influence the ecology of coastal marine ecosystems under pollution stresses. By using geochemical indicators, deep shotgun metagenomics sequencing and COI metabarcoding sequencing, we show that various pollution levels and sources have an effect on microbial populations and functional output, which in turn affect the sessile benthic metazoans.

Overall, our analysis of the shotgun metagenomics data generated using marine sediments under pollution stress suggests that diverse pollutants resulting from human activities in the coastal marine zone of the rapidly developing city impact benthic metazoans by shaping the composition and diversity of marine sediment microbiomes. Redundancy analysis revealed that pollution indices explained >50% of the variance of microbial communities. Our study's explained variation percentage closely resembles a previous finding that found water environmental stress accounted for ~56% of the variation in the seawater microbiome [68]. Arsenic, total sulfide, ammonia nitrogen, and pyrene are the most influential elements shaping microbial communities. Due to their small size, microbes have a high surface-to-volume ratio and a broad contact area with the environment, which enhances their ability to absorb [69, 70], metabolize, and interact with environmental pollutants [71–74]. At the same time, we observed significant changes in microbial functional pathways under pollutant stress, likely resulting in a unique metabolic landscape that may impact organisms in higher trophic levels. Similar to the redundancy studies based on the microbiome community, >50% of the variance of the benthic metazoans was explained by pollution indices. In addition, we identified both pollution-sensitive and pollution-tolerant microbial species, as well as benthic metazoans, that showed significant changes in abundance in response to increased pollution.

Interestingly, alpha and beta diversity analysis suggested strong positive correlations between the microorganisms and benthic metazoans, whereas mediation analyses revealed a
<figure><img src="image_11.png" /><figcaption>**Figure 5.** Interactions between microbiome and benthic metazoans under pollution stress. (A) Significant positive associations between microbiome community richness, at both species and function level, and the benthic metazoans richness. (B) Significant associations between the benthic metazoans and microbiome community dissimilarity at genus and species level and functional dissimilarity. (C and D) Mediation linkages illustrating the mediator role of microbial species and microbial functions in the associations between pollution and the benthic metazoans at the level of (C) community and (D) specific microbial species.</figcaption></figure>
possible role of the microbiome community between pollution and the benthic metazoans. In particular, from the individual-level mediation analysis of pollution-associated features, we identified 142 and 20 linkages through microbes and microbial functions to mediate the pollution impact on the benthic metazoans. We found that microbes from the *Planctomycetes* phylum are involved in more than half of all mediation linkages, significantly outnumbering those mediated by *Proteobacteria*. This is notable given that *Proteobacteria* is the most abundant phylum and contains the highest number of species significantly related to pollution. *Planctomycetes*, as a known host-associated bacteria, often live in association with other eukaryotic organisms, including giant tiger prawn *Penaeus monodon* [75], Algae, and natural sponges [76]. However, the role of *Planctomycetes* in mediating the effects of pollution on benthic metazoans remains largely unexplored. Our study also suggests that microbes from the *Nitrospinae* phylum may reduce the growth of several pollution-sensitive benthic metazoans from the *Annelida* and *Nematoda* phyla under pollution stress. Microbial mechanisms associated with energy synthesis and conservation, such as Na⁺-translocating NADH dehydrogenase, FoF1-type ATP synthase and pyrimidine salvage [77], play an active role in mediating the impact of pollution on benthic metazoans. Marine microbial communities offer enormous functional resilience [78] and play vital roles in ocean ecology and planet health [79]. Overall, our mediation analyses on chemical pollution indices, microbiome community, and benthic metazoans imply that the compositional and functional changes of the microbiome induced by the pollution might significantly influence the biodiversity and ecological dynamics of benthic metazoans.

However, our study has also limitations. At sites TPC and BI, there is a time difference of up to 2 months between the collection of pollution data and microbial/metazoan samples, which may have introduced minor temporal variability. However, given the persistence of the marine sediment pollution indices [80, 81] we studied, we expect this time difference to have a minimal effect on our results. Additionally, despite a range of pollution indices and environmental variables included in our study, not all possible environmental factors that potentially affect microbiome and benthic metazoans community were considered. A more comprehensive characterization of sediment environmental factors in future research could enable more precise identification of pollution-associated signatures. Our research is primarily correlational and, therefore, does not establish causality, and experimental validation in well-controlled artificial environments is thus needed. An example of an artificial environment is demonstrated in a previous study [82] on benthic annelids in polluted marine sediment. Following the in vitro model used for mediation validation in a prior study, the effects of treatments on both the mediator (microbes) and the responder (benthic communities) should first be evaluated independently. To assess the mediation effect of microbes, future experiments should implement sterilized marine sediment [83, 84] and seawater [85] in artificial environments, both with and without microbial presence to examine the influence of microbes on benthic metazoans responses. Furthermore, while it is widely accepted that the sediment microbiome serves as a source for the microbiome of benthic metazoans, characterization of the microbiome of each individual benthic metazoan was not performed in our study.

In summary, integrating pollution indices, the microbiome, and the sessile benthic metazoans has revealed a potential interaction between the bacterial community and benthic metazoans under pollution stress. Most importantly, through our quantitative

analysis, we discovered communication channels between the benthic metazoans and microbiome communities in marine sediments under pollution stress. Microorganisms account for the vast majority of the biomass in ocean environments and are accountable for most of the major biogeochemical cycles. Our study demonstrates not only how human activities directly affect the microbiome and benthic metazoans, but also how an altered microbiome composition may affect benthic metazoans in coastal marine environments. These findings also offer novel perspectives on the complex ecological dynamics of coastal marine ecosystems in urbanized coastal metropolises and how could be modified for protecting these fragile ecosystems.

## Author contributions

Lin-Lin Xu (Conceptualization, Investigation, Methodology, Writing—original draft, Writing—review & editing), Shelby E. McIlroy (Investigation, Writing—review & editing), Yueqiong Ni (Investigation, Writing—review & editing), Isis Guibert (Investigation, Writing—review & editing), Jiarui Chen (Investigation, Writing—review & editing), Ulisses Rocha (Investigation, Writing—review & editing), David M. Baker (Conceptualization, Resources, Supervision, Writing—review & editing), and Gianni Panagiotou (Conceptualization, Resources, Supervision, Writing—original draft, Writing—review & editing).

## Supplementary material

Supplementary material is available at ISME Communications online.

## Conflicts of interest

The authors declare no competing interests.

## Funding

This work was supported by Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy—EXC 2051—Project ID 390713860 (GP); Environment and Conservation Fund (ECF 2016/67); Research Grants Council of Hong Kong (CRF C7013-19G).

## Data availability

The raw Illumina sequence data of metagenomic data have been deposited in the sequence read archive at NCBI under Bioproject (accession ID: PRJNA1035294). All remaining data generated or analyzed during this study are included in this article and its supplementary information files. The scripts for all statistical analyses are available in the following repository: https://github.com/Lin-LinXu/Sediment-microbiome-mediates-the-effect-of-pollution-on-benthic-macro-organisms/.

## References

1. Halpern BS, Walbridge S, Selkoe KA et al. A global map of human impact on marine ecosystems. *Science* 2008;**319**:948–52. https://doi.org/10.1126/science.1149345

2. Selig ER, Hole DG, Allison EH et al. Mapping global human dependence on marine ecosystems. *Conserv Lett* 2019;**12**:e12617. https://doi.org/10.1111/conl.12617
3. Abdullah MI, Royle LG, Morris AW. Heavy metal concentration in coastal waters. *Nature* 1972;235:158-60. https://doi.org/10.1038/235158a0

4. Zhang M, Sun X, Xu J. Heavy metal pollution in the East China Sea: a review. *Mar Pollut Bull* 2020;159:111473. https://doi.org/10.1016/j.marpolbul.2020.111473

5. Rakib MRJ, Rahman MA, Onyena AP et al. A comprehensive review of heavy metal pollution in the coastal areas of Bangladesh: abundance, bioaccumulation, health implications, and challenges. *Environ Sci Pollut Res* 2022;29:67532-58. https://doi.org/10.1007/s11356-022-22122-9

6. Zhang Y, Yin J, Qv Z et al. Deriving freshwater sediment quality guidelines of polycyclic aromatic hydrocarbons using method of species sensitivity distribution and application for risk assessment. *Water Res* 2022;225:119139. https://doi.org/10.1016/j.watres.2022.119139

7. Shen J, Liao Y, Hopper JL et al. Dependence of cancer risk from environmental exposures on underlying genetic susceptibility: an illustration with polycyclic aromatic hydrocarbons and breast cancer. *Br J Cancer* 2017;116:1229-33. https://doi.org/10.1038/bjc.2017.81

8. Falkowski PG, Barber RT, Smetacek V. Biogeochemical controls and feedbacks on ocean primary production. *Science* 1998;281:200-6. https://doi.org/10.1126/science.281.5374.200

9. Parkes RJ, Cragg B, Roussel E et al. A review of prokaryotic populations and processes in sub-seafloor sediments, including biosphere-geosphere interactions. *Mar Geol* 2014;352:409-25. https://doi.org/10.1016/j.margeo.2014.02.009

10. Baker BJ, Appler KE, Gong X. New microbial biodiversity in marine sediments. *Annu Rev Mar Sci* 2021;13:161-75. https://doi.org/10.1146/annurev-marine-032020-014552

11. Hadfield MG. Biofilms and marine invertebrate larvae: what bacteria produce that larvae use to choose settlement sites. *Annu Rev Mar Sci* 2011;3:453-70. https://doi.org/10.1146/annurev-marine-120709-142753

12. Seeley ME, Song B, Passie R et al. Microplastics affect sedimentary microbial communities and nitrogen cycling. *Nat Commun* 2020;11:2372. https://doi.org/10.1038/s41467-020-16235-3

13. Chen J, McIlroy SE, Archana A et al. A pollution gradient contributes to the taxonomic, functional, and resistome diversity of microbial communities in marine sediments. *Microbiome* 2019;7:104. https://doi.org/10.1186/s40168-019-0714-6

14. Danovaro R, Cocozza di Montanara A, Corinaldesi C et al. Bioaccumulation and biomagnification of heavy metals in marine micro-predators. *Commun Biol* 2023;6:1-12. https://doi.org/10.1038/s42003-023-05539-x

15. Callens M, Watanabe H, Kato Y et al. Microbiota inoculum composition affects holobiont assembly and host growth in Daphnia. *Microbiome* 2018;6:56. https://doi.org/10.1186/s40168-018-0444-1

16. Shimabukuro M, Zeppilli D, Leduc D et al. Intra- and inter-spatial variability of meiofauna in halal trenches is linked to microbial activity and food availability. *Sci Rep* 2022;12:4338. https://doi.org/10.1038/s41598-022-08088-1

17. Simon J-C, Marchesi JR, Mougel C et al. Host-microbiota interactions: from holobiont theory to analysis. *Microbiome* 2019;7:5. https://doi.org/10.1186/s40168-019-0619-4

18. Leray M, Knowlton N. DNA barcoding and metabarcoding of standardized samples reveal patterns of marine benthic diversity. *Proc Natl Acad Sci* 2015;112:2076-81. https://doi.org/10.1073/pnas.1424997112

19. Yang H-S, Kim T, Lee K-T et al. Use of autonomous reef monitoring structures to monitor changes in the marine environment in Jeju, South Korea: a brief review. *Ocean Sci J* 2023;58:17. https://doi.org/10.1007/s12601-023-00108-z

20. McIlroy SE, Guibert I, Archana A et al. Life goes on: spatial heterogeneity promotes biodiversity in an urbanized coastal marine ecosystem. *Glob Change Biol* 2024;30:e17248. https://doi.org/10.1111/gcb.17248

21. Geller J, Meyer C, Parker M et al. Redesign of PCR primers for mitochondrial cytochrome c oxidase subunit I for marine invertebrates and application in all-taxa biotic surveys. *Mol Ecol Resour* 2013;13:851-61. https://doi.org/10.1111/1755-0998.12138

22. Martin M. Cutadapt removes adapter sequences from high-throughput sequencing reads. *EMBnetjournal* 2011;17:10-2. https://doi.org/10.14806/ej.17.1.200

23. Callahan BJ, McMurdie PJ, Rosen MJ et al. DADA2: high-resolution sample inference from Illumina amplicon data. *Nat Methods* 2016;13:581-3. https://doi.org/10.1038/nmeth.3869

24. Ranwez V, Harispe S, Delsuc F et al. MACSE: Multiple Alignment of Coding SEquences accounting for frameshifts and stop codons. *PLoS One* 2011;6:e22594. https://doi.org/10.1371/journal.pone.0022594

25. Rognes T, Flouri T, Nichols B et al. VSEARCH: a versatile open source tool for metagenomics. *PeerJ* 2016;4:e2584. https://doi.org/10.7717/peerj.2584

26. Heller P, Casaletto J, Ruiz G et al. A database of metazoan cytochrome c oxidase subunit I gene sequences derived from GenBank with CO-ARBitrator. *Sci Data* 2018;5:180156. https://doi.org/10.1038/sdata.2018.156

27. Deck J, Gaither MR, Ewing R et al. The Genomic Observatories Metadata Database (GeOMe): a new repository for field and sampling event metadata associated with genetic samples. *PLoS Biol* 2017;15:e2002925. https://doi.org/10.1371/journal.pbio.2002925

28. Bolyen E, Rideout JR, Dillon MR et al. Reproducible, interactive, scalable and extensible microbiome data science using QIIME 2. *Nat Biotechnol* 2019;37:852-7. https://doi.org/10.1038/s41587-019-0209-9

29. Li J, Sung CYJ, Lee N et al. Probiotics modulated gut microbiota suppresses hepatocellular carcinoma growth in mice. *Proc Natl Acad Sci USA* 2016;113:E1306-15. https://doi.org/10.1073/pnas.1518189113

30. Buchfink B, Reuter K, Drost H-G. Sensitive protein alignments at tree-of-life scale using DIAMOND. *Nat Methods* 2021;18:366-8. https://doi.org/10.1038/s41592-021-01101-x

31. Huson DH, Beier S, Flade I et al. MEGAN Community edition—interactive exploration and analysis of large-scale microbiome sequencing data. *PLoS Comput Biol* 2016;12:e1004957. https://doi.org/10.1371/journal.pcbi.1004957

32. Dixon P. VEGAN, a package of R functions for community ecology. *J Veg Sci* 2003;14:927-30. https://doi.org/10.1111/j.1654-1103.2003.tb02228.x

33. Uritskiy GV, DiRuggiero J, Taylor J. MetaWRAP—a flexible pipeline for genome-resolved metagenomic data analysis. *Microbiome* 2018;6:158. https://doi.org/10.1186/s40168-018-0541-1

34. West PT, Probst AJ, Grigoriev IV et al. Genome-reconstruction for eukaryotes from complex natural microbial communities. *Genome Res* 2018;28:569-80. https://doi.org/10.1101/gr.228429.117

35. Zhu W, Lomsadze A, Borodovsky M. Ab initio gene identification in metagenomic sequences. *Nucleic Acids Res* 2010;38:e132. https://doi.org/10.1093/nar/gkq275

36. Galperin MY, Wolf YI, Makarova KS et al. COG database update: focus on microbial diversity, model organisms, and widespread pathogens. *Nucleic Acids Res* 2021;49:D274-81. https://doi.org/10.1093/nar/gkaa1018
37. Kanehisa M, Furumichi M, Tanabe M et al. KEGG: new perspectives on genomes, pathways, diseases and drugs. *Nucleic Acids Res* 2016; **45**:D353–61. https://doi.org/10.1093/nar/gkw1092

38. Xie C, Mao X, Huang J et al. KOBAS 2.0: a web server for annotation and identification of enriched pathways and diseases. *Nucleic Acids Res* 2011; **39**:W316–22. https://doi.org/10.1093/nar/gkr483

39. Ayala-Muñoz D, Burgos WD, Sánchez-España J et al. Metagenomic and metatranscriptomic study of microbial metal resistance in an acidic pit lake. *Microorganisms* 2020; **8**:1350. https://doi.org/10.3390/microorganisms8091350

40. Pal C, Bengtsson-Palme J, Rensing C et al. BacMet: antibacterial biocide and metal resistance genes database. *Nucleic Acids Res* 2014; **42**:D737–43. https://doi.org/10.1093/nar/gkt1252

41. Mallick H, Rahnavard A, McIver LJ et al. Multivariable association discovery in population-scale meta-omics studies. *PLoS Comput Biol* 2021; **17**:e1009442. https://doi.org/10.1371/journal.pcbi.1009442

42. Hamidi B, Wallace K, Alekseyenko AV. MODIMA, a method for multivariate omnibus distance mediation analysis, allows for integration of multivariate exposure–mediator–response relationships. *Genes* 2019; **10**:524. https://doi.org/10.3390/genes10070524

43. Yue Y, Hu Y-J. A new approach to testing mediation of the microbiome at both the community and individual taxon levels. *Bioinformatics* 2022; **38**:3173–80. https://doi.org/10.1093/bioinformatics/btc310

44. Tingley D, Yamamoto T, Hirose K et al. Mediation: R package for causal mediation analysis. *J Stat Softw* 2014; **59**:1–38. https://doi.org/10.18637/jss.v059.i05

45. Zhou W, Yin K, Harrison PJ et al. The influence of late summer typhoons and high river discharge on water quality in Hong Kong waters. *Estuar Coast Shelf Sci* 2012; **111**:35–47. https://doi.org/10.1016/j.ecss.2012.06.004

46. Hua FL, Tsang YF, Chua H. Progress of water pollution control in Hong Kong. *Aquat Ecosyst Health Manag* 2008; **11**:225–9. https://doi.org/10.1080/14634980802100717

47. Mellbye BL, Spieck E, Bottomley PJ et al. Acyl-homoserine lactone production in nitrifying bacteria of the genera Nitrosospira, Nitrobacter, and Nitrospira identified via a survey of putative quorum-sensing genes. *Appl Environ Microbiol* 2017; **83**:e01540–17. https://doi.org/10.1128/AEM.01540-17

48. Epelde L, Lanzén A, Blanco F et al. Adaptation of soil microbial community structure and function to chronic metal contamination at an abandoned Pb-Zn mine. *FEMS Microbiol Ecol* 2015; **91**:1–11. https://doi.org/10.1093/femsec/fiu007

49. Swan BK, Martinez-Garcia M, Preston CM et al. Potential for Chemolithoautotrophy among ubiquitous bacteria lineages in the Dark Ocean. *Science* 2011; **333**:1296–300. https://doi.org/10.1126/science.1203690

50. Arístegui J, Gasol JM, Duarte CM et al. Microbial oceanography of the dark ocean's pelagic realm. *Limnol Oceanogr* 2009; **54**:1501–29. https://doi.org/10.4319/lo.2009.54.5.1501

51. Imhoff JF, Rahn T, Künzel S et al. Photosynthesis is widely distributed among proteobacteria as demonstrated by the phylogeny of PufLM reaction center proteins. *Front Microbiol* 2018; **8**:2679. https://doi.org/10.3389/fmicb.2017.02679

52. Chen Y, Jiang Y, Huang H et al. Long-term and high-concentration heavy-metal contamination strongly influences the microbiome and functional genes in Yellow River sediments. *Sci Total Environ* 2018; **637-638**:1400–12. https://doi.org/10.1016/j.scitotenv.2018.05.109

53. Deckers-Hebestreit G, Altendorf K. THE F0F1-TYPE ATP SYNTHASES OF BACTERIA: structure and function of the F0 complex. *Ann Rev Microbiol* 1996; **50**:791–824. https://doi.org/10.1146/annurev.micro.50.1.791

54. Bonaglia S, Nascimento FJA, Bartoli M et al. Meiofauna increases bacterial denitrification in marine sediments. *Nat Commun* 2014; **5**:5133. https://doi.org/10.1038/ncomms6133

55. Atienza S, Guardiola M, Præbel K et al. DNA metabarcoding of deep-sea sediment communities using COI: community assessment, spatio-temporal patterns and comparison with 18S rDNA. *Diversity 2020*; **12**:123. https://doi.org/10.3390/d12040123

56. Andrew-Priestley M, Newton K, Platell ME et al. Benthic infanual assemblages adjacent to an ocean outfall in Australian marine waters: impact assessment and identification of indicator taxa. *Mar Pollut Bull* 2022; **174**:113229. https://doi.org/10.1016/j.marpolbul.2021.113229

57. Todaro MA, Sibaja-Cordero JA, Segura-Bermúdez OA et al. An introduction to the study of Gastrotricha, with a taxonomic key to families and genera of the group. *Diversity 2019*; **11**:117. https://doi.org/10.3390/d11070117

58. Herranz M, Stiller J, Worsaae K et al. Phylogenomic analyses of mud dragons (Kinorhyncha). *Mol Phylogenet Evol* 2022; **168**:107375. https://doi.org/10.1016/j.ympev.2021.107375

59. Shenkar N, Swalla BJ. Global diversity of Asciidiacea. *PLoS One* 2011; **6**:e20657. https://doi.org/10.1371/journal.pone.0020657

60. Boos K, Gutow L, Mundry R et al. Sediment preference and burrowing behaviour in the sympatric brittlestars *Ophiura albida* Forbes, 1839 and *Ophiura ophiura* (Linnaeus, 1758) (*Ophiuroidea*, Echinodermata). *J Exp Mar Biol Ecol* 2010; **393**:176–81. https://doi.org/10.1016/j.jembe.2010.07.021

61. Uthicke S, Karez R. Sediment patch selectivity in tropical sea cucumbers (Holothurioida: Aspidochirotida) analysed with multiple choice experiments. *J Exp Mar Biol Ecol* 1999; **236**:69–87. https://doi.org/10.1016/S0022-0981(98)00190-7

62. Klunder L, Lavaleye MSS, Filippidi A et al. Impact of an artificial structure on the benthic community composition in the southern North Sea: assessed by a morphological and molecular approach. *ICES J Mar Sci* 2020; **77**:1167–77. https://doi.org/10.1093/icesjms/fsy114

63. Focardi A, Moore LR, Raina J-B et al. Plastic leachates impair picophytoplankton and dramatically reshape the marine microbiome. *Microbiome* 2022; **10**:179. https://doi.org/10.1186/s40168-022-01369-x

64. Di Cesare A, Pjevac P, Eckert E et al. The role of metal contamination in shaping microbial communities in heavily polluted marine sediments. *Environ Pollut* 2020; **265**:114823. https://doi.org/10.1016/j.envpol.2020.114823

65. Martinez-Varela A, Cerro-Gálvez E, Auladell A et al. Bacterial responses to background organic pollutants in the north-east subarctic Pacific Ocean. *Environ Microbiol* 2021; **23**:4532–46. https://doi.org/10.1111/1462-2920.15646

66. Cuellar-Gempeler C, Leibold MA. Key colonist pools and habitat filters mediate the composition of fiddler crab-associated bacterial communities. *Ecology* 2019; **100**:e02628. https://doi.org/10.1002/ecy.2628

67. Minich JJ, Härer A, Vechinski J et al. Host biology, ecology and the environment influence microbial biomass and diversity in 101 marine fish species. *Nat Commun* 2022; **13**:6978. https://doi.org/10.1038/s41467-022-34557-2

68. Glasl B, Bourne DG, Frade PR et al. Microbial indicators of environmental perturbations in coral reef ecosystems. *Microbiome* 2019; **7**:94. https://doi.org/10.1186/s40168-019-0705-7
69. Young KD. Bacterial morphology: why have different shapes? *Curr Opin Microbiol* 2007;**10**:596-600. https://doi.org/10.1016/j.cimb.2007.09.009

70. Broman D, Näf C, Axelman J et al. Significance of bacteria in marine waters for the distribution of hydrophobic organic contaminants. *Environ Sci Technol* 1996;**30**:1238-41. https://doi.org/10.1021/es9504505

71. Páez-Espino D, Tamames J, de Lorenzo V et al. Microbial responses to environmental arsenic. *Biometals* 2009;**22**:117-30. https://doi.org/10.1007/s10534-008-9195-y

72. Francis CA, Beman JM, Kuypers MMM. New processes and players in the nitrogen cycle: the microbial ecology of anaerobic and archaeal ammonia oxidation. *ISME J* 2007;**1**:19-27. https://doi.org/10.1038/ismej.2007.8

73. Zhang L, Yi M, Lu P. Effects of pyrene on the structure and metabolic function of soil microbial communities. *Environ Pollut* 2022;**305**:119301. https://doi.org/10.1016/j.envpol.2022.119301

74. Freitag TE, Klenke T, Krumbein WE et al. Effect of anoxia and high sulphide concentrations on heterotrophic microbial communities in reduced surface sediments (black spots) in sandy intertidal flats of the German Wadden Sea. *FEMS Microbiol Ecol* 2003;**44**:291-301. https://doi.org/10.1016/S0168-6496(03)00076-X

75. Fuerst JA, Gwilliam HG, Lindsay M et al. Isolation and molecular identification of planctomycete bacteria from postlarvae of the giant tiger prawn, Penaeus monodon. *Appl Environ Microbiol* 1997;**63**:254-62. https://doi.org/10.1128/aem.63.1.254-262. 1997

76. Kaboré OD, Godreuil S, Drancourt M. Planctomycetes as host-associated bacteria: a perspective that holds promise for their future isolations, by mimicking their native environmental niches in clinical microbiology laboratories. *Front Cell Infect Microbiol* 2020;**10**:519301. https://doi.org/10.3389/fcimb.2020.519301

77. Nygaard P. Purine and pyrimidine salvage pathways. In: Abrahams S, James H, Richard L (ed.), *Bacillus subtilis and Other Gram-Positive Bacteria: Biochemistry, Physiology, and Molecular Genetics*. Washington, D.C: American Society for Microbiology, 1993, 359-78.

78. Hutchins DA, Fu F. Microorganisms and ocean global change. *Nat Microbiol* 2017;**2**:1-11.

79. Abreu A, Bourgois E, Gristwood A et al. Priorities for ocean microbiome research. *Nat Microbiol* 2022;**7**:937-47. https://doi.org/10.1038/s41564-022-01145-5

80. Fukue M, Nakamura T, Kato Y et al. Degree of pollution for marine sediments. *Eng Geol* 1999;**53**:131-7. https://doi.org/10.1016/S0013-7952(99)00026-5

81. Shin PKS, Lam NWY, Wu RSS et al. Spatio-temporal changes of marine macrobenthic community in sub-tropical waters upon recovery from eutrophication. I. Sediment quality and community structure. *Mar Pollut Bull* 2008;**56**:282-96. https://doi.org/10.1016/j.marpolbul.2007.10.023

82. Ito M, Ito K, Ohta K et al. Evaluation of bioremediation potential of three benthic annelids in organically polluted marine sediment. *Chemosphere* 2016;**163**:392-9. https://doi.org/10.1016/j.chemosphere.2016.08.046

83. Wang B, Teng Y, Xu Y et al. Effect of mixed soil microbiomes on pyrene removal and the response of the soil microorganisms. *Sci Total Environ* 2018;**640-641**:9-17. https://doi.org/10.1016/j.scitotenv.2018.05.290

84. Wu J, Palasantzas V, Andreu-Sánchez S et al. Epigenome-wide association study on the plasma metabolome suggests self-regulation of the glycine and serine pathway through DNA methylation. *Clin Epigenetics* 2024;**16**:104. https://doi.org/10.1186/s13148-024-01718-7

85. Little B, Gerchakov S, Udey L. A method for sterilization of natural seawater. *J Microbiol Methods* 1987;**7**:193-200. https://doi.org/10.1016/0167-7012(87)90040-6