==== Front Environ Microbiol Rep Environ Microbiol Rep 10.1111/(ISSN)1758-2229 EMI4 Environmental Microbiology Reports 1758-2229 John Wiley & Sons, Inc. Hoboken, USA 36992633 10.1111/1758-2229.13148 EMI413148 Research Article Research Articles Using metatranscriptomics to better understand the role of microbial nitrogen cycling in coastal sediment benthic flux denitrification efficiency METATRANSCRIPTOMICS IN COASTAL N‐CYCLING Marshall et al. Marshall Alexis J. https://orcid.org/0000-0003-1990-6411 1 2 alexis.marshall@waikato.ac.nz Phillips Lori 2 6 Longmore Andrew 3 Hayden Helen L. https://orcid.org/0000-0002-5443-538X 2 4 Tang Caixian 1 Heidelberg Karla B. https://orcid.org/0000-0002-2645-8269 5 Mele Pauline https://orcid.org/0000-0002-1917-2553 1 2 7 1 La Trobe University AgriBio Centre for AgriBiosciences Bundoora Australia 2 Department of Jobs, Precincts and Regions AgriBio, Centre for AgriBiosciences Bundoora Australia 3 Centre for Aquatic Pollution Identification and Management Melbourne University Parkville Australia 4 School of Agriculture and Food, Faculty of Veterinary and Agricultural Sciences The University of Melbourne Parkville Victoria Australia 5 Department of Biology The University of Southern California Los Angeles California USA 6 Present address: Agriculture and AgriFood Canada Harrow Ontario Canada 7 Present address: Biomes Services Melbourne Victoria Australia * Correspondence Alexis J. Marshall, University of Waikato, Private Bag 3105, Hamilton, New Zealand 3240. Email: alexis.marshall@waikato.ac.nz 29 3 2023 8 2023 15 4 10.1111/emi4.v15.4 308323 10 10 2022 27 1 2023 © 2023 The Authors. Environmental Microbiology Reports published by Applied Microbiology International and John Wiley & Sons Ltd. https://creativecommons.org/licenses/by/4.0/ This is an open access article under the terms of the http://creativecommons.org/licenses/by/4.0/ License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. Abstract Spatial and temporal variability in benthic flux denitrification efficiency occurs across Port Phillip Bay, Australia. Here, we assess the capacity for untargeted metatranscriptomics to resolve spatiotemporal differences in the microbial contribution to benthic nitrogen cycling. The most abundant sediment transcripts assembled were associated with the archaeal nitrifier Nitrosopumilus. In sediments close to external inputs of organic nitrogen, the dominant transcripts were associated with Nitrosopumilus nitric oxide nitrite reduction (nirK). The environmental conditions close to organic nitrogen inputs that select for increased transcription in Nitrosopumilus (amoCAB, nirK, nirS, nmo, hcp) additionally selected for increased transcription of bacterial nitrite reduction (nxrB) and transcripts associated with anammox (hzo) but not denitrification (bacterial nirS/nirk). In sediments that are more isolated from external inputs of organic nitrogen dominant transcripts were associated with nitrous oxide reduction (nosZ) and changes in nosZ transcript abundance were uncoupled from transcriptional profiles associated with archaeal nitrification. Coordinated transcription of coupled community‐level nitrification–denitrification was not well supported by metatranscriptomics. In comparison, the abundance of archaeal nirK transcripts were site‐ and season‐specific. This study indicates that the transcription of archaeal nirK in response to changing environmental conditions may be an important and overlooked feature of coastal sediment nitrogen cycling. Department of Environment, Land, Water and Planning, State Government of Victoria 10.13039/100009675 Melbourne Water 10.13039/501100004170 source-schema-version-number2.0 cover-dateAugust 2023 details-of-publishers-convertorConverter:WILEY_ML3GV2_TO_JATSPMC version:6.3.0 mode:remove_FC converted:03.07.2023 Marshall, A.J. , Phillips, L. , Longmore, A. , Hayden, H.L. , Tang, C. , Heidelberg, K.B. et al. (2023) Using metatranscriptomics to better understand the role of microbial nitrogen cycling in coastal sediment benthic flux denitrification efficiency. Environmental Microbiology Reports, 15 (4 ), 308–323. Available from: 10.1111/1758-2229.13148 36992633 ==== Body pmcINTRODUCTION Untargeted sequencing‐based approaches have greatly enriched our understanding of microbial nitrogen cycling (e.g. the discovery of archaeal nitrifiers; Venter et al., 2004). We now recognize that few microbial taxa contain complete gene sets that enable full denitrification (Graf et al., 2014; Kuypers et al., 2018), and the complete microbial cycling of ammonia to nitrogen gas is generally considered a community‐level process (Anantharaman et al., 2016; Hug & Co, 2018). However, despite progress with DNA‐based approaches and the knowledge gained through metagenome assembled genomes (MAGs), functional evidence to support cooperative communities driving coupled metabolic processes is lagging. This is driven in part by a myriad of challenges associated with studying complex trophic systems and the availability of methods that study activity within an environmental context (Dar et al., 2021; Marlow et al., 2021). Metatranscriptomics is an untargeted approach that enables the study of natural microbial community responses to changing environmental conditions (Gilbert et al., 2008; Poretsky et al., 2010). This technique identifies genes that are or have recently been transcribed and require no a priori knowledge of the active microbial taxa or functions taking place (Frias‐Lopez et al., 2008; Moran et al., 2013; Poretsky et al., 2005). This technology is considered valuable for environmental monitoring because it utilizes short half‐life RNA assessments enabling the profiling of microbial responses to a range of disturbances within dynamic environmental conditions (Doney et al., 2004; Moran et al., 2013; Raes & Bork, 2008). In marine sediments, metatranscriptomics has improved our understanding of the microbial metabolic profiles of deep sea sediments under haloclines (Edgcomb et al., 2016), the interactions between the sediment and resting diatoms (Broman, Sachpazidou, et al., 2017), microbial community responses to transitions between oxic/anoxic conditions (Broman, Sachpazidou, et al., 2017; Broman, Sjöstedt, et al., 2017), the impacts of contaminated sediments (Birrer et al., 2019), the microbial response to changing redox gradients (Chen et al., 2017), and trophic interactions (Alexander et al., 2015; Dupont et al., 2015). This technique is a valuable molecular tool for monitoring microbial communities in coastal sediments as although community compositional changes have been linked to functional changes (Graham et al., 2016), community compositional changes do not always reflect changes in function (Grossmann et al., 2016). Port Phillip Bay (PPB) in south eastern Australia is globally recognized as an efficient coastal system for removing external sources of nitrogen (Berelson et al., 1998; Eyre & Ferguson, 2009). In this system, external riverine sources of organic carbon and reactive N (Nr) along with high rates of benthic recycling of Nr close to external inputs drive water column primary productivity (Harris et al., 1996). Spatial and seasonal variation in sediment nitrogen cycling efficiency within this system is well characterized through benthic flux chamber measurements (Berelson et al., 1998; Heggie et al., 1999; Marshall et al., 2021). The chambers measure the exchange of N2, NH4 +, NO3 − and NO2 − from the sediment to the water column. The denitrification efficiency (DE) is then expressed as the proportion of Nr released as N2. Within the system, two primary locations have been monitored annually in both spring and summer for ~20 years. In the centre of PPB, sediments are not directly impacted by external inputs and the DE estimates are globally high (~80%) and stable between spring and summer. In comparison, close to external riverine inputs of organic nitrogen at Hobsons Bay (HB) the DE is comparatively lower (45%–80%) and decreases between spring and summer commonly occur. The decrease in DE in HB has been attributed to increased levels of phytoplankton deposition on the sediment surface, which generates anoxic conditions leading to increased flux of Nr via the breakdown of coupled microbial nitrification–denitrification. Our previous attempts within this system to spatially and temporally couple variability in DE with changes in the sediment microbial nitrogen cycling community found no evidence to support a reduction in the capacity of the microbial community to engage in nitrification or denitrification close to external riverine inputs (Marshall et al., 2021, 2023). In contrast, higher transcript abundances associated with the archaeal nitrifier Nitrosopumilus were found closer to external inputs and unexpectedly to sediment depths of 10 cm (Marshall et al., 2023). At the same time, 16S rRNA amplicon surveys identified that the composition of the microbial community in sediments close to external inputs, but not in sediments isolated from inputs, shifted between spring and the following summer along with changes in the benthic flux DE (Marshall et al., 2021). The aim of this study was to determine the capacity for an untargeted metatranscriptomics approach to resolve site and seasonal features of microbial sediment nitrogen cycling within this well characterized coastal system. EXPERIMENTAL PROCEDURES Sample collection and RNA sequencing Surface sediment (0–1 cm) samples were collected from two locations within PPB, Victoria, Australia (Marshall et al., 2021). Site one is 24 m deep and is in the Central Port Phillip Bay (CPPB) (S38°03.495′ E144°52.242′) in a muddy sediment zone that is not directly impacted by external inputs. Site two is HB (S37°52.065′ E144°55.654′), which is 11 m deep, approximately 800 m from shore and is primarily influenced by the outflow from the Yarra River. Sediment was collected from each site at four timepoints (November 2014; March and November 2015; February 2016) in conjunction with the deployment of benthic chambers that measured the sediment DE through the flux of NH4 +, NO3 −, NO2 − and N2 (Marshall et al., 2021). Briefly, DE was consistently lower at HB than CPPB all time‐points (54–85 cf. 82%–98%) (Marshall et al., 2021). In CPPB, the DE did not decrease between spring and the following summer. In HB, DE decreased from spring 2014 (85 ± 8) to summer 2015 (54 ± 3) but not between spring 2015 (65 ± 4) to summer 2016 (66 ± 2). For these samples (total sample number = 24; CPPB = 12; HB = 12), sediment total organic carbon, total nitrogen, moisture content, pH, NH4 +—N and NO3 − + NO2 −—N, gene and transcript copy number of microbial nitrogen cycling markers, and total and active microbial community composition were determined and published by Marshall et al. (2021). Sediment cores were diver collected from each site with a stainless‐steel hand corer with a 50 mm diameter × 500 mm length cellulose acetate butyrate (CAB) plastic internal liner (Wildco). A modified polyvinyl chloride  piston was used to push sediment from the base of the CAB liner onto a sterile collection platform where the sediment was mixed prior to sampling. Sterile 5 mL syringes were used to collect surface sediment within 20 min of surfacing. Samples were immediately snap frozen in liquid nitrogen, transferred to a cryoshipper and stored at −80°C until extraction. RNA was extracted as in the study by Marshall et al. (2021). Briefly, high RNA quality (A260/280 = 1.9–2; RIN = 5.8–7.8; residual DNA digested and DNA removal quantified via QPCR) and quantity (18–334 ng/μL) (NanoDrop 2000 Thermo Fisher Scientific; Qubit 1.0 Thermo Fisher Scientific and Tape Station Agilent) was confirmed for all samples (Table 1). The RNA‐seq libraries were prepared using the Illumina TruSeq stranded mRNA sample preparation kit with RiboZero Gold ribosomal RNA depletion following manufacturer's instructions. An average fragment size of 253 bp was confirmed via TapeStation with High Sensitivity D1000 ScreenTape (Agilent). The evenness of each sample library within the total library pool was confirmed through a ‘spiked’ run on the Illumina HiSEQ 3000. The pooled library was sequenced on two separate occasions using HiSEQ 3000 with 2 × 151 bp sequencing technology at Agriculture Victoria, AgriBio Centre for AgriBiosciences, Victoria, Australia. TABLE 1 Summary of sample quality control measures and read counts of unassembled and assembled sequence reads. Sample information Extraction quality control Sequencing quality control Assembly quality control Site Date Core RNA ng ul−1 post DNase treatment A 260:280 post DNase treatment RIN Raw PE reads (M) Quality filtered PE reads (M) Non rRNA PE reads (M) Normalized non rRNA PE reads with >5× and <200× coverage (M) Non rRNA PE reads (M) that map with Bowtie against Trinity assembly Normalized non rRNA PE reads (M) that map with Bowtie against Trinity assembly CPPB Spring 2014 1 70.5 1.99 7 34.7 22.8 (65.7) 11.2 (32.3) 3.6 (10.4) 5.9 (52.2) 1.9 (52.1) 2 56.3 1.9 6.4 31.2 19.5 (62.5) 8.9 (28.5) 2.8 (9.0) 4.6 (51.4) 1.3 (46.9) 4 74.9 1.33 7.3 46.1 28.9 (62.7) 12.1 (26.2) 3.6 (7.8) 5.9 (49.0) 1.6 (45.8) Summer 2015 2 31.8 1.84 7.4 41.4 26.8 (64.7) 13.7 (33.1) 3.9 (9.4) 6.6 (47.9) 1.8 (44.9) 4 44.2 2 7.4 48.8 31.5 (64.5) 15.9 (32.6) 4.4 (9.0) 7.7 (48.3) 2.0 (46.3) 5 18.2 NA 7.8 35.1 22.6 (64.4) 12 (34.2) 4.1 (11.7) 4.5 (37.7) 1.7 (40.3) Spring 2015 1 59.7 1.67 7.5 49.7 32.4 (65.2) 15.2 (30.6) 4.7 (9.5) 7.0 (45.7) 2.2 (47.2) 2 24.8 1.83 7.5 42.1 29.4 (69.8) 21.3 (50.6) 6.8 (16.2) 8.6 (40.3) 2.8 (41.7) 4 49.8 2 6.9 48.6 31 (63.8) 14.4 (29.6) 4.7 (9.7) 5.8 (40.2) 2.0 (43.5) Summer 2016 1 57.8 1.9 6.8 44 27.5 (62.5) 13.8 (31.4) 4.3 (9.8) 5.8 (42.2) 1.8 (41.6) 2 59.6 1.95 6.9 62.5 37.9 (60.6) 16.5 (26.4) 5.2 (8.3) 7.1 (42.7) 2.2 (42.8) 3 57.6 1.97 7.3 46.8 28.7 (61.3) 12.9 (27.6) 3.9 (8.3) 5.8 (44.9) 1.8 (45.5) HB Spring 2014 1 67.9 1.93 7 43.4 25.2 (58.1) 9 (20.7) 3.5 (8.1) 3.5 (38.8) 1.4 (41.1) 3 110.3 1.97 5.9 39.9 23.2 (58.1) 8.3 (20.8) 3.4 (8.5) 3.3 (40.0) 1.4 (41.6) 4 29.6 1.67 7.3 42.1 25.7 (61.0) 8.9 (21.1) 3.5 (8.3) 3.7 (41.1) 1.4 (40.3) Summer 2015 1 333.8 2.09 6.3 42.5 26.5 (62.4) 12.1 (28.5) 5.3 (12.5) 4.9 (40.4) 2.0 (38.2) 3 49.8 1.91 7.1 47.1 29.1 (61.8) 14 (29.7) 4.3 (9.1) 6.2 (44.4) 1.8 (41.9) 4 195.2 2.04 6.4 43.7 26.4 (60.4) 9.9 (22.7) 4.8 (11.0) 3.2 (32.5) 1.6 (32.9) Spring 2015 1 97.6 2.08 6.3 45.2 25.5 (56.4) 7.3 (16.2) 3 (6.6) 2.4 (32.5) 0.9 (29.8) 3 65.7 1.92 7.3 50 28.4 (56.8) 10.6 (21.2) 4.2 (8.4) 3.3 (31.0) 1.4 (34.3) 5 112.9 2.08 6.9 51 29.1 (57.1) 8.7 (17.1) 3.5 (6.9) 2.8 (31.7) 1.2 (32.9) Summer 2016 2 80.8 2 5.8 42.9 25.3 (59.0) 11.8 (27.5) 4.6 (10.7) 3.6 (30.2) 1.6 (33.6) 3 91.7 2.02 6.9 53.8 31.3 (58.2) 10.6 (19.7) 3.9 (7.2) 3.6 (34.0) 1.5 (38.9) 4 80.8 1.95 7.1 37.8 21.7 (57.4) 7.4 (19.6) 2.9 (7.7) 2.4 (31.9) 0.9 (33.1) Note: Total RNA was extracted from 24 surface sediment samples collected from CPPB and HB in PPB, Victoria, Australia. Quality filtered paired‐end reads retained from Rcorrector, non‐rRNA paired‐end reads retained from SortMeRNA and reads retained after normalization with BBnorm. The percentage proportion of paired‐end (PE) reads in millions (M) retained from the original raw read counts are shown in parenthesis. The normalized non‐rRNA PE reads were assembled with Trinity. The impact of normalization on the mapping of data back to the Trinity assembly was assessed with Bowtie with the percentage of reads shown in parenthesis. Abbreviations: CPPB, Central Port Phillip Bay; HB, Hobsons Bay; NA, not available. Pre‐assembly quality trimming Erroneous k‐mers and low‐quality reads were identified and removed with rCorrector v 1.0.4 (Freedman, 2016; Song & Florea, 2015). Sequencing adapters, low‐quality bases (Phred <5), sequence reads with an error allowance greater than 0.1, and sequences <36 bases long were removed with the CutAdapt wrapper TrimGalore! v 0.6.4 (Krueger, 2007; Martin, 2011) following MacManes (2014). Ribosomal RNA reads were removed with SortmeRNA v 2.1 (Kopylova et al., 2012). Finally, reads with less than 5× coverage were removed from the dataset and highly represented reads were normalized to 200× coverage with BBnorm within BBTools v 38.81 (Bushnell, 2014). Sequenced reads were assessed at each quality step with FastQC v 0.11.9 (Andrews, 2012) and summarized with MultiQC v 1.9 (Ewels et al., 2016) (Table 1). Assembly, annotation and analysis The putative non‐ribosomal and quality screened reads from each sample (n = 24; CPPB = 12; HB = 12) were co‐assembled into transcript contigs with the de novo assembly software Trinity v 2.11.0 (Grabherr et al., 2011; Haas et al., 2013). Assembled transcripts were annotated via homology with the nitrogen cycle database NCycDB (Tu et al., 2019) using Blastn (v2.10) (Altschul et al., 1990). Taxonomy was inferred for those transcripts with NCycDB annotations with Blastn against the NCBI database (May 2021). All results were merged with the Trinotate sqlite database (https://www.sqlite.org/) into a single annotation data file (Supplementary datasheet S1). All computational analyses were enabled through the New Zealand eScience Infrastructure (NeSI) (https://www.nesi.org.nz). Expression values for Trinity assembled genes were estimated by RNA‐Seq Expectation Maximization (RSEM) (Li & Dewey, 2011) and bowtie (Langmead et al., 2009) taking into account strand specific data with the supplied Trinity wrapper align_and_estimate_abundance.pl (Grabherr et al., 2011; Haas et al., 2013). Trimmed mean of the M‐values (TMM) normalized read counts were calculated from the total assembly read count matrix (Supplementary datasheet S1). From this TMM matrix those transcript gene contigs that contained both an NCycDB annotation and were taxonomically identified as either bacteria or archaea against NCBI were selected from the dataset and log2 transformed for further analysis. A principal components analysis (PCA) was applied to assess the relationships among samples replicates for those transcript gene contigs that contained both an NCycDB annotation and were taxonomically identified as either bacteria or archaea with NCBI with the following settings: transformed to counts per million (CPM), minimum row sum = 10, log2 transformed and centred rows with the Trintiy supplied wrapper ‘PtR’ (Haas et al., 2013). R version 4.1.2 (R Core Team, 2021) and R studio (v 1.1.463) (RStudio Team, 2016) were used to interpret and display the data. Heatmaps were generated by summarizing and log2(× + 1) transforming the data with dplyr (v 1.0.7) (Wickham et al., 2022), prior to calculating Euclidian distance and hierarchical clustering using the complete method with heatmap3 (v 1.1.9) (Zhao et al., 2014). Heatmap gradients were displayed with viridis (v 0.6.2) (Garnier et al., 2021). Spearman correlation coefficients and adjusted Benjamini–Hochberg p values (Benjamini & Hochberg, 1995) were calculated with the core R package stats and psych (v 2.2.3) (Revelle, 2022), and significant relationships were displayed with corrplot (v0.92) (Wei & Simko, 2021). Differentially expressed (DEx) gene level transcripts were identified between locations and season sampled using EdgeR (Robinson et al., 2010). Significance was called if both the minimum false discovery rate (FDR) was <0.05 and a minimum fold change (FC) of 2 was identified between site comparisons. Differentially expressed gene level transcripts that both matched the NCycDB (Tu et al., 2019) and were taxonomically identified through Blastn (Altschul et al., 1990) as either bacteria or archaea were subsampled and then summarized in this study. RESULTS The PPB sediment metatranscriptome From 24 sediment RNA samples, a total of 1.7 billion paired‐end reads of 151 base pairs were generated with the Illumina HiSeq 3000 sequencing platform, with sequences having an average Phred score > 34 (Table 1). Of these paired‐end reads, 656 million passed k‐mer quality filtering and adapter trimming with 287 million paired‐end reads identified as non‐ribosomal (27% of the sequenced data). A greater number of reads passed these quality checks in CPPB samples (14 ± 3.1 million non‐rRNA reads) than HB (9.8 ± 2 million non‐rRNA reads) likely coupled to the lower efficiency of the sequence library rRNA depletion step to remove rRNA from HB samples (Table 1). After removing reads with <5× coverage and normalizing highly represented reads to 200× coverage the read representation between CPPB (4.3 ± 1 million non‐rRNA and normalized reads) and HB (3.9 ± 0.7 million non‐rRNA and normalized reads) was comparable. The de novo Trinity assembly was constructed with 9.3% of the initial sequenced data. The percentage of pre‐ and post‐normalized reads that mapped to the Trinity assembly with bowtie were comparable (Table 1). This assembly resulted in 357,071 unique transcript isoform contigs that clustered into 178,566 unique transcript gene level contigs. The average gene contig length was 503 bp. The number of gene contigs over 1 k bases was 26,170 with 191 gene contigs over 10 k bases and 9 gene contigs greater than 25 k bases. Annotation of nitrogen cycling transcripts within sediments of PPB Of the 178,566 unique gene level contigs, 5114 or 2.8% of the total Trinity assembled dataset were identified with homology to nitrogen cycling genes by annotation with the NCycDB (Tu et al., 2019). These NCycDB annotated transcripts have an N50 of 479 bp, median contig length of 265.5 and an average contig length of 439.17 bp. All analyses within this study were conducted at the Trinity gene level and are referred to throughout the text as transcripts. Of the 5114 NCycDB annotated transcripts Blastn annotation identified 151 archaeal and 1171 bacterial transcripts. Those transcripts containing both NCycDB and taxonomic annotations represented 0.7% (1322) of the total Trinity assembled transcripts (178,566). The vast majority (91%) of the archaeal transcripts were identified to the phylum level as Nitrososphaerota with 105 transcripts identified to the order Nitrosopumilales. Of the remaining archaeal transcripts, six were uncultured archaeon, six transcripts were associated with the Stenosarchaea group within the phylum Euryarchaeota, and a single transcript was associated with the phylum Crenarchaeota class Thermoprotei. These transcripts were functionally annotated predominantly as nitrite reductase, with 94 associated with nirk and 19 with nirS. Other dominant transcripts included 11 ammonia monooxygenase A (amoA), 8 ammonia monooxygenase B (amoB), and 8 ammonia monooxygenase C (amoC). The bacterial transcripts were represented by 875 unique species annotations: 292 transcripts were associated with Gammaproteobacteria, 178 Alphaproteobacteria, 178 Terrabacteria, 151 Deltaproteobacteria/Epsilonproteobacteria, 104 Fibrobacteres‐Chlorobi‐Bacteroidetes superphylum (FCB), 81 Betaproteobacteria, 54 Chlamydiae/Verrucomicrobia (PVC), 43 anammox/environmental group, 31 Nitrospina, 22 Acidobacteria, 14 Nitrospira, 7 Spirocheates, 5 unclassified Proteobacteria, and 11 transcripts were grouped together representing various other annotated taxa. Across all samples, the majority of the bacterial taxonomic diversity was reflected in 1055 transcripts that only had read counts in <2 of the 24 samples sequenced. The focused descriptive interpretation within this study was conducted on 236 (116 bacteria and 120 archaea) transcripts, selected from the total (178,566) assembled transcriptome (0.13% of the total transcriptome), on the basis of having both an NCycDB and Blastn taxonomic annotation and additionally contained read counts in at least 3 of the 24 sequenced samples (Supplementary datasheet S1). These transcripts have an N10 of 1813 bp and an N50 of 539 bp, with a median contig length of 365 bp, average contig length 518 bp. The RSEM mapping average for these 236 transcripts was 1781 ± 892 (CPPB) and 2698 ± 2194 (HB) with the normalized TMM mapping average of 273 ± 140 (CPPB) and 804 ± 577 (HB). These transcripts were predominantly associated with Nitrosopumilus (120), followed by Nitrospina (23), Gammaproteobacteria (23), Deltaproteobacteria/Epsilonproteobacteria (14), anammox/environmental group (13), Alphaproteobacteria (13), Terrabacteria (6), FCB group (6), Betaproteobacteria (6), Nitrospira (5), PVC group (4), Spirochaetes (1), Acidobacteria (1), and unclassified bacteria (1). Functional transcript representation was dominated by Nitrosopumilus with 82 nitrite reductase (nirK) and 9 nitrite reductase (nirS) transcripts, followed by Nitrospina with 15 nitrate reductase (narG) and 8 nitrite oxidoreductase (nxrB) transcripts assembled. Other well represented nitrogen cycling transcripts associated with various taxa included 14 glutamate dehydrogenase (gdh K15371), 10 nitrous oxide reductase (nosZ), 9 hydroxylamine reductase (hcp), and 13 nitrite reductase (7 nirK and 6 nirS) transcripts. Features of the active microbial nitrogen cycling community across sediments of PPB At both sites, transcripts were identified that are known to facilitate complete nitrification (Nitrosopumilus, Nitrospira and Nitrospina), anammox, and nitrite and nitrous oxide reduction. Nitrification was dominated by the archaea Nitrosopumilus sp. with only a single ammonia monooxygenase C (amoC) transcript identified associated with the bacterial nitrifier Nitrosospira multiformis. Copper containing nitrite reductase (nirK) was the dominant form of nitrite reductase assembled. Across the dataset, the most prevalent and abundant transcripts associated with both copper (nirk) and a cytochrome cd1‐containing (nirS) nitrite reductase were associated with Nitrosopumilus and not bacteria (Supplementary datasheet S1). Site based differences were identified at the individual transcript level by PCA (Figure 1A). When transcripts that contained TMM counts greater than 0 in at least 3 of the 12 site‐based samples were pooled (36 bacterial and 72 archaeal transcripts) and examined using Euclidean distance and complete clustering method, site‐based separation occurred coupled to variability in the abundance of Nitrosopumilus, Nitrospira, and Nitrospina transcript profiles in HB and nosZ and nrfA transcript profiles in CPPB (Figure 1B). FIGURE 1 Nitrogen cycling transcript features of sediments in Port Phillip Bay, Australia. (A) Transcript level separation of site sampled using principle coordinate analysis. (B) Dominant site‐based differences are visualized by pooling 236 TMM transformed and log2 normalized Trinity assembled gene level transcripts at the level of function and Blastx annotation. Only those transcripts that have a TMM value of 3 in at least 1 of the 24 samples are displayed. Archaeal transcripts are annotated as (A) and all other transcripts are associated with bacteria. The heatmap represents relationships calculated with Euclidean distance and complete clustering method and is scaled by row. CPPB: Central Port Phillip Bay; HB: Hobsons Bay. Ammonia monooxygenase (amoCAB), ferredoxin‐nitrite reductase (nirA), glutamate dehydrogenase (NAD(P)+) (gdh_K00261), glutamate dehydrogenase (gdh_K15371), glutamate synthase (NADPH/NADH) small chain (gs_K00266), hydrazine synthase subunit A (hzsA), hydroxylamine reductase (hcp), nitrate reductase (narG), nitrite oxidoreductase beta subunit (nxrB), nitrite reductase (cytochrome c‐552), (nrfA), nitrite reductase (NO‐forming) (nirK/S), nitrogenase iron protein (nifH), nitrous‐oxide reductase (nosZ). Within CPPB, the nitrogen cycling transcripts were represented by Nitrosopumilus (83), Gammaproteobacteria (19), Nitrospina (11), Deltaproteobacteria/Epsilonproteobacteria (10), Alphaproteobacteria (8), anammox/environmental group (8), Nitrospira (5), Terrabacteria (5), Betaproteobacteria (3), FCB group (3), PVC group (3), Acidobacteria (1), Nitrosospira (1), Spirochaetes (1), and unclassified bacteria (1) (Supplementary datasheet S1). Within CPPB sediment season and time‐point specific clustering of samples occurred. Grouping occurred in association with comparatively higher expression profiles of Nitrosopumilus (ammonia monooxygenase [amoAB] and nitrite reductase [nirk/nirS]) and Nitrospina (nxrB) in summer 2016 with a second cluster containing expression profiles associated with Nitrosopumilus (ammonia monooxygenase [amoC] and nitrite reductase [nirK]), Nitrospira (Candidatus Nitronauta litoralis; narG and nxrB) and Acanthopleuribacteraceae (nosZ) in spring 2015. Other samples collected in summer 2015 and spring 2014 clustered associated with higher abundances of Nitrosarchaeum (nirK) Planctomycetia bacterium and Geobacter sp. M18 nitrite reductase (cytochrome c‐552) (nrfA) in spring (Figure 2A). Spearman rank correlations with p‐adjusted (Benjamini–Hochberg) identified significant positive relationships between the transcript profiles of Nitrosopumilus amoA and amoB; and the anammox transcripts hzsA and hzsC (data not shown) in these sediments. FIGURE 2 Nitrogen cycling transcript features of (A) Central Port Phillip Bay (CPPB) and (B) Hobsons Bay (HB), Australia. Seasonal differences in transcript abundance at each site are visualized by pooling TMM and log2 normalized Trinity assembled gene level transcripts at the level of function and Blastx taxonomic annotation. At each site, only those transcripts that have a TMM value of 3 in at least 1 of the 12 site‐based samples are displayed (27 transcripts in CPPB and 25 in HB). Archaeal transcripts are annotated as (A) and all other transcripts are associated with bacteria. Each heatmap represents relationships calculated with Euclidean distance and complete clustering method and scaled by row. Ammonia monooxygenase subunit (amoCAB), ferredoxin‐nitrite reductase (nirA), glutamate dehydrogenase (NAD(P)+) (gdh_K00261), glutamate dehydrogenase (gdh_K15371), glutamate synthase (NADPH/NADH) small chain (gs_K00266), hydrazine synthase subunit A (hzsA), hydroxylamine reductase (hcp), nitrate reductase (narG), nitrite oxidoreductase beta subunit (nxrB), nitrite reductase (cytochrome c‐552) (nrfA), nitrite reductase (NO‐forming) (nirK/S), nitrogenase iron protein (nifH), nitrous‐oxide reductase (nosZ). In HB, the nitrogen cycling transcripts were represented by Nitrosopumilus (101), Nitrospina (20), anammox/environmental group (8), Deltaproteobacteria/Epsilonproteobacteria (8), Gammaproteobacteria (8), FCB group (5), Alphaproteobacteria (4), Betaproteobacteria (4), Terrabacteria (3), Nitrosospira (1) and Nitrospira (1). Within these sediments, sample specific clustering occurred in association with comparatively higher expression profiles of Nitrosopumilus sp. (ammonia monooxygenase [amoCAB] and nitrite reductase [nirk]), and those transcripts identified as Nitrospina and Candidatus Nitronauta litoralis (nitrite reductase [nxrB] and nitrate reductase [narG]) (Figure 2B). Spearman rank correlations (p‐adjusted BH) identified multiple positive relationships across nitrogen cycling taxa in HB, with increased transcript activity correlated between archaeal (amoCAB, nmo, hcp, nirK) and bacterial (nxrB and narG) nitrifiers and anammox (hzo) (Figure 3). FIGURE 3 Relationships between pooled TMM and log2 transformed nitrogen cycling functional genes from sediments collected in Hobsons Bay, Australia. Archaeal transcripts are annotated as (A) and all other transcripts are associated with bacteria. Matrix of Spearman's correlation coefficients, displaying only the transcripts identified as significant at FDR‐adjusted p < 0.05 (Benjamini & Hochberg, 1995). Ammonia monooxygenase subunit (amoCAB), ferredoxin‐nitrite reductase (nirA), glutamate dehydrogenase (NAD(P)+) (gdh_K00261), hydrazine oxidoreductase (hzo), hydroxylamine reductase (hcp), nitrate reductase (narG), nitronate monooxygenase (nmo), nitrite oxidoreductase beta subunit (nxrB), protein NrfC (nrfC), nitrite reductase (NO‐forming) (nirK). At each time point (spring 2014, summer 2015, spring 2015, summer 2016) using the total gene level transcript database (178,566), site‐based differences between CPPB (n = 3) and HB (n = 3) were identified with EdgeR (Table 2; log2 fold change >2 and FDR‐adjusted p < 0.05). Only transcripts that were identified by EdgeR as different between sites and able to be annotated for both the NCycDB and were taxonomically annotated with Blastn were selected for further analysis. The greatest number of site‐based differences were identified in spring 2015 (transcript total = 36; CPPB = 13, HB = 23), followed by summer 2016 (transcript total = 35; CPPB = 11, HB = 24), summer 2015 (transcript total = 23; CPPB = 3, HB = 20), and spring 2014 (transcript total = 4; CPPB = 0, HB = 4). Transcripts that were consistently greater in CPPB include those with homology to bacterial dissimilatory nitrate reduction (nrfA) and nitrous oxide reductase (nosZ) (Table 2). In HB, transcripts with homology to archaeal ammonia oxidation (amoC and nirK) and bacterial nitrite (nxrB) and nitrate (narG) reduction were consistently greater than in CPPB (Table 2). Across all time‐points, 33 distinct Nitrosopumilus nirK transcripts were identified as different between sites. Of these 33 distinct Nitrosopumilus nirK transcripts, 3 were consistently identified as different in all time points, 3 were different at 3 time points and 6 were different at 2 time points, the remainder were only identified once (Table 2). TABLE 2 The number of differentially expressed (EdgeR p < 0.05 and fold change greater than 2) and annotated transcripts between the sediment metatranscriptome of Central Port Phillip Bay (CPPB) and Hobsons Bay (HB). Pathway Gene (sub) families Annotation Blastn_Species Spring 2014 Summer 2015 Spring 2015 Summer 2016 CPPB HB CPPB HB CPPB HB CPPB HB Anammox nirS Nitrite reductase (NO‐forming) uncultured anaerobic ammonium oxidizing bacterium 1 Assimilatory nitrate reduction nirA Ferredoxin‐nitrite reductase Microbulbifer sp. ALW1 1 Biodegradation gdh_K00262 Glutamate dehydrogenase (NADP+) Citrobacter portucalensis 1 gdh_K15371 Glutamate dehydrogenase Burkholderia diffusa 1 Lutibacter profundi 1 Microbacterium sp. 1.5R 1 Myxococcus xanthus 1 ureB Urease subunit beta Stenotrophomonas maltophilia 1 1 Biosynthesis gs_K00266 Glutamate synthase (NADPH/NADH) small chain Moraxella osloensis 1 1 Denitrification nirK Nitrite reductase (NO‐forming) Stigmatella aurantiaca DW4/3‐1 1 nosZ Nitrous‐oxide reductase Anaeromyxobacter sp. Fw109‐5 1 1 1 nosZ Nitrous‐oxide reductase Desulfosarcina ovata subsp. sediminis 1 Dissimilatory nitrate reduction nrfA Nitrite reductase (cytochrome c‐552) Geobacter sp. M18 1 1 nrfA Nitrite reductase (cytochrome c‐552) Planctomycetia bacterium 1 1 Hydroxylamine reductase hcp Hydroxylamine reductase Draconibacterium orientale 1 1 Gammaproteobacteria bacterium 1 Pseudomonas resinovorans NBRC 106553 1 Nitrification amoA (A) Ammonia monooxygenase subunit A (archaea) Nitrosopumilales 3 1 amoB (A) Ammonia monooxygenase subunit B (archaea) 1 1 1 amoC (A) Ammonia monooxygenase subunit C (archaea) 1 1 1 2 hcp hydroxylamine reductase 1 nirK Nitrite reductase (NO‐forming) 3 13 5 14 4 14 nirS nitrite reductase (NO‐forming) 1 nxrB Nitrite oxidoreductase, beta subunit Candidatus Nitronauta litoralis 1 1 1 Nitrification/denitrification narG Nitrate reductase Candidatus Nitronauta litoralis 1 3 Note: Comparison was achieved at each time point sampled with only those transcripts with both NCycDB and Blastn annotation summarized. Comparisons within a site by time point sampled (e.g. spring 2014 vs. summer 2015) or season sampled (e.g. spring vs. summer) identified that the abundance of a single transcript with homology to nosZ Winogradskyella helgolandensis was higher in HB sediments in spring 2015 when compared to both summer 2015 and summer 2016. No other significant differences were identified between dates sampled (e.g. spring to the following summer) or when season (e.g. spring vs. summer) were combined at either site. Transcripts associated with Nitrosopumilus in PPB sediments The dominant nitrogen cycling transcripts within sediments of PPB were associated with Nitrosopumilus sp. In both CPPB and HB, Nitrosopumilus sp. activity profiles clustered samples into season sampled (Figure 4). In CPPB (Figure 4A,C), strong co‐expression profiles between ammonia monooxygenase (amoCAB), and nitrite reductase (nirk) were identified. In HB (Figure 4B,D), strong co‐expression profiles of ammonia monooxygenase (amoCAB), hydroxylamine reductase (hcp), nitronate monooxygenase (nmo) and nitrite reductase (nirk) were identified. At both sites, transcripts annotated as Nitrosopumilus sp. nosZ were identified with sample specific read mapping support (Supplementary datasheet S1). FIGURE 4 Relationships between transcripts with homology to Nitrosopumilus sp. in sediments of Central Port Phillip Bay (A, C) and Hobsons Bay (B, D), Australia. Seasonal differences are visualized by pooling TMM and log2 normalized Trinity assembled gene level transcripts at the level of function and phyla. Each heatmap (A, B) represents relationships calculated with Euclidean distance and complete clustering method. The heatmap is scaled by row. Matrix of Spearman's correlation coefficients (C, D), displaying only relationships with FDR‐adjusted p < 0.05 (Benjamini & Hochberg, 1995). Ammonia monooxygenase subunit (amoCAB), glutamate dehydrogenase (NAD(P)+) (gdh_K00261), hydroxylamine reductase (hcp), nitronate monooxygenase (nmo), nitrite reductase (NO‐forming) (nirK/S), nitrous‐oxide reductase (nosZ). A large number (81) of the 236 transcripts investigated in this study were annotated as Nitrosopumilus sp. nitrite reductase (nirK). Of these 81 Nitrosopumilus sp. nirK transcripts, 39 had TMM read support values greater than 3 in at least one sample across all 24 samples. These 39 nirk transcripts displayed site‐specific clustering using heatmap analysis. Within CPPB higher abundances of nirK transcripts additionally clustered samples within‐site by season sampled (e.g. spring with summer) (Figure 5). FIGURE 5 Relationships between transcripts with homology to Nitrosopumilus sp. nitrite reductase (nirk) transcript profiles of sediments collected from Central Port Phillip Bay and Hobsons Bay, Australia. Read counts with a minimum TMM value of 3 in at least 1 of the 24 samples are displayed as log2 normalized values. The heatmap is scaled by row. DISCUSSION In coastal sediments, four microbial pathways have been identified with the capacity to generate N2 or N2O and contribute to Nr‐loss. These are denitrification, anammox, nitrifier‐denitrification and N2 or N2O as a potential by‐product of metabolic oxygen production under anoxic conditions by the archaeal nitrifier Nitrosopumilus (Kraft et al., 2022). Denitrification is commonly identified as the dominant Nr‐loss pathway in coastal sediments, yet it is the only microbial pathway that requires community‐level coordination. In this study, transcripts with homology to all four of these Nr‐loss pathways were detected within sediments of PPB. However, the prevalence of transcripts associated with Nitrosopumilus activity, particularly NO‐forming nitrite reductase (nirK), was substantially greater than those associated with denitrification (nirS, nirK, nosZ) or anammox (hzs, hzo, nirS). There was also no evidence to support nitrifier‐denitrification within the single bacterial nitrifier Nitrosospira. As in our previous work (Marshall et al., 2021), no clear spring to summer variation in nitrification (ammonia monooxygenase [amoCAB]) transcript profiles were identified in this study in association with decreases in benthic flux DE. However, site‐specific signatures of Nitrosopumilus sp. nirK transcript abundance highlights that how Nitrosopumilus sp. respond to variable environmental conditions may be an important and overlooked feature of coastal sediment nitrogen cycling. The application of untargeted metatranscriptomics resolved key site‐specific features of the coastal sediment microbial nitrogen cycle, identifying that a few highly active nitrogen cycling taxa were detected within each site. We present in Figure 6 a summary of the proposed active pathways and key taxa involved in nitrogen cycling for each site studied in PPB. FIGURE 6 Graphical summary of site specific features of sediment nitrogen cycling transcripts within Central Port Phillip Bay (A) and Hobsons Bay (B), Australia. In the sediments of CPPB, the dominant nitrogen cycling transcript was associated with nitrous oxide reduction (nosZ) with homology to the swarming predatory bacterial taxa Myxobacteria (Anaeromyxobacter sp. Fw109‐5). These nosZ transcripts were ~13× higher in sediments within CPPB than HB and were consistently transcribed uncoupled from variability in transcripts with homology to the archaeal nitrifier Nitrosopumilus (amoCAB, nirk). Other prevalent transcripts within these sediments were associated with Nitrosopumilus (amoCAB, nirk) and bacterial dissimilatory nitrite reduction to ammonium (DNRA) (nrfA). We have previously demonstrated with 16S rRNA amplicon sequencing that within the active microbial community composition sequences associated with Nitrosopumilus and Myxobacteria both had a higher relative abundance in CPPB sediments compared to HB (Marshall et al., 2021). Anaeromyxobacter sp. are facultative anaerobes known to reduce various metals, dechlorinate aromatic compounds, and perform dissimilatory nitrate reduction to ammonium (Sanford et al., 2002). They have also been identified as playing a key role in the transformation of N2O to N2 in soils through the chemical transformation of NO2 − to NO via iron (Fe2 +) oxidation (Onley et al., 2018; Sanford et al., 2012). It is plausible that the activity of Myxobacteria generates the globally high denitrification efficiencies (>80%) characteristic of these central basin sediments. This could be achieved uncoupled from other microbial nitrogen cycling processes (e.g. nitrification) through predatory or scavenging behaviours, as these sediments receive a high deposition rate of phytoplankton. The subsequent cell death and lysis of intracellular phytoplankton derived NO2 − could support Myxobacteria N2 production. Alternatively, Nitrosopumilus has been demonstrated to generate NO2 −, and NO under aerobic conditions (Martens‐Habbena et al., 2015) and N2O, and N2 under anoxic conditions (Ji et al., 2018; Kraft et al., 2022; Santoro et al., 2011, 2021). It is possible that Nitrosopumilus additionally contributes N2 to the DE flux within CPPB, although in this study Nitrosopumilus activity profiles were comparatively lower at this site than HB. Our previous research additionally demonstrated that markers of archaeal nitrifier activity (amoA and nirK‐a) were restricted to the surface sediment (0–1 cm) in CPPB and at comparatively lower abundances than in HB (Marshall et al., 2023). This study indicates that the individual activity of Myxobacteria is the most probable primary driver of high denitrification efficiencies within these central zone sediments, and interactions between the production of NO2 − by Nitrosopumilus and the reduction of NO2 − by Myxobacteria may further support coordinated transcription of community‐level coupled nitrification–denitrification. Activity profiles associated with bacterial dissimilatory nitrite reduction to ammonium (DNRA) (nrfA) were detected in CPPB sediments but not in HB. The abundance of bacterial nrfA transcripts in these central zone sediments were 6× higher than the most abundant bacterial nitrite reductase (nirS) transcript associated with denitrification. In estuarine and lake sediments the selection of Nr retention via DNRA over Nr loss via denitrification has been coupled with high concentrations of Fe2 + (>258 μM) and limiting NO3 − concentrations (Cojean et al., 2020; Kessler et al., 2018; Robertson et al., 2016; Robertson & Thamdrup, 2017). Further investigation is required to determine the role that Fe2 + may play in selecting the sediment nitrogen cycling community in the large central basin of PPB, as comparable pore water concentrations of NO3 − + NO2 − were determined for these samples (Marshall et al., 2021) and historical records found no differences in Fe% in the surface sediment between the central muds and HB (Fabris et al., 1999). HB sediments are impacted by the external inflow of organic nitrogen from the Yarra River and are characterized as having lower and more seasonally variable denitrification efficiencies than CPPB (Berelson et al., 1998; Marshall et al., 2021). This seasonal variability is hypothesised to be coupled with increased phytoplankton deposition in summer months, which decreases oxygen availability at the sediment surface uncoupling microbial nitrification–denitrification. Within this study, transcriptional signatures of nitrifier activity in HB were associated with Nitrosopumilus, and activity profiles were positively coupled with increases in bacterial nitrate (nxrB) and nitrite (narG) reduction profiles, and signatures of anammox activity (hzo) but not with bacterial nitrite reduction (bacterial nirS). No clear evidence was identified in these sediments to support a decrease in nitrifier activity between spring and the following summer in association with variability in the benthic flux DE (Marshall et al., 2021). However, our interpretation may be impacted as the DE remained low in spring 2015 and did not decrease as expected between spring 2015 (65 ± 4) to summer 2016 (66 ± 2) (Marshall et al., 2021). Weak evidence to support bacterial denitrification via nosZ transcripts were detected in spring 2014 associated with Anaeromyxobacter sp. Fw109‐5 and in spring 2015 and summer 2016 associated with Winogradskyella helgolandensis. Despite being significantly less abundant than in the Central Port Phillip Bay sediments the presence of nosZ in Hobsons Bay provides some support that N2O to N2 transformation via denitrification is occurring at this site, and at higher levels in spring 2014 than summer months when the benthic flux DE was comparatively more efficient (spring 2014; 85 ± 8) (Marshall et al., 2021). These findings provide transcriptional support within HB sediments for coupled interactions between archaeal and bacterial nitrifiers and little evidence to support coupled nitrification–denitrification. Within HB, the dominant nitrite reduction activity profiles were associated with archaeal NO‐forming nitrite reductase (nirK). Kraft et al. (2022) reported that axenic cultures of N. maritimus persist with ammonia oxidation under anoxia to produce both N2 and O2, whilst generating the intermediates NO and the greenhouse gas N2O. The mechanisms that Nitrosopumilus uses to produce metabolic oxygen are still unresolved yet support is growing for a role for nirK. In addition to Kraft et al. (2022), expression of nirK by the representative N. maritimus strain SCM1 was strongly coupled to both increased ammonium and copper availability (Qin et al., 2018), and in our previous work (Marshall et al., 2021, 2023) nirK‐a activity profiles were high in HB to sediment depths of 10 cm and strongly correlated with amoA. Within this study, Nitrosopumilus transcripts were also associated with homologues to hydroxylamine reductase (hcp) and nitronate monooxygenase (nmo). This study applied a homology‐based approach using the NCycDB (Tu et al., 2019) and NCBI taxonomy and a role for hcp and nmo in nitrification within Nitrosopumilus is unresolved. The high abundance of archaeal nitrifiers in anoxic coastal sediments have historically been coupled to oxygen availability at the exchange layer between the sediment–water interface and to physical and biological processes that resuspend and re‐oxygenate sediment (Beman et al., 2012; Lipsewers et al., 2014). Here within an environmental context, we report transcriptional evidence that diverse homologues for Nitrosopumilus nirk are abundant and transcribed across seasons in anoxic sediments impacted by external sources of organic nitrogen. This finding is in conflict with the long‐standing hypothesis that archaeal ammonia oxidation requires environmental oxygen and is inhibited under environmentally induced anoxia. Here we highlight that differences in how Nitrosopumilus engages in nitrogen cycling, especially the environmental drivers that lead to nirK transcription, may be an overlooked and important feature of coastal sediment nitrogen cycling. A finding that is relevant to the management of coastal regions as nirK transcription has been coupled to the production of the greenhouse gas N2O (Kraft et al., 2022; Qin et al., 2018). Genomics based approaches have associated a community‐based theory to nitrogen cycling (Anantharaman et al., 2016; Hug & Co, 2018). Many taxa contain the genomic repertoire to engage in nitrogen cycling, especially via denitrification (Graf et al., 2014; Kuypers et al., 2018), leading to a consortia of microbial taxa required to collectively engage in Nr‐loss. In this study, we identified that there is strong evidence for coordinated transcription between archaeal and bacterial nitrifier activity profiles and weak evidence to support a coordinated community‐level microbial nitrification–denitrification transcriptional response. As in our previous work (Marshall et al., 2021, 2023), we found no evidence to couple community‐level functional transcript abundances with benthic flux DE measures. In addition, the untargeted investigation of individual nitrogen cycling transcripts identified that at each location active nitrogen cycling genes were restricted to a few taxa that reflected site‐specific functional profiles. Transcriptionally active multi‐species community‐level nitrogen cycling was not overly supported with this untargeted approach, highlighting that we may have potentially overestimated community‐level reliance to facilitate high levels of coastal Nr‐loss. Incorporating quantitative DNA‐based microbial trait‐based data into ecosystem‐level biogeochemical models is challenging (Raes et al., 2020; Zhang et al., 2018). Our findings suggest that biogeochemical models that look to estimate the retention and loss of Nr within coastal systems may be improved by first identifying the specific nitrogen cycling taxa that are active through untargeted approaches and then selectively targeting temporal and spatial variability using high‐throughput quantitative methods (e.g. RT‐QPCR). Many environmental studies take the approach of sequencing DNA, assembling genomic scaffolds and grading high quality MAGs prior to mapping RNA to functionally annotated genes. By assembling genes into genomes, this approach strengthens the association of a functional gene to an annotated MAG strengthening the functional credibility of the mapped RNA profiles. However, support is building for the application of untargeted metatranscriptomics without DNA as a proxy for studying transcriptionally active organisms within environmental research (Söllinger et al., 2018; Täumer et al., 2022). One of the constraints for the combined DNA and RNA approaches is that substantial sequencing effort is required to resolve MAGS in diverse environments and commonly only dominant taxa with relatively small genome sizes assemble to meet high‐quality thresholds. This generates a bias to DNA centric approaches as there is no guarantee that the most prevalent taxa will be the most active taxa. Although it is not without its own technical assumptions, using a de novo assembly approach without constraining the sediment activity profiles to genomic features is a viable alternative to circumvent this challenge in complex samples. In this study, this untargeted and unconstrained approach identified a key feature of Nr‐loss in CPPB via nosZ with homology to Anaeromyxobacter. Although the relative abundance of this taxa was higher in CPPB sediments than HB (Marshall et al., 2021), the total relative abundance of the dominant Myxobacteria Amplicon Sequence Variant (ASV) was associated with the family Polyangiaceae and represented <0.2% of the active community composition. The potential role for Anaeromyxobacter was overlooked with 16S rRNA amplicon sequencing as the single ASV identified was removed during quality control due to low relative abundance (<0.005%)  (Marshall et al., 2021). The low relative abundance of Anaeromyxobacter suggests that a relatively large Myxobacteria associated genome would be challenging to resolve via untargeted environmental metagenomics, identifying that substantial sequencing effort would be required to confirm these results. Profiling this community through untargeted metatranscriptomics is only the first step in identifying the role of Myxobacteria within the sediments of PPB. Selective culturing and targeted long‐read genome sequencing will be required to identify the role that Myxobacteria play in nitrogen cycling within the large central basin of PPB. In conclusion, we find that untargeted environmental metatranscriptomics can resolve differences in how microbial communities engage in sediment nitrogen cycling. In contrast to studies that apply gene‐based quantifications with universal primer targets (Marshall et al., 2023) untargeted metatranscriptomics can resolve both the functional gene and the taxonomic identity. By applying this approach, we identified that there is strong evidence for coordinated community‐level transcription of archaeal and bacterial nitrifier activity, and weak evidence to support a coordinated community‐level microbial nitrification–denitrification transcriptional responses within HB sediments. Therefore, seasonal changes in DE at this site may not be driven by impacts to coupled community‐level nitrification–denitrification as it is a weakly supported microbial interaction. Alternatively, site‐specific environmental selection of metabolic processes that determine how Nitrosopumilus engages in nitrogen cycling via pathways that utilize nirK, and spatial drivers that select for Anaeromyxobacter sp. Fw109‐5 nosZ transcription may be overlooked yet important features of Nr‐loss in this coastal system. Future targeted research that qualifies the impact of increased nirK transcription on archaeal nitrification is required to support coastal management and resolve the contribution of archaeal nitrifiers in the production of the greenhouse gas N2O. AUTHOR CONTRIBUTIONS Alexis Marshall: Conceptualization (supporting); data curation (lead); formal analysis (lead); methodology (lead); visualization (lead); writing – original draft (lead); writing – review and editing (equal). Lori A Phillips: Conceptualization (equal); funding acquisition (equal); methodology (supporting); project administration (equal); supervision (equal); writing – review and editing (equal). Helen L Hayden: Formal analysis (supporting); methodology (supporting); supervision (supporting); writing – review and editing (equal). Andrew Longmore: Conceptualization (equal); funding acquisition (equal); project administration (equal); resources (equal); writing – review and editing (supporting). Caixian Tang: Project administration (supporting); supervision (supporting); writing – review and editing (supporting). Karla Heidelberg: Conceptualization (equal); funding acquisition (equal); methodology (supporting); project administration (equal); supervision (equal); writing – review and editing (equal). Pauline Mele: Conceptualization (equal); funding acquisition (equal); project administration (equal); resources (lead); supervision (lead); writing – review and editing (equal). CONFLICT OF INTEREST STATEMENT The authors declare no conflict of Interest. Supporting information Data S1: Supporting Information   emi413148‐sup‐0001‐supinfo.csv Click here for additional data file. ACKNOWLEDGEMENTS The authors acknowledge funding from Department of Environment, Land, Water and Planning, State Government of Victoria and Melbourne Water for carrying out the experimentation. Open access publishing facilitated by The University of Waikato, as part of the Wiley ‐ The University of Waikato agreement via the Council of Australian University Librarians. DATA AVAILABILITY STATEMENT Sequence data is available in Supporting Information Data S1. ==== Refs REFERENCES Alexander, H. , Rouco, M. , Haley, S.T. , Wilson, S.T. , Karl, D.M. & Dyhrman, S.T. (2015) Functional group‐specific traits drive phytoplankton dynamics in the oligotrophic ocean. Proceedings of the National Academy of Sciences, 112 , E5972–E5979. Available from: 10.1073/pnas.1518165112 Altschul, S.F. , Gish, W. , Miller, W. , Myers, E.W. & Lipman, D.J. (1990) Basic local alignment search tool. Journal of Molecular Biology, 215 , 403–410. Available from: 10.1016/S0022-2836(05)80360-2 2231712 Anantharaman, K. , Brown, C.T. , Hug, L.A. , Sharon, I. , Castelle, C.J. , Probst, A.J. et al. (2016) Thousands of microbial genomes shed light on interconnected biogeochemical processes in an aquifer system. Nature Communications, 7 , 13219. Available from: 10.1038/ncomms13219 Andrews, S. (2012) FastQC: a quality control tool for high throughput sequence data. http://www.bioinformatics.babraham.ac.uk/projects/fastqc/ Beman, J.M. , Bertics, V.J. , Braunschweiler, T. & Wilson, J. (2012) Quantification of ammonia oxidation rates and the distribution of ammonia‐oxidizing archaea and bacteria in marine sediment depth profiles from Catalina Island, California. Frontiers in Microbiology, 3 , 263. Available from: 10.3389/fmicb.2012.00263 22837756 Benjamini, Y. & Hochberg, Y. (1995) Controlling the false discovery rate: a practical and powerful approach to multiple testing. Journal of the Royal Statistical Society: Series B: Methodological, 57 , 289–300. Berelson, W.M. , Heggie, D. , Longmore, A. , Kilgore, T. , Nicholson, G. & Skyring, G. (1998) Benthic nutrient recycling in Port Phillip Bay Australia. Estuarine, Coastal and Shelf Science, 46 , 917–934. Available from: 10.1006/ecss.1998.0328 Birrer, S.C. , Dafforn, K.A. , Sun, M.Y. , Williams, R.B.H. , Potts, J. , Scanes, P. et al. (2019) Using meta‐omics of contaminated sediments to monitor changes in pathways relevant to climate regulation. Environmental Microbiology, 21 , 389–401. Available from: 10.1111/1462-2920.14470 30411468 Broman, E. , Sachpazidou, V. , Dopson, M. & Hylander, S. (2017) Diatoms dominate the eukaryotic metatranscriptome during spring in coastal ‘dead zone’ sediments. Proceedings of the Royal Society B, 284 , 20171617. Available from: 10.1098/rspb.2017.1617 28978732 Broman, E. , Sachpazidou, V. , Pinhassi, J. & Dopson, M. (2017) Oxygenation of hypoxic coastal Baltic Sea sediments impacts on chemistry, microbial community composition, and metabolism. Frontiers in Microbiology, 8 , 2453. Available from: 10.3389/fmicb.2017.02453 Broman, E. , Sjöstedt, J. , Pinhassi, J. & Dopson, M. (2017) Shifts in coastal sediment oxygenation cause pronounced changes in microbial community composition and associated metabolism. Microbiome, 5 , 96. Available from: 10.1186/s40168-017-0311-5 28793929 Bushnell, B. (2014) BBTools software package. Chen, J. , Hanke, A. , Tegetmeyer, H.E. , Kattelmann, I. , Sharma, R. , Hamann, E. et al. (2017) Impacts of chemical gradients on microbial community structure. The ISME Journal, 11 , 920–931. Available from: 10.1038/ismej.2016.175 28094795 Cojean, A.N.Y. , Lehmann, M.F. , Robertson, E.K. , Thamdrup, B. & Zopfi, J. (2020) Controls of H2S, Fe2 +, and Mn2 + on microbial NO3 − reducing processes in sediments of an eutrophic lake. Frontiers in Microbiology, 11 , 1158. Dar, D. , Dar, N. , Cai, L. & Newman, D.K. (2021) Spatial transcriptomics of planktonic and sessile bacterial populations at single‐cell resolution. Science, 373 , eabi4882. Doney, S.C. , Abbott, M.R. , Cullen, J.J. , Karl, D.M. & Rothstein, L. (2004) From genes to ecosystems: the ocean's new frontier. Frontiers in Ecology and the Environment, 2 , 457–468. Available from: 10.1890/1540-9295(2004)002[0457:FGTETO]2.0.CO;2 Dupont, C.L. , McCrow, J.P. , Valas, R. , Moustafa, A. , Walworth, N. , Goodenough, U. et al. (2015) Genomes and gene expression across light and productivity gradients in eastern subtropical Pacific microbial communities. The ISME Journal, 9 , 1076–1092. Available from: 10.1038/ismej.2014.198 25333462 Edgcomb, V.P. , Pachiadaki, M.G. , Mara, P. , Kormas, K.A. , Leadbetter, E.R. & Bernhard, J.M. (2016) Gene expression profiling of microbial activities and interactions in sediments under haloclines of E. Mediterranean deep hypersaline anoxic basins. The ISME Journal, 10 , 2643–2657. Available from: 10.1038/ismej.2016.58 27093045 Ewels, P. , Magnusson, M. , Lundin, S. & Käller, M. (2016) MultiQC: summarize analysis results for multiple tools and samples in a single report. Bioinformatics, 32 , 3047–3048. Available from: 10.1093/bioinformatics/btw354 27312411 Eyre, B.D. & Ferguson, A.J.P. (2009) Denitrification efficiency for defining critical loads of carbon in shallow coastal ecosystems. Hydrobiologia, 629 , 137–146. Available from: 10.1007/s10750-009-9765-1 Fabris, G. , Monahan, C. & Batley, G. (1999) Heavy metals in waters and sediments of Port Phillip Bay Australia. Marine and Freshwater Research, 50 , 503. Available from: 10.1071/MF98032 Freedman, A. (2016) Filter uncorrectable FE fasta. github.com/harvardinformatics/TranscriptomeAssemblyTools/blob/master/FilterUncorrectablePEfastq.py Frias‐Lopez, J. , Shi, Y. , Tyson, G.W. , Coleman, M.L. , Schuster, S.C. , Chisholm, S.W. et al. (2008) Microbial community gene expression in ocean surface waters. Proceedings of the National Academy of Sciences of the United States of America, 105 , 3805–3810. Available from: 10.1073/pnas.0708897105 18316740 Garnier, S. , Ross, N. , Rudis, B.B. , Filipovic‐Pierucci, A. , Galili, T. , Timelyportfolio Greenwell, B. et al. (2021) viridis ‐ Colorblind‐Friendly Color Maps for R. Zenodo. 10.5281/ZENODO.4679424 Gilbert, J.A. , Field, D. , Huang, Y. , Edwards, R. , Li, W. , Gilna, P. et al. (2008) Detection of large numbers of novel sequences in the metatranscriptomes of complex marine microbial communities. PLoS One, 3 , e3042. Available from: 10.1371/journal.pone.0003042 18725995 Grabherr, M.G. , Haas, B.J. , Yassour, M. , Levin, J.Z. , Thompson, D.A. , Amit, I. et al. (2011) Full‐length transcriptome assembly from RNA‐Seq data without a reference genome. Nature Biotechnology, 29 , 644–652. Available from: 10.1038/nbt.1883 Graf, D.R.H. , Jones, C.M. & Hallin, S. (2014) Intergenomic comparisons highlight modularity of the denitrification pathway and underpin the importance of community structure for N2O emissions. PLoS One, 9 , e114118. Available from: 10.1371/journal.pone.0114118 25436772 Graham, E.B. , Knelman, J.E. , Schindlbacher, A. , Siciliano, S. , Breulmann, M. , Yannarell, A. et al. (2016) Microbes as engines of ecosystem function: when does community structure enhance predictions of ecosystem processes? Frontiers in Microbiology, 7 , 214. Available from: 10.3389/fmicb.2016.00214 26941732 Grossmann, L. , Beisser, D. , Bock, C. , Chatzinotas, A. , Jensen, M. , Preisfeld, A. et al. (2016) Trade‐off between taxon diversity and functional diversity in European lake ecosystems. Molecular Ecology, 25 , 5876–5888. Available from: 10.1111/mec.13878 27747959 Haas, B.J. , Papanicolaou, A. , Yassour, M. , Grabherr, M. , Blood, P.D. , Bowden, J. et al. (2013) De novo transcript sequence reconstruction from RNA‐seq using the Trinity platform for reference generation and analysis. Nature Protocols, 8 , 1494–1512. Available from: 10.1038/nprot.2013.084 23845962 Harris, G. , Batley, G. , Fox, D. , Hall, D. , Jernakoff, R. , Murray, A. et al. (1996) Port Phillip environmental study final report. Dickson: CSIRO. Heggie, D.T. , Skyring, G.W. , Orchardo, J. , Longmore, A.R. , Nicholson, G.J. & Berelson, W.M. (1999) Denitrification and denitrifying efficiencies in sediments of Port Phillip Bay: direct determinations of biogenic N2 and N‐metabolite fluxes with implications for water quality. Marine and Freshwater Research, 50 , 589–596. Available from: 10.1071/MF98054 Hug, L.A. & Co, R. (2018) It takes a village: microbial communities thrive through interactions and metabolic handoffs. mSystems, 3 , e00152–e00117. Available from: 10.1128/mSystems.00152-17 29556533 Ji, Q. , Buitenhuis, E. , Suntharalingam, P. , Sarmiento, J.L. & Ward, B.B. (2018) Global nitrous oxide production determined by oxygen sensitivity of nitrification and denitrification. Global Biogeochemical Cycles, 32 , 1790–1802. Available from: 10.1029/2018GB005887 Kessler, A.J. , Roberts, K.L. , Bissett, A. & Cook, P.L.M. (2018) Biogeochemical controls on the relative importance of denitrification and dissimilatory nitrate reduction to ammonium in estuaries. Global Biogeochemical Cycles, 32 , 1045–1057. Available from: 10.1029/2018GB005908 Kopylova, E. , Noé, L. & Touzet, H. (2012) SortMeRNA: fast and accurate filtering of ribosomal RNAs in metatranscriptomic data. Bioinformatics, 28 , 3211–3217. Available from: 10.1093/bioinformatics/bts611 23071270 Kraft, B. , Jehmlich, N. , Larsen, M. , Bristow, L.A. , Könneke, M. , Thamdrup, B. et al. (2022) Oxygen and nitrogen production by an ammonia‐oxidizing archaeon. Science, 375 , 97–100. Available from: 10.1126/science.abe6733 34990242 Krueger, F. (2007) TrimGalore! https://www.bioinformatics.babraham.ac.uk/projects/trim_galore/ Kuypers, M.M.M. , Marchant, H.K. & Kartal, B. (2018) The microbial nitrogen‐cycling network. Nature Reviews. Microbiology, 16 , 263–276. Available from: 10.1038/nrmicro.2018.9 29398704 Langmead, B. , Trapnell, C. , Pop, M. & Salzberg, S.L. (2009) Ultrafast and memory‐efficient alignment of short DNA sequences to the human genome. Genome Biology, 10 , R25. Available from: 10.1186/gb-2009-10-3-r25 19261174 Li, B. & Dewey, C.N. (2011) RSEM: accurate transcript quantification from RNA‐Seq data with or without a reference genome. BMC Bioinformatics, 12 , 323. Available from: 10.1186/1471-2105-12-323 21816040 Lipsewers, Y.A. , Bale, N.J. , Hopmans, E.C. , Schouten, S. , Sinninghe Damsté, J.S. & Villanueva, L. (2014) Seasonality and depth distribution of the abundance and activity of ammonia oxidizing microorganisms in marine coastal sediments (North Sea). Frontiers in Microbiology, 5 , 472. Available from: 10.3389/fmicb.2014.00472 25250020 MacManes, M.D. (2014) On the optimal trimming of high‐throughput mRNA sequence data. Frontiers in Genetics, 5 , 13. Available from: 10.3389/fgene.2014.00013 Marlow, J. , Spietz, R. , Kim, K.‐Y. , Ellisman, M. , Girguis, P. & Hatzenpichler, R. (2021) Spatially resolved correlative microscopy and microbial identification reveal dynamic depth‐ and mineral‐dependent anabolic activity in salt marsh sediment. Environmental Microbiology, 23 , 4756–4777. Available from: 10.1111/1462-2920.15667 34346142 Marshall, A. , Longmore, A. , Phillips, L. , Tang, C. , Hayden, H. , Heidelberg, K. et al. (2021) Nitrogen cycling in coastal sediment microbial communities with seasonally variable benthic nutrient fluxes. Aquatic Microbial Ecology, 86 , 1–19. Available from: 10.3354/ame01954 Marshall, A.J. , Phillips, L. , Longmore, A. , Hayden, H.L. , Heidelberg, K.B. , Tang, C. et al. (2023) Temporal profiling resolves the drivers of microbial nitrogen cycling variability in coastal sediments. Science of the Total Environment, 856 , 159057. Available from: 10.1016/j.scitotenv.2022.159057 36174701 Martens‐Habbena, W. , Qin, W. , Horak, R.E.A. , Urakawa, H. , Schauer, A.J. , Moffett, J.W. et al. (2015) The production of nitric oxide by marine ammonia‐oxidizing archaea and inhibition of archaeal ammonia oxidation by a nitric oxide scavenger. Environmental Microbiology, 17 , 2261–2274. Available from: 10.1111/1462-2920.12677 25420929 Martin, M. (2011) Cutadapt removes adapter sequences from high‐throughput sequencing reads. EMBnet.journal, 17 , 10–12. Available from: 10.14806/ej.17.1.200 Moran, M.A. , Satinsky, B. , Gifford, S.M. , Luo, H. , Rivers, A. , Chan, L.‐K. et al. (2013) Sizing up metatranscriptomics. The ISME Journal, 7 , 237–243. Available from: 10.1038/ismej.2012.94 22931831 Onley, J.R. , Ahsan, S. , Sanford, R.A. & Löffler, F.E. (2018) Denitrification by Anaeromyxobacter dehalogenans, a common soil bacterium lacking the nitrite reductase genes nirS and nirK . Applied and Environmental Microbiology, 84 , e01985–e01917. Available from: 10.1128/AEM.01985-17 29196287 Poretsky, R.S. , Bano, N. , Buchan, A. , LeCleir, G. , Kleikemper, J. , Pickering, M. et al. (2005) Analysis of microbial gene transcripts in environmental samples. Applied and Environmental Microbiology, 71 , 4121–4126. Available from: 10.1128/AEM.71.7.4121-4126.2005 16000831 Poretsky, R.S. , Sun, S. , Mou, X. & Moran, M.A. (2010) Transporter genes expressed by coastal bacterioplankton in response to dissolved organic carbon. Environmental Microbiology, 12 , 616–627. Available from: 10.1111/j.1462-2920.2009.02102.x 19930445 Qin, W. , Amin, S.A. , Lundeen, R.A. , Heal, K.R. , Martens‐Habbena, W. , Turkarslan, S. et al. (2018) Stress response of a marine ammonia‐oxidizing archaeon informs physiological status of environmental populations. The ISME Journal, 12 , 508–519. Available from: 10.1038/ismej.2017.186 29053148 R Core Team . (2021) R: a language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing. Raes, E.J. , Karsh, K. , Kessler, A.J. , Cook, P.L.M. , Holmes, B.H. , van de Kamp, J. et al. (2020) Can we use functional genetics to predict the fate of nitrogen in estuaries? Frontiers in Microbiology, 11, 1261. Raes, J. & Bork, P. (2008) Molecular eco‐systems biology: towards an understanding of community function. Nature Reviews. Microbiology, 6 , 693–699. Available from: 10.1038/nrmicro1935 18587409 Revelle, W. (2022) psych: procedures for psychological, psychometric, and personality research. Evanston, Illinois: Northwestern University. Robertson, E.K. , Roberts, K.L. , Burdorf, L.D.W. , Cook, P. & Thamdrup, B. (2016) Dissimilatory nitrate reduction to ammonium coupled to Fe(II) oxidation in sediments of a periodically hypoxic estuary. Limnology and Oceanography, 61 , 365–381. Available from: 10.1002/lno.10220 Robertson, E.K. & Thamdrup, B. (2017) The fate of nitrogen is linked to iron(II) availability in a freshwater lake sediment. Geochimica et Cosmochimica Acta, 205 , 84–99. Available from: 10.1016/j.gca.2017.02.014 Robinson, M.D. , McCarthy, D.J. & Smyth, G.K. (2010) edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics, 26 , 139–140. Available from: 10.1093/bioinformatics/btp616 19910308 RStudio Team . (2016) RStudio: integrated development for R. Boston, MA. http://www.rstudio.com/ Sanford, R.A. , Cole, J.R. & Tiedje, J.M. (2002) Characterization and description of Anaeromyxobacter dehalogenans gen. nov., sp. nov., an aryl‐halorespiring facultative anaerobic myxobacterium. Applied and Environmental Microbiology, 68 , 893–900. Available from: 10.1128/AEM.68.2.893-900.2002 11823233 Sanford, R.A. , Wagner, D.D. , Wu, Q. , Chee‐Sanford, J.C. , Thomas, S.H. , Cruz‐García, C. et al. (2012) Unexpected nondenitrifier nitrous oxide reductase gene diversity and abundance in soils. Proceedings of the National Academy of Sciences, 109 , 19709–19714. Available from: 10.1073/pnas.1211238109 Santoro, A.E. , Buchwald, C. , Knapp, A.N. , Berelson, W.M. , Capone, D.G. & Casciotti, K.L. (2021) Nitrification and nitrous oxide production in the offshore waters of the Eastern Tropical South Pacific. Global Biogeochemical Cycles, 35 , e2020GB006716. Available from: 10.1029/2020GB006716 Santoro, A.E. , Buchwald, C. , McIlvin, M.R. & Casciotti, K.L. (2011) Isotopic signature of N2O produced by marine ammonia‐oxidizing archaea. Science, 333 , 1282–1285. Available from: 10.1126/science.1208239 21798895 Söllinger, A. , Tveit, A.T. , Poulsen, M. , Noel, S.J. , Bengtsson, M. , Bernhardt, J. et al. (2018) Holistic assessment of rumen microbiome dynamics through quantitative metatranscriptomics reveals multifunctional redundancy during key steps of anaerobic feed degradation. mSystems, 3 , e00038–18. Available from: 10.1128/mSystems.00038-18 Song, L. & Florea, L. (2015) Rcorrector: efficient and accurate error correction for Illumina RNA‐seq reads. GigaScience, 4 , 48. Available from: 10.1186/s13742-015-0089-y 26500767 Täumer, J. , Marhan, S. , Groß, V. , Jensen, C. , Kuss, A.W. , Kolb, S. et al. (2022) Linking transcriptional dynamics of CH4‐cycling grassland soil microbiomes to seasonal gas fluxes. The ISME Journal, 1–10 , 1788–1797. Available from: 10.1038/s41396-022-01229-4 Tu, Q. , Lin, L. , Cheng, L. , Deng, Y. & He, Z. (2019) NCycDB: a curated integrative database for fast and accurate metagenomic profiling of nitrogen cycling genes. Bioinformatics, 35 , 1040–1048. Available from: 10.1093/bioinformatics/bty741 30165481 Venter, J.C. , Remington, K. , Heidelberg, J.F. , Halpern, A.L. , Rusch, D. , Eisen, J.A. et al. (2004) Environmental genome shotgun sequencing of the Sargasso Sea. Science, 304 , 66–74. Available from: 10.1126/science.1093857 15001713 Wei, T. & Simko, V. (2021) R package “corrplot”: visualization of a correlation matrix. Wickham, H. , François, R. , Henry, L. & Müller, K. (2022) dplyr: a Grammar of data manipulation. Zhang, X. , Zhang, Q. , Yang, A. , Hou, L. , Zheng, Y. , Zhai, W. et al. (2018) Incorporation of microbial functional traits in biogeochemistry models provides better estimations of benthic denitrification and anammox rates in coastal oceans. Journal of Geophysical Research‐Biogeosciences, 123 , 3331–3352. Available from: 10.1029/2018JG004682 Zhao, S. , Guo, Y. , Sheng, Q. & Shyr, Y. (2014) Heatmap3: an improved heatmap package with more powerful and convenient features. BMC Bioinformatics, 15 , P16. Available from: 10.1186/1471-2105-15-S10-P16