Metabolism in a deep hypertrophic aquatic ecosystem with high water-level fluctuations: a decade of records confirms sustained net heterotrophy

Environmental variables

Samples were collected and incubated ∼monthly from July 2006 to December 2015. Environmental parameters were determined before and after oxygen incubation, and then averaged. Temperature, dissolved oxygen concentration and its saturation were determined at 1 m vertical intervals down the full water column by means of a multi-parametric probe (Yellow Springs Instruments model 6600). Secchi depth was measured with a standard disk. Chl-a concentrations were determined in water samples collected at the 1, 2, 4, 6 and 8 m depths in a Niskin bottle. Samples for Chl-a were filtered with 0.45 µm membranes (Millipore, Burlington, MA, USA), extracted with 90% acetone and determined with a spectrophotometer.

Oxygen dynamics

Oxygen evolution was quantified by incubations of light and dark bottles in situ following Wetzel & Likens (1991) and Valdespino-Castillo et al. (2014). The oxygen incubations approach to metabolism measurement has some important limitations, mainly that it is very labor intensive because it requires the deployment and recovery of the incubations, it can cause containment artifacts (mainly if the required incubation time is long), and issues related to upscaling if samples are not representative (Staehr et al., 2012). In contrast, it has the advantage of involving a direct measurement of the process, and it does not require a quantification of the air-water oxygen flux (Staehr et al., 2012). Because of the strong winds, incubations are better suited for monitoring metabolism in VB than free-water oxygen monitoring because of the impossibility of calculating the strong air-water exchange (Merino-Ibarra et al., 2008). This was verified in the analysis of the first year of monitoring (Valdespino-Castillo et al., 2014), during which the metabolic rates obtained were one order of magnitude higher than the errors of the method.

Samples and measurements were taken monthly at a central station (cf. Valeriano-Riveros et al., 2014 for the exact location); this is supported by the verification of horizontal homogeneity in the mixed layer of the reservoir (Nandini, Merino-Ibarra & Sarma, 2008) due to intense daytime wind (Merino-Ibarra et al., 2008). The bottles were incubated for four to six light-hours at depths of 0, 1, 2, 4, 8, 12, 20, 24 m, and as close to the bottom as possible, depending on the depth of the reservoir at the moment of sampling. At each depth, nine oxygen bottles were filled (three for initial oxygen determination, three for light incubation and three for dark incubation). During the sampling, extreme precautions were taken to completely avoid bubbling that could alter the oxygen, as recommended by Valdespino-Castillo et al. (2014). Dissolved oxygen concentration in each bottle was determined in the laboratory, in triplicate for each sample bottle to minimize and assess error (cf. Valdespino-Castillo et al., 2014 for further details on the method).

Metabolism calculations

Gross primary production (GPP), net primary production (NPP), and community respiration (R) were calculated using the oxygen change rate in the light and dark bottles, respectively, following Wetzel & Likens (1991), therefore dividing the differences between initial and final oxygen concentrations by the specific incubation time of each set of bottles. Following Valeriano-Riveros et al. (2014), the depth of the production layer (Z_pl; assumed to end where GPP = 0) was calculated using Secchi depth (SD). To do this, we analyzed the correlation between SD and the GPP = 0 depth in those production profiles in which this point was found in one of the sampled depths. The data fitted better when stratification and circulation data were separated, therefore two correlations (stratification and circulation) were obtained: a coefficient of 4.149 (n = 26, R² = 0.85) for stratification, and 3.295 (n = 29, R² = 0.90) for circulation. The corresponding coefficient was multiplied by SD to obtain Z_pl in those samplings in which it could not be easily identified in the vertical profile, because GPP = 0 did not concur with one of the sampled depths.

To obtain area-based rates, volumetric rates were integrated over the production layer: the rate from each of the depths sampled was multiplied by the height of the water layer it represented. The same procedure was used when integrating R rates below the production layer (R_bpl), which were integrated from Z_pl to the bottom of the reservoir during circulation, or to the depth where total anoxia was found during stratification. Both respirations (R_pl, respiration in the production layer) were added to obtain the total respiration (R_Total = R_pl + R_bpl) for the full water column.

To calculate diel rates, the hourly production rates were multiplied by the corresponding photoperiod for each sampling date at this latitude. In the case of diel respiration, the respiration during the night in the production layer was calculated by multiplying the nighttime hours by the dark respiration rate, conservatively estimated as 10% of the GPP rate following Geider & Osborne (1989). Below the production layer, diel respiration was assumed to be constant and was therefore estimated by multiplying R_bpl by 24 h. The proportion GPP: R_Total in carbon units was used to assess the metabolic balance of the ecosystem as a whole.

Conversion of oxygen rates to carbon rates used the theoretical and most widely used conversion values PQ = 1.3 and RQ = 1.0 (Gazeau et al., 2005). The f-ratio (NPP/GPP), sensu Falkowski et al. (2003), was used to assess the fraction of the production that could potentially be exported from the production layer through the sinking of biomass).

Data analyses

Nutrient data from Sacristán-Ramírez (2017) and Barjau-Aguilar (2018) were used for preliminary assessment of the relations between nutrient availability and metabolic fluxes. To assess the effect of water-level fluctuations, we calculated the relative lake-level fluctuation index (RLLF) as proposed by Kolding & Van Zwieten (2012, Eq. (1)). The RLLF is a simple empirical indicator defined as the mean amplitude of the annual or seasonal fluctuations in lake level divided by the mean depth of the lake or reservoir, times 100. We calculated the RLLF both by year and by period (stratification and circulation) to identify the time scale in which water-level fluctuations exert the greatest influence on metabolism of the reservoir. (1)${\displaystyle \text{RLLF}=\text{mean lake level amplitude/mean depth}\ast 100.}$ Contour diagrams were constructed in Surfer 11.0.642 (Golden Software, Inc.). The Kriging algorithm with the default linear variogram was used. For the temperature and oxygen diagrams (Figs. 2 and 3) a grid of 2,632 ×  108 was used in order to obtain point spacing of approximately one month on the horizontal axis and 1.1 m on the vertical axis. Vertical variation contours of the metabolic rates (Fig. 4) were obtained from monthly averages from the surface to 4 m for an average year, using a search ellipse of four months and 2 m.

Estimation of error propagation

Errors and their propagation were determined according to classical error propagation theory (Ku, 1966). To obtain confidence intervals (CI) around our calculated variables, the standard error (SE) of the mean was converted directly to CI by multiplying with the t value for the desired alpha and specific degrees of freedom (df) for each case (Lehrter & Cebrian, 2010).

Since our calculations required only sums, products and quotients, we used the simplified equations proposed by those authors for simple step-by-step calculations, where SE_W is the SE of a function W calculated from the means of two or more variables (i.e., $W=f\left(\overline{U},\overline{V}\right)$). For the sum we used Eqs. (2) and (3) for products and quotients: (2)${\displaystyle {\mathrm{SE}}_W\sqrt{{\left(a\mathrm{SE}\overline{U}\right)}^2+{\left(b{\mathrm{SE}}_{\overline{V}}\right)}^2}}$ (3)${\displaystyle {\mathrm{SE}}_W\sqrt{\left(\frac{{\mathrm{SE}}_{\overline{U}}}{\overline{U}}\right)+{\left(\frac{{\mathrm{SE}}_{\overline{V}}}{\overline{V}}\right)}^2}.}$ These equations assume that the means involved in the calculation follow a normal distribution, which was checked with the data. Propagation of the degrees of freedom (df) from the metabolic rates used Eq. (4) adapted by Lehrter & Cebrian (2010) from the Welch-Satterthwaite formula (Ku, 1966). (4)${\displaystyle d{f}_W\frac{{\mathrm{SE}}_W}{\sum _{i=1}^J\frac{c_i^4{\mathrm{SE}}_{\frac{4}{W_i}}}{d{f}_{{\overline{W}}_i}}}.}$

Environmental parameters

Temperature in VB ranged during 2006–2015 between 17.7 and 23.9 °C. Thermal stratification occurred every year, with averages of 22.3 °C during stratification and 19.9 °C in the circulation (Fig. 2). Mean epilimnetic temperatures were always above 22.5 °C during the summer (Jun–Sep, well-established stratification). The lowest temperatures (17.3–19.0 °C) were found in January and February, during the circulation period. These minimal temperatures varied considerably among the years sampled, reaching a minimum during the circulation of 2010–2011 (mean 18.77 °C), and a maximum in the 2013–2014 circulation period (mean 20.07 °C) (Table 1).

The vertical distribution of oxygen found in VB during 2006–2015 (Fig. 3) sharply outlined the monomictical behavior of VB. During the stratification, dissolved oxygen concentrations were close to saturation (∼8 mg L⁻¹) from the surface down to ∼6 m and thence declined abruptly down to ∼10 m depth, below which an anoxic hypolimnion was found. In contrast, during circulation, oxygen concentrations were nearly homogeneous throughout the water column but remained under saturated (50–80 %). Because of this sharp contrast between stratification and circulation, we used the oxycline -and in particular the 1.0 mg L⁻¹ oxygen isoline- as a proxy for the depth of the functional mixing layer (Z_mix) in VB (Fig. 3), following Catalan & Rondón (2016). Z_mix ranged between 8 and 13 m during the stratification periods and deepened to the bottom of the reservoir during most of the circulation periods. Secchi disk transparency varied generally from ∼1 m during stratification to ∼2 m during circulation (Table 1), although it extraordinarily reached 6.3 m (December 2010). The depth of the production layer (Z_pl), also shown in Fig. 3, averaged 5.9 m during stratification and 6.8 m during circulation. Chl-a concentrations averaged 11.6 µg L⁻¹ overall in VB during the period sampled. The mean concentration during the stratification periods (12.2 µg L⁻¹) was higher than during the circulation periods sampled (10.7 µg L⁻¹, Table 1), but the difference was not statistically significant because variability among periods was high.

Vertical variations of the metabolic rates (gO₂ m⁻³ h⁻¹)

As expected, metabolic rates exhibited strong vertical gradients. Figure 4 summarizes these gradients and their annual variations by plotting the mean monthly rates at each depth to obtain the average annual pattern for 2006–2015. Average production (both net and gross) rates were maximal around the 1 m depth (0.09–0.14 gO₂ m⁻³ h⁻¹ for NPP, and 0.15–0.20 gO₂ m⁻³ h⁻¹ for GPP) throughout the year, and diminished rapidly and relatively regularly below this depth, clearly in relation to the decrease in light availability. The maximum was higher during the stratification months, from April to July, and showed a secondary maximum around the onset of the circulation period, between September and November. The maximum was smaller and less evident during the central months of the circulation period, December and January, when the vertical variation of production rates was also smaller. The highest variability of production through time was observed in the surface (Fig. 4), where averages ranged from <0.06 to >0.14 gO₂ m⁻³ h⁻¹ for NPP, and from <0.11 to >0.19 gO₂ m⁻³ h⁻¹ for GPP. Respiration also had a vertical gradient, and was maximal (around −0.06 gO₂ m⁻³ h⁻¹, Fig. 4B) at the surface and down to 1 m, presumably as a result of photic stress and photorespiration.

Vertically integrated metabolic rates (gO₂ m⁻² h⁻¹)

Once vertically integrated for the whole production layer, metabolic rates could be graphically depicted for the full period sampled to show the temporal heterogeneity of production and consumption of oxygen (Fig. 5A). The depth of the production layer (Z_pl), to which the rates were integrated, ranged from 3 to 21 m (Fig. 3). Hourly GPP rates ranged from 0.15 to 1.26 gO₂ m⁻² h⁻¹, oxygen consumption from −0.13 to −0.83 gO₂ m⁻² h⁻¹ and NPP from −0.36 to 0.66 gO ₂m⁻² h⁻¹ (Fig. 5A). The three rates showed variability between samplings, but a relative long-term stability.

Oxygen consumption below the production layer (R_bpl) ranged from −0.23 to −1.38 gO₂ m⁻² h⁻¹. R_bpl rates during the stratification periods were similar to the production layer values, and the highest rates occurred during the circulation periods (Fig. 5B), when they were integrated from Z_pl to the bottom of the reservoir. When the respiration of both layers was considered, the resulting total respiration (R_Total) for the water column of VB ranged from −0.28 to −2.11 gO₂ m⁻² h⁻¹, and the highest rates also occurred during the circulation periods.

Metabolic balance and carbon flux

To assess the metabolic balance and carbon flows, daily rates (24 h-integrated) were considered (Table 2). Although during the daytime the production layer was most of the time autotrophic, when the nighttime respiration was also considered the production layer was sometimes heterotrophic on a daily basis, particularly during the circulation samplings. Moreover, when the respiration from the layer below the production layer (R_bpl) was also considered, VB as a whole was net heterotrophic throughout the studied period. In fact, the GPP: R_Total ratio (in C units) of the reservoir ranged from 0.07 to 0.93 and averaged 0.52 for the full 2006–2015 period (Fig. 6). Although it indicated heterotrophy throughout the full period, the GPP: R_Total ratio also showed a very slight long-term trend with an increase of 0.02 y⁻¹ (Mann–Kendall trend test, p < 0.01).

In terms of carbon fluxes, GPP showed VB to be a very productive system, with an average C fixation of 3.60 gC m⁻² d⁻¹ during the 10 years sampled (Table 2). The NPP averaged 1.64 gC m⁻² d⁻¹, so the potential carbon exportation from the production layer through biomass sinking (f = NPP∕GPP) of VB averaged 45% (Table 3) and was higher during the stratification periods (48%) than during the circulation periods (40%). On the annual scale, our measurements of NPP ranged from 258 to 892 g C m⁻² y⁻¹ and averaged 588 g C m⁻² y⁻¹. When compared to the more recent carbon burial rates reported for VB (250 g C m⁻² y⁻¹; Carnero-Bravo et al., 2014), the estimated recycling efficiency in the water column of VB would average 58% of the annual NPP and would have ranged from 3% to 72% in the years 2006 to 2015. This means that on the average, only 19% of the total annual 1,344 g C m⁻² y⁻¹ fixed in VB would be sequestered in the sediments, and most of the potential carbon burial would be remineralized back to the water column. Overall, the system’s total respiration rate was higher than GPP or NPP, averaging −6.89 gC m⁻² d⁻¹, which in steady state would imply a potential net carbon release of −3.29 gC m⁻² d⁻¹ to the atmosphere on average during 2006–2015.

Correlations with environmental drivers

Overall, although highly variable among samplings, metabolic rates showed significant differences between the stratification and circulation periods. Total respiration (R_Total) rates were significantly (p < 0.01) higher (1.11 gO₂ m⁻² h⁻¹) during circulation than during stratification (0.73 gO₂ m⁻² h⁻¹). R_bpl rates were also significantly (p < 0.01) higher (0.77 gO₂ m⁻² h⁻¹) during circulation than during stratification (0.42 gO₂ m⁻² h⁻¹), while R_pl rates showed a similar pattern, but the differences were not statistically significant (p < 0.33).

In contrast, production rates showed the opposite pattern. NPP rates were higher (0.31 gO₂ m⁻² h⁻¹) during the stratification than during the circulation (0.26 gO₂ m⁻² h⁻¹) and this was also a statistically significant difference (p < 0.04). GPP rates were also slightly higher (0.63 gO₂ m⁻² h⁻¹) during the stratification than during the circulation (0.60 gO₂ m⁻² h⁻¹), but this difference was not statistically significant (p < 0.20)

Significant correlations among the rates and the environmental drivers considered were also obtained for each of the two limnological conditions of the reservoir. GPP and NPP did not correlate significantly with temperature, but respiration did correlate positively for both periods. For the stratification periods, R_bpl showed a high and significant positive correlation (r² = 0.452, p < 0.05) with temperature. For the circulation periods, when respiration occurs throughout the water column (Fig. 7), R_Total also correlated positively (r² = 0.314, p < 0.12) with the mean water column temperature, which varied between years depending on each winter’s conditions, as described earlier.

For the stratification periods, significant correlations were also found with the water level index RLLF (higher values of RLLF mean wider water-level fluctuation and lower levels below reservoir capacity). In particular, Secchi depth showed a significant inverse correlation with RLLF (r² = 0.569, p < 0.04). GPP also showed an inverse correlation (r² = 0.393, p < 0.07), although the correlations were much less significant during circulation or at the annual scale. Metabolic rates showed significant correlations with nutrients only during the stratification periods. GPP correlated inversely with total P (r² = 0.078, p < 0.02) and with total N (r² = 0.124, p < 0.004). NPP only correlated with total N (r² = 0.129, p < 0.003). R_bpl did not correlate with either P or N.

Limnological processes

The vertical distribution of temperature in VB found during 2006–2015 confirms the continuity of a monomictic behavior of the reservoir, as found before 2006 (Merino-Ibarra et al., 2008; Ramírez-Zierold et al., 2010) in spite of the drastic water-level decreases that occurred during this period. Functional thermal stratification occurred every year, although the thermocline gradient was small (≤ 0.4 °C m⁻¹), as expected for tropical lakes (Catalan & Rondón, 2016). Hypolimnetic temperatures increased in VB during each stratification period during 2006–2015, confirming the persistence of this trend in VB during the stratification, as reported for previous years (Merino-Ibarra et al., 2008; Ramírez-Zierold et al., 2010). Maximum temperatures in the hypolimnion occurred during the stratification periods when the water level of the reservoir was lowest (i.e., in 2008, 2009 and 2013, Fig. 2), supporting the suggestion (Ramírez-Zierold et al., 2015) that the rate of hypolimnetic warming is related with the magnitude of the water level decrease. Additionally, the 2006–2015 time series here described also shows the interannual variability of temperature in VB, particularly during the circulation periods, when the full water column is affected by heat exchange through the surface, and the water temperature is likely to be directly related to the coldness of each winter.

Oxygen vertical distribution was consistent with the temperature distribution and the monomictic behavior confirmed for VB during the decade sampled. As suggested for tropical eutrophic systems (Catalan & Rondón, 2016), we find that oxygen is a better descriptor of the limnological cycle of VB than is temperature. The mainly anoxic hypolimnions (Fig. 3) offer a clear representation of the system’s monomictic behavior and its sustained eutrophic condition. Therefore, we concur with Catalan & Rondón (2016) that this possibility should be considered for other tropical eutrophic systems, where vertical variations of oxygen are sharper than vertical variations of temperature.

Secchi depth was also a simple but very useful tool to extend the assessment of the vertical range of production (Z_pl) throughout the full sampled period (Fig. 3). In particular, correlations between SD and GPP = 0 allowed the depth where GPP equals zero to be modeled; this was frequently a depth where experimental incubations were not performed, in spite of the effort we applied to resolve vertical variations through incubating at up to 8 different depths. Hence, the use of Secchi depth to estimate Z_pl can be a useful tool to manage large datasets containing vertical assessments of GPP and R, and to improve their yield of integrated rates. It is notable that among the metabolic studies in tropical lakes or reservoirs summarized in Table 4, fewer than half report on the measurement of metabolism at multiple water depths, and only one in four reports on the determination of the depth of the production layer. This demonstrates the need for studies that deal in detail with vertical variations of the metabolic rates and therefore offer reliable integrated rates.

Effects of temperature and water-level fluctuations on metabolic rates

Although temperature is the ultimate driver of metabolic rates at organism level, at the ecosystem scale its effect can be hindered by the many other changes that can occur simultaneously. In particular, (unlike marine systems, cf. Downing, 2014), monomictic inland systems exhibit contrasting physicochemical conditions coupled to the seasonal variations of temperature, which are mainly dependent on the variable intensity and spatial extent of mixing processes. These conditions include: the formation of an anoxic hypolimnion where aerobic respiration cannot occur; the transport of plankton across light boundaries during full circulation; and the exposure of reduced compounds to oxidizing conditions and/or the transport of nutrients associated with boundary mixing and hypolimnetic entrainment (Ramírez-Zierold et al., 2010; Ramírez-Zierold et al., 2015).

All of these processes and conditions can affect production and respiration rates, hindering the effect of concomitant temperature changes. This apparently was the case in the initial assessment of Valdespino-Castillo et al. (2014), who did not find a significant positive correlation between respiration and temperature in tests of 2006–2007 stratification data, or even of annual data for VB. Valdespino-Castillo et al. (2014) attributed this to the dominance of the effect of mixing on respiration over that of temperature when stratification and circulation data were assessed together. In contrast, the long-term data now reported reveal a positive correlation of respiration with temperature when circulation and stratification periods were analyzed separately. In the case of the circulations, having data from a series of years allowed enough thermal variability to reveal the expected direct effect of temperature on metabolic rates, and particularly on respiration. Additionally, this expected correlation could also be identified in the case of the stratifications within the lower part of the epilimnion (Fig. 7), where the effects of daily wind mixing would be smaller than in the rest of the epilimnion.

Because of the wide water-level fluctuations in VB during 2006–2015, this data set allowed assessment of the effects of these variations, parameterized here as RLLF. Water level decreases may increase the frequency and intensity of boundary mixing events and hypolimnetic entrainment during stratification in VB (Ramírez-Zierold et al., 2010; Ramírez-Zierold et al., 2015). In turn, it has also been shown by Valeriano-Riveros et al. (2014) that increased boundary mixing favors diatoms over noxious cyanobacteria (Nostocales), owing both to mixing itself and to the nutrient inputs involved. This would mean that water-level fluctuations could affect the food web through changes in mixing. In fact, Jiménez-Contreras (2009) also found significant variation in the composition of the zooplankton community, observing that the dominance of rotifers over cladocerans inverted during low-level periods, changing from a short microbial loop to a longer food web. Hence, the impact of water-level fluctuations on metabolic rates likely involves its effect over multiple processes, including mixing itself, changes in the planktonic food web and an increase in nutrient availability to the surface layer.

In our 2006–2015 metabolic data set, where ten stratifications with different water levels could now be compared, both Secchi depth and GPP decreased significantly as a function of the RLLF. The decrease of Secchi depth could be due to the increase in nutrient supply to the surface layer through hypolimnetic entrainment and boundary mixing, where nutrients might be limiting the expansion of phytoplankton. That there may be a certain degree of nitrogen limitation during the stratification in VB was suggested by Valeriano-Riveros et al. (2014) and this view is also supported by the low but significant positive correlation we obtained between NPP and TN.

Another process that could be causing decreased Secchi depth and GPP during periods of lower water-levels could be resuspension of fine sediment along the new shorelines; sediments that had settled under conditions of higher water levels would be increasingly exposed as water levels fell. The energy needed for this is available at VB, where strong winds blow daily, but the relative abundance of the sediments suspended in the water column under different water-level conditions still needs to be measured to verify the importance of this possibility.

In either case, these results are consistent with those previously found by Valdespino-Castillo et al. (2014) who concluded that sharp water level decreases may shift the community metabolism from autotrophy towards heterotrophy. Altogether, this pattern is also important to direct the needed reassessment of the contribution of reservoirs and lakes to the global carbon cycle (Cole et al., 2007; Tranvik et al., 2009; Lewis, 2011; Raymond et al., 2013; Downing, 2014), because it shows that lowering of the level would cause a decrease in the net carbon sequestration of deep stratified systems.

Metabolic balance, spatial and temporal long-term trends

The high GPP found in VB throughout the decade confirms that it has remained a very productive system, in spite of the wide water-level fluctuations and other changes. In fact, its C fixation rates are higher that those found during the initial assessment of 2006–2007 and are now within the range (2.7–5.0 gC m⁻² d⁻¹) of other tropical epicontinental water bodies considered hypertrophic (e.g., Oloiden in Kenya, Lake Conway in Florida, McIlwaine in Zimbabwe, Pao Caniche in Venezuela, and Shahidullah Hall in Bangladesh, Table 4), but not yet as high as, for example, the highly hypertrophic Lake Xolotlán in Nicaragua. Hence, our results in terms of primary production indicate that VB would now be more accurately classified as hypertrophic than as eutrophic.

Although productivity in VB is high, and nearly half of it is net production -which can be exported from the production layer to the sediments and become sequestered C- our respiration data and the comparison with net C burial assessed through independent radiometric methods of Carnero-Bravo et al. (2014) indicate that most of this carbon is remineralized and only about 40% of the NPP is permanently buried in the sediments of this ecosystem.

Furthermore, because the total respiration in this ecosystem is also very high -nearly doubling GPP rates- the system has a net heterotrophic metabolism. These results are consistent with the findings of Gupta et al. (2008), Almeida et al. (2016) and Räsänen et al. (2018) and confirm—now with a long-term data set—that tropical hypertrophic reservoirs can be highly productive and simultaneously be important sources of atmospheric emissions; these emissions may be larger than the burial of organic carbon in their sediments, because of high rates of mineralization in the water column and sediments, as occurs in VB. These findings indicate that the pattern found for temperate systems by Duarte & Prairie (2005) and Hoellein, Bruesewitz & Richardson (2013) may also take place in eutrophic reservoirs. More metabolic studies on tropical systems that include detailed measurements of vertical variation and respiration are needed, to reassess the contribution of epicontinental water bodies to global carbon balance (Cole et al., 2007; Alin & Johnson, 2007; Tranvik et al., 2009; Lewis, 2011; Raymond et al., 2013; Downing, 2014).

It is expected that long-term data will be key to assessing metabolic variability (Staehr et al., 2010; Sarmento, 2012; Solomon et al., 2013) in the possible scenarios of climate change (Kosten et al., 2010). Our results show that a decade of metabolic records can be enough to start identifying trends, as recently found by Agusti et al. (2017). Furthermore, the correlations found here with environmental drivers allow the exercise of simple predictions for VB that can inform global expectations. For instance, our results indicate that for each 1 °C increase in the lake temperature its respiration and C emission could further increase by 0.4–0.9 gC m⁻²d⁻¹; this is consistent with the findings of Kosten et al. (2010) on CO₂ emission from lakes, and suggests one of the negative feedbacks between climate change and eutrophication that could be expected in heterotrophic tropical systems. Similarly, our results on the inverse relationship between GPP and water level decrease (RLLF) also demonstrate that the water deficiency expected for the latitude of VB—if coupled with the intensification of its use as a source of fresh water—could enhance its role as a carbon source, a process that may also occur in the numerous water reservoirs that will be more intensively used throughout similar latitudes.