Alternation of heterotrophic bacterial and archaeal production along nitrogen and salinity gradients in coastal wetlands

Coastal wetlands are valuable ecosystems with high biological productivity and diversity, which provide ecosystem services such as a reduction in the inputs of nitrogen into coastal waters, and storage of organic carbon, thus, acting as net carbon sinks. The rise of sea level as a consequence of climatic warming will salinize many coastal wetlands, but there is considerable uncertainty about how 20 salinization will affect microbial communities and biogeochemical processes. We analyzed prokaryotic abundance and heterotrophic bacterial and archaeal production in 112 ponds within nine coastal wetlands from the western Mediterranean coast. We determined the main drivers of prokaryotic abundance and production in these wetlands using generalized linear models (GLMs). The best GLM, including all the coastal wetlands, indicated that the concentration of total dissolved nitrogen (TDN) 25 positively affected the abundance of heterotrophic prokaryotes and heterotrophic archaeal production. In contrast, heterotrophic bacterial production was negatively related to TDN. This negative relationship appeared to be mediated by salinity and virus abundance. Heterotrophic bacterial production declined as salinity, and virus abundance, increased. We observed a switch from heterotrophic bacterial production towards heterotrophic archaeal production as salinity and virus abundance increased. Our results imply 30 https://doi.org/10.5194/bg-2020-60 Preprint. Discussion started: 3 April 2020 c © Author(s) 2020. CC BY 4.0 License.

The discovery that archaea are not exclusively extremophiles (DeLong 1992;Fuhrman et 65 al. 1992;Karner et al. 2001) led to studies showing that archaea are abundant, ubiquitous and diverse, with essential roles in carbon and nitrogen cycling (Herndl et al. 2005;Justice et al. 2012;Offre et al. 2013), including mineralizing organic matter under oxic conditions in the ocean (Herndl et al. 2005;Ingalls et al. 2006). The organo-heterotrophic nature of some archaeal groups has been shown for specific carbohydrates (Lazar et al. 2016) and dissolved proteins (Orsi et al. 2016). Archaea also 70 oxidize ammonia, and perform denitrification under anoxic or suboxic conditions (Francis et al. 2007;Offre et al. 2013) with potential implications for nitrogen removal in wetlands (You et al. 2009).
However, the prevalence of heterotrophic metabolism in archaea and its biogeochemical relevance are poorly studied. As well, archaeal heterotrophic production patterns along gradients of dissolved organic carbon and nitrogen, and salinity, remain completely unexplored in the literature. 75 Currently, most wetlands are salinizing as a consequence of climatic warming and freshwater extraction (Herbert et al. 2015;Jeppesen et al. 2015). Coastal wetlands are particularly affected since the rise of the sea level is introducing seawater high up into estuaries (Herbert et al. 2015;Schuerch et al. 2018). Salinization alters water and sediment chemistry and changes biogeochemical reactions (Herbert et al. 2015). Although there is still considerable uncertainty about how salinization affects 80 wetland biogeochemistry, it promotes the release of inorganic nitrogen (NH 4 + ) and phosphorus (PO 4 3-), with implications for internal wetland eutrophication (Ardón et al. 2013;Weston et al. 2006). Besides, salinization increases the generation of toxic sulfides (Lamers et al. 1998), which play an essential role in wetland N cycling, inhibiting the final steps of denitrification, and resulting in the emission of nitrogen oxides (Brunet andGarcía-Gil 1996, Laverman et al. 2007). Overall, salinization disrupts 85 https://doi.org/10.5194/bg-2020-60 Preprint. Discussion started: 3 April 2020 c Author(s) 2020. CC BY 4.0 License. many ecosystem services provided by wetlands, such as their capacity to store organic carbon (Weston et al. 2011;Luo et al. 2017) and remove nitrogen from water (Herbert et al. 2015;Franklin et al. 2017).
The projected scenario of climatic warming includes both wetland salinization and more nitrogen inputs. These environmental risks interact, with increased nitrogen boosting organic carbon 90 mineralization and reducing carbon storage in coastal marshes (Deegan et al. 2012), while salinization can affect microbial nitrogen removal (Franklin et al. 2017). Our work describes patterns in microbial abundance and heterotrophic bacterial and archaeal production in a broad set of semi-natural and constructed wetlands covering gradients of salinity, dissolved organic carbon, and nitrogen along the Western Mediterranean. Furthermore, we identified the main drivers of these patterns and discussed 95 potential changes in prokaryotic heterotrophic production under future scenarios.

Study sites
We sampled 112 saline ponds from nine coastal wetlands during late spring-midsummer of   Table S1). The coastal wetlands were the Odiel marshes (OdielM) at the 100 mouth of the Odiel river, the Veta la Palma (VPalma) fishponds at the mouth of the Guadalquivir river, the Cabo de Gata (CGata) salt pans, Santa Pola (SPola) and El Hondo (Hondo) at the mouth of the Vinalopó river, and the Ebro Delta (EbroD) ponds at the mouth of the Ebro river in Spain, Salin-de-Giraud and Saintes-Maries-de-la-Mer at the mouth of the Rhône (Camargue, France), the Molentargius, Santa Guilla, and Santa Caterina ponds (Sardinia, Italy), and Thyna solar salterns (Sfax, Tunisia) 105 ( Fig.1a). All the study ponds are located in semiarid or arid areas, under typical Mediterranean climatic conditions, covering large gradients of salinity from hypersaline (solar salt pans) to oligo-mesohaline (brackish) waters (Fig.1b), and microbial productivity (Fig.1c). Solar salterns are often multi-pond systems connected for commercial salt production, which creates a strong salinity gradient from evaporation ponds to crystallizer ponds (Fig.1b). They are natural laboratories across which microbial 110 communities can be compared (Fig.1c). The brackish wetlands that we studied are also crucial for https://doi.org/10.5194/bg-2020-60 Preprint. Discussion started: 3 April 2020 c Author(s) 2020. CC BY 4.0 License. sustainable aquaculture and waterfowl conservation (Britton & Johnson 1987;Walton et al. 2015).
Complete location, sampling dates, and physical-chemical and microbiological details for each pond are in Supplementary Table S1.

Chemical analyses 115
We recorded salinity with a multi-parameter probe (HANNA HI 9828). Water samples with salinity higher than 70 ppt were diluted with Milli-Q water until they were in the operating range of the probe.
We measured total nutrient concentrations in unfiltered samples, while we filtered the samples through pre-combusted 0.7-µm pore-size Whatman GF/F glass-fiber filters for dissolved nutrient analysis. We measured total phosphorus (TP), and total dissolved phosphorus (TDP) using the molybdenum blue 120 method (Murphy & Riley, 1962) after persulfate digestion (30 min, 120ºC). We analyzed total nitrogen (TN) and total dissolved nitrogen (TDN) by high-temperature catalytic oxidation (Álvarez-Salgado and Miller 1998) using a total nitrogen analyzer (Shimadzu TNM-1), and potassium-nitrate standards. We analyzed dissolved organic carbon (DOC) concentration as non-purgeable organic carbon by hightemperature catalytic oxidation using a total organic carbon analyzer (Shimadzu TOC-V CSH). DOC 125 samples were pre-filtered through pre-combusted Whatman GF/F filters (2 h at 500ºC) and acidified with H 3 PO 4 (final pH<2). We calibrated the instrument using a four-point standard curve of potassium hydrogen phthalate. We purged the samples with phosphoric acid for 20 min, and we set up three to five injections for each sample.

Biological analyses 130
We determined chlorophyll-a (Chl-a) concentration by extracting Whatman GF/F filters in 95% methanol in the dark for 24 h at 4 ºC (APHA 1992) through which 50 to 2000 ml of water was filtered.
Absorbance was measured at 665 nm and 750 nm using a Perkin Elmer Lambda 40 spectrophotometer.
Samples for determining the abundances of prokaryotes (PA) and cyanobacteria (CyA) by flow cytometry (Gasol and del Giorgio, 2000) were collected in cryovials, fixed with 1% paraformaldehyde 135 and 0.05% glutaraldehyde in the dark (30 min at 4 ºC), then frozen in liquid nitrogen and stored at -80 ºC until analysis. Once defrosted, the samples were diluted (≥10-fold) with Milli-Q water to avoid the https://doi.org/10.5194/bg-2020-60 Preprint. Discussion started: 3 April 2020 c Author(s) 2020. CC BY 4.0 License. electronic coincidence of the prokaryotic cells. We counted the PA and CyA in a FACScalibur flow cytometer with a laser emitting at 448nm and a suspension of yellow-green 1-µm latex beads (Polysciences) per sample as an internal standard. The PA samples were stained in the dark (10 min) 140 with a 10 µM DMSO solution of Sybr Green I stain (Molecular Probes), run at low speed for 2 min, and detected in bivariate plots of side scatter (SSC) vs. FL1 (Green fluorescence). The CyA samples were run at high speed for 4 min and detected in bivariate plots of side scatter (SSC) vs. FL3 (Red fluorescence). We obtained the heterotrophic prokaryotic abundance (HPA) by subtracting the cyanobacteria abundance (CyA) from prokaryotic abundance (PA). 145 Samples for virus abundance (VA) were collected in cryovials, fixed with glutaraldehyde at a final concentration of 0.5 % (15 min at 4ºC) in the dark, frozen in liquid nitrogen, and stored at -80ºC until analysis by flow cytometry (Brussard et al., 2010). Before analysis, we diluted the samples ≥ 100fold with TE-buffer pH 8.0 (10 mM Trishydroxymethyl-aminomethane, 1 mM ethylenediaminetetraacetic acid) to avoid the electronic coincidence in virus particle counts. We stained 150 the VA samples with a working solution (1:200) of SYBR Green I (10,000X concentrate in DMSO, Molecular Probes) for 10 min in the dark and then kept them at -80º until counting. We added fluorescent microspheres (FluoSpheres carboxylate modified yellow-green fluorescent microspheres; 1.0 µm diameter) as an internal standard. We acquired the data in log mode and detected their signature in bivariate plots of side scatter (SSC) vs. FL1 (Green fluorescence). We processed the flow cytometry 155 data using BD CellQuest Pro software.
Heterotrophic bacterial production (BP) and archaeal production (ArP) were measured by 3 H-Leucine incorporation into proteins (Smith and Azam, 1992). The heterotrophic archaeal activity was measured using erythromycin to inhibit bacterial activity, as proposed by Yokokawa et al. (2012).
Although the specificity of erythromycin to selectively inhibit bacteria has recently been questioned for 160 the open ocean (Frank et al. 2016), it appears to have better efficiencies (ca. 80%) in water of higher salinity (Oren et al. 1990;Pedrós-Alió et al. 2000) and for specific functional groups as nitrifiers, particularly Firmicutes (Du et al. 2016). For each pond, two sets of three replicate (1.5 ml) and two trichloroacetic-acid (TCA, 50%)-killed (final concentration10%) samples were incubated at in situ temperature with 54.6 nM or 58.4 nM leucine (1:2 hot: cold v/v) for 2 to 5 h. One of each set of 5 samples also contained 10 µg ml -1 of erythromycin (final concentration) to inhibit bacterial production, and obtain heterotrophic production attributed mostly to archaea. We added TCA at a final concentration of 10% to end the incubations. In the laboratory, the samples were centrifuged (15366 g for 10 min), rinsed with TCA (5%), vortexed, and centrifuged again. Finally, we added 1.5 ml of liquid scintillation cocktail (Ecoscint A) to each sample and determined the incorporated radioactivity using an 170 auto-calibrated scintillation counter (Beckman LS 6000 TA).

Statistical analyses
To determine the main drivers of the observed microbial patterns, we carried out four sets of generalized linear models (GLMs). In the first set of GLMs, the dependent microbial variables were heterotrophic prokaryotic abundance (HPA), cyanobacterial abundance (CyA), heterotrophic bacterial 175 production (BP), and heterotrophic archaeal production (ArP). The predictor variables selected were site (i.e., each wetland complex as a categorical variable), and as continuous variables, salinity, total dissolved nitrogen (TDN), and total dissolved phosphorus (TDP). When necessary, log and square-root transformations were applied to the dependent and predictor variables to improve the fit to a normal distribution, and to avoid heteroscedasticity. In the second set of GLMs, we also included the 180 concentration of DOC as a predictor variable but did not include the Camargue site since this variable was not measured there. In the third set of GLMs, we determined the best predictors of virus abundance (dependent variable) in the six sites (Odiel M, VPalma, CGata, SPola, Hondo, and EbroD) where this variable was available. Nutrients were not included as predictors because viruses are mostly dependent on host (bacterial or archaeal) density. We considered the study site as a categorical variable, and 185 salinity, PA, and CyA as continuous variables. To determine the potential effect of viruses on the other microbial components, we performed the fourth set of GLMs in which VA is considered as a predictor variable and BP and ArP as dependent variables. We selected the model with the smallest value of the Akaike Information Criterion (AIC) as the best model for each dependent variable considered. More details of alternative models with ∆AIC < 2.0 are in the Supplementary Tables S2-S6. All analyses were 190 performed using Statistica 7.0 software. https://doi.org/10.5194/bg-2020-60 Preprint. Discussion started: 3 April 2020 c Author(s) 2020. CC BY 4.0 License.

Results
The study ponds represented a wide range of environmental conditions. Salinity ranged more than four orders of magnitude, from 0.22 to 343 ppt (Table 1), spanning oligo-to hyper-saline conditions; the highest median value was in EbroD ponds, and the lowest median value in El Hondo ponds (Fig. 2a). 195 DOC concentrations ranged from 0.24 to 5.76 mmol-C l -1 (Table1), with the highest median concentration in EbroD ponds (as for salinity) and the lowest in CGata ponds (Fig. 2b). TDN concentrations ranged from 0.02 to 0.62 mmol-N l -1 (Table 1), with EbroD ponds also showing the highest median TDN concentration and SPola salt pans the lowest values (Fig. 2c). TDP concentrations ranged from 0.50 to 40.09 µmol-P l -1 (Table 1), with the highest median value at OdielM wetlands and 200 the lowest median at CGata salt pans (Fig. 2d).
Biological variables also ranged widely. Prokaryotic heterotrophic abundances (PHA) ranged almost three orders of magnitude, from 0.35 x10 6 to 252 x 10 6 cells ml -1 (Supplementary Table S1), with the highest median value in Tyna salt pans and the lowest in CGata salt pans (Fig. 3a). The range in cyanobacterial abundances was even higher, ranging from 0.01 to 38105 (x10 3 ) cells ml -1 205 (Supplementary Table S1; Fig. 3b). Heterotrophic bacterial production ranged from 19 to 3007 pmoles of leucine l -1 h -1 (Supplementary Table S1; Fig. 3c) with the highest median value at SPola salt pans and the lowest at Thyna salt pans. Heterotrophic archaeal production ranged from 56 to 2290 pmoles of leucine l -1 h -1 with the highest median value at EbroD ponds and the lowest at Camargue ponds. Table S1; Fig. 3d). Chlorophyll-a concentrations ranged more than 10000-fold from 210 0.04 to 617.41 µg l -1 (Table1), with the highest median value at Thyna salt pans and the lowest at EbroD ponds ( Fig. 3e).

(Supplementary
In the first set of GLMs, the best model for heterotrophic prokaryotic abundance included the TDN concentration and the wetland site (categorical) as predictors ( Table 2). The residuals of heterotrophic prokaryotic abundance (once the site was controlled for) were significantly and positively 215 related to TDN (Fig. 4 a). Alternative (less explanatory) models also included salinity or TDP (Table   supplementary S2). The best GLM for cyanobacterial abundance included the site as a categorical variable and TDN and TDP as continuous variables, with a significant positive partial relationship with TDN and negative partial relationship with TDP (Table 2). In the case of heterotrophic bacterial https://doi.org/10.5194/bg-2020-60 Preprint. Discussion started: 3 April 2020 c Author(s) 2020. CC BY 4.0 License. production, the best GLM included TDN and the wetland site as predictors ( Table 2). The residuals of 220 heterotrophic bacterial production (once the site was controlled for) were significantly and negatively related to TDN (Fig. 4 b). Indeed, the higher the TDN concentration, the lower the bacterial production, irrespectively of the site. Alternative (less explanatory) models also included salinity or TDP (Supplementary Table S2). For archaeal heterotrophic production, the best GLM also included TDN and site as predictors (Table 2). In contrast to heterotrophic bacterial production, the residuals of 225 heterotrophic archaeal production (once the site was controlled for) were significantly and positively related to TDN concentration (Fig. 4 c). Indeed, the higher the TDN, the higher the archaeal production, irrespective of site. Alternative (less explanatory) models also included salinity or TDP (Supplementary   Table S2).
In the second set of GLMs, we also included the concentration of DOC as a predictor variable 230 but excluded the Camargue site since these data were not available. This predictor variable only affected the results previously exposed to the heterotrophic prokaryotic abundance and the heterotrophic bacterial production (Table supplementary S3). In this second analysis, heterotrophic prokaryotic abundance was affected negatively by salinity and positively by DOC concentration, while depending on the site, DOC concentration negatively affected bacterial production. 235 In the third set of GLMs, we determined the best GLM for the virus abundance, excluding nutrients as predictors (Table 3). This model included salinity as a continuous predictor and site as a categorical predictor. We observed a positive relationship between the residuals of virus abundance (once the site was controlled for) and salinity ( Figure 5). Alternative (less explanatory) models also included heterotrophic prokaryotic and cyanobacterial abundances (Table supplementary S4). 240 Finally, in the fourth set of GLMs, we included viral abundance as a predictor of heterotrophic bacterial and archaeal production (Table supplementary S5), but only considered the following sites: Odiel M, VPalma, CGata, SPola, Hondo, and EbroD. The best GLM for heterotrophic bacterial production included salinity and virus abundance as predictors (Table 4). Once the site effect was controlled for, we observed a significant negative relationship between heterotrophic bacterial 245 production and salinity (Figure 6a), as well as viral abundance (Figure 6b). On the other hand, the best GLM for heterotrophic archaeal production was coherent with the first set of GLMs, which included all https://doi.org/10.5194/bg-2020-60 Preprint. Discussion started: 3 April 2020 c Author(s) 2020. CC BY 4.0 License. the saline wetlands; the best GLM included TDN and site as predictors (Table supplementary S5). This result confirms the robust relationship between TDN and heterotrophic archaeal production in the study wetlands. 250

Discussion
The best generalized linear model (GLM), including all the coastal wetlands, indicated that the concentration of TDN positively affected the abundance of heterotrophic prokaryotes and heterotrophic archaeal production. TDN was consistently associated with greater archaeal heterotrophic production, even considering the alternative GLMs with other predictor variables, such as DOC concentration and 255 viral abundance (VA). In contrast, TDN affected negatively heterotrophic bacterial production, but this result changed when DOC and VA were also considered as predictor variables in alternative GLMs.
The negative relationship between TDN and heterotrophic bacterial production appeared to be mediated by salinity and VA. Heterotrophic bacterial production declined as salinity and viruses increased; therefore, there was a switch from heterotrophic bacterial production to heterotrophic archaeal 260 production as salinity and VA increased. Archaea appeared to be the main microorganisms processing TDN in the most saline wetlands.
In this across-wetlands study, salinity did not directly affect heterotrophic archaeal production, whereas total dissolved nitrogen was consistently the best predictor, irrespectively of the variables considered, in the generalized linear models. Historically, archaea were associated with extreme 265 environments, such as salterns, where salt concentrations are close to saturation (Oren, 1994(Oren, , 2011Antón et al. 1999). Now, we know that archaea occur in a wide range of environments and perform diverse functions in N cycling, including nitrification and denitrification (Francis et al. 2007, Offre et al. 2013). In our study, ammonia oxidation by archaea during nitrification likely is not a significant process due to the high concentrations of dissolved nitrogen in most wetlands, particularly in the Ebro Delta 270 (Fig. 2). Ammonia-oxidizing archaea (AOA) are abundant in marine systems (Könneke et al. 2005, Wuchter et al. 2006, Lam et al. 2007), some lakes (Jiang et al. 2009), and wetlands (Sims et al. 2010).
The relative contribution of AOA to ammonium oxidation is still uncertain, but they appear to dominate On the other hand, our results emphasize that heterotrophic archaeal activity, measured as the uptake of 3H-leucine, was coupled to total dissolved nitrogen concentration (Fig. 4c). Denitrification is primarily a heterotrophic, microaerophilic, or anoxic process widespread along salinity gradients, which can be done by archaea (Zumft, 1997). Coastal wetlands provide optimal conditions for archaeal 280 denitrification due to their low oxygen concentrations, as a result of high salinities, and high concentrations of organic matter and nitrate (Seitzinger, 1988;Seitzinger et al. 2006). Furthermore, Orsi et al. (2016 demonstrated the organo-heterotrophic nature of some archaeal groups, which take up dissolved proteins. Therefore, heterotrophic archaeal metabolism via denitrification and protein assimilation affects the nitrogen cycle ameliorating eutrophication of recipient coastal systems. 285 However, a more detailed analysis of functional genes of N processing, beyond heterotrophic production, is needed to quantitatively assess the contribution of archaea to nitrogen cycling in coastal wetlands. Early works emphasized the importance of organic carbon derived from phytoplankton as a critical driver of heterotrophic prokaryotic (bacterial + archaeal) production in diverse aquatic 290 ecosystems (e.g., Cole et al. 1988;Baines and Pace 1991;Ortega-Retuerta et al. 2008). However, in this study, we did not find difference between the best GLMs for heterotrophic archaeal production with and without the concentration of dissolved organic carbon (DOC) (Table S3). Indeed, the total pool of DOC concentration was not a predictive variable at least for heterotrophic archaeal production. In saline wetlands, extreme evapoconcentration (Batanero et al. 2017) and high loads of organic matter from 295 freshwater inflows (Walton et al. 2015) result in high concentrations of DOC in oligo-to eusaline wetlands, likely ensuring plenty DOC for heterotrophic metabolism. Similarly, TDP was not included as a predictor variable in the best GLMs for the different microbial variables, except for the abundance of cyanobacteria (CyA). Bacterioplankton growth is assumed to be limited by the availability of phosphorus in many aquatic ecosystems (Cotner et al. 1997;Smith and Prairie 2004;Ortega-Retuerta et 300 al. 2007). However, more recently, Cotner et al. (2010) demonstrated that bacteria exhibit high flexibility in their P content and stoichiometry, questioning the ranges of P-limitation. Our data suggest https://doi.org/10.5194/bg-2020-60 Preprint. Discussion started: 3 April 2020 c Author(s) 2020. CC BY 4.0 License. that neither phosphorus nor dissolved organic carbon appear to affect significantly the prokaryotic activity in the full range of saline wetlands studied.
In the study wetlands, heterotrophic bacterial production (BP) appeared to be controlled by viral 305 abundance (VA) and salinity (Fig. 6), providing an indirect explanation for the negative correlation between total dissolved nitrogen and BP (Fig. 4b). The results showed that, within a given wetland complex, VA increased significantly with salinity (Table 3, Fig. 5), and negatively affected the BP (Table 4, Fig. 6b), but did not have a significant effect on heterotrophic archaeal production (Table   supplementary S5). These findings suggest that most of the viruses were likely bacteriophages that 310 caused significant bacterial mortality in these extremely saline wetlands. Thus, external factors such as evaporation would trigger an increase in VA, which in turn would increase bacterial mortality, while archaeal communities would be more resilient to such changes. Guixa-Boixereu et al. (1996) observed the highest prokaryotic losses by viral lysis at the highest salinities in solar salterns. Another complementary explanation is that as the salinity increases, there is a coupled decline in oxygen 315 solubility, promoting suboxic conditions. High salinity and low oxygen may be more energetically favourable conditions for some halophilic archaeal groups (Oren, 2001).
The projected scenario of climatic warming includes wetland salinization due to sea-level rise and decreases of freshwater discharges from rivers, which will allow marine waters further into the estuaries, salinizing the surrounding wetlands (Herbert et al. 2015;Schuerch et al. 2018). Human 320 population growth and crop production are leading to more nitrogen inputs (Mulholland et al. 2008;Mekonnen & Hoekstra, 2015). These two environmental risks can interact, affecting the potential services of these coastal wetlands. Eutrophication can boost the mineralization of organic carbon, reducing the carbon storage capacity of coastal marshes (Deegan et al. 2012), and salinization can affect the bacterial capacity to remove nitrogen (Franklin et al. 2017). Our results imply that microbial activity 325 will change from bacterial dominated processes to archaeal-dominated processes along with the increases of nitrogen and salinity. More studies linking the processing rates of dissolved nitrogen and organic carbon with specific archaeal taxa will allow more accurate predictions in the next scenarios of global change.

Conclusions
The concentration of total dissolved nitrogen appeared to determine the abundance of heterotrophic prokaryotes and, mainly, the heterotrophic archaeal production in the suite of coastal wetlands studied.
In contrast, total dissolved nitrogen affected negatively to heterotrophic bacterial production. We think that this negative relationship between total dissolved nitrogen and heterotrophic bacterial production 335 appeared to be mediated by salinity and virus abundance. Heterotrophic bacterial production declined as salinity and viruses increased; therefore, there was a switch from heterotrophic bacterial production in wetlands with lower nitrogen content and salinity to heterotrophic archaeal production in wetlands with higher nitrogen and salinity. Archaea appeared to be the main prokaryotes processing nitrogen in the most saline wetlands. Therefore, changes in nitrogen inputs and salinities associated with global change 340 can alter prokaryotic functional groups and, consequently, the ecosystem services provided by coastal wetlands.

Data availability
Additional figures and tables can be found in the supplementary information.

Competing interests
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.      production and (c) archaea production and TDN, based on the best models shown in Table 2. In each case the Y variable represents the residuals taken from the selected model after fitting the other predictor variables. https://doi.org/10.5194/bg-2020-60 Preprint. Discussion started: 3 April 2020 c Author(s) 2020. CC BY 4.0 License. Figure 5. Partial effect for relationship between virus abundance and salinity based on the best models shown in Table 3. https://doi.org/10.5194/bg-2020-60 Preprint. Discussion started: 3 April 2020 c Author(s) 2020. CC BY 4.0 License. Figure 6. Linear regressions between bacterial production and their biological and physico-chemical predictors. Partial effects are shown for relationships between (a) salinity (b) virus abundance, based on the best models shown in Table 4. In each case the Y variable represents the residuals taken from the selected model after fitting the other predictor variables.   https://doi.org/10.5194/bg-2020-60 Preprint. Discussion started: 3 April 2020 c Author(s) 2020. CC BY 4.0 License. Table 3. Summary of the best-generalized linear model (according to AIC) selected for the virus abundance in the sites: Odiel M, VPalma, CGata, SPola, Hondo, and EbroD (see Figure 1 and Table 1). 730 Columns show the estimates, the standard errors, the Wald statistics, and the p-values for the predictor variables. The complete set of alternative (less explanatory) models with δAIC< 2.0 are shown in Supplementary Table S4