Introduction

Wetlands constitute major biodiversity hotspots and provide multiple ecosystem services (Davidson et al. 2018; Finlayson et al. 2018). However, since the beginning of the twentieth century, wetland loss and degradation has intensified due to increasing anthropogenic pressures (Davidson 2014). Eutrophication is a major cause of wetland degradation worldwide (Ramsar Convention on Wetlands 2018), mainly driven by agricultural, urban and industrial activities (Kingsford et al. 2016; Verhoeven et al. 2006), and causes multiple negative effects such as alteration of the biogeochemical cycling, toxic cyanobacterial blooms, anoxia and collapse of fish and submerged plant communities (Ansari et al. 2011; Sánchez-Carrillo et al. 2011). Eutrophication may also lead to negative socio-economic impacts, due to declining quality and increased disease risk of drinking water supplies, decline of fish stocks, or decreasing tourism and leisure activities (e.g. due to bad odors or loss of water transparency) (Schuyt 2005). Eutrophication impacts are further enhanced due to synergy with climate change effects such as warming, soil erosion and droughts (Cramer et al. 2018; Geijzendorffer et al. 2019; Green et al. 2017). Thus, maintaining good water quality status of wetlands is essential for human well-being, biodiversity maintenance, and to enhance resilience to climate change (Green et al. 2017).

Implementation of broad-scale measures (e.g. European Union Directives) has led to major improvements of water quality through reductions of anthropogenic nutrient inputs in much of northern and central Europe (EEA 2015). However, less information is available about southern Europe. Countries of the Mediterranean region have greater challenges due to lower precipitation and rapid increases in human population, particularly around coastal wetlands (MWO2 2018). In the most important remaining Mediterranean wetlands, it is vital to increase knowledge on spatio-temporal nutrient inputs at the catchment scale, and plan effective actions for preventing and controlling eutrophication.

Although wetlands provide a vital ecosystem service by reducing nutrient loads and enhancing water quality (Johnston 1991; Fisher and Acreman 2004), high water quality should be ensured at the catchment level before surface water reaches wetlands, so as to guarantee their functioning and biodiversity value. This action may be accomplished by conserving natural vegetation in the catchment area, tertiary water treatment, or use of constructed wetlands. Green infrastructure, such as constructed wetlands and riparian buffer zones, has been used to reduce nutrient loadings into protected natural wetlands, such as the iconic Everglades (Chimney and Goforth 2006; Tonderski et al. 2017).

Doñana in SW Spain is an iconic wetland complex (Scheffer et al. 2015) recognized as one of the most biodiverse areas in Europe and the Mediterranean region (Green et al. 2018). However, the core area protected as a National Park and a UNESCO World Heritage Site (WHS), which contains an extensive marshland, is threatened by impacts on water quantity and quality, especially due to intensive agriculture and municipal wastewaters in the catchment (Camacho-Muñoz et al. 2010; Green et al. 2017; Paredes et al. 2019). Nevertheless, there is a lack of detailed studies of nutrient status and eutrophication both inside and outside the boundaries of Doñana National Park (Serrano et al. 2006).

In this study, our objective was to quantify water quality in the Doñana marsh and its entry streams by determining nutrient concentrations and phytoplankton abundance in the water column, and to compare our findings with accepted water quality standards to assess whether poor water quality is a threat to biodiversity in the area. We studied the spatio-temporal variability of standard eutrophication indicators (nitrogen (N), phosphorus (P) and chlorophyll-a (chla) concentrations) across Doñana (both the marshland and the streams that supply it) during four consecutive years (2013 to 2016). We evaluated long-term trends in intensive agriculture as a likely cause of eutrophication by quantifying the changes in surface area and distribution of greenhouses in the catchment between 1995 and 2018, since these contain berry farms relying on extraction of groundwater and intensive chemical use (Green et al. 2017). We hypothesized that: (1) streams would show higher nutrient and chla concentrations than the marsh because they are located closer to pollution sources (agriculture and urban wastewaters); (2) streams located in sub-catchments most directly affected by greenhouse expansion and urban effluents would show higher nutrient and chla concentrations; (3) variation in nutrient and chla concentrations within the marsh would be explained by distance to the stream entry points, as well as by depth variation and evaporation processes within the marsh; (4) total nutrient and chla concentrations in the marsh and the streams would be lower in the cool, wet months than the warm, dry months. Our study provides insight into the eutrophication processes in Doñana, which may be common to other complex Mediterranean aquatic ecosystems affected by increasing anthropogenic pressures.

Materials and methods

Study area

Doñana marsh

The Doñana wetland complex occupies a seasonal brackish floodplain located in the estuary of the Guadalquivir River on the Atlantic coast in Southwestern Spain (37°0′N 6°37′W) (Fig. 1). The natural marsh is of international importance for biodiversity conservation, covering a large area (270 km2) of impermeable clay in Doñana National Park (543 km2), which was declared in 1969 and later designated as a Wetland of International Importance under the Ramsar Convention in 1982, EU-Special Protection Area for birds in 1988, and WHS in 1994 (Green et al. 2018). The entire wetland complex is included in a UNESCO Biosphere Reserve declared in 1980 (Fig. 1). In the past, the natural marsh covered a much larger area, but about 80% was transformed into other land uses (mainly agriculture) during the twentieth century (Zorrilla-Miras et al. 2014a).

Fig. 1
figure 1

Study area, indicating the boundaries of protected areas and sub-catchments, and the location in Spain as an insert. Orange dots show the sites used for surface water sampling. All geographical coordinates refer to the study area. Points A and B within the marsh were sampled more often than most, and are referred to in the results section

Full size image

Doñana is in a region with a subhumid Mediterranean climate with an Atlantic influence. Mean annual temperature is 17.3 °C, ranging from 10.1 °C in January to 25.5 °C in August (mean values calculated for the period 2000–2015. Source: El Rocío Metheorological Station). The mean annual precipitation is 550 mm, with high interannual variation ranging between 170 and 1000 mm (see Fig. S1 and Díaz-Delgado et al. 2016). Precipitation is mainly concentrated between October and April with a hot, dry season from May to August. Hydrological years begin in September, when reflooding may begin, and end in August by which time the marsh is always dry. Mean annual precipitation for the four hydrological years (from September 2012 to August 2016) relevant to our study was 480.6 ± 90.3 mm, with 85.5% of precipitation falling between October and April. Direct precipitation and surface flow from streams constitute the main water inputs into the marsh, which have been estimated to contribute 99–144 hm3/year (www.chguadalquivir.es) and 15–25 hm3/year (https://sig.mapama.gob.es/redes-seguimiento/index.html?herramienta=Aforos), during the period Oct.2013 to Aug.2016, respectively (surface flow corresponds to the supply of Rocina and Partido streams. No data were available for Guadiamar). Additional inputs to the marsh come from direct groundwater inputs and tidal inputs from the Guadalquivir estuary. There is high interannual variation in the extent and depth of flooding in the marsh (see Fig. S2 and Díaz-Delgado et al. 2016). Owing to flat, shallow topography, hydrodynamics and water distribution in the marsh are strongly affected by wind, which may produce changes in depth due to accumulation of water downwind (Ramos-Fuertes et al. 2014). In addition, water slowly moves by gravity in a southerly direction until flowing through sluices into the estuary of the Guadalquivir River. Evaporation meanwhile leads to a spatial salinity gradient with higher salinity towards the south, away from entry streams (Fig. S5; see also Espinar and Serrano 2009; Grillas et al. 1993).

Hydrological connectivity in the marsh depends on the flooding level and microtopography. In general, below a flooding level of about 1.3 masl, the marshes are highly fragmented with surface waters divided into different units (Fig. S2). The highest flooding levels are reached between November and March, and the marsh dries out completely in July–August (Díaz-Delgado et al. 2016). The mean annual maximum water column at the deepest point in the marsh is 0.51 m. Evaporation in the marsh is fairly constant between October and February (2 mm day−1), then increases about 5% per day from March until reaching 8 to 10 mm day−1 in June/July (Ramos-Fuertes, 2012). The marsh is dominated by Bolboschoenus maritimus and other helophytic vegetation (Lumbierres et al. 2017; Paredes et al. 2019).

Entry streams

The most important entry streams draining into the Doñana marsh are located within four sub-catchments (Guadiamar (1300 km2), Rocina (416 km2), Partido (291 km2), Sotos (60 km2)) (Fig. 1). Apart from direct precipitation and runoff, all these entry streams are also partly fed by groundwater discharge, which has been estimated at 34 hm3/year for the Rocina and 11 hm3/year for the Partido, although flows have since been reduced due to groundwater extraction (Arambarri et al. 1996; Guardiola-Albert et al. 2011). Flow rates for February to May 2016 are provided in Table S3. All four streams are now intermittent throughout the year. Streams of the Rocina, Partido and Sotos catchments enter the marsh in its north-west corner (Fig. 1). The Rocina and Sotos catchments are impacted by similar agricultural activities without any influence of wastewater treatment plants, and owing to the small number of sampling points for Sotos, we combined these adjacent catchments (i.e. Rocina/Sotos) to facilitate statistical analyses. Since a restoration project completed in October 2014, water from the Guadiamar catchment flows again into the marsh through a channel at the north east corner (Fig. 1), after previously being disconnected by dykes for decades due to drainage of adjacent areas for agriculture. Since it was affected by a spill of mine waste in 1998, the Guadiamar catchment has been partly restored, with an increase in riparian vegetation (Ontiveros et al. 2013).

Field sampling

We collected surface water samples over a network of 59 points across the Doñana marsh and entry streams (Fig. 1) between January 2013 and June 2016, covering part or all of the four consecutive hydrological years from September 2012 to August 2016. From here on, these hydrological years are referred to as 2013, 2014, 2015 and 2016 respectively (e.g. September-December 2012 is part of “hydrological year” 2013). The samples we analyzed were mainly collected between October and early June since surface water is absent in most streams and the marsh from late June to mid September due to scarce precipitation and high temperatures, and because groundwater extraction (mainly for greenhouses) has severely reduced discharges into the streams (Guardiola-Albert and Jackson, 2011). Sampling time, location and number of water samples collected varied between months and years due to variation in the spatial distribution of surface water, and only a subset of the network of points was sampled in any given month. During 2014, flooded areas were particularly scarce across the marsh due to low levels of precipitation (Fig. S2).

We collected 1-2L of surface water to measure N, P and chla concentrations. At the end of each field day, we filtered part of our samples through FILTER-LAB MFV5047 glass-fiber filters (0.45 µm pore size) using a low-pressure vacuum pump. We stored all filtered samples in the freezer (-20 °C) prior to analysis of dissolved nutrient concentrations. To determine hydrogen isotopic ratios (δ2H) in water, we filtered 2 ml of water through CHROMAFIL® Xtra PET-20/25 0.20 µm and transferred it into a glass vial with a septum cap. We stored these samples at 5 °C prior to analysis. At the spot where the water samples were taken, we recorded in situ electrical conductivity using a WTW (Weilheim, Germany) Multi-340i handheld meter and water column depth using a measuring stick. Some parameters (e.g. chla or δ2H) were only measured in later years, thus leading to differences in sample sizes.

Laboratory analyses

Nitrogen and phosphorus

We used standard colorimetric methods to measure the concentration of three dissolved inorganic nitrogen (DIN) species (nitrate NO3−, nitrite NO2− and ammonium NH4+) and phosphate PO4−3 (ISO 13395:1996 for nitrate and nitrite; ISO 11732:2005 for ammonium; ISO 15681–2:2003 for phosphate) on a multi-channel SEAL Analytical AA3 AutoAnalyzer (Norderstedt, Germany). We used unfiltered water samples to analyze total nitrogen (TN) by digestion with potassium persulfate as described by Nydahl (1978), and total phosphorus (TP) by digestion with a mixture of potassium persulphate and sulphuric acid following the colorimetric method (Murphy and Riley 1962).We performed all the nutrient analyses at the Aquatic Ecology Laboratory (LEA), Estación Biológica de Doñana, (EBD-CSIC, Seville, Spain).

Chlorophyll-a

We determined chla concentrations using acetone extraction (UNESCO, 1966). We first filtered water samples through FILTER-LAB MFV5047 glass-fiber filters (0.45 µm pore size) and then soaked the filters in acetone (90%) overnight at 4 °C in the dark. The acetone samples were then filtered again. Finally, we read absorbance values at 750, 663, 645 and 630 nm using a UV VIS spectrophotometer (AnalytiK Jena AG, model SPECORD ® 205), from which we calculated chla concentrations (the 750 nm reading was used to correct for turbidity).

Water isotopic analyses

To obtain hydrogen isotopic values (2H/1H, δ2H, ‰), we analyzed the water samples by cavity ring-down spectroscopy (CRDS) using a L2130-i spectrometer (Picarro Inc., 480 Oakmean Parkway, Sunnyvale, California, 94085, USA; www.picarro.com) in the Stable Isotopes Laboratory (LIE-EBD) at the EBD-CSIC.

Marsh flooding masks

Flooding masks were generated by the Remote Sensing and GIS Laboratory (LAST-EBD) using mid-infrared band 5 (1.55–1.75 µm, TM and ETM +) and band 4 (0.8–1.1 µm, MSS) to produce final inundation masks based on 30 × 30 m pixels from Landsat images (see details in Bustamante et al. 2009; Díaz-Delgado et al. 2016). These masks were used to illustrate the variation in flooding patterns during our sampling campaigns.

Spatio-temporal evolution of greenhouses

In recent decades, greenhouses dedicated to strawberry, raspberry, blueberry and blackberry cultivation have expanded in the surroundings of Doñana, and now occupy more than half of the total cultivated area of the Rocina, Partido and Sotos catchments, based on Aeolian sands well suited to these crops (WWF 2018). These crops depend on intensive groundwater extraction and agrochemical use, and so contaminate the entry streams with nutrients originating from fertilizers (Manzano et al. 2009). Hence we quantified their expansion within the Doñana catchment between 1995 and 2018 via remote sensing, applying an ensemble of target detection methods on radiometrically normalized Landsat images (Díaz-Delgado et al. 2016) to determine the annual greenhouse surface area in the Rocina, Partido and Sotos catchments (see details in “Methods S1” section). This work was carried out at the Remote Sensing and GIS Laboratory LAST-EBD.

Data analyses

Based on the natural limits of the stream catchments and the marsh, we defined four waterbody categories: (1) Partido, (2) Rocina/Sotos, (3) Guadiamar and (4) marsh. To visualize the spatio-temporal patterns of nutrients and chla across these different waterbodies, we mapped the concentrations of DIN, TN, PO4−3, TP and chla from samples collected during December 2014 and 2015, February 2015 and 2016, and May of 2015 and 2016, and laid them over the flooding masks available for the closest dates. These were the months with the most extensive spatial coverage of sampling points.

Principal component analysis (PCA)

We used PCA to determine whether there was a clear spatial pattern among our four waterbodies. Our dataset included 338 samples collected from 56 different sites during 19 different sampling periods between January 2013 and June 2016. Each observation contained information on nine different variables: TP, TN, PO4−3, NH4+, NO3−, NO2−, conductivity, depth and δ2H.

We first standardized the variables to make the magnitude of the different variables comparable to each other (Gotelli and Ellison 2004). Then we constructed a correlation matrix and extracted the eigenvalues and eigenvectors which inform about the significance of PCs and the importance of each original variable to a particular component. We only considered as significant those PCs with eigenvalues of 1 or higher, since they retain the highest proportion of the variance (Le et al. 2017). We used the ‘factoextra’ (Kassambara and Mundt, 2017) and ‘ggplot2’ (Wickham 2016) packages from CRAN-R software to perform the PCA and to plot the PCA output and the scree plots, respectively.

Two-way ANOVA and post-hoc pairwise comparisons

To test for significant differences over time and space between our four waterbodies, we used two-way ANOVA with TP, TN, PO4−3, NH4+, NO3−, NO2− and chla as the dependent variables. We used visual (‘Normal Q-Q’ and ‘Residuals vs. fitted values’ plots) as well as statistical tools (global validation test of linear model assumptions: ‘gvlma’ package from CRAN-R) to validate model assumptions. We applied log-transformation to all variables to meet assumptions of normality and homoscedasticity. ‘Waterbody’ was included as a categorical explanatory variable with four levels, and ‘sampling period’ as one of 16 levels. When there were significant effects of waterbody, we determined pairwise differences with Tukey’s “honestly significant difference” (HSD) post–hoc tests. Data used in these analyses were collected from 59 sampling sites between Jan 2013 and May 2016 (for nutrients) and from 55 sampling sites between Dec 2014 and May 2016 (for chla) (Fig. S3). Data from June, July and August were excluded because lack of surface water severely reduced the number of available sampling sites. We performed the analyses using the ‘stats’ package in CRAN-R software.

Analyses of variation in water quality within the marsh

To analyze variation in nutrient and chla concentrations within the marsh itself, we used multiple linear regression models and hydrological predictor variables (conductivity, depth and δ2H). The marsh received relatively low conductivity water (Table S1) from multiple stream entry points, and conductivity served as a proxy for the distance of a given point within the marsh to such entries. Understanding of how water moves and mixes within the marsh is limited, and the contribution of different sources to the water present at each sampling point is unclear. The dataset for these analyses included nutrient data from 26 sampling sites (Fig. S4) within four different periods (March–April 2013, December 2014, February 2015 and May 2015) when there was more extensive flooding across the marsh and we were able to sample a greater number of localities. In the case of chla, we included data from 24 sampling sites (Fig. S4) during February 2015 and May 2015.

To run the models and obtain the regression parameters and partitioning of variance we used the functions ‘lm’ and ‘aov’, respectively, from the ‘stats’ package in CRAN-R. To validate model assumptions we used visual (‘Normal Q-Q’ and ‘Residuals vs. fitted values’ plots) as well as statistical tools (global validation test of linear model assumptions: ‘gvlma’ package from CRAN-R). We applied log or square root transformations to the variables when necessary to meet assumptions of normality and homoscedasticity. When transformations failed to meet assumptions, we then ranked the dependent variable and performed the models again (Seaman et al. 1994). For each response variable and sampling period, we retained the model with the smallest Akaike Information Criteria (AIC) value (from hereon “best model”). We obtained best models using the ‘step’ function (where “direction = both”) from the ‘stats’ package in CRAN-R software 3.4.1. (https://www.R-project.org).

Comparing our data to water quality reference values

European countries, Spain among them, carry out the monitoring and evaluation of surface water quality status in compliance with the Water Framework Directive (WFD 2000/60/EC). In Spain, reference values for those WFD indicators used in each surface waterbody type (i.e. rivers, lakes, coastal water and transitional water) are described in the Royal Decree 817/2015. However, for the group “lakes” in which waterbodies within the Guadiamar, Marsh and Sotos catchments are classified within the Royal Decree 817/2015, information about nutrient reference values is lacking (Table S4). Moreover, the “rivers” group in which Rocina, Partido and one stream of the Sotos catchment are classified, only provides limits for the very good and good status based on NH4+, PO43− and NO3−, but information is lacking about other nutrient indicators such as NO2−, TP and TN, and no limits are established for worse water quality status. Thus, to provide a more complete assessment of the current surface water quality status in both the Doñana marsh and entry streams, we used nutrient reference values suggested by the OECD (2007) for five different use classes. We considered the OECD index as a good proxy to assess waterbodies included within this study because it is based on general European legislation for surface waterbodies in which Spanish surface waters are included. Class I represents an undisturbed, natural aquatic system, and may be considered equivalent to the “higher status” classification of the WFD. Classes I and II are the only ones supporting complete biodiversity functioning, whereas no fish are expected to survive in classes IV and V. We calculated the percentage of samples classified within each class for each nutrient parameter and waterbody. We presented percentages with and without data from June to September), given that the low water level and high evaporation rate during these months were expected to increase nutrient concentrations.

For chla concentrations, we used reference values based on trophic classification (OECD 1982), such that values < 2.5 µg L−1 can be considered oligotrophic, 2.5–8 µg L−1 mesotrophic, 8–25 µg L−1 eutrophic and > 25 µg L−1 hypereutrophic.

Results

Greenhouse expansion in the Doñana catchment

The surface area covered by greenhouses in the Partido and Rocina/Sotos catchments increased by 487% in 21 years from 939 ha in 1995 to a peak of 5510 ha in 2016, but then declined to 3204 ha by 2018 (Fig. 2). During our study period for water quality (2013 to 2016), the mean area (± SD) of greenhouses was 1506 ± 338 ha in the Partido catchment and 3144 ± 501 ha in the Rocina/Sotos catchment. This represented 5.1% and 7.2% of the total subcatchment surface area respectively. There was relatively more expansion in the lower part of the Partido catchment and the upper part of the Rocina catchment (Fig. 2).

Fig. 2
figure 2

Surface area covered by greenhouses (mainly for berries) between 1995 and 2018 in the Partido and Rocina/Sotos catchments (a), and the spatial distribution of these greenhouses in 1995, 2008 and 2018 (b). There were no greenhouses in the Doñana marsh and almost none (17 ha in 2018) in the Guadiamar catchment

Full size image

Spatial heterogeneity in water quality parameters

General differences in water quality between waterbodies

Within a PCA of the set of water quality variables, PC1 (eigenvalue PC1 = 3.84) explained 42.7% of the total variance and was highly correlated with PO4−3, TP, NO3−, TN and NO2−. PC2 (eigenvalue PC2 = 1.70) explained 18.9% of the total variance and was highly correlated with conductivity and δ2H. The resulting biplot of PC1 vs. PC2 (Fig. 3) reveals major variation between waterbodies (Partido, Rocina/Sotos, Guadiamar streams and the marsh). Samples from the Partido and some samples from the Rocina/Sotos streams had higher PC1 values, reflecting higher nutrient concentrations. Samples from the marsh generally had higher PC2 values, reflecting higher conductivity and δ2H.

Fig. 3
figure 3

PCA of surface water quality from the Doñana marsh and three main entry streams

Full size image

The marsh showed the highest conductivity values (mean ± s.d. = 4330 ± 5737 µS cm−1), followed by the Guadiamar (2114 ± 1908 µS cm−1), the Partido (mean ± s.d. = 1191 ± 766 µS cm−1) and the Rocina/Sotos (mean ± s.d. = 544.2 ± 195.8 µS cm−1). Mean water depth for each waterbody ranged from 26.8 cm in the marsh to 67.1 cm in the Guadiamar. The maximum depth recorded was in the Partido stream in December 2014 (480 cm). Mean (± s.d.) δ2H values for Partido (− 26.9 ± 5.5‰), Rocina/Sotos (− 26.5 ± 7.3‰) and Guadiamar (− 23.8 ± 13.9‰) were very similar, and lower than the marsh (− 1.5 ± 25.9‰). The marsh had the largest range of values (− 52‰ to 68.2‰) (see details in Table S1).

Comparing nutrient and chla concentrations between waterbodies

Maps showing the nutrient concentrations (dissolved inorganic nitrogen [DIN, the sum of NH4+, NO3− and NO2−], TN, PO4−3 and TP) in different waterbodies for six major sampling events (December 2014, February 2015, May 2015, December 2015, February 2016 and May 2016) confirmed that the streams had generally higher values compared to the marsh, a difference that was consistent between years and seasons (Figs. 4, 5, 6 and 7). In contrast, within a given waterbody there were no consistent differences in the concentrations recorded between cool, wet months (December or February) or May. The Partido and the lower part of the Rocina/Sotos catchments were the areas with the highest nutrient concentrations. These spatial differences were more marked for dissolved inorganic nutrients than for Total N and Total P.

Fig. 4
figure 4

Concentration values of dissolved inorganic nitrogen (DIN) in 2015 and 2016 in the Doñana marsh and entry streams, and the extent of flooding in the marsh at the time of sampling

Full size image
Fig. 5
figure 5

Concentration values of phosphate (PO4−3) in 2015 and 2016 in the Doñana marsh and entry streams, and the extent of flooding in the marsh at the time of sampling

Full size image
Fig. 6
figure 6

Concentration values of total nitrogen (TN) in 2015 and 2016 in the Doñana marsh and entry streams, and the extent of flooding in the marsh at the time of sampling

Full size image
Fig. 7
figure 7

Concentration values of total phosphorus (TP) in 2015 and 2016 in the Doñana marsh and entry streams, and the extent of flooding in the marsh at the time of sampling

Full size image

DIN concentrations during these six sampling events ranged from < 0.1 to 9 mg N L−1, and clear separation between confidence intervals showed they were much higher in streams (Geometric Mean = 1.16 mg N L−1; 95% Confidence intervals (CIs) 0.85–1.58; N = 132) than in the marsh (GM = 0.04 mg N L−1; 95% CIs 0.03–0.05; N = 92). Total N concentrations ranged from 0.4 to 12.5 mg N L−1, with a GM of 4.84 mg N L−1 (95% CIs 4.28–5.48) in streams compared to 2.20 mg N L−1 (95% CIs 1.96–2.47) in the marsh.

PO4−3 concentrations ranged from < d.l. to 2.8 mg P L−1, and were much higher in streams (GM = 0.21 mg N L−1; 95% CIs 0.16–0.29) than in the marsh (GM = 0.02 mg N L−1; 95% c0.01–0.04). Total P concentrations ranged from 0.03 to 3 mg P L−1, with a GM of 0.48 mg P L−1 (95% CIs 0.41–0.57) in streams and of 0.17 mg P L−1 (95% CIs 0.14–0.21) in the marsh. At one southern location in the marsh adjacent to a large heronry, PO4−3 and TP concentrations were particularly high in May 2015 (Fig. 5 and 7).

Differences between waterbodies were less consistent for chla, although the lower part of the Rocina and Partido catchments often had the highest concentrations (Fig. 8), and the concentrations in streams were sometimes clearly higher than in the marsh (e.g. in February 2015). Chla concentrations ranged from 0.5 to 88 µg L−1 (n = 172), and values in streams (GM = 6.78 µg L−1; 95% CIs 5.66–8.13; N = 109) were generally higher than in the marsh (GM = 4.27 µg L−1; 95% CIs 3.34–5.45; N = 63), with no overlap between confidence intervals.

Fig. 8
figure 8

Concentration values of chlorophyll-a in 2015 and 2016 in the Doñana marsh and entry streams, and the extent of flooding in the marsh at the time of sampling. Thresholds used for the legend are based on the OECD (1982) classification

Full size image

Two-way ANOVAs confirmed that there were highly significant differences between ‘waterbodies’ (P < 0.0001) for all dependent variables (nutrients and chla) (Table S2), where the Partido, Rocina/Sotos and Guadiamar showed consistently higher concentrations than the marsh as shown by model estimates (Fig. 9). There were also highly significant differences between sampling events, although these differences were not consistent between seasons (Table S2). Concentrations in the marsh were significantly lower than in all three streams for all nutrients except NH4+ and PO4−3. The Partido had significantly higher concentrations than other streams for all nutrients except NO3−. There were no significant differences in nutrient concentrations between Rocina/Sotos and Guadiamar. Rocina/Sotos and Guadiamar both had significantly higher chla concentrations than the faster flowing Partido and the marsh (Fig. 9).

Fig. 9
figure 9

Comparison of the model estimates from two-way Analysis of Variance (ANOVA) showing the relative differences between the three streams (Partido, Rocina/Sotos and Guadiamar) and the marsh (M; model estimates are zero for this waterbody in all cases). Higher model estimates indicate higher average concentration values, and all response variables were log-transformed prior to analysis. See Table S2 for further details

Full size image

Explaining spatial variation of nutrient and chla concentrations within the marsh

We used multiple regression to explore the relationship between nutrient concentrations and conductivity, depth and isotopic variation δ2H within the marsh. For each nutrient response variable, best models included at least one physicochemical variable for at least two of the sampling periods, and the proportion of variation explained by the models was relatively high for TN, TP and NO2− (Table 1). In the final models, conductivity and depth repeatedly showed significant negative partial effects on nutrient concentrations. δ2H had a positive effect on nutrients when this was the only explanatory variable in a final model, but in contrast had a negative partial effect when the final model also included depth (Table 1).

Table 1 Best regression models of variation in water quality within the Doñana marsh in different sampling periods. Predictors were conductivity, depth and δ2H, plus dissolved nutrients for Chla. Symbols in brackets indicate positive (+) or negative (−) partial effects of predictor variables. Blank cells indicate that null models were the best. Significance range for p-value: * < 0.05; ** < 0.01; *** < 0.001
Full size table

As expected given the negative correlations with conductivity, nutrient concentrations in the marsh were highest in those sampling sites that are closest to the entry points of the Partido, Rocina and Sotos streams (Figs. 4, 5, 6, 7). This is further demonstrated when taking into account the additional data collected outside the sampling periods represented in Figs. 4, 5, 6, and 7. For example, the GM of NO3− is 0.22 mg N L−1 (95% CIs 0.10–0.49, N = 9, mean conductivity = 438.1 µS cm−1) at point A in Fig. 1 (37°4′18″N, -6°27′17″W) situated at 6 km from the nearest upstream greenhouse area within the Rocina/Sotos, compared to only 0.002 mg N L−1 (95% CIs 0.0004–0.01; N = 8; mean conductivity = 1487.4 µS cm−1) at point B (37°2′21″N, − 6°24′52″W) 5 km further downstream within the marsh.

We used multiple regression to explore the relationship between chla concentrations in the marsh, with the above physicochemical variables as well as the six nutrient variables as predictors. Best models included at least one nitrogen variable together with conductivity for two sampling periods (Table 1). Conductivity had a significant positive partial effect on chla in May, yet an insignificant negative partial effect in February (Table 1).

Surface water quality status

Water quality status based on nutrient concentrations

Comparing all our nutrient data with the water quality classification of OECD (2007) (Table 2), we observed extensive differences between waterbodies in water quality. The Partido stream had the lowest quality, and showed the highest proportion of concentration values (> 63.8%) classified within the three worst water quality categories. In particular, 93.4% of PO4−3 and TP values, 71.4% of TN values and 69.2% of NO2− values in the Partido were within the two worst classes (IV and V) not suitable for fish life.

Table 2 Percentage of spot samples (ntotal = 337) for each habitat (Partido (n = 91), Rocina/Sotos (n = 78), Guadiamar (n = 32) and Marsh (i = 136) corresponding to five use classes defined by OECD (2007). Each class defines which water uses are supported for a given water quality (uses: ecosystem functioning, fish breeding/protection, drinking water supply, bathing/recreation, irrigation, industrial water use, power use, mineral extraction and transportation). Class I and II are the only ones considered compatible with biodiversity functioning. Class I may be considered equivalent to “high status” under the WFD, i.e. a virtually undisturbed, natural system. Class V is equivalent to the WFD’s “bad status”. No fish are expected to survive Classes IV and V. Figures in brackets are %s when excluding data from June to September (hot, dry months), so that ntotal = 304, Partido (n = 81), Rocina/Sotos (n = 72), Guadiamar (n = 29) and Marsh (n = 122)
Full size table

Water quality was consistently higher in the marsh, followed by the Guadiamar stream. In the Rocina/Sotos between 14.1% and 48.7% of samples were in the worst two classes (IV and V) for four parameters (Table 2). In the marsh, for all parameters except TP, more than 80% of all samples were classified within the best two water quality categories (class I and II), and  < 9% in the worst two classes (IV and V). For TP, 16.8% of values were in the worst two classes, these being recorded close to stream entry points and in the western area of the marsh (Fig. 7). General results were the same when data from June to September were excluded (Table 2).

Water quality status based on chla concentrations

Following the chla-based classification of trophic level from OECD (1982), 20.3% of all samples (63% of which were from streams; nstreams = 109; nmarsh = 63) were oligotrophic, 44.2% were mesotrophic, 26.8% were eutrophic and 8.7% were hypereutrophic. The majority (73%) of the 15 samples classified as hypereutrophic (chla > 25 µg L−1) were found in streams (Fig. 8). Similarly, 76% of 45 samples classified as eutrophic (chla > 8–25 µg L−1) were recorded in streams (nstreams = 34 vs. nmarsh = 11). Most cases of hypereutrophy were recorded during February or May (Fig. 8).

Discussion

In this study we addressed the physico-chemical quality of surface waters and evidence for anthropogenic eutrophication in Doñana, the largest wetland complex in western Europe and one of the most important in the Mediterranean region (Green et al. 2018). We recorded persistent poor water quality in the major entry streams, despite the protection of the core area of Doñana as a WHS since 1992, and the earlier protection of the entire catchment within a Biosphere reserve. Indeed, concentrations of nitrites, ammonia and other nutrients are often above toxic thresholds for aquatic life. This is likely to be related to the five-fold expansion of agriculture in greenhouses observed over two decades since 1995, the increasing human population in the region and the poor treatment of urban wastewaters. These anthropogenic impacts may also explain the differences we observed in nutrient concentrations between stream subcatchments, which were consistent in different seasons and years, despite major fluctuations in precipitation and surface water distribution. The more strictly protected marsh has better water quality than the less protected streams, but the edge of the marsh is clearly impacted by anthropogenic nutrient inputs and is undergoing eutrophication (see also Paredes et al. 2019). This may compromise the biodiversity conservation objectives of this WHS, and suggests a need for policies aimed at reducing nutrient concentration in the catchment itself (Chimney and Goforth 2006; Tonderski et al. 2017). All our initial hypotheses were supported, except that we did not identify consistent seasonal patterns in nutrient and chla concentrations. This is likely related to the complex hydrology of the study area, and the influence of high temporal variation in precipitation and irregular drought periods on flow rates and nutrient cycling. Most of our data were collected during hydrological years 2015 and 2016, between which the timing and extent of flooding varied greatly (Figs. 4, 5, 6, 7, 8), as did the precipitation pattern (Fig. S1). Stream flow rates at a given sampling point varied by more than an order of magnitude, even within the space of 10 days (Table S3). The influence of any predictable seasonal trends (e.g. in temperature) on water quality was therefore hidden under the greater influence of the variable timing and extent of flooding events.

Unfortunately, the poor water quality and negative trends in Doñana may reflect trends at a broader scale across southern Spain and much of the Mediterranean region. In many areas of Europe, water quality has improved in recent decades in major rivers with regard to BOD (Biochemical Oxygen Demand), nitrates, phosphates, and ammonium, although a quarter of European groundwaters have poor status, mainly due to nitrate concentrations (European Environment Agency 2015). In contrast, Spain stands out at a global level as one of the countries most affected by anthropogenic phosphorus inputs to freshwaters, mainly due to agriculture (Mekonnen and Hoekstra 2018). The south of Spain is also generally heavily affected by anthropogenic nitrogen pollution, but does not stand out at a global scale as much as it does for phosphorus (Mekonnen and Hoekstra 2015).

Overall, Mediterranean wetlands have experienced rapid degradation and biodiversity loss in the past 50 years due to socio-economic pressures and human population growth (Perennou et al. 2020). Like much of the Mediterranean Basin, Doñana is located in a region of high water stress caused by the synergies between climatic conditions (e.g. scarce precipitation and long dry periods) and exploitation of a high proportion of available water resources (Hofste et al. 2019). This inevitably puts pressure on natural freshwater ecosystems dependent on these same resources. Information about water quality in other Mediterranean wetlands is currently limited, because little monitoring of pollutants other than nutrients is carried out, and few data are publicly available (Haener 2008). The limited data available (e.g. UNEP 2016; MEDPOL 2015) suggest a general decline in water quality in other wetlands across the Mediterranean Basin in recent decades, a disturbing trend which is predicted to continue (Veolia and IFPRI 2015), and which matches our findings in Doñana.

Drivers of eutrophication in Doñana: agriculture and urban development

Agricultural intensification and urban development have greatly increased in the surroundings of Doñana since the 1960s, owing to rapid economic growth (Ales et al. 1992). In response to the imminent threat from land conversion, road construction and economic development, protection of Doñana began in 1963 when the World Wildlife Fund (WWF) and the Spanish government purchased and protected 6800 ha of what since has become a much larger National Park (Green et al. 2018). However, less attention has been paid to limiting development in the catchment area, where land conversion into intensive agrosystems has since occurred, with clear impacts on water quality and quantity (Zorrilla-Miras et al. 2014b). Our results for greenhouse expansion reflect broader trends in recent decades across the Guadalquivir River Basin, with a shift into more intensive, high value irrigated crops largely for export to other European countries (Expósito and Berbel 2017; Rodríguez and De Stefano 2012). In the Doñana catchment, this growth has partly been driven by the destruction of native forests and other vegetation, and partly by the conversion of vineyards and other crops with lower impacts on water quality and quantity (Junta de Andalucía 2005). In the case of the Guadiamar catchment, this trend has been partly reversed since 1998 (see below), and this is likely to explain the higher water quality we recorded in this stream.

Agricultural intensification, increasing car ownership and the expansion of tourism (Fernández-Delgado 2017) have led to urban expansion and human population increase in the Doñana catchment since 1990, generating more urban wastewaters which are often poorly treated. On the one hand, as for many waste water treatment plants (WWTPs) in Andalusia, many of those around Doñana do not comply with the EU waste-water directive (https://uwwtd.eu/Spain/). On the other hand, the many agricultural workers in greenhouses produce waste which is not sent to WWTPs. Furthermore, the town and shrine of El Rocío on the very edge of the Doñana marsh at the outflow of Rocina stream has a WWTP designed for only c.5000 people, yet the town receives over a million visitors during the week-long annual pilgrimage in May or June.

As a consequence of increasing anthropogenic nutrient loading in the catchment, Doñana National Park was declared a sensitive area to eutrophication in May 1998 under the European Urban Waste-water Treatment Directive 91/271/CEE (EEC 1991). Then in 2008, part of the marshland and catchment area were declared vulnerable to nitrate from agricultural origin by the Andalusian Government (Decree 36/2008), according to the European Nitrate Directive 91/676/EEC (EEC 1991b). These and other relevant Directives such as the WFD (EC 2000) have been insufficient to prevent further eutrophication in the catchment, as indicated by the 20 fold increase in phosphate concentration from 1996 to 2010 at a handful of sampling points reported by Espinar et al. (2015), which mirrors the expansion of greenhouses we documented in the Partido, Rocina and Sotos catchments (see also WWF 2019a).

Heterogeneity within the catchment: explaining water quality differences between streams

As predicted, we found that sampling points closer to greenhouses or WWTPs (Rocina/Sotos, Partido) showed higher nutrient concentrations than points located further away from any intensive anthropogenic activity within the marsh. These results were in agreement with those previously reported for nitrogen pollution in some of the same sampling points over a shorter period (Paredes et al. 2019). Stable isotope analysis of dissolved nitrates from the Rocina, Sotos and Partido suggested that the major origin is fertilizers, followed by wastewaters (Paredes et al. 2020). Samples collected within the Guadiamar catchment (not studied by Paredes et al. 2019) had intermediate nutrient levels.

The marked reduction in the flow rate of the Rocina/Soto streams in recent decades due to groundwater extraction for the greenhouses (Guardiola-Albert et al. 2011) may have contributed to high nutrient concentrations owing to reduced dilution of pollutants. Moreover, sampling points in the Partido are the most directly affected by point sources of pollution due to their shorter distance to upstream WWTP effluents and relatively high flow rates (see Table S3). The latter reduces the time for nutrient retention and transformation by in-stream biogeochemical processes (Peterson et al. 2001). The Guadiamar and the Rocina/Sotos streams are mainly affected by diffuse sources (i.e. fertilizers and decentralized waste from agricultural workers spread around the area) since WWTPs are either located further away from the sampling sites than in the Partido (Guadiamar) or they are not present at all (Rocina/Sotos). The influence of wastewaters in the Partido catchment has also been demonstrated by analyzing stable isotopes in emergent plants growing along the stream banks (Paredes et al. 2019). High concentrations of pharmaceuticals in the Partido stream and areas of marsh close to El Rocío village, some of which exceed toxic thresholds for aquatic invertebrates by an order of magnitude, have also demonstrated the influence of urban wastewaters (Camacho-Muñoz et al. 2013, 2010).

Riparian vegetation also acts as a buffer zone along the stream margins by reducing nutrient inputs into the streams, stimulating in-stream nutrient retention and attenuating flood runoff (Pinay et al. 2018; Weigelhofer et al. 2012; Wilkinson et al. 2019). In the Partido, riparian zones are considerably more degraded than in the Rocina/Sotos and Guadiamar catchments, thus increasing the vulnerability of surface waters to nutrient inputs by runoff (see SIOSE 2013 for detailed land use information). The Guadiamar catchment has benefitted from restoration since the 1998 mine spill which directly affected this area (Madejón et al. 2018). Actions such as abandonment of surrounding agricultural areas in the surroundings and the creation of a “green corridor” with extensive Riparian vegetation (Ontiveros et al. 2013) have likely reduced nutrient inputs into the Guadiamar River.

In general, the nitrate concentrations recorded are not particularly high and would not alone indicate “bad status” of surface waters, even in the Partido stream (Table 2). This probably reflects high rates of denitrification within the groundwater before discharge into streams, as indicated by microbial activity, greenhouse gas emissions (Tortosa et al. 2011) and N and O isotopes in dissolved nitrate (Paredes et al. 2020). In contrast, nitrite concentrations are frequently well above levels lethal to fish life in the Partido and Rocina/Sotos streams, and ammonia concentrations are also high (a classic indication of waste-water influence). For example, the Freshwater Fish Directive 78/659/EEC (EEC, 1978) identified a guidance limit of 0.009 mg L−1 N-NO2− and 0.15 mg L−1 N-NH4+ for cyprinid waters, thresholds greatly exceeded in the Doñana streams (Table 2).

Chla concentrations were generally higher in the streams than in the marsh, but the differences were much less marked than for nutrient concentrations. There was much spatial variation in chla concentration within and between streams, largely connected with varying depth. Towards the mouth of the Rocina, the stream widens and a dam creates a deeper pool where residence time is longer and phytoplankton can grow better and accumulate. This is not the case in the Partido for example, where chla concentration in the water column was particularly low. Measuring chla in benthic algae (i.e. periphyton) rather than in planktonic ones may have provided results that were more strongly related to nutrient concentrations in streams, since benthic algae often represent an important fraction of the primary productivity in shallow streams and rivers (Jarvie et al. 2003). In the Partido, lower chla concentrations may also be due to reduced amount of phytoplankton due to lower tolerance to continuous high NH4+, NO3−or NO2− concentrations in the water column (Domingues et al. 2011; Glibert et al. 2016). We found evidence that high NH4+may reduce Chla concentration in the marsh (Table 1).

Temporal and spatial variation in water quality within the WHS

In general, we found considerably lower nutrient concentrations within the marsh protected within the WHS, than in catchment streams. This difference may be partly explained by high nutrient removal capacity of the marsh due to biogeochemical processes such as plant uptake, P adsorption to sediments or nitrogen removal by denitrification (Fisher and Acreman 2004), but also due to higher distance from the polluting sources in the catchment. Only a few sampling sites in the marsh showed high N and P concentrations, particularly in the north-west under the influence of catchment-derived nutrients entering from adjacent agricultural and urban areas (Fig. 4, 5, 6 and 7).

According to our PCA and linear regression models, conductivity, depth and stable isotopes (δ2H) play a key role in explaining the variation and nutrient concentrations between sites within the Doñana marsh. We found a negative effect of depth and a positive effect of δ2H on nutrient concentrations consistent with water evaporation effects that cause decrease in depth, and increase in δ2H and solute concentrations (Fellman et al. 2011; Gat, 2010). Although we might therefore assume a positive effect of conductivity on nutrient concentrations due to evapoconcentration of solutes, we repeatedly found a negative effect because conductivity is a proxy for distance to stream entry points. There is a previously documented conductivity gradient in the marsh, where the lowest conductivity values are found in the western and north-western areas of the marsh due to groundwater discharges and stream inputs, respectively, and the highest conductivity values are found in the eastern and south-eastern areas (Fig. S5). This negative correlation between nutrients and conductivity is also consistent with a nutrient purification effect as surface water flows further into the marsh, indicating its major ecosystem service as a “green filter” to purify contaminated waters from the catchment. The marsh is dominated by emergent vegetation (Lumbierres et al. 2017), as are constructed wetlands used for water purification. However, part of the reduction in nutrient concentrations towards the centre of the marsh can be explained by dilution by the water reaching the marsh by direct precipitation. Unfortunately, it is not possible to estimate the contribution of this dilution effect (but see 2.1 Study area for available information). Despite its capacity for purification, the marsh has experienced a long-term increase in phosphate concentration since 1995 which is thought to have facilitated the invasion of the nitrogen-fixing alien floating fern Azolla filiculoides since 2001 (Espinar et al. 2015). The expansion of this and other floating plants during eutrophication is likely to promote anoxia and reduce aquatic biodiversity in Doñana, as observed in many other wetlands (Green et al. 2017).

Not all spatial variation in nutrient concentrations within the marsh can be explained by the influence of streams, groundwater and evaporation. For example, particularly high PO4−3 and TP concentrations were recorded in the marsh in May 2015 far from any entry area, in the proximity to a large heronry, suggesting the occurrence of in situ animal-derived nutrient loadings (i.e. guanotrophy) rather than catchment-derived nutrients. Guanotrophy has often been reported in freshwater systems as a direct consequence of waterbird excreta where birds concentrate for breeding or roosting (Hahn et al. 2007; Martín‐Vélez et al. 2019). Indeed, guanotrophication is responsible for the decline of nesting oak trees in part of Doñana National Park, and expansion of heronries causes intense nutrient inputs into specific areas of the Doñana marsh where colonies are located (Fedriani et al. 2017; Ramo et al. 2013).

Furthermore, the Doñana marsh has high concentrations of wild and domestic ungulates (Junta de Andalucía 2015) which are controversial given evidence of negative impacts on insect communities, aquatic vegetation and nesting birds (Sharps et al. 2015; Verdú et al. 2018). Ungulates promote eutrophication through treading and their faeces (Declerck et al. 2006), and this may partly explain the relatively high phosphorous concentrations in the water column of some sites in the marsh. Nevertheless, our study suggests that any eutrophication problems in Doñana caused by bird or ungulate populations are of far less overall significance than the anthropogenic eutrophication of the entry streams.

Eutrophication management: gaps, failures and future perspectives

Lack of nutrient reference values for Doñana surface waters

We used water quality classes which were derived for wetlands in general in the Mediterranean region (OECD 2007), because the national reference system available that has been developed to assess waterbodies such as the Doñana shallow seasonal Mediterranean marshes is lacking key nutrient reference values.

Based on the OECD (2007) classification, the marsh showed considerably higher water quality than most entry streams. However, the status for phosphorus in the marsh is of concern, and only a third of samples were equivalent to a “high status” for Total P (Table 2). These results are a strong indication that eutrophication is ongoing, and are consistent with the reported increase in PO4−3 concentration since 1995 (Espinar et al. 2015). However, Serrano et al. (2017) argue that high Total P does not necessarily mean “poor” water quality status in naturally eutrophic shallow aquatic systems such as the Doñana marsh. Phosphate concentrations in the water column of shallow systems are not only the result of external inputs, and the chemical equilibrium between P- adsorption and release in the sediment plays a key role in the overall dissolved inorganic P availability in the water (Golterman 2004). Thus, our results demonstrate the need for developing specific eutrophication indices for shallow, temporary Mediterranean systems which may better distinguish between external (“anthopogenic”) and internal (“natural”) P loadings.

Failures in the past

Doñana is an iconic European wetland in which important conservation efforts have been made at national and international levels to protect the natural values against the impacts of anthropogenic land-use changes that have dominated the European landscape over the past century. Nevertheless, these measures have not been sufficient to protect the area from eutrophication. Concern has been expressed for decades about nutrient pollution in the marsh and its tributary streams due to changing land use and an increasing human population in the catchment (Serrano et al. 2006; Junta de Andalucía 1992; Novo 1994). Spanish administrations have historically focused on engineering based solutions, and projects to modify water quantity through earthworks have been favored ahead of projects focusing on improving water quality (Méndez et al. 2012). For example, in an extensive restoration project in Doñana that followed a toxic mine spill in 1998, various engineering works increased the supply of water into the marsh by restoring surface water inputs diverted during past transformations for agriculture (García-Novo et al. 2007; García-Novo and Marín 2006), despite concerns about the quality of some of those sources (Serrano et al. 2006) and without clear criteria for acceptable quality.

The eutrophication of the Doñana entry streams and those areas of marsh surrounding the stream mouths has accelerated in recent decades, and now often reaches levels incompatible with biodiversity conservation. Our study of the growth of greenhouses demonstrates one of the most likely causes of this increase. As has also occurred in many other protected wetlands (Peterson et al. 2001; Pusch et al. 1998), in the case of Doñana a basin has been strictly protected whilst the catchment area has not, leading to increases in pollution sources and land-use changes within the catchment that translate into anthropogenic nutrient export into the protected wetlands.

How is it possible that such eutrophication has been permitted to occur in arguably Europe’s most iconic wetland complex? One explanation lies in the huge economic interests involved in part of what has become Europe’s largest strawberry culture area (in Huelva province). Another likely explanation is that conservation interests since the foundation of the National Park have largely focused on “heroic megafauna” such as the Iberian lynx and the Spanish imperial eagle (Fernández-Delgado 2017), which do not act as indicators of water quality. Remaining fish biomass has long been dominated by alien species such as the carp Cyprinus carpio and mosquitofish Gambusia holbrooki, and relatively little attention has been paid to aquatic fauna other than waterbirds, for which Doñana is particularly important (Ramo et al. 2013; Rendón et al. 2008). However, waterbirds have high plasticity and are relatively insensitive to eutrophication (Amat and Green, 2010; Almeida et al. 2019). Nevertheless, those wintering waterbird species with negative population trends in Doñana, such as herbivorous ducks, are those most likely to be impacted by eutrophication in the marsh (Rendón et al. 2008). Direct monitoring of eutrophication has previously been neglected in Doñana, hence the importance of our study.

Conservation actions for the future

Despite the existence of the Biosphere reserve, since its declaration the Doñana area has been managed under a conservation versus development paradigm in which the strictly protected National Park (including the marsh) has been managed for conservation purposes. In contrast, the surrounding area (including the entry streams) has been largely earmarked for economic development (Zorrilla-Miras et al. 2014), hence the expansion of greenhouses. A similar paradigm has led to protection of many shallow lakes in Andalusia without effective protection of their catchment areas (Rodríguez-Rodríguez et al. 2012).

To mitigate the excess of downstream anthropogenic nutrient transport into the Doñana WHS, a radical shift in land and water management is required towards a more holistic approach. Water management should be based on multi-scale governance systems capable of meeting the basic demands of a variety of stakeholders and beneficiaries on local and larger scales, without compromising biodiversity conservation (Defries and Nagendraz 2017). Water purification of anthropogenic inputs is necessary before water reaches the marshland, instead of relying on the marsh itself to provide these services. This requires specific actions such as restoring and increasing the area of riparian buffer zones along the affected streams, reversing agricultural encroachment up to the edge of streams, and other measures to mitigate soil erosion. Such riparian buffers can be highly effective (Weigelhofer et al. 2012; Teufl et al. 2013) and are particularly required in the Partido and Rocina/Sotos catchments. Creation of a network of upstream constructed wetlands before the water enters the Doñana marsh, and the improvement of the existing wastewater collection and treatment would considerably reduce downstream nutrient transport. Furthermore, a higher spatio-temporal resolution of water quality monitoring is crucial, particularly to identify point sources and to understand the impact of events such as runoff from agricultural lands during heavy rains, or El Rocío pilgrimage celebrations. In addition, the extensive illegal expansion of greenhouses should be reversed, including closure of illegal wells, acquisition of water rights by the Government, or declaring the underlying aquifer as overexploited (WWF 2019b). Given the synergy between eutrophication and climate change, it is essential that plans to reduce nutrient inputs to Doñana are ambitious, and this synergy should be recognized and accounted for by adopting a “safe operating system” approach to management of the entire Doñana catchment (Green et al. 2017).

Conclusions

In conclusion, we found strong evidence that eutrophication is a major and ongoing problem in Doñana, particularly in the streams flowing into the Doñana WHS and the areas of marsh around the mouth of streams. Nutrient concentrations were consistently lower within the marsh partly due to its “purification capacity”, and the greater distances to agricultural and urban sources of nutrients.

The current situation is partly the consequence of decades of polarized management, in which the strictly protected National Park (including the marsh) has been managed for conservation whereas the catchment area (including the streams) has been prioritized for economic development, leading to agricultural intensification and human population growth. Thus, the high amount of nutrients entering the catchment, together with decreasing water quantity due to groundwater abstraction for agriculture, have seriously compromised surface water quality, and threaten the ecological value of this emblematic wetland. These conflicts are likely to be common in other Mediterranean wetlands. An urgent shift to more holistic, integrated management is necessary to improve water quality, and so ensure the long-term conservation of the Doñana wetland complex and its resilience against climate change.