Patterns of microbial abundance and heterotrophic activity along nitrogen and salinity gradients in coastal wetlands

Coastal wetlands are valuable aquatic ecosystems with high biological productivity, which provide services such as a reduction in nitrogen loading into coastal waters and storage of organic carbon acting as carbon dioxide sinks. The predicted rise of sea level or freshwater extractions, particularly in the arid Mediterranean biome, will salinize many coastal wetlands. However, there is considerable uncertainty about how salinization will affect microbial communities and biogeochemical processes. We determined the abundance of total prokaryotes, cyanobacteria, and viruses and quantified the heterotrophic production of prokaryotes sensitive- (predominantly Bacteria) and resistant- (predominantly Archaea) to erythromycin in 112 ponds from nine coastal wetlands. We explored the main drivers of prokaryotic abundance and heterotrophic production using generalized linear models (GLMs). The best GLM, including all the wetlands, indicated that the concentration of total dissolved nitrogen (TDN) positively affected the total abundance of prokaryotes and the heterotrophic erythromycin-resistant (ery-R) production. In contrast, heterotrophic erythromycin-sensitive (ery-S) production was negatively related to TDN. This negative relationship appeared to be mediated by salinity and virus abundance. Heterotrophic ery-S production declined as salinity and virus abundance increased. Consequently, we observed a switch from heterotrophic ery-S production towards ery-R production as salinity and virus abundance increased. Our results imply that microbial activity will change from heterotrophic bacterial-dominated processes to archaeal-dominated processes with anthropogenic nitrogen and salinization increases. However, more studies are required to link the mineralization rates of dissolved nitrogen and organic carbon with specific archaeal taxa to enable more accurate predictions on future scenarios in coastal wetlands.


Introduction
Wetlands have high biological productivity and functional diversity (Mitsch and Gosselink 2000;Mitsch et al. 2015). Among the ecosystem services they provide, wetlands reduce nutrient loading from freshwater inflow (Verhoeven et al. 2006) and have been described as "the kidneys of the landscape" (Mitsch et al. 2015). Coastal wetlands, in particular, are considered global reservoirs of organic carbon because of their elevated primary productivity and high sedimentation rates that, in the long-term, minimize carbon dioxide emissions to the atmosphere, acting as carbon sinks (Roehm 2005;Chmura et al. 2003;Geertz-Hansen et al. 2011). Therefore, microbial activities such as denitrification and organic carbon production and mineralization mediate these wetland functions and services. These aquatic Aquatic Sciences * Isabel Reche ireche@ugr.es 1 Departamento de Ecología e Instituto del Agua, Universidad de Granada, 18071 Granada, Spain ecosystems are currently facing two main threats: eutrophication and salinization. The production of nitrogen fertilizer has enabled increased food production but has drastically disrupted the nitrogen cycle by doubled its inputs to the Earth´s surface leading to eutrophication of many aquatic ecosystems (Gruber and Galloway 2008;Canfield et al. 2010). This change likely exceeds all other human effects on nutrient cycles (Gruber and Galloway 2008;Schlesinger 2009) but has received relatively less attention than the carbon cycle (Battye et al. 2017). A substantial fraction of this nitrogen, mainly derived from agricultural lands, is transported to rivers and groundwater through runoff and is lost through emissions of ammonia and other nitrogen compounds by denitrification (Schlesinger 2009;Battye et al. 2017). In many estuaries, natural and constructed wetlands modify the nitrogen loading into coastal zones, alleviating eutrophication (Walton et al. 2015). These transitional waters include saline wetlands such as natural marshes and multi-pond solar salterns (Razinkovas-Baziukas and Povilanskas 2012).
Multi-pond solar salterns are present along the coasts of the Mediterranean and Portugal since the Phoenicians, who used evaporation procedures dating from the Emperor Huang era about 2500 years B.C. to obtain salt (Baas-Becking 1931;García-Vargas and Martínez-Maganto 2006). They are excellent model systems to study microbial ecology along salinity gradients from that of seawater to salt saturation (Ventosa et al. 2014). Coastal saline wetlands also have great importance for waterbird conservation and sustainable aquaculture (Britton and Johnson 1987;Sánchez et al. 2006;Athearn et al. 2009;Walton et al. 2015). Most studies of microbial ecology in coastal wetlands have focused on extremophile microorganisms and their bioenergetic and diversity constraints (Antón et al. 1999;Pedrós-Alió et al. 2000;Oren 2001;Casamayor et al. 2002;Gasol et al. 2004;Ghai et al. 2011;Ventosa et al. 2014); but much less is known about moderately halophilic bacteria and archaea in oligo-to eusaline waters (Hahn, 2006;Herbert et al. 2015).
The discovery that archaea are not exclusively extremophiles (DeLong 1992;Fuhrman et al. 1992;Karner et al. 2001) led to studies showing that archaea are abundant, ubiquitous, and diverse, with essential roles in both 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 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. Archaeal heterotrophic production patterns along gradients of dissolved organic carbon and nitrogen, and salinity, remain completely unexplored in the literature.
Currently, most wetlands are salinizing due to climatic warming or freshwater extraction (Herbert et al. 2015;Morrissey et al. 2014;Jeppesen et al. 2015). Coastal wetlands are particularly affected since the sea-level rise introduces 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 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 (Brunet and García-Gil 1996;Laverman et al. 2007). Overall, salinization disrupts many ecosystem services provided by wetlands, such as their capacity to store organic carbon (Weston et al. 2011;Luo et al. 2019) and remove nitrogen from water (Herbert et al. 2015;Franklin et al. 2017).
The projected scenario of global change for wetlands includes wetland salinization and more nitrogen inputs. These environmental risks can interact. Nitrogen inputs can boost organic carbon mineralization and reduce 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 abundances (total prokaryotes, cyanobacteria, and viruses) and heterotrophic production of prokaryotes sensitive (predominantly Bacteria domain) and resistant (predominantly Archaea domain) to erythromycin 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 potential changes in prokaryotic heterotrophic production under future scenarios associated with global change.

Study sites
We sampled 112 saline ponds from nine coastal wetlands during the late spring-midsummer of 2011 Table S1). The coastal wetlands from west to east were: (1) the Odiel marshes (OdielM) at the mouth of the Odiel river that included natural marshes and salterns, (2) the Veta la Palma (VPalma) fishponds at the mouth of the Guadalquivir river, (3) the Cabo de Gata (CGata) salterns, (4) El Hondo (Hondo) wetlands and (5) Santa Pola (SPola) salterns at the mouth of the Vinalopó river, and (6) the Ebro Delta (EbroD) marshes and salterns at the mouth of the Ebro river in Spain, (7) the Salin-de-Giraud and Saintes-Maries-de-la-Mer marshes and salterns at the mouth of the Rhône river (Camargue) in France, (8) the Molentargius, Santa Guilla, and Santa Caterina marshes and salterns (Sardinia) in Italy, and (9) Thyna solar salterns (Sfax) in Tunisia (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 (crystallizers in solar saltpans) to oligo-mesohaline (brackish and evaporation ponds) waters (Fig. 1b), and microbial heterogeneity (Fig. 1c). Multi-pond solar salterns are excellent model systems to study microbial ecology along salinity gradients from that of seawater to salt saturation (Ventosa et al. 2014). They are natural laboratories across which microbial communities can be compared. The brackish wetlands that we studied are also crucial for sustainable aquaculture and waterfowl conservation (Britton and Johnson 1987;Walton et al. 2015). Coordinate location, sampling dates, and physical-chemical and microbiological data for each pond are in Supplementary Table S1.

Chemical analyses
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 method (Murphy and 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 high-temperature catalytic oxidation using a total organic carbon analyzer (Shimadzu TOC-V CSH). DOC samples were pre-filtered through precombusted Whatman GF/F filters (2 h at 500 °C) and acidified with H3PO4 (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
We determined chlorophyll-a (Chl-a) concentration by filtering 50-2000 ml of water through Whatman GF/F filters and extracting the pigments using 95% methanol in the dark for 24 h at 4 °C (APHA 1992). We measured the absorbance at 665 nm and 750 nm using a Perkin Elmer Lambda 40 spectrophotometer.
Samples for determining the abundances of total prokaryotes (PA) and cyanobacteria (CyA) by flow cytometry (Gasol and del Giorgio, 2000) were collected in cryovials, Orthophoto of the Odiel marshes including evaporation and crystallization ponds; Orthophoto credit: Junta de Andalucía. Consejería de Agricultura, Ganadería, Pesca y Desarrollo Sostenible (ámbito de Desarrollo Sostenible). http:// www. junta deand alucia. es/ medio ambie nte/ site/ rediam/. c Picture of GF/F glass-fiber filters to illustrate the pigment diversity. Two top filters correspond to crystallizers and the bottom filter corresponds to evaporation ponds from water samples collected at the Odiel marshes 22 Page 4 of 13 fixed with 1% paraformaldehyde and 0.05% glutaraldehyde in the dark (30 min at 4 °C), and frozen in liquid nitrogen and stored at − 80 °C until analysis. Once defrosted, the samples were diluted (≥ tenfold) with Milli-Q water to avoid the electronic coincidence of the prokaryotic cells. We counted the PA and CyA in a FACScalibur flow cytometer with a laser emitting at 448 nm and a suspension of yellow-green 1-μm latex beads (Polysciences) per sample as an internal standard. We stained the PA samples with a 10 μM DMSO solution of Sybr Green I stain (Molecular Probes) in the dark (10 min), 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).
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 ≥ 100-fold 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 the VA samples with a working solution (1:200) of SYBR Green I (10,000 × 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 data using BD CellQuest Pro software.
We measured the heterotrophic prokaryotic production using the rate of 3 H-Leucine incorporation into proteins (Smith and Azam, 1992). [ 3 H]-Leucine incorporation by microbial eukaryotes appears irrelevant in short-term incubations even with high leucine concentrations (Buesing and Gessner 2003; Miranda et al. 2007). We assumed that, using erythromycin, we inhibit most bacterial production and determine the production mostly attributable to Archaea following the procedure proposed by Yokokawa et al. (2012). Although some authors have questioned the specificity of erythromycin to inhibit some bacteria phylum in the open ocean (Frank et al. 2016), this procedure appears to have high efficiencies (ca. 80%) in waters with salinities higher than the seawater (Oren et al. 1990;Pedrós-Alió et al. 2000) and for specific groups of bacteria (Du et al. 2016). For each pond, two sets of five samples (1.5 ml) (three replicate 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-5 h. In one of the sets of 5 samples, we also added 10 μg ml −1 of erythromycin (final concentration) to inhibit most of the bacterial production (hereafter, heterotrophic erythromycin-sensitive (ery-S) production). To obtain the production attributed mainly to archaeal and the few bacterial groups that were erythromycin-resistant (hereafter, heterotrophic erythromycinresistant (ery-R) production), we subtracted from the total production (first set of samples) the ery-S production. We added TCA at a final concentration of 10% to end the incubations. In the laboratory, the samples were centrifuged (15,366g 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 autocalibrated 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 prokaryotic abundance (PA), cyanobacterial abundance (CyA), heterotrophic erythromycin-sensitive (ery-S) production, and heterotrophic erythromycin-resistant (ery-R) production. 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 avoid heteroscedasticity. In the second set of GLMs, we also included the 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 determined. Nutrients were not included as predictors because viruses are mostly dependent on the host (bacterial or archaeal) density. We considered the study site a categorical variable and 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 a predictor variable and ery-S production and ery-R production 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. We performed all analyses using Statistica 7.0 software.

Results
The study ponds represent a wide range of environmental conditions. Salinity ranged from 0.22 to 343 ppt (Table 1), spanning from oligo-to hyper-saline conditions; the highest median value was in EbroD ponds, and the lowest median value in El Hondo ponds (Fig. 2a). 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. 2c). 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 saltpans the lowest values (Fig. 2b). TDP concentrations ranged from 0.50 to 40.09 μmol-P l −1 (Table 1), with the highest median value at OdielM wetlands and the lowest median at CGata saltpans (Fig. 2d).
In the first set of GLMs, the best model for prokaryotic abundance included the TDN concentration and the wetland site (categorical) as predictors ( Table 2). The residuals of prokaryotic abundance (once the site was controlled for) were significantly and positively related to TDN (Fig. 4a). Alternative (less explanatory) models also included salinity or TDP (Table supplementary S2). The best GLM for the 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 a negative partial relationship with TDP (Table 2). In the case of heterotrophic ery-S Table 1 Ranges of basic physicochemical variables, total dissolved nitrogen (TDN), total dissolved phosphorus (TDP), dissolved organic carbon (DOC) and concentration of chlorophyll a (Chl a) and virus abundance for each set of ponds in the study wetlands The complete data set is in the Supplementary Table S1.  production, the best GLM included TDN and the wetland site as predictors ( Table 2). The residuals of heterotrophic ery-S production (once the site was controlled for) were significantly and negatively related to TDN (Fig. 4b).
Indeed, the higher the TDN concentration, the lower the heterotrophic ery-S production, irrespectively of the site. Alternative (less explanatory) models also included salinity or TDP (Supplementary Table S2). For heterotrophic ery-R production, the best GLM also included TDN and site as predictors (Table 2). In contrast to heterotrophic ery-S production, the residuals of heterotrophic ery-R production (once the site was controlled for) were significantly and positively related to TDN concentration (Fig. 4c). Indeed, the higher the TDN, the higher the production associated to archaea and bacterial groups resistant to erythromycin, 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, but excluded the Camargue site since these data were not available for that wetland. This predictor variable only affected the results previously exposed to the prokaryotic abundance and the heterotrophic ery-S production (Table supplementary S3). In this second analysis, prokaryotic abundance was affected negatively by salinity and positively by DOC concentration, while depending on the site; DOC concentration negatively affected heterotrophic ery-S production.
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 (Fig. 5). Alternative (less explanatory) models also included prokaryotic and cyanobacterial abundances (Table supplementary S4).
Finally, in the fourth set of GLMs, we included viral abundance as a predictor of heterotrophic ery-S and ery-R production (Table supplementary S5), but only considered the sites where these data were available: Odiel M, VPalma, CGata, SPola, Hondo, and EbroD. The best GLM for heterotrophic ery-S production included salinity and virus abundance as predictors (Table 4). Once the site effect was controlled, we observed a significant negative relationship between heterotrophic ery-S production and salinity (Fig. 6a), as well as viral abundance (Fig. 6b). On the other hand, the best GLM for heterotrophic ery-R production was coherent with the first set of GLMs, including all the saline wetlands; the best GLM included TDN and site as predictors (Table supplementary S5 confirms the robust relationship between TDN and heterotrophic ery-R production in the study wetlands.

Discussion
The best-generalized linear model (GLM), including all the coastal wetlands, indicated that the concentration of TDN positively affected the abundance of prokaryotes and heterotrophic ery-R production. TDN was consistently associated with higher ery-R heterotrophic production, even considering all the alternative GLMs with other predictor variables, such as DOC concentration and viral abundance (VA). In contrast, TDN negatively affected heterotrophic ery-S 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 ery-S production appeared to be mediated by salinity a b c d e f Fig. 3 Violin plots to visualize the data distribution of the abundance of prokaryotes (a), cyanobacteria (b), viruses (c), heterotrophic erythromycin-sensitive (predominantly bacteria) production (d), heterotrophic erythromycin-resistant (predominantly archaea) production (e) and chlorophyll concentration (f) for each wetland site studied. White dots are the data, black lines are the median values, and the grey shadows are the data distributions for each wetland. Note panels (a-c), and (f) have log scales and VA. Heterotrophic ery-S production declined as salinity and viruses increased; therefore, we observed a switch from heterotrophic ery-S production (predominantly bacterial production) to heterotrophic ery-R production (predominantly achaeal production) as salinity and VA increased. Archaea could be the main microorganisms processing dissolved nitrogen in the most saline and eutrophic wetlands.
In this across-wetlands study, salinity did not directly affect heterotrophic ery-R production (predominantly achaeal production), whereas total dissolved nitrogen was consistently the best predictor, regardless of the variables considered, in all the alternative generalized linear models. Historically, archaea were associated with extreme 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 broader variety of environments and perform several functions in N cycling, including nitrification and denitrification (Francis et al. 2007, Offre et al. 2013). On the other hand, our results emphasize that heterotrophic ery-R production, measured as the uptake of 3 H-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 carried out by archaea (Zumft 1997). Coastal wetlands provide optimal conditions for archaeal 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. 2007). 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 coastal eutrophication. However, a more detailed analysis of N functional genes, beyond leucine incorporation, is needed to assess the actual 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 ecosystems (e.g., Cole et al. 1988;Baines and Pace 1991;Ortega-Retuerta et al. 2008). However, in this study, we did not find a difference between the best GLMs for heterotrophic ery-R 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 ery-R production. In saline wetlands, extreme evapoconcentration (Batanero et al. 2017) and high loads of organic matter from freshwater inflows (Walton et al. 2015) result in high concentrations of DOC in oligo-to eusaline wetlands, likely ensuring plenty of 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. 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 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 results suggest that neither phosphorus nor dissolved organic carbon appears to significantly affect the prokaryotic production in the studied saline wetlands.
In the study wetlands, heterotrophic ery-S production appeared to be controlled by viral abundance (VA) and salinity (Fig. 6), providing an indirect explanation for the negative correlation between total dissolved nitrogen and heterotrophic ery-S production (Fig. 4b). The results showed that, within a given wetland complex, VA increased significantly with salinity (Table 3, Fig. 5) and negatively affected the heterotrophic ery-S production (Table 4, Fig. 6b) but did not have a significant effect on heterotrophic heterotrophic ery-R production (Table supplementary S5). These findings suggest that most of the viruses were likely bacteriophages that caused significant bacterial mortality in these extremely saline wetlands. Thus, external factors such as evaporation would trigger an increase in VA, which 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 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). On the other hand, human population growth and crop production lead to more nitrogen inputs (Mulholland et al. 2008;Mekonnen and Hoekstra, 2015). These two environmental risks can interact, affecting the potential services of 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 will change from bacterial-dominated processes to archaealdominated processes along with the increases of nitrogen and salinity. More studies linking the processing rates of dissolved nitrogen and organic carbon with specific archaeal and bacterial taxa will allow more accurate predictions in under future scenarios of global change.

Conclusions
The concentration of total dissolved nitrogen appeared to determine the abundance of heterotrophic prokaryotes and, mainly, the heterotrophic ery-R (predominantly archaea) production in the suite of coastal wetlands studied. In a b c Fig. 4 Patterns of prokaryotic abundance and production and the concentration of total dissolved nitrogen. Partial effects for relationships between the residuals of prokaryotic abundance (a), the residuals of heterotrophic erythromycin-sensitive production (b) and the residuals of heterotrophic erythromycin-resistant production (c) and the concentration of total dissolved nitrogen, 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 contrast, total dissolved nitrogen negatively affected the heterotrophic ery-S (predominantly bacterial) production. This negative relationship between total dissolved nitrogen and heterotrophic ery-S production appeared to be mediated by salinity and virus abundance. Heterotrophic ery-S production declined as salinity and viruses increased; therefore, there was a switch from heterotrophic ery-S production in wetlands with lower nitrogen content and salinity to heterotrophic ery-R production in wetlands with higher nitrogen and salinity. Archaeal taxa could be the main prokaryotes processing nitrogen in the most saline wetlands. Therefore, changes in nitrogen inputs and salinities associated with global change can alter prokaryotic functional groups and, consequently, the ecosystem services provided by coastal wetlands. Table 3 Summary of the best-generalized linear model (according to AIC) selected for the virus abundance in the sites: Odiel M, VPalma, CGata, Hondo, SPola, and EbroD (see Fig. 1 and Table 1) 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

Fig. 5
Partial effect for relationship between the residuals of virus abundance and salinity based on the best models shown in Table 3 Table  Fig. 6 Patterns between heterotrophic erythromycin-sensitive production and their chemical and biological predictors. Partial effects for relationships between salinity (a) and virus abundance (b) and the residuals of heterotrophic erythromycin-sensitive production, 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