Multiple Evidence for Climate Patterns Influencing Ecosystem Productivity across Spatial Gradients in the Venice Lagoon

Abstract

Effects of climatic changes in transitional ecosystems are often not linear, with some areas likely experiencing faster or more intense responses, which something important to consider in the perspective of climate forecasting. In this study of the Venice lagoon, time series of the past decade were used, and primary productivity was estimated from hourly oxygen data using a published model. Temporal and spatial patterns of water temperature, salinity and productivity time series were identified by applying clustering analysis. Phytoplankton and nutrient data from long-term surveys were correlated to primary productivity model outputs. pmax, the maximum oxygen production rate in a given day, was found to positively correlate with plankton variables measured in surveys. Clustering analysis showed the occurrence of summer heatwaves in 2008, 2013, 2015 and 2018 and three warm prolonged summers (2012, 2017, 2019) coincided with lower summer pmax values. Spatial effects in terms of temperature were found with segregation between confined and open areas, although the patterns varied from year to year. Production and respiration differences showed that the lagoon, despite seasonality, was overall heterotrophic, with internal water bodies having greater values of heterotrophy. Warm, dry years with high salinity had lower degrees of summer autotrophy.

1. Introduction

Current global trends in greenhouse gas emissions will lead to at least 1.5 °C of atmospheric warming by 2030–2050 and at least 3 °C by 2100 [1,2], and half the world’s population already experiences temperatures 1.5 °C warmer than preindustrial [3]. The rise of sea surface temperature (SST) is tightly connected with these changes. According to observation data beginning in 1880, the global ocean has warmed by 0.005 °C year⁻¹ [4,5]. However, in recent years (1986–2005), considerably stronger increasing tendencies for SST have been estimated (0.02 °C year⁻¹ [6]). This increase in SST is not homogeneous, with estimates for 1985–2006 in the Mediterranean of circa 0.03 °C (peaking up to increases of 0.16 °C in summer months [7]). In aquatic environments, coastal areas may be more vulnerable than deeper waters. Transitional coastal areas are also of a heterogeneous nature, meaning that climatic impacts in these areas are context-specific and may be difficult to pinpoint [8] and to generalize depending upon local ecosystem characteristics such as depth and flushing times [9]. According to [9], lagoons were found to be the fastest-warming ecosystem types, falling within the “warm spot” of shallow average depths and short–medium flushing times.

Climate change can have important consequences for coastal ecosystems, including on nutrient dynamics (e.g., shifting DIN: DIP seasonality [10]) and oxygen depletion, due to the effects of climate on solubility and stratification [11]. The effects of climate change can increase risks of eutrophication due to light reduction at the bottom that may lead to the decay of benthic primary producers [12], further aggravating the situation. More than that, plankton blooms are directly influenced by temperature in terms of bloom magnitude, timing and phenology of species blooming [13]. All of this can influence primary productivity, an essential part of ecosystem functioning necessary for supporting higher trophic levels and a large player in aquatic metabolism with a key role in carbon cycling [14].

The Venice lagoon, located in the northwest Adriatic Sea, is a very peculiar habitat compared to the rest of the Mediterranean, being a large lagoon (550 km²; 632 × 10⁶ m³) with a tidal regime unusual for the Mediterranean (tidal range of circa 0.84 m) [15]. The lagoon is considered to be a transitional coastal system due to the heterogeneity in salinity and water residence time [16,17], receiving an annual average of about 35 m³ s^–1 of freshwater into the lagoon from twelve tributaries [18], and its depth ranges from 1.52 to 39 m [19]. The heterogeneity of the lagoon is also derived from the heterogeneity of localized anthropogenic pressures (e.g., nutrient discharges, chemical industry, boating and navigational channel maintenance, aquaculture, etc. [18,20,21]). This lagoon is also highly studied, and long term data from surveys and multiparametric sensors are available [22], making it the optimal study system to understand the recent climatic evolution and its effects on productivity. Through analysis of environmental data collected in the last decade at multiple stations in the Venice lagoon, the aims of this study were to (1) understand temperature dynamics in the past decade (2) understand ecological consequences in terms of ecosystem productivity and (3) consider spatial heterogeneity of these responses.

2. Methods

2.1. Study Site

The Venice lagoon is a semi-enclosed bay in the Adriatic Sea northeast of Italy, and it is the largest wetland in the Mediterranean Sea. It has a surface area of around 550 km². About 11% is permanently covered by open water, while around 80% consists of mud flats, tidal shallows and salt marshes. It is connected to the Adriatic Sea by three inlets: Lido, Malamocco and Chioggia. Based upon its characteristics, mainly salinity and water residency time, the lagoon has been subdivided into multiple water bodies (see

Table 1. Water body classification, where Polyhaline equals 18–30, Euhaline > 30.

Type	Code	Name	Surface (km2)
Polyhaline confined	PC1	Dese	18
Polyhaline confined	PC2	Millecampi Teneri	37
Polyhaline confined	PC3	Val di Brenta	7
Polyhaline confined	PC4	Teneri	10
Euhaline confined	EC	Palude Maggiore	40
Euhaline not confined	ENC1	Centro Sud	106
Euhaline not confined	ENC2	Lido	10
Euhaline not confined	ENC3	Chioggia	3
Euhaline not confined	ENC4	Sacca Sessola	24
Polyhaline not confined	PNC1	Marghera	28
Polyhaline not confined	PNC2	Tessera	25

).

2.2. Environmental Data

The “SAMANET” network consists of a network of ten, multi-parametric sensors deployed across the whole lagoon (see Figure 1). These sensors collect half-hourly data. Data have been made available for the 2008–2019 period for water temperature (°C), or T_w, oxygen saturation, or DO (%) and Salinity, S. From this dataset we extracted a time series of hourly mean values for each water body where buoys were present. Water bodies with no buoys were not considered for the purpose of this study. For the purpose of the analyses of this study, the water bodies PC1 and PNC2 were considered in this study to be a single water body because the ‘Ve-7′ buoy lay on the border of the two water bodies (see map, Figure 1). When multiple buoys were present for each water body, a mean value was taken, which allowed us to deal with some missing data present in the timeseries.

2.3. Climate Effects

To assess similarities between the seasonal evolutions of water temperature in the 2008–2019 decade, data for all buoys were pooled together to avoid missing entire periods due to missing data, and a median water temperature value for each fortnight was calculated to input into K-means clustering algorithms. The same was done for salinity values. To understand whether spatial differences in climate occurred for each year, separate median water temperatures for each water body where data were available in each season were calculated to input into the algorithm. The choice of seasonal medians was made to avoid losing a whole basin in case of one missing data point, which happened if data were averaged fortnightly or monthly due to the large number of missing data.

K-means clustering [23] was implemented in R. The optimal number of clusters was chosen based on the elbow method [24] while also ensuring that there was no overlapping between clusters when plotted in the multidimensional space.

2.4. Ecological Responses

Production and Respiration in the System

Using the model developed by Ciavatta et al. (2008), based upon Equation (1) of oxygen mass balance and following exactly the protocol described in the paper, the oxygen (DO) data for each water body were used to calculate respiration (r) based on autoregressions of nighttime hourly data and production (pmax) and based on daytime hourly data on DO. Days characterized by either negative values of r and/or above the 0.95 quantile were removed from the dataset. For the r calculation, based on linear model regression, the standard error of the regression was also outputted.

$$\frac{dDO}{dt}=pmaxfe-r+\mathrm{k}({\mathrm{DO}}_{\mathrm{sat}}-\mathrm{DO})$$

(1)

Input data were taken from the SAMANET network. DO values, which were provided as SAT (%), were firstly converted to mg L⁻¹ by dividing the values by the oxygen saturation (DO_sat) which was calculated following Equation (2), where T_w is water temperature and S is salinity:

DO_sat (t) = 14.6244 − 0.367134T_w(t) + 0.0044972T_w(t)² − 0.0966S(t) + 0.00205S(t)T_w(t) + 0.0002739S(t)²

(2)

The model also requires irradiance data in kJ h⁻¹, which are needed also to calculate the function fe, i.e., the Smith–Talling formulation of the production–irradiance relationship. These data were obtained from ARPAV as hourly data which were made available until July 2018, and thus, for this part of the research, the timeframe 2008–2017 was considered.

Means of r and pmax were then estimated for each fortnight, allowing us to deal with some of the missing data. The relationship between pmax and temperature was also modeled with ordinary least squares regressions, finding the best fitting curve between a linear fit and polynomial regressions of various degrees using a forward step approach, stopping when the first non-significant term was found, and comparing the models based on the Akaike Information Criteria (AIC).

The correlations among productivity and respiration, i.e., pmax and r, and inorganic and biotic variables were investigated using Pearson’s coefficients. The following variables were taken into account: total number of phytoplankton cells, concentration of phaeopigments and chlorophyll-a data (https://www.arpa.veneto.it/dati-ambientali/open-data/idrosfera/acque-di-transizione/acque-di-transizione-laguna-di-venezia-monitoraggio-ecologico, accessed on 1 February 2021), inorganic nutrients (nitrogen (NH₄⁺ + NH₃⁺ + NO₃⁻ + NO₂⁻), phosphorus (PO₄³⁻)) (https://www.arpa.veneto.it/dati-ambientali/open-data/file-e-allegati/soaml/laguna-di-venezia/tss-nutr-acqua, accessed on 1 February 2021). All variables were all log transformed to obtain a Gaussian-like distribution, which was visually checked with histograms. Field measurements were only available from 2011 and were collected seasonally (in February (2013–2015)/March (2016–2018), May, August and November); thus, only data from the model belonging to these months and years were used in the correlation analysis.

The same procedure of clustering to identify temporal and spatial effects using the k-means approach described in the “Climate” paragraph was applied to the clustering of r and pmax.

Finally, the model estimates relative to productivity for each fortnight (pmax and r) were transformed from mg O l⁻¹ h⁻¹ to mmol C day⁻¹ m⁻² (to P and R) in order to quantify net ecosystem primary productivity following Ciavatta et al. (2008)’s procedures, where R is the sum of the mean r value in the fortnight for all hours in a day (24), whereas P is computed following Equation (3) where fe is the average of the function for the fortnight and dt refers to the average length of daylight on a given fortnight.

$$\int pmax fe dt$$

(3)

From the average daily values, the annual values of NEPP (Net Ecosystem Primary Production, calculated as P − R) were then calculated and transformed from mmol C to t CO₂ in order to estimate whether and to what degree the pelagic compartment of the Venice lagoon acted as a source or sink. To obtain comparable data, results will be presented both by unit area (ha) and for whole water bodies and lagoons.

3. Results

3.1. Climate: Temporal Effects

For temperature, the k-means clustering identified four temporal clusters (Figure 2):

Cluster 1: 2014—milder year

Cluster 2: 2009, 2010, 2011, 2016—”normal” seasons

Cluster 3: 2012, 2017, 2019—colder in winter and warmer long summer (end of June to end of August)

Cluster 4: 2008, 2013, 2015, 2018—warm peak at the end of summer (August)

With regards to salinity, the k-means clustering identified three temporal clusters (Figure 3):

Cluster 1: 2008, 2012, 2016, 2017—high salinity between 28–33

Cluster 2: 2010, 2014, 2019—medium salinity values with low values down to 25 in spring (May) and winter (November–December)

Cluster 3: 2009, 2011, 2013, 2015, 2018—medium salinity values, with lower values averaging around 27 in early spring (March–April) and upper values around 32 in late summer–autumn (August–September) and winter (December).

3.2. Climate: Spatial Effects

Maps of temperatures as identified with k-means are presented in the Supplementary Materials. Lagoon temperatures were spatially heterogeneous, but the spatial differences based on k-mean clustering varied between years and there was no consistently warmest or coldest water body for the whole decade. The area near Marghera (PNC1) was the warmest in spring and summer 2008–2009 and coldest in the same winters, while in autumn, it was warmest in 2008 but coldest in 2009. In 2011, 2012, 2014 and 2016 the area of Millecampi (PC2) was the warmest in summer; in 2010 and 2011 in summer and spring; and in 2012, 2014, 2015 and 2016 it was the warmest area in autumn, while in 2011, 2014, 2015 and 2016 it was cooler in winter. The northernmost water body, Palude Maggiore, (EC) tended to be the coolest in summer 2012, 2013 and 2014, but warmer than the rest in 2015 and 2016, and together with other northern water bodies, it became the warmest in the latest two years (2018, 2019). In winter 2011, 2014, 2015 and 2016 the unit made of center-north water bodies, Dese and Tessera (PNC2PC1), was the coldest (full maps in Supplementary Materials). Aside from 2012 and 2015, where the difference in summer mean temperature between basins was up to 4 °C, in the other years, a difference of up to 2 °C could be observed. In winter, spring and autumn the maximum intra-annual variability observed was up to 2 °C.

3.3. Production & Respiration

Mean (± SD) measured values of oxygen saturation across the whole dataset varied from 82.6 ± 29.4% (12 ± 4.3 mg L⁻¹) in August to 99.9 ± 17.4% (10.9 ± 2.1 mg L⁻¹) in March (values for each year and basin in the Supplementary Materials).

The calculation of r using nighttime data had a larger error during summer months (

Table 2. Estimates and residual standard errors in the estimation of r (mg O l⁻¹ h⁻¹) and estimates of pmax (mg O l⁻¹ h⁻¹). Values reported are monthly means for all water bodies and years pooled ± the standard deviations.

Month	January	February	March	April	May	June	July	August	September	October	November	December
Estimate (r) (mean ± SD)	0.729 ± 0.527	0.703 ± 0.505	0.663 ± 0.492	0.643 ± 0.440	0.588 ± 0.389	0.508 ± 0.355	0.437 ± 0.311	0.421 ± 0.323	0.483 ± 0.349	0.555 ± 0.408	0.707 ± 0.491	0.733 ± 0.502
RSE (r) (mean ± SD)	0.225 ± 0.194	0.297 ± 0.221	0.389 ± 0.315	0.594 ± 0.6	0.726 ± 0.520	0.950 ± 0.651	1.02 ± 0.609	0.977 ± 0.540	0.7 ± 0.502	0.503 ± 0.483	0.275 ± 0.212	0.219 ± 0.145
Estimate (pmax) (mean ± SD)	0.323 ± 0.199	0.381 ± 0.223	0.674 ± 0.292	1.041 ± 0.363	1.292 ± 0.435	1.380 ± 0.465	1.291 ± 0.451	1.170 ± 0.384	1.071 ± 0.334	0.831 ± 0.353	0.619 ± 0.258	0.336 ± 0.197

). In this period, less dark hours were available for the regression, and thus estimates were less accurate. Therefore, despite an apparent decrease in estimated respiration means from winter to summer, this was not statistically significant.

There was a slight positive correlation between pmax and values related to plankton abundance (number of cells, chlorophyll and phaeopigments) and phosphorous, while there was a negative correlation with values of nitrogen (Figure 4).

The best fitting temperature–pmax model was a third-degree polynomial (4th degree = N.S., AIC = 310; 3rd degree, p < 0.01, RSS = 62.1, AIC = 308; 2nd degree, p < 0.001, RSS = 63.5, AIC = 312; 1st degree, RSS = 67.7, AIC = 310).

Where: pmax = 0.01 + 0.001T+ 0.006 T² − 0.0001 T³ (Figure 5).

K-means grouping for pmax identified 5 groups (Figure 6). Cluster 3 made by the years 2009, 2010, 2011 and 2016 represented a typical curve with an increased production in spring that plateaued from the beginning to the end of summer and declined in autumn through winter. Cluster 4 formed by years 2008, 2015 and 2017 was represented by a similar spring dynamic but was followed by a decline in summer. The other clusters were made by single, anomalous years—2012, 2013 and 2014, which all reached higher production values but at different times and with different dynamics: 2012 had a second delayed peak towards the end of summer, 2013 had a high production peak in mid-summer, whereas 2014 had high productivity until June, followed by a steep decline.

No clear grouping could be identified for respiration.

When total production (P) and respiration (R) were taken in account, it could be seen that in winter months, respiration always exceeded production, while in summer, production exceeded respiration (see Figure 7). In terms of annual differences in production and respiration, the lagoon results to be always heterotrophic (

Table 3. Differences between production and respiration (P-R) expressed in tons of CO₂eq per unit area for each water body and each year. NA values represent missing data. The last column represents median values of the differences scaled up to the whole area of the lagoon based only on the water bodies for which data were available in this study. The last row represents the median value of the differences scaled up to the whole basin area. The last column represent the median value for the whole period considered for the whole lagoon. Negative values represent higher respiration. N.B. median values are calculated excluding NAs.

Year	EC (t CO2eq ha−1 Year−1)	ENC1 (t CO2eq ha−1 Year−1)	ENC2 (t CO2eq ha−1 Year−1)	ENC4 (t CO2eq ha−1 Year−1)	PC2 (t CO2eq ha−1 Year−1)	PNC1 (t CO2eq ha−1 Year−1)	PNC2PC1 (t CO2eq ha−1 Year−1)	Median Whole Lagoon Area (t CO2eq)
2008	−40.8	−19.9	−47.6	−15.5	−43.7	−26.5	−61	−152,296
2009	−22.2	−15.9	−37.5	−8.4	−32.5	−31.2	−32.3	−116,376
2010	−16.9	−21.4	−50.4	−0.8	−38.3	−27.3	−63.2	−101,833
2011	−42.2	−11.9	−30.9	0.6	−33.4	−24.1	−43.3	−115,482
2012	−10.4	−7.9	−29.4	−2.7	−39.3	−11.9	−26.2	−44,667
2013	−31.5	−12.3	−30.5	−5.5	NA	−25.8	−28.9	−102,033
2014	−38.2	−12.3	−20.2	−1.5	NA	−30.5	−58.5	−94,704
2015	−35.5	−14.4	NA	NA	NA	−30.6	−46.2	−122,725
2016	−23.6	−19	NA	NA	NA	−38.9	−18.9	−79,661
2017	−25.2	−28.9	NA	NA	NA	−30.9	−33.7	−111,570
Median whole water body (t CO2eq)	−9053	−20,426	−5557	−1472	−19,335	−9261	−22,784	−107,850

). Spatially, water bodies ENC4 and ENC1 had the lowest difference between yearly production and respiration per unit area (see Table 3).

4. Discussion

This study shows that in the Venice lagoon in the latest decade, there have been years where high water temperatures were recorded for long periods over the summer while winters had colder values than other years (2012, 2017, 2019) as well as a cluster of years when extremely high summer water temperatures were observed (2008, 2013, 2015, 2018). Hourly photosynthetic production interannual clustering for the Venice lagoon resembled the clustering for the temperature, with a mild year (2014) showing higher values of planktonic maximum oxygen production rates in spring that were at their lowest in late summers of years classified as warm (2008, 2015, 2017) whereas it followed a “typical” curve in years considered to have normal seasonality (2009, 2010, 2011, 2016). This study also showed that the median values of annual production (246 mmol C day ⁻¹ m²) were lower than the values of annual respiration (430 mmol C day ⁻¹ m²), in line with results from [25], and that the pelagic compartment of the Venice lagoon was overall heterotrophic [26]. A large variation around these median values was found that was due to both spatial and interannual difference, which will be further discussed below.

Climate

The fast warming of the Venice lagoon was already observed by [27] in data up to 2014, and in this study, the use of clustering techniques allowed identification of shifts in patterns of seasonality, manifesting in longer summers or more extreme summers and winters with steep changes between seasons. Data and model simulations with data from 2001 showed lower averages of summer temperatures with simulations showing summer temperatures of up to 27 °C [28], while in the warm cluster years here presented, summer values could be above this value for long periods. Furthermore, the clusters identified in this study were in agreement with the descriptions of local weather (ARPAV, “commenti meteoclimatici” (https://www.arpa.veneto.it/temi-ambientali/climatologia/dati/commenti-meteoclimatici, accessed on 20 February 2021). The isolated year in Cluster 1, 2014, presented with a relatively cool summer with precipitation above average. For Cluster 2 summers, 2009 was described as normal, 2010 was generally warm with some perturbation, 2011 was described as summer-like with some instabilities and 2016 was described as the norm. For the years in Cluster 3, 2012 was defined as one of the warmest summers, with an elevated number of days at high temperatures from June, and for its winter, January showed a cold spell and spring approached rapidly. 2017 was also described as an anomalously warm summer with no discrete peaks, with January presenting temperatures often lower than average climatological means. 2019 was described as very warm summer with two heat waves in June and July and again, a cold January. For Cluster 4, 2015 and 2018 were classified as two years with heat waves in the summer, 2008 had a start of summer with cool periods, which then significantly warmed up from the third decade of July with a sultry August. 2013 was described as generally fluctuating, with an intense and lengthy heatwave between the end of July and August, and generally, summer was 1 °C above average. 2015 was classified as generally warm, especially from July, with less precipitation. Finally, 2018′s summer was defined as one to be divided in two, with a relatively cool beginning and a warm end.

The spatial patterns of temperature revealed some slight heterogeneity in intra-annual temperatures observed, which was only significantly marked (up to 4 °C) in summer 2012 and 2015, two warm years. In the other cases, mean seasonal temperature had only a difference of up to 2 °C. Unfortunately, due to some missing data, it was not possible to look at spatial heterogeneity at a finer scale, and therefore, only average seasonal temperatures could be compared. This component of the study nonetheless adds some information to the variability that can exist in such habitats, which are highly dependent upon exchanges with seawater and to some extent, inputs from riverine waters. It should be noted that the differences shown in this study refer to the water column in areas where the stations are located and thus do not consider the variation that can occur in shallower areas and that are shown by [29] to be up to 8 °C at certain times of day depending upon the state of the tide.

The clustering of pmax values revealed somewhat similar patterns to the SST clusters, with 2009, 2010, 2011 and 2016 grouped together showing a “standard” curve, while warm years 2008, 2015 and 2017 showed reduced production. The relation of pmax and temperature showed a positive relation until 20 °C and then stabilized and decreased in line with phytoplankton growth and performance curves for temperate species (e.g., [30,31]), suggesting that summers that reach extreme high temperatures could be deleterious for primary production. In this study, we found significant correlations between production estimates and values from the survey, but the latter were sampled at a low frequency (trimestral). We would recommend continuing monitoring planktonic abundances and species composition, especially in the spring–summer months. However, we demonstrated that continuous monitoring of oxygen together with temperature could be useful to follow these relationships over a long time, and large spatial scales could aid understanding of the changes of production in conjunction with climate. Furthermore, we demonstrated the usefulness of reanalysis of available, long-term time-series from local monitoring schemes. In this study, with a large amount of data from the past decade, a direct relationship between annual NEPP and climate (intended simply as water temperatures) could not be identified, while there seemed to be some deleterious effects from high salinity years. This was likely due to the complex interaction between production, respiration and nutrients, which may vary according to precipitation patterns [32]. Furthermore, given the heterogeneity of the lagoon in terms of its salinity and the recognition of some more and less stable water-basins [16] and given the link between salinity, precipitation and discharge, it would be important to include the changes in this variable: Identifying climate patterns also based upon precipitation patterns could be a useful addition to the work presented in this study in order to further detangle the relationships with productivity.

The calculated production and respiration was in line with previous estimates (e.g., [25,33] and references within). The Venice Lagoon as a whole, in terms of its water column and thus phytoplankton compartment, was always heterotrophic when the whole year was considered. This is in line with [34], showing how the Mediterranean Sea area together with the North Atlantic Subtropical Gyre have a prevalence of heterotrophic communities. Spatially, some water bodies acted more heterotrophic than others. The fine balance between autotrophy and heterotrophy in estuarine systems was found to be impacted by nutrients and trace metals, with nutrient addition alone shifting the system towards autotrophy but with the addition of trace metals found to reduce this effect, and if trace metals were added alone, heterotrophic conditions were consistently met [35]. In the study presented here, the confined areas and the water body ENC2 (where the multiparametric buoy was located just outside of the island city of Venice) had greater degrees of heterotrophic behavior. These areas had long residence times and higher inputs from terrestrial sources. Nitrogen and phosphorous values recorded were higher in the more internal “polyhaline confined” areas, coinciding with these deriving from terrestrial inputs. Nonetheless, nitrogen was lower than previously recorded values (0.25–0.55 mg/L [36,37]), while phosphorus remained similar (0.002–0.012 mg/L [36,37]). Inorganic nitrogen appeared to be negatively correlated with pmax and with plankton variables, which could be due to the uptake of these nutrients by the planktonic communities during growth, leading to a nutrient limitation [38]. This would result in a lag between nutrient peaks and pmax peaks, visible with higher frequency nutrient samples which were not available here. On the other hand, respiration results positively correlated with inorganic nitrogen, coinciding with greater respiration in more internal areas where less water exchanges due to longer water residence time may have given rise to microbial proliferations and higher remineralization rates, a process involving the consumption of dissolved oxygen and leading to the release of inorganic nutrients back to the water column (e.g., [39]).

Other communities beside the planktonic one could play a role in the spatial differences observed in production estimates. Benthic primary producers such as seagrasses account for a relevant fraction of the primary production of the Venice lagoon [40] and are mostly located in water bodies ENC1 and ENC4 (http://cigno.atlantedellalaguna.it/maps/6/view, accessed on 20 February 2021), in proximity of two of the parametric buoys (Ve-3 and Ve-6), with ENC4 showing a seagrass community mostly dominated by Z. marina instead of C. nodosa, which appeared instead as the predominant species in the ENC1 water body in proximity of the Ve-3 buoy in the latest maps (dated 2010). Z. marina has a productivity circa 5 folds greater than C. nodosa [41]), which may explain why this water body acts as the least heterotrophic. Furthermore, the water body ENC1 is also an area of bivalve production, which may contribute to higher respiration rates from the biomass of mussels and clams and lower production if plankton is ingested, which may contribute to explaining some of the observed differences with ENC4. According to back- and forecasting scenarios by [42], the lagoon shifted from a scenario where production exceeded respiration (beginning of 1900) to one where respiration was the predominant player (current and future), and this could be attributed to the presence and activity of Manilla clam farming. The role of other species such as macrophytes, likely to become more abundant with rising temperatures (e.g., Ulvaceae [43]), should also be considered with ad-hoc abundance surveys for the whole lagoon, which may become more accessible with the use of new technologies (e.g., [44]).

Temporally, there were seasonal effects wherein from spring (May) to summer (August), production exceeded respiration in the water column. This switching between heterotrophic and autotrophic is therefore in line with seasonal phytoplankton abundances, peaking in spring and summer, as seen already for other temperate lagoons (e.g., [45]). In some years (2008, 2015, 2017) this effect was minor, with production exceeding respiration by less than 50 mmol C m⁻² day⁻¹, or below average (2016), with production exceeding respiration by 80 mmol C m⁻² day⁻¹ instead of in excess of 200 mmol C m⁻² day⁻¹ in all other years. These were warm years with high salinity, suggestive of an interaction between warming, low precipitation (high salinity) and a reduction in net primary productivity.

5. Conclusions

With this study, we demonstrated warming occurrences in recent years in the Venice lagoon and some concomitant effects on pmax. This study provides further evidence of the usefulness of using a production model that requires oxygen, temperature and salinity as inputs, things commonly measured at high frequencies in routine monitoring schemes. It would be useful to replicate this in other areas where such schemes exists and for such monitoring schemes to continue in order to identify changes as they occur. Furthermore, this study confirmed that the water column in the Venice lagoon is mostly heterotrophic and therefore calls for caution in considering pelagic coastal zones as “carbon sinks”. Using a combination of statistical and process-based models such as that presented in this study can provide very relevant information, not just for understanding the effects of long term climatic impacts but also short term variability and changes driven both by climate and local pressures (warm water industrial discharges, Mo.S.E (Modulo Sperimentale Elettromeccanico) barrier closures [27,46]). These datasets should be integrated with ad hoc studies and further installation of sensors, e.g., in areas with macrophytes and macroalgae and closer to the benthos, to better quantify CO₂ fluxes and ecosystem services of lagoons for climate change mitigation strategies.