Organic matter and metal loadings influence the spatial gradient of the benthic bacterial community in a temperate estuary

Keywords: 16S rRNA gene, bacteria, bacterial eDNA, environmental monitoring, functional diversity, metals, organic matter, temporal change, temperate estuary, urban.

Estuaries provide valuable ecological, economic and cultural services (Gillanders et al. 2011). These services include, among others, coastal protection, nursery habitats for commercially and ecologically important fish species and cultural, spiritual experiences (Martin et al. 2020). Estuaries are spatially and temporally dynamic systems, and although estuaries vary in size, geomorphology and energy exchange between the land and oceans, many of them have one key thing in common, namely, a growing human population (Kennish 2002). Since 2000, the global human population has increased by 21%, and the expansion of urbanisation along coast lines throughout the world is leading to an array of negative environmental impacts (United Nations 2019). Urbanisation is directly and indirectly affecting all the aforementioned services that estuaries provide by habitat loss and fragmentation, a decrease in biodiversity, harmful algal blooms, over-exploitation of resources, hypoxia, increases in water temperatures, contamination by sewage (both biological and chemical), pesticides, heavy metals and other organic and inorganic pollutants (Diaz and Rosenberg 2008; Greaver et al. 2016; Scanes et al. 2020). Building a baseline understanding of natural variation and seasonal patterns is, therefore, paramount if we, as a society, want to disentangle natural from anthropogenic pressures in urbanised estuaries.

Management strategies, protection and restoration plans are crucial to mitigate current and future impacts on estuaries and coastal zones in general. In Australia, indicators of water-column quality are typically measured and managed at a local (council) and state government level, with little data centralisation at a national level (Raes et al. 2021). To indicate water quality and toxicity, local monitoring programs generally survey abiotic and biotic variables, such as concentrations of nitrate, nitrite and phosphate, turbidity, total suspended and dissolved solids, O₂ concentrations, pH values, and chlorophyll-a (Chl-a) concentrations (Coughanowr et al. 2009; Visby and Weller-Wong 2020). However, there are opportunities in using molecular approaches such as genomic bacteria DNA from environmental samples to improve ecosystem surveillance strategies, monitoring programs and, ultimately, environmental management decisions (Thomsen and Willerslev 2015; Clarke et al. 2020; Raes et al. 2021). Studies utilising metabarcoding and quantitative PCR to investigate sediment and water column bacterial communities (bacterial eDNA) provide both qualitative and quantitative data. For example, bacterial eDNA can be used to investigate bacterial removal and retention of nitrogen (Francis et al. 2013; Raes et al. 2020) to identify the footprint of aquaculture pens (Li et al. 2011; Stoeck et al. 2018), reveal global trends in anti-microbial resistance (Hendriksen et al. 2019) and even the evolutionary responses of bacteria to long-term mercury exposure (Ruuskanen et al. 2020). However, routinely integrating metabarcoding data in environmental monitoring programs remains challenging.

Here, we present the findings from a 2.5-year monitoring survey from sediment samples collected monthly and analysed for genomic bacterial DNA, and physical and biochemical metadata, including sediment total nitrogen (TN), total organic carbon (TOC), analyses of the δ¹⁵N and δ¹³C and metal concentrations in the Derwent estuary in Tasmania, Australia. The Derwent estuary is a microtidal-influenced system connected to the Southern Ocean. Non-biological monitoring programs covering the past 15 years have shown a strong metal pollution gradient in the Derwent estuary, with zinc (Zn), arsenic (As), cadmium (Cd), copper (Cu) and lead (Pb) all exceeding the Australian quality guideline values at sites close to the zinc refinery (Simpson et al. 2007; Coughanowr et al. 2009). Earlier work on sediment bacterial communities in the Derwent estuary showed strong spatial patterns of the bacterial community composition (Abell et al. 2013). More specifically, bacterial nitrifiers, archaeal ammonia oxidisers, denitrifiers and total bacterial and archaeal biomass were positively correlated with C:N ratios and the isotope ratios of N (δ¹⁵N) and C (δ¹³C) of the organic matter (OM) in the sediment.

On the basis of the previous work from Abell et al. (2013) and the results from several studies that have highlighted that bacteria exhibit substrate preferences for the size and age of OM (Findlay 2003; Ding et al. 2015; Wang et al. 2015), we hypothesise that the variability of the bacterial community would be best explained by the source (indicated by the isotopic δ¹³C values; see Peterson (1999)) and quantity of OM. Second, we hypothesise that, although the quantity and source of the OM would be the main explanatory variable for the bacterial community composition on an estuary scale, the sediment metal concentrations would significantly determine the bacterial community composition at the most metal-affected sites near the zinc refinery, but not on an estuary scale. Last, we selected several independent environmental parameters that may have some influence over the bacterial community composition, and thus help determine its variability, including the inferred metabolic potentials of these communities. As a result of this research, we provide a baseline from which we hope to further increase our understanding and the relationships between the sediment geo-chemistry, bacterial assemblages and the ecosystem functions that the latter perform.

We present trends in physico-biochemical water column characteristics at 11 sites in the Derwent estuary, Tasmania, since 2010 (Fig. 1, 2). Monthly water profiles (every metre from the surface to 0.5 m from the bottom) for temperature, salinity and dissolved oxygen were collected (using a YSI EXO 3 multiparameter sonde and a Hydrolab 4a DataSonde) by the following three organisations: The Derwent Estuary Program (DEP), Norske Skog Boyer, and Nyrstar Hobart. Dissolved measurements for inorganic nitrate (NO₃^–), nitrite (NO₂^–) and Chl-a concentrations were collected using Niskin bottles at the surface and bottom by these agencies. Samples were collected 0.5 m from the bottom and the depths ranged from 3 to 26 m (see also Supplementary Table S1). All samples were analysed at Analytical Services Tasmania (AST; https://analyticalservices.tas.gov.au/). All biochemical and physical water column data were requested through the DEP (https://www.derwentestuary.org.au/contact/); the data were received on 20 August 2020. The data and detailed methodologies are available in the Supplementary material and the State of the Derwent Update 2020 (Visby and Weller-Wong 2020); https://www.derwentestuary.org.au/state-of-the-derwent/.

Fig. 1.  Sampling locations in (a) Tasmania and (b) the Derwent estuary. Top-right photo shows the local conditions at the upper-estuary, the middle photo the Nyrstar Hobart zinc smelter (hereafter called zinc refinery) and the bottom photo the open-ocean influenced lower-estuary. Photo credit: Eric J. Raes.

Fig. 2.  Physical and chemical bottom water measurements collected over a decade for 11 sites in the Derwent estuary, Tasmania (Fig. 1a, b), from January 2010 until April 2020. Bottom water (a) temperature (°C; n = 3102), (b) concentrations for dissolved oxygen (DO, mg L^–1; n = 2982), (c) nitrate and nitrate (NO_x, µg L^–1; n = 1998), (d) inorganic phosphorus (PO₄^3–, µg L^–1; n = 1336), (e) ammonium (NH₄⁺, µg L^–1; n = 1336) and (f) zinc concentrations (Zn, µg L^–1; n = 1286). The warmer months are colour coded red (i.e. January, February and March), the colder months blue (i.e. June, July, August) and the transition months yellow (i.e. April, May, September, October, November, December). Note: this temporal colour coding is also used in Fig. 5c and shown in detail in Supplementary Fig. S2.

Monthly sediment surveys were conducted between May 2018 and February 2020 at the same 11 sites where the DEP collected physico-biochemical water-column data. Over the span of 2.5 years, 160 sediment samples were collected for abiotic and genomic bacterial DNA analyses from the lower to the upper basin of the estuary. At each site, three sediment cores were collected using a triangular configured sediment corer with polyethylene sample tubes (4.5-cm internal diameter) located at each corner. The top 1 cm of three sediment cores were homogenously mixed in glass jars, and subsamples were taken for genomic DNA (gDNA) extractions, TN, TOC, including analyses of the δ¹⁵N and δ¹³C, total phosphorus (P) as phosphate from sediment digests (referred to as total PO₄^3– dry weight) and metals, which included As, Cd, Cu, iron (Fe), Pb, Zn and mercury (Hg).

Samples for elemental, organic matter and stable isotope analyses in the sediments were collected between May 2018 and February 2020. Sediment metal concentrations and total phosphorus were analysed by inductively coupled plasma–atomic emission spectroscopy (ICP-AES; Method 2301). Sediment Hg analyses were performed by cold-vapour atomic fluorescence spectroscopy (CV-AFS; Method 2304; please see the ‘Metals and stable isotope analyses’ section in the Supplementary material for more details). Metal concentrations and total P are expressed in milligrams per kilogram on a dry matter basis (DMB). Stable isotope analysis was used to determine the fractionation and potential source of C and N in the estuary. The isotopic composition of the total N and organic C in the sediments (δ¹⁵N and δ¹³C) were analysed at the CSIRO laboratories in Hobart, Tasmania. A Carlo Erba NA1500 CNS analyser was interfaced with a Conflo IV to a Thermo Scientific Delta V Plus isotope ratio mass spectrometer and operated in the continuous-flow mode during sample analyses. Results are presented in the standard δ notation (please see Supplementary material for more details).

Sediment samples for genomics were collected between May 2018 and February 2020 and stored at –80°C until gDNA extraction. Approximately 0.25 g of wet sediment was weighed into powerbead tubes and gDNA was extracted using the QIAGEN DNeasy PowerSoil Total DNA Kits (QIAGEN; Cat. No. 12888-100) according to the manufacturer’s instructions. Nucleic acids were quantified on a QuBIT 2.0 fluorometer (Supplementary Table S1).

Bacterial diversity was investigated by tag sequencing targeting the V1–V3 region of the 16S rRNA gene with the bacterial forward 27F and reverse 519R primer sets (Lane et al. 1985; Weisburg et al. 1991), by using the Illumina MiSeq platform. Amplicons targeting the 16S rRNA gene were amplified from environmental sediment gDNA extracts and sequenced at the Ramaciotti Centre for Genomics (UNSW, Sydney, NSW, Australia). Nextera XT barcode incorporation, purification, library generation and sequencing using the Illumina MiSeq platform (Illumina, Inc., San Diego, CA, USA), with 300-bp paired reads, were performed according to the manufacturer’s directions.

Amplicon sequence variant (ASV) tables were prepared after Bissett et al. (2016) (see also; https://www.australianmicrobiome.com/protocols/16sanalysisworkflow/). Briefly, FLASH (ver. 1.2.11, see https://ccb.jhu.edu/software/FLASH/; Magoč and Salzberg 2011) was used to merge paired-end reads and unique sequences were denoised into ASVs with the UNOISE3 algorithm (Edgar and Flyvbjerg 2015), using USEARCH 64 bit (ver. 8.0.1517, see https://www.drive5.com/usearch/download.html; Edgar 2010). Abundance profiles per sample were constructed by mapping all the reads to the unique ASVs by using the USEARCH ‘otutab’ command. The SILVA database (ver. 138, see https://www.arb-silva.de/; Quast et al. 2012; Yilmaz et al. 2014) was used to derive taxonomy by using the naive Bayesian classifier method of Wang et al. (2007) with the classify.seqs command in Mothur (settings: cutoff = 60, probs = TRUE; Schloss et al. 2009). ASVs were classified as follows: Domains: ‘unknown’, ‘Archaea’ and ‘Eukaryota’; Phylum: ‘Bacteria_unclassified’; Family: ‘Mitochondria’; Order: ‘Chloroplast’, and likely erroneous sequences were removed before the analyses. The final 16S rRNA gene library had 13 803 ASVs and the sequencing depth ranged between 4121 and 40 253 reads per sample. Genomic data are available at the NCBI under the Bioproject accession number PRJNA611638.

PICRUSt2 (Langille et al. 2013; Douglas et al. 2020) was used to infer an approximate phylogenetically inferred metabolic potential of the bacterial communities in the Derwent estuary. To reduce run time, and because PICRUSt2 estimates function from the nearest ‘ancestor’, we clustered the ASVs at the 97% similarity threshold (97% of operational taxonomic units (OTUs) were generated using the USEARCH ‘cluster_fast’ function with -id 0.97; see Supplementary material for the full workflow). The PICRUSt2 pipeline (v. picrust/2.3.0b) was run with default settings, and sequences with Nearest Sequenced Taxon Index (NSTI) scores of >2 were removed. The average NSTI score, based on 163 samples, was 0.52 ± 0.28 (±s.d.). Only 3 of the 13 803 OTUs were above the maximum NSTI cut-off of 2 (these ASVs had less than 31 reads). Pathways with fewer than 10 reads were removed. MetaCyc pathway abundances are the main high-level prediction output by PICRUSt2. MetaCyc is a curated database (Caspi et al. 2016) of experimentally deduced metabolic pathways from all domains of life and an open-source alternative to the Kyoto Encyclopaedia of Genes and Genomes (KEGG; Kanehisa et al. 2007). A superclass or parental class contains a set of pathways that accomplish roughly the same biological function, such as degradation of a given starting material, or biosynthesis of an end product.

The Phyloseq package (ver. 1.32.0, see https://joey711.github.io/phyloseq/; McMurdie and Holmes 2013) was used to analyse, visualise and plot the microbiome and physical and bio-chemical metadata using R (ver. 4.0.2, R Foundation for Statistical Computing, Vienna, Austria) in RStudio (ver. 1.1.442, RStudio Inc., Boston, MA, USA, see https://www.rstudio.com/). Statistical tests were conducted using the Vegan package (ver. 2.5-6, J. Oksanen, R. Kindt, P. Legendre, B. O’Hara, M. H. H. Stevens, M. J. Oksanen and M. Suggests, see https://cran.r-project.org/web/packages/vegan/). Analyses of similarities (ANOSIM) were used to test whether we could identify statistical differences among the 11 estuarine sites and the 4 estuarine zones using the PRIMER-e (ver. 7.0.17, see https://www.primer-e.com/; Clarke and Gorley 2006). The plyr package (ver. 1.8.4, see https://cran.r-project.org/web/packages/plyr/index.html) was used to calculate means and summarise the data (Wickham 2011).

Differential abundance analyses were used to identify changes in taxa and functions between the four different estuarine zones, by using ANOVA-like differential expression (ALDEx2 ver. 1.20.0, see https://bioconductor.org/packages/release/bioc/html/ALDEx2.html; Gloor et al. 2017). Prior to the analyses, ASVs were agglomerated at the family level (on average and across all sites <5% of the ASVs at the phyla were assigned as ‘uncultured’). To account for the compositionality of the data, both the ASV count data and the inferred metabolic pathway (PICRUSt2) data were transformed using the centre log ratio (CLR) for all statistical analyses. The CLR transformations were performed using the microbiome package (ver. 1.16.0, see https://microbiome.github.io/) with a pseudo-count of min(relative abundance) ÷ 2 to exact zero relative abundance entries. The differential abundances were generated by 128 Monte Carlo samples sourced from a Dirichlet distribution (Fernandes et al. 2014).

Unconstrained (principal component analyses; PCA) and constrained (redundancy analyses; RDA) were constructed from Aitchison (Euclidean) distance matrices generated from CLR-transformed ASVs. Environmental parameters were standardised with the ‘standardise’ function (variables were scaled to zero mean and unit variance) using decostand from the Vegan package (J. Oksanen et al., see https://cran.r-project.org/web/packages/vegan/) and significant (P < 0.05) environmental parameters were derived using the ‘envfit’ function in Vegan and overlaid as vectors to identify multiple explanatory variables among the estuarine zones (J. Oksanen et al., see https://cran.r-project.org/web/packages/vegan/).

During the past decade, annual bottom water temperatures ranged between 6.5 and 21.3°C. Bottom water temperatures were significantly different among the winter, transition and summer months (Fig. 2a; Wilcoxon test, P < 0.05). Dissolved oxygen (DO) concentrations displayed similar significant changes among the winter, transition and summer months, with highest DO concentrations being recorded in the colder months (>7.5 mg L^–1), and significantly lower concentrations in summer (Wilcoxon test, P < 0.05; Fig. 2b). Dissolved inorganic NO₃^– + NO₂^– (NO_x) concentrations showed a clear annual trend, with significant differences among the winter, transition and summer months, with highest concentrations in the colder and lowest concentrations in the warmer months (Wilcoxon test, P < 0.05; Fig. 2c). Dissolved inorganic PO₄^3–, NH₄⁺ and Zn concentrations showed a general increase from the lower to the upper estuary (Fig. 2d–f). Dissolved inorganic PO₄^3– was significantly higher in winter, whereas NH₄⁺ concentrations were significantly higher in the transition months (Fig. 2e; Wilcoxon test, P < 0.05). Bottom water Zn concentrations were significantly higher in the warmer months (Wilcoxon test, P < 0.05). Salinity measurements from the bottom waters showed an increasing trend from fresh water at the upper estuary (i.e. site U12), up to 32 PSU in the lower estuary (i.e. B1 and B3 sites; Supplementary Fig. S2). Across the estuary, all bottom water measurements showed clear seasonal cycles across the entire estuary (Supplementary Fig. S2a–f).

Concentrations of TOC and TN in the first centimetres of the sediments showed similar, and significant, trends across the estuary, with the lowest concentrations occurring at the lower estuarine sites (i.e. B1, B3 and RBN sites; ~5 µg C g^–1 and ~0.5 µg N g^–1), an increase at the mid-estuary sites (i.e. E, G2 and KB sites; ~50 µg C g^–1 and ~4 µg N g^–1), and a declining trend for the upper estuary sites (i.e. U7 and U12 sites; ~20 µg C g^–1 and ~1 µg N g^–1; Fig. 1b, 3a, b; Wilcoxon test, P < 0.05). The isotopic composition of the TOC (δ¹³C) showed a significant trend to more negative values from the lower to the upper estuary (Wilcoxon test, P < 0.05). The highest variability for δ¹⁵N was measured at the sites with the lowest TOC and TN concentrations (Fig. 3c, d). Besides at Site U12 (where seasonal differences were noted for TOC, TN and δ¹⁵N between the transition and warm months; Wilcoxon tests, P > 0.05), no seasonal effect was noted at the other sites for the TOC, TN and isotopic measurements (Wilcoxon tests, P > 0.05). At most sites, metal (As, Cd, Cu, Hg, Pb and Zn) concentrations in the sediments exceeded the Australian and New Zealand upper sediment quality guideline values (SQGV), with the highest concentrations at the sites near the zinc refinery (Sites U2 and NTB5; Supplementary Table S2). Overall highest metal concentrations were always recorded at Site NTB5 (Fig. 4, Supplementary Table S2). Similar trends along the estuary were noted for all other metal concentrations (Supplementary Table S2), and no seasonal differences were found at each site (Wilcoxon tests, P > 0.05).

Fig. 3.  Organic sediment data in the Derwent estuary, Tasmania. Site numbers present a gradient upstream from the lower to the upper estuarine sites. (a) Total organic carbon concentrations (TOC, µg g^–1); (b) total nitrogen concentrations (TN, µg g^–1); (c) the isotopic values for the TOC (δ¹³C); and (d) the isotopic values for the TN (δ¹⁵N). Sample numbers for all sites are shown on top of d. TOC, TON and stable isotope samples were collected between May 2018 and February 2020.

An unconstrained PCA analysis showed that, similar to the bio-geochemical sediment data, the bacterial communities displayed a clear spatial pattern along the Derwent estuary, with sites at the lower estuary, sites in the mid-estuary, sites close to the zinc refinery and sites in the upper estuary clustering closer to each other respectively (Fig. 5a). Analysis of similarities (ANOSIM) complemented these findings and showed a high and significant separation among the 11 sites and the different estuary zones (ANOSIM statistic R: 0.893 and 0.840 respectively; significance: 0.001; pairwise combinations are shown in Supplementary Tables S3, S4).

Redundancy analyses (RDAs) showed that a variety of sediment and bottom water parameters were significantly correlated with the bacterial community composition (Supplementary Table S5). More specifically, regression fitting of single explanatory variables including, TN, δ¹³C, TOC, Fe, bottom water salinity, dissolved organic C, sediment PO₄^3–, Hg, the C:N ratio of the organic matter in the sediment, Pb concentrations and changes in the δ¹⁵N each explained more than 20% of the variance in the prokaryotic community composition (R²: 0.65, 0.60, 0.58, 0.58, 0.47, 0.41, 0.31, 0.29, 0.28, 0.27, 0.26 respectively, P < 0.001 for all; Supplementary Table S5). Bottom water dissolved O₂, NH₄⁺, Cu, PO₄^3–, As, Cd and surface Chl-a concentrations also showed significant correlations with the bacterial community composition, but each parameter contributed <20% of the variance (R²: 0.19, 0.19, 0.19, 0.12, 0.10, 0.06 and 0.12 respectively, P < 0.001 for all; Supplementary Table S5). All metals showed significant correlations with the bacterial community composition, but Fe (mg kg^–1 DMB, which showed a trend similar to the TOC and TN concentrations; Fig. 4c, d, 5b) showed the strongest effect on the community composition (R²: 0.58 Supplementary Table S5). A corresponding redundancy analysis constrained by seasonality showed that a very low percentage of the variance (~1%) could be explained by seasonal changes; Fig. 5c). Most of the measured sediment and bottom water measurements significantly (P < 0.001) correlated with the bacterial community data, except for bottom water temperatures, bottom water TN, bottom water NO_x and sediment Cd concentrations (Supplementary Table S5).

Fig. 4.  Sediment metal concentrations in milligrams per kilogram on a dry matter basis (DMB) for (a) zinc (Zn), (b) mercury (Hg), (c) iron and (d) phosphorus in the Derwent estuary. Red lines indicate the upper quality guideline values (SQGV; see Supplementary Table S2). Note that there are no SQGVs for Fe nor for PO₄^3–. Sample numbers for all sites are shown on top of d. Sediment metal samples were collected between May 2018 and February 2020.

Fig. 5.  Changes in bacterial β diversity in the Derwent estuary, Tasmania. Colours show sites and symbols estuary zones. (a) Principal coordinates analysis (PCA) of the bacterial community; data were CLR transformed with an Aitchison distance matrix. (b) Redundancy analyses (RDA) to characterise the specific abiotic variables in the sediment which exerted the largest influence on the bacterial community composition and RDA, and (c) constrained by warmer (January, February and March), colder (June, July, August) and transition months (the other months; see Fig. 2a).

Fig. 6.  Differential relative abundance analyses at a phylum level between the lower and the other estuary zones. All estuary zones were significantly different from each other; here we show only the comparisons between the lower and other estuary zones. The lower estuary zone is the least metal affected zone and can be seen as a ‘pristine’ zone compared with the others. (a) The lower v. the mid-estuary, (b) the lower v. zinc refinery, and (c) the lower v. upper estuary sediment communities. The plots show ASVs binned in their respective phyla for pairwise comparisons. The log₂ standardised difference (the effect size) is coloured red (effect size of ≥1) or blue (effect size of ≤1) for the comparisons. Grey points are ASVs that did not show an effect >1 or <1.

No significant differences were found within each site between the sediment bacterial communities in the warmer, transition and colder months (PERMANOVA P > 0.05). These results confirmed the strong differences in the bacterial communities between the sites and estuary zones (see ANOSIM analyses; Supplementary Tables S3, S4), along with the low percentage of the variance (1.6%) explained by seasonal changes in bottom water (RDA results; Fig. 5c). Differential relative abundance analyses were then used to identify the main phyla that contributed to the variations in the sediment communities between the estuary zones (as observed in the PCA analysis, Fig. 5). Results showed that the lower estuarine sediments were enriched with the phyla Proteobacteria, Nitrospinota, Gemmatimonadota and Desulfobacterota when compared with mid-estuarine samples. The zinc refinery sites were shown to have higher relative abundances of the phyla Zixibacteria, SAR324 clade, Proteobacteria, Planctomycetota and Myxococcota than in sediments in the lower estuary. The upper sediments could mostly be characterised by an increase in the phyla Zixibacteria, Verrucomicrobiota, SAR324 clade, Myxococcota and Firbobacterota, compared with the lower estuary (Fig. 6a, c).

The 16S rRNA gene-based metabolic reconstructions resulted in 369 inferred MetaCyc pathways and could be collapsed into 40 superclasses. Once again, site was a highly explanatory variable for bacterial functional community composition (ANOSIM; R: 0.430; significance: 0.001; with significant differences for all pairwise combinations except between B1–B3 and U7–U4 sites; see Supplementary Table S6). Estuarine zone showed a statistically significant effect albeit with a lower explanatory R value than for the site comparisons (ANOSIM; R: 0.406; significance: 0.001; see Supplementary Table S7 for pairwise combinations and Fig. 7a). RDAs showed that again a variety of sediment and bottom water parameters were significantly correlated with the phylogenetically inferred functional community composition, but only three parameters (δ¹³C, TOC and salinity in bottom waters) contributed >20% of the variance (Fig. 7a, Supplementary Table S8). Differential relative abundance analysis identified a higher number of metabolic superclasses related to biosynthesis and methanogenesis pathways in the lower estuary. The sediments near the zinc refinery were associated with a higher relative abundance or more diverse and complex degradation pathways and siderophore biosynthesis pathways (Fig. 7b).

Fig. 7.  (a) RDA of the phylogenetically inferred functional community with the parameters that explained >20% of the variance. Data were CLR transformed with an Aitchison distance matrix. (b) Differential abundance analysis of the phylogenetically inferred functional sediment communities between the lower estuary and zinc refinery sites. The log₂ standardised difference (the effect size) is coloured red (effect size of ≥1) or blue (effect size of ≤1) for the comparisons. Grey points are Metacyc pathways which did not show an effect >1 or <1.

Our study complements the previous results from Abell et al. (2013) in the Derwent estuary, who found that both the abundance and community structure of bacterial nitrifiers, archaeal ammonia oxidisers, and denitrifiers are largely influenced by the sediment organic matter composition. However, the results from the current study extend these findings from those bacteria involved in the N-cycle to the whole bacterial community. Abell et al. (2013) already noted that the differences in OM sources are potentially affecting the rate of mineralisation by heterotrophic microorganisms. The trend to more negative δ¹³C values of the OM, from the lower to the upper estuary, indicates an increase of more terrestrial C sources (Stevens et al. 2021). By contrast, the low OM concentrations (both TN and TOC), the low C:N ratios, and the less negative δ¹³C values in the sediments of the lower estuary indicate a greater contribution of marine phytoplankton (Peterson 1999). The study from Abell et al. (2013) suggested that the OM in the lower estuary is more labile and, hence, more easily remineralised to NH₄⁺, which was indirectly confirmed by observing a higher biomass of heterotrophic organisms (the authors reported higher 16S rRNA gene copy numbers for both bacteria and archaea in the lower estuary). In this study, we have also provided evidence that sites near the zinc refinery, with higher organic matter quantity, had a greater diversity of degradation pathways, including the presence of more complex compound degradation pathways. These results suggest that although the organic matter is less labile (from a C:N ratio perspective) in the mid-estuary and zinc refinery sites, the larger abundance of different organic substrates also seems to provide niche diversification for the bacterial community.

Although metal concentrations in the sediments exceeded the Australian and New Zealand upper sediment quality guideline values at most sites, and always at the sites near the zinc refinery, they did not show a large effect on the bacterial community composition on an estuary-wide basis (Supplementary Table S5). These results suggest that metal-tolerant species have adapted to the high metal concentrations in the estuary, and highlight the need to have a baseline from which we can untangle the natural variation owing to past and future pressures. These findings are also supported by the recent study from Hillyer et al. (2022) who showed that changes in the metabolite profiles for three macroinvertebrate groups were not significantly correlated with Zn concentrations on an estuary basis in the Derwent. However, zinc sediment concentrations were significantly and positively correlated with changes in metabolite profiles at the most metal-affected sites near the zinc refinery. These results, including those from our study, confirmed our second hypothesis that metal concentrations will influence site-specific spatial variation but not on a whole of estuary scale, because metals did not exceed the upper sediment quality guideline values at every site. In general, these results also align with those presented by Chen et al. (2018), who showed that the microbiome is able to adapt to long-term heavy-metal pollution exposure, and that sediment bacterial diversity, therefore, can be best examined among similar affected sites and not on an estuary-wide scale.

Water column data are the most commonly collected metadata, and we believed it would be a reasonable proxy to explain variability in the sediment bacterial communities. However, our study, and the earlier results from Abell et al. (2013) showed that the bacterial communities had a high site-specific community structure with distinct communities over time. These results are interesting, given the significant changes in environmental conditions in the estuary between seasons (e.g. temperature differences of 8°C, 3-fold increases for NO_x and NH₄⁺ and doubling of PO₄^3– concentrations in the bottom waters between summer and winter months). Neither bottom water temperature or dissolved NO_x concentrations showed significant correlations with the sediment microbiomes, suggesting that bacterial groups are more tightly controlled by local sediment abiotic and biotic parameters rather than large fluctuations in the overlying bottom water. The disconnect between the bottom water and the sediment has important ramifications, because it has been suggested that dynamic water column changes, including nutrient bursts, may have little effect on the compositional dynamics of sediment microbes (Orland et al. 2020). These findings suggest that rather than an indication of current state, perhaps our data give some indication of trajectory or temporal stability. Sediment microbes may instead respond stronger to fine‐scale variations in the composition (quantity, quality and source, e.g. marine and terrestrial) of organic matter (Fagervold et al. 2014; Ding et al. 2015; Wang et al. 2015), and, consequently, the physicochemical gradients and geochemistry of their surroundings. These findings are consistent with observational evidence that the abundance rather than presence of bacterial taxa is limited by resource availability in sediments (Louca et al. 2018; Orland et al. 2019). Our study, therefore, supports the necessity to sample multiple sites and consider sediments differently from the water column when bacterial processes are upscaled to entire estuary basins (Orland et al. 2020).

The cumulative effects of climate change and increasing pressures from urbanisation are affecting coastal ecosystems worldwide (United Nations 2019). Over the past 12 years in Australia, estuarine temperatures have increased on average by 2.16°C, with waters acidifying at a rate of 0.09 pH units and freshening at 0.086 PSU year^–1 (Scanes et al. 2020). The benefits of protecting and managing estuaries is directly related to the fact that these systems provide valuable ecological, cultural and economic services, including aquaculture industry, nursery grounds for commercial and wild fisheries, N removal, C storage, and recreational tourism (Evans et al. 2017; Fulford et al. 2020). Identifying reliable and inexpensive methods to track the status of estuaries and to predict the effects of specific stressors or management remains a major challenge (United Nations 2019). The scientific community are increasingly incorporating genomic surveys into their methodologies for biomonitoring because of its high accuracy, taxonomically holistic lens and ease of deployment (Berry et al. 2021). Although the integration of genomic data comes with its own challenges (Porter and Hajibabaei 2018; Quinn et al. 2018), it is becoming a core technology for ecosystem and water-quality assessment (Glasl et al. 2017; Berry et al. 2021). However, the potential to include genomic surveys into environmental monitoring programs can be achieved only once a baseline (the aim of this study) is established; from there, change can be estimated and quantified.

Omics-based monitoring technologies using bacterial marker genes have been shown to reveal valuable mechanistic insights into the functioning of ecosystems and the metabolic potential of bio-geochemical cycles (Ruuskanen et al. 2020). The underlying reason is that microbes are seen as useful, reliable and responsive indicators of ecosystem health because of their high abundance, the diversity of nutrient cycling functions they perform, and because their short life-cycle makes them able to respond quickly to environmental change. Integrating bacterial eDNA data in environmental monitoring programs remains a challenge. Our study did provide insights into the inferred metabolic potential of the sediment bacterial community, by showing relatively larger proportions of pathways associated with the biosynthesis of organic molecules in the sites from the lower estuary (referred to as the more ‘pristine’ sites compared with those close to the zinc refinery in terms of metal concentrations). These results also highlighted significant differences in the metabolic potential between these pristine and metal-polluted sites and open the possibility to get insights into the functional bacterial diversity. The opportunity to infer and measure this functional bacterial diversity, or metabolic potential, may provide new meeting points for microbial ecologists, biogeochemists and ecosystem modellers. However, one major gap in many omics surveys including ours is the lack of rate or metabolite data to truly allow extrapolation of marker gene predictions to metabolic processes. The congruent data collection for marker gene data, metabolic rate data (in situ or from experiments) and metabolomic data (e.g. Hillyer et al. 2022) are highly needed to fully unravel the interactions among different bacterial taxa within a community and their influence on the environment. These insights will be key to appreciate the pivotal role microbes have in an ecosystem.

How ecosystem functioning changes with bacterial communities remains an open question in many natural ecosystems, and merits investigation; hence, the main purpose of this study was to provide a descriptive analysis of the sediment bacterial communities of the Derwent estuary. Overall, our study supports evidence that small-scale environmental conditions (such as organic sediment composition) have a larger influence over the bacterial community assemblages in the Derwent estuary than do large‐scale environmental conditions (such as seasonal bottom water changes). We hope that these data will increase our understanding and ability to predict or infer the metabolic state of coastal ecosystems in a holistic way, and one day may improve environmental management decisions.

All biochemical and physical water-column data are available through the DEP (https://www.derwentestuary.org.au/contact/). Genomic data are available at the NCBI under the Bioproject accession number PRJNA611638.

The authors declare that they have no conflicts of interest.

This research was funded by CSIRO Environomics Future Science Platform.

E. J. Raes, K. Karsh, B. H. Holmes, K. E. Hillyer, M. Green and A. Bissett collected monthly sediment data for genomic DNA, metal and organic matter analyses. A. T. Revill and E. A. Brewer analysed and quality controlled the organic matter analyses. S. Whitehead, B. Proemse, U. Taylor and A. Weller-Wong collected, analysed and provided long-term metadata. L. Bodrossy designed the study. A. Bissett, J. van de Kamp and E. J. Raes analysed the sequence data. E. J. Raes provided ecological interpretation and wrote the manuscript. All authors contributed to the editing of the paper.

Supplementary material is available online.