Using Stable Carbon and Nitrogen Isotopes to Investigate the Impact of Desalination Brine Discharge on Marine Food Webs

Stable isotope ratios were used to trace the impact of anthropogenically derived brine from desalination plants on organisms at different trophic levels (primary producers and consumers) along the highly urbanized ultra-oligotrophic Israeli coast (southeast Mediterranean). Primary producer and consumer organisms were collected from two sampling stations at two desalination plants sites: an “Impacted station”, near the brine discharge outlets, and a “Control station” situated further offshore to the impacted zone. δ(_^13)C and δ(_^15)N values of both producers and consumers displayed minor variations between the impacted and control stations, indicating little effect of brine discharge on the coastal trophic structure. The coastal δ(_^15)N values were generally higher than those of similar pelagic communities of the southeastern Mediterranean. These were particularly high in benthic invertebrates and benthic carnivores (fish) from the southern site, where high anthropogenically N loads from ground water amelioration are discarded alongside the brine. The observed differences in the δ(_^15)N of the benthic components between the two study sites suggest that brine derived density plumes from desalination plants are a possible vector of nutrients to benthic communities. The results indicate that the benthic components were the most sensitive group to anthropogenic derived N pollution, and provide insight into site-specific processes.


INTRODUCTION
Ecosystems structure and functioning are undergoing profound changes due to the increase in anthropogenic impacts on natural environments. One example of anthropogenic influence is the discharge of pollutants from infrastructures along the shoreline, such as desalination plants and groundwater amelioration facilities. The increasing need for potable water, in conjunction with technological advances, has transformed seawater desalination into an exponentially growing industry comprising 59% of global desalination capacity (WWAP, 2012). The main by-product of all seawater reverse osmosis (SWRO) desalination processes is large quantities of concentrated brine, discharged back to the marine coastal environment (Lattemann and Höpner, 2008). SWRO plants discharge brine with nearly twice the salt concentration of the ambient seawater, as well as chemicals used during the SWRO process, e.g., coagulants, biocides, neutralizers, antiscalants, cleaning solutions, and pH and hardness adjustors (NRC, 2008;UNEP, 2008;Spiritos and Lipchin, 2013). The temporal and spatial dispersion pattern of the brine differs among sites and seasons, in accordance with the discharge technology of the plant, local current regime, and annual physicalchemical characteristics of the water column (Dawoud, 2012). Furthermore, the brine flow as a density plume that sink and disperse near the bottom (Biton et al., 2008). The environmental and ecological impacts of desalination plants on the marine ecosystem have been poorly documented (Roberts et al., 2010). Yet, the elevated salinity of the brine at discharge locations may impact marine life and water quality (Belkin et al., 2015;Frank et al., 2017), and lead to alterations in biodiversity and succession, as well as impact food web and ecosystem functioning and structure (Zingone and Oksfeldt Enevoldsen, 2000). Limited information has been published to date, either locally or globally, on the impact of desalination discharge on coastal microbial communities, phyto-and bacterioplankton, comprising the foundation of the aquatic food webs (Drami et al., 2011;Van der Merwe et al., 2014;Belkin et al., 2015Belkin et al., , 2017Frank et al., 2017). Studies on brine impact were performed on benthos (Sánchez-Lizaso et al., 2008;Riera et al., 2012;Begher Nabavi et al., 2013), as well as fish (Iso et al., 1994). However, the changes occurring to the structure and functioning of an entire food web under desalination stresses have not been studied.
Another source of anthropogenic stress on coastal marine ecosystems is the discharge of groundwater amelioration products. These products may affect coastal food webs in multiple ways (Lee and Kim, 2007;Zhang and Mandal, 2012;Amato et al., 2016). Groundwater amelioration that contains concentrates of high-level nutrients, rather than salt, can potentially impact coastal food webs (Amato et al., 2016). Submarine groundwater discharge is a significant source of N to coastal waters, regardless of the level of human impact on the adjacent land (e.g., Amato et al., 2016;Cho et al., 2018) and is a particularly important source of nutrients in oligotrophic regions (Weinstein et al., 2011).
Food webs dynamics include the exchange of organic matter between organisms of an ecosystem, and energy flow from primary producers to top predators (Krumins et al., 2013). Stable isotopes (SI) analysis has been commonly used to study trophic pathways in various ecosystems (Layman et al., 2012;Bongiorni et al., 2015), and is based on the fractionation of stable isotopes in animal tissue due to metabolic processes. Furthermore, SI analyses are used to quantitatively determine the origin and distribution of anthropogenic discharges to the sea, and utilized to understand the community-wide aspects of coastal trophic structure under such stressors (e.g., Layman et al., 2007Layman et al., , 2012Turner et al., 2010;Mancinelli and Vizzini, 2015). The stable carbon ratio is commonly used to study the carbon source in the diet composition, and carbon flow to consumers (Post, 2002). This measurement is useful in coastal waters which receive a complex mix of organic sources of both terrestrial and marine origin in particular (Bongiorni et al., 2015). The stable nitrogen isotope ratio is increasingly employed as a sensitive indicator of N sources in many ecosystems (McClelland et al., 1997;Vizzini and Mazzola, 2006;Orlandi et al., 2014;Graniero et al., 2016). Natural and anthropogenic N sources (the latter including fertilizer, sewage, aquaculture and manure) differ in terms of their δ 15 N value, and, consequently, values of benthic communities can reflect the relative contribution of these different sources (McClelland et al., 1997;Vizzini and Mazzola, 2006;Orlandi et al., 2017).
Five SWRO desalination plants are currently operating along Israel's 190 km long coastline with additional plants planned for the near future. Three of the five plants are among the largest in the world. The plants are located adjacent to power stations and dispose of the brine at the shoreline along with the power plant's cooling waters. Besides limited compliance monitoring and a number of studies focusing on specific taxa (Drami et al., 2011;Belkin et al., 2015Belkin et al., , 2017Frank et al., 2017), to date, no study has been conducted on the influence of the brine discharge, from the desalination plants, on the entire local food web, and only limited studies have been reported worldwide (e.g., Ozair et al., 2018).
In this study, we use stable isotopic characterization to address the impact of the brine discharge from desalination plants on the food web. We hypothesized that the concentrated discharge of brine, in the form of density plumes, mostly impacts the benthic communities, as well as the lower trophic levels. This is the first attempt to characterize the coastal food web structure of the southeastern (SE) Levantine coast via N and C stable isotopes. The regional expansion of freshwater shortage, due to climate change, and the resulting increase of desalination production and discharge of brine into coastal water can impose multiple stresses on coastal food webs. To that end, the ultra-oligotrophic coastal waters of Israel (Krom et al., 1993;Kress et al., 2014), in which severe N deficiency was recently recognized to limit primary production (Rahav et al., 2018) offers unique insight to the effect of brine discharge in a highly sensitive ecosystem.

Study Sites
In this study, we focused on two desalination plants along the Israeli coast employing SWRO methods. Ashqelon, in southern Israel (Figure 1, 31.63 • N, 34.51 • E), and Hadera, located midway along the Israeli coast (Figure 1, 32.46 • N, 34.88 • E), nearby a river outlet and thus, affected by terrestrial input (Careddu et al., 2015). The facility in Ashqelon, located near a power station, has been operating since 2006, and has a capacity of 330,000 m 3 of water production per day. The plant produces around 13% of the country's domestic consumer demand, at one of the world's lowest cost for desalinated water. Groundwater amelioration products from groundwater drilling, located near the Ashqelon plant, are discharged alongside the brine from the desalination plant. NO − 3 concentrations of the by-products from the groundwater amelioration are hundred times higher than typical coastal levels (Israel Electric Company [IEC], 2017). The desalination plant in Hadera has been operational since 2009 and is situated inside the Orot Rabin Power Station compound. Its designed production is 350,000 m 3 per day, and it is the largest SWRO-based desalination plant in the world. Both plants discharge the brine at the shoreline, at sea level. FIGURE 1 | Study sites, with the two desalination plants, Ashqelon in circles, and Hadera in triangles, solid marker for "impacted" station, and hollow marker for "control" station. Bathymetry is taken from Hall et al. (2017).

Data Collection
The majority of the samples analyzed for this study were collected during cruises conducted offshore of the Ashqelon and Hadera desalination plants. Four cruises were conducted in January, April and September 2017, and January 2018, aboard EcoOcean's R/V Mediterranean Explorer. Data were collected at both the discharge point (brine plume -"impacted station"), and at an ambient environment ("control station") outside the plume which was delineated by salinity measurements using a Hydrolab MS5 Multiprobe and located as far as ∼5 km from the brine discharge point ("stations, " Figure 1). We collected specimens from five different guild groups: zooplankton, phytoplankton, fish, benthic invertebrates, and benthic primary producers. Zooplankton were collected during the night using a 200 µm zooplankton Bongo-net, towed horizontally for approximately 15 min. Two tows were conducted at the "control" station at the surface (at a depth 1-2 m) and at a depth of 8 m, and one tow was conducted at the "impacted" station at surface depth only. Specimens were isolated and preserved in 96% ethanol. Phytoplankton sampling was performed by filtering 10 L of surface seawater collected during the cruises at each station. The water was filtered on glass fiber filters (GF/D -2.7 µm mesh size for micro-, and GF/F -0.7 µm mesh size for nanophytoplankton) pre-combusted at 550 • C for 4 h. Following rinsing of the samples with distilled water, the filters were frozen for future analysis. Blank filters were analyzed for each cruise. Various species of fish with different feeding habits were captured with the help of fishermen near the Ashqelon site, in addition to one herbivore fish from Hadera, and were frozen after their collection. Those samples were not acidified (10% HCL), as recommended by Brodie et al. (2011), since it is shown to significantly alter δ 15 N values. Due to fishing limitation at the study sites, all fish were collected from the "control" stations only. Benthic primary producers and benthic invertebrates were collected by diving with the collaboration of the Israel Electric Company (IEC) marine biology team, from both Hadera and Ashqelon Coal Jetty pillars, at the "impacted" stations (maximum depth of 5 m). The "control" stations from which the benthos were retrieved, were located one to two km from the "impact" stations at depths of 12 and 26 m, respectively. These "control" stations are therefore, situated within the general impact of the brine plume (Israel Electric Company [IEC], 2017). The soft tissues of benthic primary producers and benthic invertebrates were separated manually, but were not acidified. Their C/N ratios and δ 13 C values may therefore, include carbonates. We note that even minor amounts of carbonates could lead to a significantly enriched δ 13 C in such samples, but not in δ 15 N values, which tend to be negatively biased (Jacob et al., 2005). All fish and benthos taxa data can be found in the Supplementary Information. All samples were washed with distilled water and dried over night at 60 • C, sealed, and sent to analysis. In order to evaluate the effect of lipids on the carbon isotopic values of zooplankton, fish, and most benthos taxa (Supplementary Material) we used the relationship between C/N ratios and δ 13 C ratio . Values of %C and %N required for calculating C/N ratios were measured during the SI analysis and were used to normalize the raw δ 13 C values for lipid concentration. Several benthic invertebrates, namely bryozoan samples from the genus Watersipora sp. and the stony coral Oculina patagonica were not normalized to lipid concentration, due to their high C/N ratio and δ 13 C values. Those samples are suspected to include carbonates (were not acid washed prior to the SI analysis).

Stable Isotope (SI) Analyses
The stable isotope composition of C and N provides indication of the carbon source and trophic level, respectively, from the individual to the community level. Dry samples of organic solids were weighed into tin capsules. None of the samples were acid washed to remove carbonates prior to SI analysis. Calibrated internal standards were prepared with every batch of samples for normalization of the data. The isotopic composition of organic carbon (and nitrogen) was determined by the analysis of CO 2 (and N 2 ) continuous-flow produced by combustion on a Carlo Erba NC2500 connected on-line to a DeltaV isotope ratio mass spectrometer coupled with a ConFlo III interface. For detailed methodology, see Brodie et al. (2011). All SI analyses were performed at Cornell University Stable Isotope Laboratory.

Trophic Level Calculation
Measured isotope ratios are reported in the δ-notation, i.e., as the deviation in per mill ( ) from the internationally accepted standard: where, R represents the 15 N/ 14 N or 13 C/ 12 C ratio. The analytical precision for the in-house standard was ± 0.04 [1σ] for both δ 13 C and δ 15 N. Trophic position (TP) for a given species, or guild group, was estimated, according to the equation proposed by Post (2002): where, is the mean enrichment factor for each trophic position increase, and λ is the trophic level based on the study site (Post, 2002). We used a mean enrichment factor of 3.4 and the trophic level base of 2 (Post, 2002) for all the calculations which are supported by our results. For the grazer gastropods Monodonta turbinata a mean δ 15 N value of 2.8 ± 0.8 was used as a baseline for TP calculation in the study area (McCutchan et al., 2003;Fanelli et al., 2015). The only grazer gastropod that was measured during this study, the Littorina puncatata yielded a mean δ 15 N value of 4.1 ± 0.1 . A similar value was measured for a single omnivore fish, Siganus sp. We, thus, use the mean δ 15 N value of 4.1 ± 0.1 as a baseline for our TP calculation.

Statistical and Multivariate Analyses
For the statistical analyses, we pooled all our collected taxa into five guilds: fish, benthic invertebrates, phytoplankton, zooplankton, and benthic primary producers. We sub-divided our zooplankton guild into several sub-groups: planktonic crustacean, chaetognaths, planktonic invertebrates (e.g., annelids, mollusks, invertebrate larvae, etc.). Other collected zooplankton (jellyfish and fish larva) were not included in the statistical analyses. Zooplankton from both deep and shallow tows of the "control" stations were pooled. The phytoplankton guild is composed of two sub-groups: nano-and microphytoplankton. All statistical and multivariate analyses were performed under permutation concepts and non-parametric tests that analyze quantified data that does not satisfy statistical assumptions underlying traditional parametric tests (Collingridge, 2013).
To compare the effect of Station (impacted vs. control), Site (Hadera vs. Ashqelon), and Guild (different guild groups), a multivariate analysis was performed on both δ 13 C and δ 15 N values. Prior to analysis, the data was arcsin-square root (on absolute fraction) transformed. A three-way PerMANOVA (factors "Station, " "Site, " and "Guild, " all fixed and crossed within each other) was performed on a Euclidian dissimilarity matrix, with 1 × 10 4 permutations. A post hoc pairwise PerMANOVA followed the analysis with a Holm correction. PerMANOVA was performed also for phyto-and zooplankton sub-groups, as well as for seasonality patterns in these guilds. In the seasonality comparisons, the factor "Guild" was replaced with "Season." TP comparison for the same factors was analyzed with a univariate three-way permutation ANOVA, followed by pairwise permutation test. All values are presented at a confidence interval of 95%.
Trophic position was calculated for consumers as the TP of producers is considered as "1." Therefore, when comparing TP, only consumers were compared. We also present a comparison of the TP ranges of zooplankton subgroups (chaetognaths, planktonic crustacean and planktonic invertebrates). All statistical and multivariate analyses were performed with R i386 3.3.3 (R Core Team, 2014) using vegan, ggplot2, lmPerm, rcompanion and boxplotdou packages; isospaces (δ 13 C vs. δ 15 N) were plotted with Matlab.

General Overview
A total of 104 samples, including primary producers, invertebrates and fish, were analyzed in this study (Table 1 and Supplementary Table S1). Generally, there was no significant differences in δ 13 C and δ 15 N between sites (Ashqelon and Hadera) or between the impacted and control stations (Table 2 and Figure 2). The only significant difference was found among guilds (three-way PerMANOVA, p < 1 × 10 −4 , Table 2), where post hoc pairwise of all combinations of guilds had p-adjusted < 0.05, excluding the combinations algae-fish, benthic invertebrate-fish and algae-benthic invertebrates (padjusted > 0.11, Supplementary Table S2). Similarly, when comparing consumer guilds TP values, a significant difference was found only among guilds (three-way permutation ANOVA, p < 0.005, Supplementary Table S3), and all post hoc pairwise consumer guild combinations showed p-adjusted < 0.01 (Supplementary Table S4). Likewise, there was no significant difference in TP between sites and stations ( Figure 6C).

Benthic Components
Results from benthic algae and benthic invertebrates showed the highest diversity, with large differences in δ 13 C and δ 15 N values occurring between the two study sites and between the different species (Figure 3). However, there were no statistical differences in the δ 13 C and δ 15 N values of benthic invertebrates between the impacted and control stations (Figures 6A,B, and Supplementary Table S1). Considering all benthos species (algae and benthic invertebrates), Ashqelon δ 15 N values ranged from 4.2 to 9.1 , while Hadera δ 15 N values ranged from -1.1 to 6.7 ( Table 1, Supplementary Table S1, and Figure 3). A maximum 15 N range of 3.6 was measured in the Bryozoan samples from the genus Watersipora sp. between Ashqelon (5.8 ) and Hadera (2.2-2.5 ) allowing a better discrimination between the two sites (Supplementary Table S1). δ 15 N enrichment at the southern site was also detected in the filter and suspension feeders, such as the ascidian Phallusia nigra and barnacles, and the stony coral O. patagonica (Table 1 and Figure 3). Comparison of δ 15 N, δ 13 C values between Hadera and Ashqelon, and "impacted" and "control" among the examined guilds, trophic position (TP), and C:N ratio (mean ± SD). TP is calculated for consumers only. Values without standard deviation represent one sample. For pairwise multivariate statistics and extra details, see Supplementary Table S1.
Likewise, the δ 13 C data of the benthic invertebrate group also displayed the highest variability, with a maximum 13 C range of 22.6 (between -2.5 and -25.2 ) between the different species (Figures 3, 4, and Supplementary Table S1). However, the high δ 13 C values (>-10.0 ) that were measured in some benthos species, like: Watersipora sp. and O. patagonica indicate that those samples are likely biased, due to the presence of carbonates (see methods). However, similar δ 13 C values, ranging between -8.5 and -6 , were measured in some benthic invertebrate species in the Lebanese coast (Fanelli et al., 2015). Those samples are characterized by marine C/N ratio (3.8-4), thus, not likely biased due to the carbonate presence.

Planktonic Components and Seasonality
Unlike the benthic invertebrate data, the phytoplankton and zooplankton groups didn't exhibit marked spatial variability in δ 13 C and δ 15 N (Figures 5A,C). The δ 13 C and δ 15 N of the two fractions of phytoplankton varied by 9.8 and 4.5 , yielding mean values of -20.7 ± 2.9 and 3.7 ± 1.1 , respectively. The difference between nano-and phytoplankton was significant (three-way PerMANOVA, pseudo F = 5.29, p < 0.02, Supplementary Table S5), however, like in the general analysis, there was no difference between the sites and control stations. Marked seasonality in δ 13 C and δ 15 N observed in phytoplankton from the Hadera site (Figures 5B,D, three-way PerMANOVA, pseudo F = 8.48, p < 0.001, Supplementary Table S6), which is located nearby a river outlet.
crustacean and chaetognaths (pairwise PerMANOVA, pseudo F = 17.09, p-adjusted = 0.009). However, no statistical difference was observed between the impacted and the control samples at both sites or between sites ( Figure 6D and Supplementary  Table S1). Seasonality in zooplankton communities was not significant (p > 0.1, Supplementary Table S8).

Fish Component and Trophic Enrichment
Trophic positions (TP, see details in Supplementary Table S1) of predatory fish species (TP = 3.98 ± 0.45, range between 3.17 and 4.49) were one level higher than herbivores (Siganus sp., TP range between 2.03 and 3.94) (Supplementary Table S1).
Overall δ 13 C values ranged from -19.2 to -16.1 for predatory fish species, indicating feeding on zooplankton communities, while the herbivores δ 13 C ranged between -16.0 and -9.1 and was consistent with the range of macroalgae (-15.5 to -7.6 ). Sparus aurata and Diplodus sp. samples yielded the highest δ 15 N values of 12.57 and 12.16 , respectively. The average 15 N of phytoplankton-zooplankton and zooplankton-fish yielded enrichment factors was 2.5 and 3.9 , respectively.

DISCUSSION
The impact of human-related processes on coastal environments is characterized by multiple sources of pollution that superimpose their effects on coastal food webs (e.g., Mancinelli and Vizzini, 2015). It is challenging to link changes in food web structure, or functioning, to a specific pollutant (both natural and anthropogenic), as it intertwines pressures, acting both directly and indirectly, on different components of the ecosystem (Raffaelli, 2005;Graniero et al., 2016). The basic approach of using stable isotopes to discriminate natural and anthropogenic factors, and their effect on food webs, is based on species-specific isotopic signatures of an impacted site that are distinct from those collected at a reference site (e.g., Vizzini and Mazzola, 2006;Orlandi et al., 2014). Several studies using this approach have shown that wastewater nutrients and organic matter (animal waste and sewage disposal) derived from onshore and offshore human activities are typically characterized by enriched δ 15 N signatures compared to those naturally occurring in marine environments (e.g., McClelland et al., 1997;Mazzola, 2004, 2006;Orlandi et al., 2014;Bongiorni et al., 2015;Graniero et al., 2016;Rossi et al., 2018). The potential impact of exposures to desalination brines on the community structure, were studied through experiments in the field and laboratory (Roberts et al., FIGURE 2 | δ 15 N and δ 13 C isospace of all samples. Circles -Ashqelon, triangles -Hadera, filled marker -sample taken from "impacted" brine, empty markersample taken from ambient "control" area. (A) All data; (B) Data by functional ecological guilds (black -benthic invertebrates, orange -benthic algae, greenphytoplankton, red -zooplankton, blue -fish).
FIGURE 4 | δ 15 N and δ 13 C isotopic ranges of all samples, presented by ecological guilds. The 2-axes boxplot shows on each axis the interquartile range (IQR) between the 1st and the 3rd quartiles, whiskers (minimum and maximum excluding outliers) which cross the IQR at the other axis' median, and outliers in circles. 2010). Studies from the Israeli coast have demonstrated relatively small-scale alterations in microbial community composition and production in elevated salinities locations of brine discharge compared to control sites (Belkin et al., 2015;Frank et al., 2017). Some of these differences were related to brine discharge byproducts, such as phosphonates and Fe-hydroxides, reducing microbial production and biomass (Belkin et al., 2015(Belkin et al., , 2017. Such effects of brine discharge on microbial production can potentiality affect coastal food webs. This study attempted to trace the impact of brine discharge on coastal food webs via SI analysis.

The East Mediterranean Food Web Structure and the Direct Effect of the Desalination Facilities
The isotopic trophic structure described here, both trophic level and carbon sources, were found similar to the only study conducted in the region, in Lebanon (Fanelli et al., 2015). Although Fanelli et al. (2015) studied mostly native and invasive fish species, their results were compared to a food web model of the Israeli EEZ by Corrales et al. (2017), and a good fit was found (R = 0.9). The enrichment between guilds was found to match other food webs worldwide (Post, 2002;McCutchan et al., 2003). Our average 15 N enrichment factors for the three main guilds: phytoplankton-zooplankton and zooplankton-fish of 2.5 and 3.9 , respectively, correspond well with the mean factor of 3.4 suggested by Post (2002) and 2.3 considered by McCutchan et al. (2003). Corrales et al. (2017) found in the Israeli EEZ, via a food web model, a similar enrichment structure, where the enrichment between zooplankton and fish groups is ∼1 trophic level, similar to our observation (3.9 , equivalent to 1.11 trophic levels based on an enrichment of 3.4 per trophic level). Some trophic levels calculated by Corrales et al. (2017) were comparable to our results (e.g., cephalopods 3.50/3.30, jellyfish 3.20/2.90, shrimps 2.80/2.90 in this study /Corrales et al. (2017), respectively), and some were not (rocky fish 4.70/3.02, benthic invertebrates 3.10/2.15 in this study /Corrales et al. (2017), respectively). They, however, modeled the entire Israeli continental shelf while our measurements were local.
We found no significant difference in the food web structure between the impacted and control stations. This observation does not mean that the brine itself is not affecting the fauna and flora as was found in some studies (Chesher, 1971;Gacia et al., 2007;Frank et al., 2017), however, if such effect exists, it is not reflected through the isotopic composition of the organisms. The lack of a significant differences between impacted and control stations is supported by several studies conducted along the Israeli coast on planktonic microbial communities (Belkin et al., 2017(Belkin et al., , 2018, in Saudi Arabia on photosynthesis efficiency of coral symbionts , and of the Al-Jubail desalination plant, in Saudi Arabia as well, on phyto-and zooplankton (Abdul-Azis et al., 2003). When distinguishing between the pelagic vs. the benthic food webs, however, evidence of brine effects in the form of N pollution are detectable in the benthic system (see next subsection). We found no evidence of brine effects on the zooplankton guild subgroups. On average, the δ 15 N of both the phytoplankton (3.6 ± 1.0 ) and zooplankton (6.4 ± 1.0 ) groups were equal to, or higher than, the δ 15 N range between 1.0 and 4.0 in zooplankton communities (caught with a mesh grid of 100 and 333 µm) of the surface open waters of the southeastern Mediterranean (site M15, south east of Cyprus, Koppelmann et al., 2009). Our results were also significantly higher than δ 15 N values of mesozooplankton communities in the northeast Aegean Sea and surrounding the coast of Cyprus (1-4 , Hannides et al., 2015), but agree well with values of zooplankton (grazers) off the Lebanese coast (5.8 ± 0.9 , Fanelli et al., 2015).

Differences Between Facilities and Benthic Anthropogenic Nitrogen Pollution
Enriched δ 15 N signatures were observed in our study. This was especially true for the benthic invertebrate group of the southern site, where ground water amelioration, rich in nutrients, is released into the coastal water as part of the discharge of brine (Israel Electric Company [IEC], 2017). The difference in δ 15 N values of Bryozoan tissue between the Ashqelon and Hadera sites was 3.3 , equivalent to an increase in trophic position. This difference suggests different N sources between Ashqelon and Hadera. The Hadera site exhibited values which are equal to, or lower than, natural marine N sources, such as POM (0-2 ; Koppelmann et al., 2009). The values at Ashqelon, however, more likely represent a higher contribution of anthropogenic sources of N. Mean measured NO − 3 and NH + 4 levels were found to be significantly higher in both impacted stations, but particularly in the southern site, compared to the values routinely measured in the pelagic surface water of the eastern Mediterranean (NO − 3 = 0.2-0.3 µM, Kress et al., 2014;Raveh et al., 2015;Rahav et al., 2018). In Ashqelon NO − 3 surface values ranged between 0.2 and 12.2 µM, while in Hadera the values ranged between 0.1 and 3.0 µM (Israel Electric Company [IEC], 2012). Furthermore, in Ashqelon, surface NO − 3 and salinity have been shown to be highly and positively correlated, annually, indicating that the brine outflow is the source for high N loads at this site. These differences in N loads between the two sites reflect the combination of desalinization brine discharge with ground water amelioration rich in NO − 3 and NH + 4 , at the Ashqelon site, is likely the cause for the differences in δ 15 N of the benthic invertebrate groups between the two sites, as the amelioration products sink along with the heavy brine.
The most notable results were found in the analysis of benthic groups. In addition to the differences in bryozoan (Watersipora sp.) tissue there was also a difference of one trophic level between the two sites in the tunicates (P. nigra), as well as a generally higher δ 15 N in the benthos group of the southern site (Figure 3, Table 1, and Supplementary Table S1). Benthic species have been shown to be the most suitable species for distinguishing between unpolluted and disturbed locations (e.g., Vizzini and Mazzola, 2006;Orlandi et al., 2014;Vizzini et al., 2017). Utilization of photosynthetic carbon from either algae and/or phytoplankton, or from particulate and Frontiers in Marine Science | www.frontiersin.org sedimentary organic matter, as well as the incorporation of HCO − 3 could explain the wide range of δ 13 C values of the benthos species in this study and the lack of site-specific δ 13 C differences unlike the δ 15 N data. However, the high δ 13 C values of bryozoan and coral may likely result from their carbonate skeletons, since those samples were not pretreated with HCl. This may also be the case for the high δ 13 C values measured in the Littorina sp. and worm (unidentified species). Nevertheless, our results support spatial and species-specific variations in dietary sources of nitrogen and carbon in the benthos group.
Fish are situated high in the food web hierarchy, and therefore carnivores do not generally mirror the "isotopic state" of a single location, but show a "spatially integrated" isotopic signal (e.g., Vizzini and Mazzola, 2006). Unlike the carnivores, herbivorous fish feed directly on primary producers such as macroalgae and seaweeds, and hence their δ 15 N ranges reflects the ranges of the local flora. The Siganus sp. is considered as an obligate herbivore fish, resulting in a wide δ 15 N range between 4.2 and 10.7 , which is consistent with the range of macroalgae in the study area (between 2.2 and 9.0 ). However, studies on Siganus sp. feces show that it also feeds on fauna such as mollusks, annelids and more (Supplementary Table S1 in Guy-Haim et al., 2017), which also supports the wide range of observed δ 15 N and high TP measured for the Hadera specimen. A similar δ 15 N range between 5.9 and 9.7 was measured in Siganus sp. off the coast of Lebanon (Fanelli et al., 2015). In contrast to the Siganus sp., both S. aurata and Diplodus sp. that feed on benthic bivalves, gastropods and crustaceans (Tancioni et al., 2003) provided the highest δ 15 N values of 12.57 and 12.16 , respectively. For comparison, Fanelli et al. (2015) found a value of 9.2 ± 0.4 in δ 15 N for Diplodus sargus off the coast of Lebanon, ∼3 lower than our results. This δ 15 N enrichment may reflect the penetration of the anthropogenically N pollutants up the food webs of the Israeli coast. For the S. aurata specimen, such high values could also indicate that the individual escaped from an offshore fish farm (Serrano et al., 2007(Serrano et al., , 2008, located approximately 10 km north of the sampling site. However, in the case of escaped S. aurata from Mediterranean farms the enriched δ 15 N values coincided with depleted δ 13 C values (Arechavala-Lopez et al., 2013). This was not the case in our samples that yielded δ 13 C values comparable with wild fish, further supporting δ 15 N enrichment from other anthropogenic sources. Since only a limited number of fish were analyzed, almost entirely from the southern "control" site, the extent of our conclusions are restricted. Nevertheless, these results point to the potential use of the δ 15 N ratio of benthic carnivores to trace anthropogenic N pollution, as seen in Ashqelon.

Brine Density Plume as Pollution Vector
The dispersal range and trajectories of discharged material may vary spatially, due to both physical and biological differences. Site-specific dispersal characteristics depend mostly on the type and size of particles, advection processes and resuspension events after deposition (Dawoud, 2012). Therefore, it is difficult to estimate a general rule of brine/pollution dispersal. Modelbased estimates of fish-farm-derived particle transport as a function of resuspension rate indicated trajectories of several hundred meters' from the outlet (Cromey et al., 2002). A longerdistance trajectory was predicted for particles leaving fish farms along the Israeli coast based on a Lagrangian tracking model (Grossowicz et al., 2017). It is worth noting that in Ashqelon a wide dispersal range of varying directions has been observed for the brine (Israel Electric Company [IEC], 2017). The heavy brine is also characterized as a density plume, flowing under the force of gravity down to the bottom (Biton et al., 2008). For example, in Ashqelon, NO − 3 levels, following amelioration and dilution with the brine, were around 150 µM at the brine outlet (Israel Electric Company [IEC], 2017). These high nutrient loads can be transported to the benthic communities via density plumes that sink and disperse near the bottom (Biton et al., 2008).
Ground water in the Israeli coastal aquifer is typically enriched in nitrogen pollutants in the form of synthetic fertilizers from agriculture (Kass et al., 2005) with a δ 15 N range between +4.7 and +11.4 (Kaplan and Magaritz, 1986). Thus, it is not surprising that the benthic invertebrate guild of Ashqelon yielded the most enriched δ 15 N, which agrees well with an anthropogenic source in the form of groundwater. Brine derived density plumes may therefore, act as a vector transporting high loads of nutrients from groundwater amelioration to benthic communities. It is also worth noting that another source of anthropogenic N loads may be the Gaza Strip sewage outlet located 17 km south of the impacted location (Nassar and Afifi, 2006;El-Nahhal et al., 2014). However, we lack information regarding the dispersal range and trajectories of this source as well as its δ 15 N signatures and therefore, can't estimate its general impact on food webs at the southern site.
This is a first of its kind study, where carbon and nitrogen stable isotope ratios of primary producers, invertebrates and fish, covering most trophic levels, collected from the highly urbanized Israel coastal ecosystem were analyzed in an attempt to explore the impact of brine discharge on coastal food webs. Samples of all examined guilds from the brine plume of two desalination plants outlets, were analyzed and compared with similar samples from control stations. Statistically, we found no differences in the δ 13 C and δ 15 N signatures between the impacted and control stations of both sites. Thus, the immediate effect of the brine on coastal food webs was undetectable using the SI method. However, a significant difference in δ 15 N was observed in the benthic invertebrates between the two sites, with a general enrichment in δ 15 N in the southern site, where ground water amelioration products are discarded alongside to the brine. Brine derived density plumes rich in nutrients from groundwater amelioration is suggested here as a possible vector, delivering anthropogenic nitrogen to benthic communities. High sensitivity of benthic components to anthropogenic derived N pollution via density plumes allowed insight into site-specific processes. Benthic carnivores, such as S. aurata and Diplodus sp. which feed on benthic components are therefore the most likely predators to be affected by this N pollution. These results from the warm, and ultra-oligotrophic, N-limited coastal waters of the southeast Mediterranean presents important insights into the possible future effects of the growing industry of desalination on coastal food webs.