The sampled estuaries were assessed against the Australian and New Zealand Guidelines for Fresh and Marine Water Quality in sediments , using the collected physicochemical data (See supplementary information). Overall, the four estuaries sampled are considered healthy and below the trigger values for concern stated in the guidelines for tropical estuaries in Australia. This supports the most recent aquatic ecosystem report card published by the Queensland Government that stated the region as being “good quality” in terms of water quality, pollution, and introduced flora and fauna species.
It is noted that some water quality and sediment parameters were at, or above, the guideline trigger values, including chlorophyll-a, total nitrogen (TN), nitric oxide (NOx), total phosphorus (TP), and turbidity (Supplementary Tables S1-S2). Furthermore, measured concentrations of iron and aluminium were also high compared to the marine sediment guideline trigger values across the four estuaries (Supplementary Table S1). However, this was to be expected based on the geological composition of the surrounding soil/substrate and, for the region, is considered naturally occurring. It should be noted that bauxite mining does occur out of the Embley-Hey estuary. Again, the concentrations of iron and aluminum were observed consistently high in all estuaries irrespective of localized mining activities. Sediment particle size analysis was also consistent for the region (Supplementary Figure S2 and Table S3). The distinct wet-dry seasonality of the semi-arid Cape York estuaries has been shown to trap and redistribute sediments over tidal and seasonal cycles, leading to an extreme, but highly variable, turbidity [50, 51]. As such, turbidity above 20 NTU for all estuary sites was expected and may account for elevated chlorophyll-a and primary nutrients. The analyzed sediments were free from any pesticide residues screened as part of the Agilent Pesticide PCDL (data not shown).
3.1 Community metabolomics analysis
In total, 1681 discovery polar metabolite features were detected, of which 727 metabolite features were identified in the sediment samples collected from across four estuaries. Of all the identified metabolites, 53 central carbon metabolism metabolites and 674 discovery polar metabolites were identified against the METLIN metabolites and Agilent PCDL. Supplementary Figure S3 presents the PCA plots for the metabolomics data and provides a preliminary overview of sample/site clustering. The results indicate no clear discrimination between estuaries. The PCA model was found to be significant, with good linearity R2X (cum) = 0.639 and predictability Q2 (cum) = 0.484. Outliers were examined using a distance of observation (DModX) analysis ; no metabolites were considered outliers at the two times DCrit threshold (p = 0.05) (Supplementary Figure S3). To explore potential for metabolic pathway activity levels in each estuarine system, the metabolomic data were subjected to pathway enrichment analysis. Pathway enrichment was used to decipher the metabolic pathways to which the metabolites are associated and resulted in 445 identified metabolites being mapped to the KEGG database. Pathway impact measured the relative significance of all metabolites of a single pathway for the overall spectrum of metabolites. Supplementary Figure S3 summarises the pathways mapped. We then performed a ChemRICH class annotation of the KEGG identified metabolites and found that 70 chemical clusters were represented. These chemical groups can be associated with genes and phenotypes. The analysis highlighted saturated and unsaturated fatty acids (FA), glycosides, aminoglycosides, deoxyadenosines, glutamates, indole alkaloids, ketones, monoterpenes, pregnanediones, and purine nucleosides as the most abundant chemical clusters (Supplementary Figure S4). Figure 2 illustrates the chemical clusters of significance (p < 0.05) contributing to metabolic differences between the estuaries (Supplementary Figure S4).
Metabolites that differentiated individual estuaries were harmine (KEGG ID C06538; an indole alkaloid found in fruit and seeds , which was less abundant in the Skardon estuary), deoxyadenosine (C00559, a deoxyadenosines linked to marine microbial natural products , which was elevated in Skardon), deoxyguanosine (C00330; a purine nucleosides found in plants, and elevated in Skardon), undecylenic acid (C13910; an unsaturated FA elevated in Skardon that is a known antifungal drug and eukaryotic metabolite produced in plants), and a glycerol 1-phosphate derivative (C04590; a ketone found in marine sediments , and less abundant in Skardon). The Wenlock samples were found to have decreased levels of biphenyl (C06588; a deoxyadenosine found in some plants and animal tissues). Other significant key chemicals were fumigaclavine C (C20438; an indole alkaloid derived from marine-based fungus ), itaconic acid (C00490; a succinate), and L-serine (C00065; an unsaturated FA), all of which were slightly decreased. Deoxyloganin (C06071; a monoterpene) and oxoglutaric acid (C00026; a glutarate) were lower in the Embley-Hey catchment, likewise 4-hydroxy-L-glutamic acid (C0449; a glutamate) was lower in the Archer-Watson estuary. It should be noted that all the significant metabolites here were identified based on p-values alone and no metabolites were observed with changes greater than 1.5-fold across all the estuary sites. This is indicative of ecosystems that are considered similar in terms of their metabolic function at the time of sampling. The specific pathways that are the same across all the estuaries are discussed in more detail below and form the basis for establishing microbial blueprints of the sampled estuaries.
3.2 Bacterial community analysis
Bacterial communities (with relative abundance ≥ 1%) from the four estuaries were taxonomically diverse, comprising 15 bacterial phyla. A total of 16,066 phylotypes (97% OTUs) across all four estuaries were obtained. The four estuaries shared 3,582 common phylotypes (Supplementary Figure S5).
Taxonomic profiling demonstrated that the most abundant phyla (70.5 – 71.6% relative abundance) were, Proteobacteria (mostly classes g-proteobacteria and a-proteobacteria; Embley-Hey > Skardon > Wenlock > Archer-Watson), Bacteroidota (mostly class Bacteroidia; Archer-Watson > Wenlock > Skardon > Embley-Hey), Chloroflexi (mostly class Anaerolineae; Wenlock > Archer-Watson > Skardon > Embley-Hey) and Desulfobacterota (mostly classes Desulfobulbia and Desulfobacteria; Skardon > Embley-Hey > Archer-Watson > Wenlock) (Supplementary Figure S5). Members of Planctomycetota (4.5 – 6.1%) and Acidobacteriota (3.9 – 4.5%) were also found to have high relative abundances in all four estuaries. While the four estuaries share a largely similar community composition at the phylum level, a distinctive yellow microbial mat on the intertidal area in the upper and mid regions of the Wenlock and the Archer-Watson estuaries was observed. This is most likely due to an increased relative abundance of cyanobacteria in these estuaries. It is worth noting that elevated turbidity was found in these estuaries (Supplementary Table S1), which might offer some competitive advantage to benthic cyanobacteria that are exposed to direct sunlight at low tide . The presence of cyanobacteria, albeit in low relative abundances (2.3% in Wenlock and 3.5% in Archer-Watson), could explain the higher chlorophyll-a content in these samples.
A significant difference in OTU diversity indices (Shannon’s diversity (H’), Simpson’s diversity (D1), and Simpson’s dominance (D2); Supplementary Table S4) was detected among the four estuaries (Kruskal–Wallis; P < 0.05) (Supplementary Table S4). A significant difference in H’, D1 and D2, were observed between Archer-Watson and Skardon estuaries (Dunn’s Test, P < 0.05). It should be noted that these estuaries are the most distant geographically and may also be the most different in terms of watershed characteristics and wind/tidal forcing. This was also illustrated by lower beta-diversity indices (CJ and Cs, Supplementary Table S4).
3.3 Overall functional gene diversity and structure of microbial communities
The predicted functional profiling of the estuarine microbiome datasets showed the presence of 7,303 KEGG IDs in the four estuaries. The most active pathways are indicated in Supplementary Figure S6. We observed a highly similar functional repertoire of the estuarine microbiome present across the four estuaries, despite the differences in their overall taxonomic profile composition. This was also illustrated by the OPLS-DA statistics of the metagenome dataset (Supplementary Figure S7).
We did however observe a significant difference in predicted gene diversity indices (H’, D1 and D2) of microbial communities among the four estuaries (Kruskal–Wallis; P < 0.05) (Supplementary Table S5). We also observed a significant difference between Archer-Watson and each of the other three estuaries for Shannon’s diversity index (H’, Dunn’s test, P < 0.05). A low Jaccard’s and Sorensen’s index also indicated that there was only a small dis-similarity between these estuaries (Supplementary Table S5).
3.4 Carbon metabolism and biosynthesis of amino acids in estuarine sediments
As a result of the significant similarities observed between these geographically distinct estuaries, the estuaries were further analyzed in terms of characterizing the microbial blueprints important for microbially-mediated processes. As such, all identified metabolites were integrated with 317 functional genes predicted from the 16S rRNA gene sequencing data that are involved in carbon, nitrogen, and sulphur cycles, amino acid metabolism and metal homeostasis and resistance. Figure 3 maps all the detected and identified metabolites with the predicted functional genes. A high relative abundance of enzymes involved in carbon metabolism was also predicted based on the 16S rRNA gene sequences (Figure 4). Central carbon metabolism includes the Embden-Meyerhof-Parnas (EMP) pathway of glycolysis, the pentose phosphate pathway (PPP), the citric acid cycle (or TCA cycle), six known carbon fixation pathways, and some pathways of methane metabolism .
High relative abundances of TCA associated enzymes were predicted in the sediments based on the predicted functional gene data. Iron (Fe), present at elevated concentrations across all four estuaries (61.4 – 95.2 mg kg-1), is known to modulate the expression of critical enzymes of the TCA cycle including aconitase (ACO), citrate synthase (CS), isocitric dehydrogenase (IDH1), and succinate dehydrogenase (SDHB) . Fe also increases the formation of reducing equivalents such as NADH by the TCA cycle and thus increased mitochondrial O2 consumption and ATP formation via oxidative phosphorylation; NAD+, ADP, and ATP were all detected in the current study and demonstrate the importance of the TCA cycle here for establishing a microbial blueprint. Several genes such as sdhB and frdB that encode enzyme iron-sulphur subunits with a role in oxidative phosphorylation were predicted from the gene sequence data. Furthermore, increased Fe concentrations can also cause repression of glycolysis . This supports the negligible relative abundance of genes involved in glycolysis being predicted in the sampled sediments.
Several intermediates of the pentose-phosphate pathway such as fructose-6-phosphate, fructose-1,6-diphosphate, xylulose-5-phosphate, ribose-5-phosphate, ribulose-1,5-diphosphate, sedoheptulose-7-phosphate and glycerate-3-phosphate were detected and identified. Glycerate-3-phosphate is oxidised to produce pyruvate (Figure 3). Prediction of high relative abundance of genes encoding pyruvate oxidation enzymes dihydrolipoamide S-Acetyltransferase (DLAT) and dihydrolipoamide dehydrogenase (DLD) coupled with the identification of acetyl-CoA in the sediments suggests this pathway is indeed active amongst the sediment microbiomes.
A very high abundance of the genes acnA (encoding aconitate hydratase), sdhB (encoding succinate dehydrogenase) and frdB (encoding fumarate reductase) indicate enhanced reductive citric acid cycle (also called the Arnon-Buchanan cycle) activity (Figure 3). The reductive citric acid cycle is found in microaerophiles and anaerobes, such as green sulphur bacteria belonging to the Bacteroidota phylum. In the current study, members of the Bacteroidota, phylum were suggested to contribute towards carbon fixation in the metagenome of all estuaries (Supplementary Figure S8).
The reductive acetyl-CoA pathway (also called the Wood-Ljungdahl pathway) is found in strictly anaerobic Proteobacteria and Planctomycetes, some of which are methane-forming. The most important bifunctional enzymes, carbon monoxide dehydrogenase (encoded by cooS) and acetyl-CoA synthase (encoded by cdhE), catalyse CO2 1 CO and CO2 to a methyl group conversion, respectively. In the current study, members of Proteobacteria phylum were suggested to contribute towards methane metabolism in the metagenome of all estuaries (Supplementary Figure S9).
Noticeably, high relative levels of genes such as accABCD, which encode various subunits of acetyl-CoA carboxylase were predicted in this study (Figure 3). These results indicate that the predicted metagenome fixes CO2 by the 3-hydroxypropionate bicycle pathway commonly found in prokaryotes. This pathway is found in some green non-sulphur bacteria of the phylum Chloroflexi. Taxa-functional analysis of the predicted metagenomes also indicated that members of Chloroflexi substantially contributed to this pathway (Supplementary Figure S10).
Biomass is generated by fluxes branching out from central carbon metabolism, originating at twelve well-known precursor substances that include glucose-6-phosphate, fructose-6-phosphate, glyceraldehyde-3-phosphate, glycerate-3-phosphate, phosphoenolpyruvate, pyruvate, acetyl-CoA, ribose-5-phosphate, erythrose-4-phosphate, a-ketoglutarate, succinyl-CoA and oxaloacetate, and glycerate-1,3-diphosphate . Most of these precursor molecules were detected in the sediments (Figure 3). Glucose-6-phosphate, fructose-6-phosphate, glyceraldehyde-3-phosphate, and acetyl-CoA produce glycogen, cell wall components, and lipids. Several fatty acids were also identified in the sediments (Figure 3). Glycerate-3-phosphate is an important precursor in the biosynthesis of cysteine, serine, and glycine. Phosphoenolpyruvate produces tyrosine and tryptophan downstream (Figure 3).
Pyruvate branches out to form alanine, isoleucine, lysine, leucine, and valine. Isoleucine and leucine were identified in the sediments while several genes encoding enzymes catalysing these conversions were predicted in the sediments (Figures 3 and 5). Leucine is also produced from acetyl-CoA. Several amino acids including arginine, proline, glutamine, and glutamate are produced from a-ketoglutarate and were identified in the sediments. A very high relative proportion of genes involved in the biosynthesis of various amino acids was predicted (Figure 5). Members of the phyla Bacteriodota, Chloroflexi, and Proteobacteria were identified as contributing to the biosynthesis of amino acids in the sediment samples (Supplementary Figures S8-S10).
Community-wide biosynthetic and degradation potential for each metabolite in each sample was predicted through an integrated 16S rRNA gene sequencing-metabolomics approach using the MIMOSA2 webtool. Metabolites whose variation is consistent with expectations based on well-predicted microbial metabolic potential were identified. From these data, a community-wide metabolic model was constructed using MIMOSA2 for each sample and a community metabolic potential (CMP) was calculated, representing the relative capacity of the predicted community enzyme content in that sample to synthesize or degrade each metabolite based on metabolic reference information that links gene predictions to their substrates and products (metabolites) . Twenty-four metabolites with predicted CMPs were identified; however, none of these metabolites were found to be significant using the multi-variate statistical approaches described above.
Figure 6 illustrates the top 5 positively correlated CMP metabolites. Higher positive CMP scores indicate metabolites for which this agreement is statistically “well-predicted.” Negatively predicted metabolites can be interpreted in several different ways. Possible reasons for a negative correlation between a metabolite’s levels and its metabolic potential include incorrectly annotated or missing reactions (i.e., missing metabolites within the identified gene-metabolite pathway). However, contributors identified for negatively correlated metabolites have the potential to represent true relationships but should be interpreted more cautiously than positively correlated metabolites. The top 5 negatively correlated CMP metabolites are presented in Supplementary Figure S11.
The CMP successfully predicted alanine à pyruvate conversion through alanine-glyoxylate transaminase activity (K00830). Pyruvate may form 2-oxoisovalerate which produces 2-isopropylmalate downstream. 2-isopropylmalate à 3-isopropylmalate, followed by 3-isopropylmalate à 2-oxoisocaproate is predicted due to the activity of 3-isopropylmalate dehydratase and 3-isopropylmalate dehydratase, respectively. 2-oxoisocaproate may lead to the biosynthesis of leucine. Pyruvate is possibly oxidized to acetyl-CoA by pyruvate dehydrogenase (K00163) or to acetolactate by acetolactate synthase (K01652, K01653, K11258). Acetolactate may further proceed downstream to produce valine or leucine. O-Succinyl-L-homoserine produces succinate (TCA cycle intermediate) and 2-ketobutanoate (the keto acid precursor for isoleucine), by O-Succinyl-L-homoserine succinate-lyase (K01739). Isocitrate was also predicted to form succinate by isocitrate lyase (K01637). Isocitrate, an intermediate of the TCA cycle, may form a-ketoglutarate which branches out to form homocitrate and glutamate, possibly leading to the production of 2-aminoadipate and g-amino butanoate. 4-guanidinobutanoate is obtained from the degradation of g-amino butanoate by glycine amidinotransferase (K00613).
Several amino acids such as methionine, histidine, serine, tryptophan, phenylalanine, and arginine with predicted CMP were indicated in the current study. Arginine degradation also produces 4-guanidinobutanoate. In the current study, arginine à ornithine by arginase (K01476) à citrulline by ornithine carbamoyl transferase (K00611) was predicted. Citrulline is also predicted to be produced from N-acetyl ornithine by acetylornithine deacetylase (K01438). Methionine is synthesized from 4-methylthio-2-oxobutanoic acid by aromatic aminotransferase (K08969, K00832). Degradation of methionine forms 5’-methylthioadenosine via S-adenosyl-L-methionine intermediate by S-adenosyl-L-methionine methylthioadenosine-lyase (K01762).
Histidinal is degraded to histidine by histidinol dehydrogenase (K00013, K14152). Tryptophan synthase (K01694, K01695, K01696, K06001) catalyzes the conversion of serine à tryptophan. Serine à pyruvate is also catalyzed by serine dehydratase (K01754). Tryptophan is further degraded to kynurenine by indoleamine 2,3-dioxygenase (K00463). Another aromatic acid phenylalanine à phenethylamine is predicted by phenylalanine carboxy-lyase (K01593). L-cystine thiocysteine-lyase (K01760) catalyzes cystine à thiocysteine.
3.5 Metal homeostasis and metal resistance in estuarine sediments
Our results indicated very high concentrations of Fe (61.4 – 95.2 mg kg-1) and Al (45.5 – 68.2 mg kg-1) in all estuarine sediments. Some metals, such as Cu, Fe, Ni, and Zn, are essential nutrients and play important roles in the regulation of gene expression, and the activity of biomolecules, including enzymes and cofactors for essential biochemical reactions . Other metals, such as Al, Cd, Pb, As and Hg, have no biological role and are considered nonessential . Unsurprisingly, 125 KEGG genes related to metal homeostasis and metal resistance. These results are summarized in Supplementary Figure S12 and indicate the presence of a very high proportion of functional genes related to Fe and Zn metabolism, mostly related to membrane transport functions. The role of Fe in transport systems for Al uptake has been previously reported by Auger, Han . Relatively higher proportions of Fe-related genes may reduce the toxic effects of Al found in the sediments. There were no metabolites detected and identified that related to these specific genes, however, additional metallomics-based approaches could be applied to overcome this analytical limitation (i.e., assessing metal efflux in marine sediment microbiomes) .
3.6 Sulphur metabolism in estuarine sediments
High relative proportions of functional genes such as dsrAB (encoding dissimilatory sulphite reductase), aprAB (encoding adenylylsulphate reductase) and sat (encoding sulphate adenylyl transferase) were predicted in the metagenomes from all four estuaries (Supplementary Figure S13). SO42- à S2- is driven by the oxidation of organic carbon, supplemented by the anaerobic oxidation of methane (CH4) at the subsurface SO42--CH4 transition . Most of the S2- is ultimately re-oxidized back to SO42-, via diverse sulphur intermediates, by geochemical or microbial reactions . The assimilatory pathway produces reduced sulphur compounds for the biosynthesis of S-containing amino acids such as cysteine.
A high proportion of genes (Supplementary Figure S13) involved in sulphur oxidation such as soxAX (encoding L-cysteine S-thiosulfotransferase), soxB (encoding S-sulfosulfanyl-L-cysteine sulfohydrolase), soxC (encoding sulfane dehydrogenase), and soxYZ (encoding sulphur-oxidizing protein) were predicted in the current study. The large amounts of Fe in the estuarine sediments could increase the sulphur oxidation. Fe acts as an oxidant for S2- in the deeper sediment layers where it partly binds to S2- to form iron sulphide (FeS) and pyrite (FeS2).
The sulphur cycle of marine/estuarine sediments is largely microbial-driven, via dissimilatory sulphate (SO42-) reduction to sulphide (S2-) . Members of the Proteobacteria (mainly classes Alphaproteobacteria and Gammaproteobacteria) and Desulfobacterota (class Desulfobacteria) phyla contribute the most to the sulphur metabolism (Supplementary Figure S13). Members of Proteobacteria have been previously shown to contribute considerably to Fe cycling . While the metabolomics data did not provide any metabolites that directly feed into these pathways, characterizing the relevant proteins that drive these processes may enable greater insight into these functions and curation of microbial blueprints of activity.
3.7 Nitrogen metabolism in estuarine sediments
Several key functional genes (Supplementary Figure S14) related to nitrogen metabolism were predicted in all four estuaries, including genes (in decreasing order of relative proportion) involved in dissimilatory nitrate (NO3-) reduction to ammonium (NH4+) (DNRA), denitrification, assimilatory NO3- reduction to NH4+ (ANRA), nitrification, and nitrogen fixation. However, the inorganic metabolomics related to these processes is not measurable by the LC-MS metabolomics-based approaches used herein. Nitrogen cycling through the environment is influenced by microbially-driven processes . The denitrification and DNRA processes are mediated by the products of several genes illustrated in Supplementary Figure S14. Predicted functional gene analysis indicated the likely presence of these genes in all four estuaries. Taxa-functional analysis indicated that nitrogen metabolism is mediated by a diverse polyphyletic group of bacteria. The highest contribution to nitrogen cycling genes was attributed to Proteobacteria, Desulfobacterota, and Bacteroidota (Supplementary Figures S8-S9, S13).
Denitrification is one of the major nitrogen loss pathways that uses multiple electron donors such as hydrogen gas, hydrogen sulphide, or other organic compounds. Microorganisms can grow using DNRA, the key reaction of which is NO3- à NH4+, by coupling it to the oxidation of electron donors, such as organic matter, methane, hydrogen, sulphur compounds, and iron. The environmental importance of DNRA is not well established ; however, DNRA appears to be favoured over denitrification in marine and lake sediments when there is an excess of electron donors relative to NO3- . The results of this study agree with this observation (Supplementary Figure S14). High DNRA in estuarine samples indicates that nitrogen is retained in the system. This is evident from amino acid production discussed earlier.
The competition for nitrogen in estuaries is mostly for the inorganic forms of ions, commonly NH4+ or NO3-. NH4+ is the reduced form of nitrogen and, therefore, requires less energy to build amino acids. As such, NH4+ is preferred over NO3- and makes the competition for NH4+ fierce in estuaries.
The assimilation of NH4+ in prokaryotes occurs mainly via ammonium transporter proteins. Primary products of assimilated ammonia are glutamine and glutamate, which constitute the central reservoir of nitrogen for many biosynthetic pathways. The first pathway involves a high-affinity ammonium transporter, encoded by amtB, that moves NH4+ into the cell, where NH4+ + glutamate à glutamine by glutamine synthetase . The functional gene amtB was predicted in all the four estuaries indicating an assimilatory process in the studied environments. The other mechanism involves a low-affinity transport system in which NH3 is passively (diffusion) transported into the cell and NH4+ + α-ketoglutarate à glutamate by glutamate dehydrogenase . Glutamate, α-ketoglutarate, and glutamine were identified in the estuarine samples. Glutamate is the principal source of nitrogen for the production of N-amine and is involved in transamination reactions at the core of amino acid metabolism.