Carbon Dioxide Concentration and Emissions along a Trophic Gradient in Tropical Karst Lakes

Abstract

Inland aquatic ecosystems are valuable sentinels of anthropic-associated changes (e.g., agriculture and tourism). Eutrophication has become of primary importance in altering aquatic ecosystem functioning. Quantifying the CO₂ emissions by inland aquatic ecosystems of different trophic statuses may provide helpful information about the role of eutrophication on greenhouse gas emissions. This study investigated diel and seasonal carbon dioxide (CO₂) concentrations and emissions in three tropical karst lakes with different trophic statuses. We measured CO₂ emissions using static floating chambers twice daily during the rainy/warm and dry/cold seasons while the lakes were thermally stratified and mixed, respectively. The CO₂ concentration was estimated by gas chromatography and photoacoustic spectroscopy. The results showed a significant seasonal variation in the dissolved CO₂ concentration (C_CO2) and the CO₂ flux (F_CO2), with the largest values in the rainy/warm season but not along the diel cycle. The C_CO2 values ranged from 13.3 to 168.6 µmol L⁻¹ averaging 41.9 ± 35.3 µmol L⁻¹ over the rainy/warm season and from 12.9 to 38.0 µmol L⁻¹ with an average of 21.0 ± 7.2 µmol L⁻¹ over the dry/cold season. The F_CO2 values ranged from 0.2 to 12.1 g CO₂ m⁻² d⁻¹ averaging 4.9 ± 4.0 g CO₂ m⁻² d⁻¹ over the rainy/warm season and from 0.1 to 1.7 g CO₂ m⁻² d⁻¹ with an average of 0.8 ± 0.5 g CO₂ m⁻² d⁻¹ over the dry/cold season. During the rainy/warm season the emission was higher in the eutrophic lake San Lorenzo (9.1 ± 1.2 g CO₂ m⁻² d⁻¹), and during the dry/cold the highest emission was recorded in the mesotrophic lake San José (1.42 ± 0.2 g CO₂ m⁻² d⁻¹). Our results indicated that eutrophication in tropical karst lakes increased CO₂ evasion rates to the atmosphere mainly due to the persistence of anoxia in most of the lake’s water column, which maintained high rates of anaerobic respiration coupled with the anaerobic oxidation of methane. Contrarily, groundwater inflows that provide rich-dissolved inorganic carbon waters sustain emissions in meso and oligotrophic karstic tropical lakes.

1. Introduction

Inland aquatic ecosystems (IAEs) play a relevant role in the global C cycle and, therefore, in the climate change modulation, although they only represent 0.02% of the water on Earth [1,2,3,4,5]. The IAEs are active sites that transport, transform, exchange, recycle, and sequester C from the atmospheric exchange, catchment, or groundwater inflow [6,7]. IAEs release significant amounts of greenhouse gases (GHGs) into the atmosphere, mainly CO₂ and CH₄, although there is some disagreement about the magnitude of these emissions [8,9]. CO₂ and CH₄ outgassing from IAEs account for around 76% of the C that terrestrial landscapes may deliver to them (3.9 Pg Cyr⁻¹; [10]), although considerable uncertainty remains about this global estimate. Although DelSontro et al. [8] estimated that GHG emissions from lakes and impoundments are equivalent to ~20% of global fossil fuel CO₂ emissions, Pollard [11] recently quantified that current emissions from freshwaters occur at an annual rate six times that of fossil fuel burning. Moreover, GHG emission data from tropical and subtropical freshwaters remain very sparse, and global assessments also experience this geographic bias [12]. Nonetheless, there is quite a lot of unanimity that GHG emissions will intensify even further with the continued eutrophication of lentic ecosystems [8,9,13], but its magnitude is not yet clear.

IAEs are sensitive to accelerated changes caused by anthropogenic activities (mainly agriculture and urban development; [14,15]). C processing in lakes is closely related to trophic status [10]. Several environmental parameters linked to eutrophication, such as chlorophyll-a (Chl-a) and nutrients concentrations, are well-known drivers of GHG fluxes in freshwaters [8,9,16]. However, the effects of increasing nutrient loading into lakes on the CO₂ emission to the atmosphere still need to be better understood. While some studies (e.g., [17,18,19]) have shown that the organic matter (OM) mineralization and CO₂ outgassing are stimulated in lakes receiving agricultural nutrient-rich effluents, other (e.g., [20,21,22]) found that the increased aquatic primary productivity promoted by high nutrient loading reduces CO₂ emissions. The shift from one state to another is still being determined. Nonetheless, several mechanisms involved in the biological processes that produce and/or consume CO₂, including those anaerobic processes such as methane oxidation [23], and abiotic such as carbonate weathering [24], could drive these observed patterns.

Karst lakes cover about 20% of the planet’s ice-free land, the primary drinking water source for hundreds of millions worldwide [25]. Groundwater is key to karst lakes’ hydrological and ecological functioning [26]. Groundwater transports large amounts of dissolved inorganic carbon (DIC = HCO₃⁻ and CO₂) in the karst landscape, which is added by: (i) the biological activity through respiration and decomposition of OM in the soil [27] and (ii) the geochemical interactions (rock weathering/dissolution) with carbonate minerals [28]. The amount of CO₂ released in karst lakes is the net budget of the C inflow through groundwater which depends on carbonate weathering, the lake metabolism, the OM fermentation, the CH₄ oxidation (aerobic and anaerobic), and the precipitation and dissolution kinetics of dissolved carbonates.

There are few studies of GHG emissions from karst lakes, but notably fewer in tropical regions [28,29,30,31]. Very little is known about which variables control GHG emissions in karst lakes and how these lakes can give feedback or modulate climate change. Evidence suggests that very alkaline lakes are more supersaturated in CO₂ than other lakes [24], but does this imply that these lakes emit more CO₂? Or does biogenic carbonate precipitation and stimulated primary production significantly reduce CO₂ outgassing to the atmosphere compared to other lakes? It is necessary to identify the main mechanisms underlying the release of CO₂ in these lakes and how they respond when productivity is stimulated by the inputs of nutrient-rich waters and then when eutrophication changes metabolic processes.

Karst landscapes cover more than 15% of Mexican territory and host some of the most important lakes in the country [32]. The “Lagunas de Montebello” National Park (LMNP), located in Chiapas State, extends over one of Mexico’s main karst lake districts. The LMNP has more than 130 karst water bodies of different morphometric and physicochemical characteristics [33,34,35]. Since these lakes can be found throughout gradients related to groundwater and surface water inputs, contaminated effluents, surrounded by agriculture or pristine forests, provide a unique opportunity to assess CO₂ emissions mechanisms together with eutrophication effects, considering a wide variety of lakes.

In this work, our main aim was to determine the principal drivers of the dissolved CO₂ in the water column and the CO₂ emissions into the atmosphere of tropical karst lakes, assessing the potential effects of eutrophication through a seasonal comparison. For this purpose, we conducted a study in three lakes at the LMNP, which displayed contrasting trophic statuses sampling during two consecutive seasons, rainy/warm and dry/cold. Our central hypotheses were (a) if the NPP driven by high nutrient loading and mirrored by high Chl-a concentrations leads to CO₂ concentrations below the atmospheric equilibrium, then a direct negative relationship between eutrophication proxies and CO₂ emission is to be expected, (b) if the increase in temperature promotes metabolic activity, then it is expected that during the rainy/warm season a more significant amount of CO₂ is emitted compared to the dry/cold season due to the high productivity of the lakes and the lower DO concentration and, (c) if DIC-rich groundwater inflow contributes to supersaturation and subsequent CO₂ emission, then tropical karst lakes are expected to be significant emitters of CO₂ to the atmosphere.

Therefore, to answer the questions mentioned above, after the main results on physical, chemical, and biological variables are shown along with the differences observed due to trophic statuses and seasons, the main drivers of the observed changes in dissolved CO₂ concentrations and emissions are elucidated using multiple linear regressions. CO₂ dynamics in these tropical karst lakes depend on a complex interplay of biotic and abiotic variables exposed to strong tropical seasonality which are deeply discussed.

2. Materials and Methods

2.1. Study Area

The LMNP (16°04′–16°10′ N, 91°37′–91°47′ W, 1500 m a.s.l.) is in southeastern Chiapas, Mexico, in the Río Grande de Comitán basin, which is part of the Hydrological-Administrative Region No. 3056 Grijalva-Usumacinta and has an area of 810 km² (Figure 1; [36]). The LMNP extends over 6425 ha between La Independencia and La Trinitaria municipalities [37]. The “Lagunas de Montebello” was a designated Natural Protected Area, Biosphere Reserve, and Ramsar site in 1959, 2011, and 2003 respectively [38]. The climate is tropical rainy, with a relatively cold and dry season (dry/cold from January to May) and a typical summer rainfall regime (rainy/warm from June to November; Cb(m)(f)ig [39]) with December as a transitional month. The average annual temperature and precipitation are 18.7 °C and 1960 mm, respectively (CONANP Automatic Meteorological Station N15DA7496, 16°06′52.5″ N, 91°43′48.2″ W). The landscape in the basin is karstic, developed under Cretaceous stratigraphic units, mainly on limestone–dolomite and limestone–shale rocks [32]. In LMNP, leptosols are predominant, with coniferous forests occupying 73% of the total area and farming (irrigated and rainfed) 9%. Water bodies extend to 16% of the total area. In the last 20 years, agriculture has increased by 24% in the NW zone [36].

Although lakes receive surface inflows from the Río Grande de Comitán, which runs through urban and agricultural areas, the primary water source is groundwater, which has generated a complex karst landscape crossed by a system of dolines, uvalas, and poljes created by the dissolution of carbonate rocks [32]. The study was carried out on three selected lakes that typify three different trophic statuses (based on previous records of concentration of Chl-a, TSS, total phosphorous, and total nitrogen in the water column) that can be founded in the lake district: oligotrophic Tziscao (TZ), mesotrophic San José (SJ), and San Lorenzo (SL). The main characteristics of the three lakes are described in

Table 1. Location and main morphological characteristics of the studied lakes [33]. Z_MAX = maximum depth, Z_M = mean depth.

Lake	Lat	Long	Alt	Volume	Area	ZMAX	ZM
	(N)	(W)	(m a.s.l.)	(km3)	(ha)	(m)	(m)
TZ	16.075	91.665	1490	0.09	306.6	86	28.9
SJ	16.106	91.738	1454	0.006	60.6	30	10.3
SL	16.126	91.753	1455	0.02	181.3	67	11.2

. The lakes are warm monomictic, circulating during the dry/cold winter, and remain stratified for the rest of the year along the rainy/warm seasons.

2.2. Physico-Chemical and Biological Characterization

Samplings were carried out in the two contrasting seasons, typical of tropical America: (1) the rainy/warm (2021) and (2) the dry/cold season (2022). All measurements were performed twice daily, in the morning (~10:00 h) and the evening (~18:00 h). Sampling consisted of in situ profiles of temperature (T, °C), dissolved oxygen (DO, mg L⁻¹), electrical conductivity (K₂₅, µS cm⁻¹), and pH using a Hydrolab DS5 multiparameter water quality sonde (at 1 m of depth resolution). We calculated the thermo- and oxyclines’ location, width, and gradient (∇). The mixing layer (Z_MIX) was estimated according to the vertical profiles of T and DO. The water transparency was recorded by Secchi disk depth (Z_SD).

Sampling was performed at three depths (epilimnion (Z_MIX), metalimnion (planar thermocline), and hypolimnion (close to the bottom)) when stratified, while a half meter below the surface, mid-water, and one meter above the lake bottom when circulating. Water samples for analyses were taken with a UWITEC bottle (5 L). Water samples for nutrient analyses were filtered in situ through cellulose acetate filters (0.22 μm pore size), placed in polyethylene bottles, and stored in a dark and frozen place (−4 °C) until analysis. The phosphorous as soluble reactive phosphorous (P-SRP), nitrogen as ammonia (N-NH₄), nitrites (N-NO₂), and nitrates (N-NO₃) were performed in a segmented flow autoanalyzer Skalar Sanplus System following standard methods [40,41,42]. Total suspended solids (TSS, seston) determination followed the gravimetric method [43,44,45].

Chlorophyll-a (Chl-a) concentration was determined, at the same three depths, following the EPA method 445.0 [46] using a Turner Designs 10-AU fluorometer. Net primary production, NPP, rates at the surface (0.5 m) were measured following the light and dark bottles method through the changes in DO concentration [44,47]; which was measured with a portable HACH model HQ40d oximeter with an LBOD10101 luminescent oxygen probe.

Carbonate precipitation was assessed through the saturation index for calcite (SIcalcite [48]) which is an approximate indicator of the degree of saturation of calcium carbonate in water. SIcalcite was calculated for each lake separately based on T, K₂₅, and pH values (this study) and HCO₃^– and Ca⁺ concentrations ([49]; Alcocer, 2022 personal communication).

2.3. DIC, Dissolved CO_2, and Flux Measurements

Total dissolved inorganic carbon (DIC) was determined after filtering 60 mL of water sample through Whatman GF/F filters which were then acidified with H₃PO₄ (85%). The DIC samples were stored in a cold and dark place for 24 h to force the equilibrium of the carbonate system. Then the CO₂ gas was extracted using the headspace technique [50] and stored in 12 mL Exetainer vials (Labco) previously evacuated until analyses. Dissolved CO₂ concentration (C_CO2) was also measured using the headspace equilibration technique, avoiding the acidification of the water sample, as described by Goldenfum [51]. Then the CO₂ was transferred into Exetainer vials (12 mL) for subsequent analysis. All samples were duplicated and analyzed by gas chromatography (Agilent model 6890N) equipped with a single-stage dual-packed column, where CO₂ was detected using a thermal conductivity detector (TCD).

Direct measurements of CO₂ total evasion fluxes (F_CO2) from the water surface were performed using static floating chambers (5 L polyethylene) similar to those described in DelSontro et al. [52]. A total of eight chambers were randomly placed in the center of each lake, which lasted for 30–45 min to monitor the changes in CO₂ concentration inside the chamber. In the rainy/warm season, CO₂ was measured at 15-min intervals after being extracted with a syringe through a butyl rubber stopper located in the upper part of the chambers. As previously cited, each extracted sample was injected into an Exetainer vial (12 mL) for subsequent analysis by gas chromatography. In the dry/cold season, the concentration inside the chamber was measured in situ in the field using a Gasera ONE PULSE based on photoacoustic spectroscopy through NDIR-PAS technology (mechanically chopped broadband IR source with optical bandpass filters). Gas sampling by the Gasera instrument was carried out automatically at 9-min intervals through 2 m Teflon tubing, recirculating the measured gas back into the chambers. Both analytical methods (chromatography and photoacoustic spectroscopy) were previously intercalibrated, and similar control values were obtained by both methods (R² = 0.987; Spectros (CO₂ ppm) = 0.99*Chromat (CO₂ ppm)).

The CO₂ flux was determined from the slope of the change in the CO₂ concentration (ppm) inside the chamber concerning the sampling time (Δt; s) by means of linear regression. The area (A_C; m²) and the volume (V_C; m³) of the chamber were also considered (Equation (1)). The concentration was adjusted for pressure and temperature according to the ideal gas law.

$${\mathrm{F}}_{\mathrm{C}\mathrm{O}2}=\frac{∆\mathrm{C}}{∆\mathrm{t}}\left(\frac{{\mathrm{V}}_{\mathrm{c}}}{{\mathrm{A}}_{\mathrm{c}}}\right)$$

(1)

2.4. Statistical Analyses

Several statistical tests were used to describe the central tendency of the measured variables in order to detect significant differences between lakes (TZ, SJ, and SL), between seasons (rainy/warm and dry/cold), and to determine the main drivers of the CO₂ concentration and emissions in these tropical karst lakes. After checking the normality of variables using the Shapiro–Wilk test [53], the non-parametric Kruskal–Wallis’ one-way analysis of variance [54] was used to comparatively describe the lakes and to assess the influence of cultural eutrophication (hypothesis one), through physical and chemical variables (T, DO, pH, K₂₅, TSS, Z_SD, and nutrients) as well biological variables (Chl-a and NPP). Seasonal effects on metabolic processes involved in CO₂ budgets and emissions (second hypothesis), were also tested by means of the Kruskal–Wallis test. Analyses were performed by the functions shapiro.test() and kruskal.test() on the “stats” package [55]. The Kruskal–Wallis analyses were performed following an H distribution, and differences at the p < 0.05 level were considered statistically significant. Subsequently, a post hoc least significant difference test was used to determine the differences between the measured groups.

Multiple stepwise linear regressions were used to assess the influence of CO₂-rich groundwater input (third hypothesis) as well as to identify which environmental variable best explained the variation in C_CO2 for each kind of lake. Multiple regression models (MRM) for each lake were run only with data of T, DO, pH, K₂₅, TSS, Chl-a, nutrients, DIC, and C_CO2 from the Z_MIX. We used the ols_step_forward_p() function from the “olsrr” package [56]. Finally, a simple linear regression analysis was applied to the data to estimate the correlations between C_CO2 and F_CO2. Autocorrelation was checked using the Durbin–Watson test [57] by using the dwtest() function on the “lmtest” package [58]. Linear functions were then adjusted through the cochrane.orcutt() function on the “orcutt” package [59].

3. Results

3.1. Physico-Chemical and Biological Characterization

There were significant differences among the three lakes with different trophic statuses for T, DO, K_25, TSS, and nutrients (p < 0.05). Due to this statistically significant variation, each lake’s environmental characterization will be described as a unit.

The oligotrophic TZ presented significant seasonal variations in the values of T, DO, K₂₅, and N-NO₂ (

Table 2. Mean values (X), standard deviation (SD), minimum (MIN), and maximum (MAX) of the physico-chemical variables of the three studied lakes (R/W: rainy/warm season, D/C: dry/cold season, TZ: Tziscao, SJ: San José, SL: San Lorenzo; T: temperature, DO: dissolved oxygen, K₂₅: electrical conductivity, Z_MIX: mixing layer, Z_SD: Secchi disk depth, TSS: total suspended solids, KW_L: Kruskal–Wallis analysis among lakes, KW_S: Kruskal–Wallis analysis among seasons, n: number of observations). * = entire water column, – = no data. Bold values indicate statistically significant differences (p < 0.05).

Variable	n	Parameter	Lake						KWL
			TZ (Oligotrophic)		SJ (Mesotrophic)		SL (Eutrophic)
			R/W	D/C	R/W	D/C	R/W	D/C
T	253	X ± SD	21.8 ± 0.2	19.2 ± 0.4	20.9 ± 0.3	21.4 ± 0.3	21.2 ± 0.1	20.8 ± 0.7	X2 = 9.6, p < 0.01
(°C)		KWS	X2 = 126.5, p < 0.01		X2 = 4.8, p < 0.05		X2 = 2.5, p = 0.11
DO	253	X ± SD	7.2 ± 0.73	5.4 ± 0.4	7.5 ± 0.3	6.0 ± 0.1	4.7 ± 1.1	6.1 ± 3.4	X2 = 7.4, p < 0.05
(mg L−1)		KWS	X2 = 111.36, p < 0.01		X2 = 34.71, p < 0.01		X2 = 4.95, p < 0.01
pH	253	X ± SD	8.5 ± 0.1	8.5 ± 0.1	9.1 ± 0.04	7.6 ± 0.04	7.7 ± 0.07	9.1 ± 0.3	X2 = 4.5, p = 0.10
		KWS	–		–		–
K25	253	X ± SD	235.4 ± 2.	241.9 ± 1.6	338.4 ± 1.7	333.2 ± 0.6	518 ± 0.1	609.7 ± 7.1	X2 = 172, p < 0.01
(µS cm−1)		KWS	X2 = 112.7, p < 0.01		X2 = 36.5, p < 0.01		X2 = 22.9, p < 0.01
ZMIX (m)	6	X	34	86 *	10	30 *	5	8
ZSD (m)	6	X	7.2	7	3.5	2.3	0.6	0.4	X2 = 4.8, p = 0.1
		KWS	–		–		–
Chl-a	84	X ± SD	0.54 ± 0.01	0.45 ± 0.01	1.0 ± 0.2	0.8 ± 0.1	46.3 ± 10.1	50.8 ± 8.5	X2 = 72.5, p < 0.01
(µg L−1)		KWS	X2 = 7.5, p < 0.01		X2 = 12.3, p < 0.01		X2 = 1.6, p = 0.2
TSS	57	X ± SD	1.2 ± 0.3	1.6 ± 1.0	4.4 ± 3.0	2.9 ± 1.3	9.4 ± 2.4	7.2 ± 1.3	X2 = 9.8, p < 0.01
(mg L−1)		KWS	X2 = 0.6, p = 0.44		X2 = 0, p = 1		X2 = 2.4, p = 0.12
P-SRP	24	X ± SD	0.1 ± 0.01	0.08 ± 0.01	0.3 ± 0.1	0.2 ± 0.03	1.0 ± 0.6	43.5 ± 15.6	X2 = 61.9, p < 0.01
(µmol L−1)		KWS	X2 = 1.9, p = 0.2		X2 = 13.7, p < 0.001		X2 = 1.8, p = 0.1
N-NO3	24	X ± SD	0.8 ± 0.13	3.4 ± 0.1	3.3 ± 1.7	5.3 ± 0.9	2.6 ± 1.5	2.0 ± 1.1	X2 = 24.3, p < 0.01
(µmol L−1)		KWS	X2 = 1.7, p = 0.1		X2 = 13.7, p < 0.001		X2 = 0.9, p = 0.4
N-NO2	24	X ± SD	0.1 ± 0.01	0.2 ± 0.03	0.4 ± 0.1	0.4 ± 0.3	0.4 ± 0.3	0.6 ± 0.4	X2 = 28.8, p < 0.01
(µmol L−1)		KWS	X2 = 4.3, p < 0.05		X2 = 0.21, p = 0.6		X2 = 5.4, p < 0.05
N-NH4	24	X ± SD	0.6 ± 0.4	1.9 ± 0.8	12.4 ± 12.2	3.7 ± 1.7	0.5 ± 0.1	0.7 ± 0.1	X2 = 34.5, p < 0.01
(µmol L−1)		KWS	X2 = 1.3, p = 0.2		X2 = 3.4, p = 0.06		X2 = 1.9, p = 0.2

). The pH, TSS concentration, N-NO₃, N-NH₄, and P-SRP did not show seasonal variations. No environmental variable showed a significant daily variation (p > 0.05). TZ showed thermal stratification during the rainy/warm season, which was also confirmed by a DO clinograde profile (Z_MIX = 34 m). During the dry/cold season, the TZ was mixed and showed homogeneous vertical profiles of T, OD, pH, and K₂₅ (Z_MIX = 86 m). The T in the Z_MIX was about 2 °C higher in the rainy/warm season than in the dry/cold one (p < 0.01). The DO concentration in the Z_MIX was higher and more variable in the rainy/warm season (7.2 ± 0.7 mg L⁻¹) than in the dry/cold season (5.4 ± 0.4 mg L⁻¹). In both seasons, the values of K₂₅, TSS, and nutrients were the lowest compared to the other two lakes under study (p < 0.05). The N-NO₃ was the most abundant inorganic nitrogen fraction in both seasons.

The mesotrophic SJ presented significant seasonal variations in the values of T, DO, K₂₅, P-PO₄, and N-NO₃ (Table 2). There were no seasonal variations in the pH, TSS concentration, N-NO₂, or N-NH₄ values. Daily variation in any environmental variable was significant (p > 0.05). A DO clinograde profile also confirmed the thermal stratification during the rainy/warm season in SJ (Z_MIX = 10 m). Over the dry/cold season, T and DO vertical profiles confirmed that the SJ lake was completely mixed (Z_MIX = 30 m). In this case, the T in the Z_MIX was about 0.5 °C higher in the dry/cold season than in the rainy/warm one (p < 0.01). In both seasons, the values of DO concentration were the highest compared to TZ and SL lakes (p < 0.05). The N-NH₄ was the most abundant inorganic nitrogen fraction during the rainy/warm season and the N-NO₃ during the dry/cold one.

The eutrophic SL presented significant seasonal variations in the values of DO, K₂₅, and N-NO₂ (Table 2). The T, pH, TSS, N-NO₃, N-NH₄, and P- SRP did not show seasonal variations. SL showed thermal stratification and anoxic conditions below 8 m in both seasons (Z_MIX = 5 m in the rainy/warm season and Z_MIX = 8 in the dry/cold season). There was also no daily variation in any measured variables (p > 0.05). The highest values of K₂₅ were recorded during the two seasons in the eutrophic SL. Significant higher concentrations of TSS and nutrients were observed in eutrophic SL (p < 0.05) than in mesotrophic SJ and oligotrophic TZ without seasonal variation (except for N-NO₂). The N-NO₃ was the most abundant inorganic nitrogen fraction in both seasons.

Chl-a concentration in the Z_MIX of the lakes ranged from 0.3 to 63.3 µg L⁻¹ with significant differences between lakes (Kruskal–Wallis: X² = 72.493, p < 0.001; Figure 2A). The eutrophic SL was the one that presented the highest average concentration (48.6 ± 3.2 µg L⁻¹) without seasonal variation (Kruskal–Wallis: X² = 1.6, p = 0.2). The Chl-a concentration in the mesotrophic SJ was about 25% higher during the rainy/warm season (1.0 ± 0.2 µg L⁻¹) compared to the dry/cold one (0.8 ± 0.1 µg L⁻¹) (Kruskal–Wallis: X² = 12.3, p < 0.01). Finally, the oligotrophic TZ presented the lowest Chl-a concentration compared to SJ and SL and showed seasonal variation (Kruskal–Wallis: X² = 7.5, p < 0.01) with the highest values during the rainy/warm season (0.54 ± 0.1 µg L⁻¹) compared to the dry/cold one (0.45 ± 0.1 µg L⁻¹). The NPP rate did not present significant differences between seasons (p = 0.5; Figure 2). However, during the warm season, the NPP rate was negative (GPP < R) in the oligotrophic TZ (−0.03 ± 3.0 mg C m⁻³ h⁻¹) and positive (GPP > R) in the mesotrophic SJ and oligotrophic SL (11.4 ± 6.8 mg C m⁻³ h⁻¹ and 33.0 ± 12.7 mg C m⁻³ h⁻¹, respectively) (p < 0.05). During the dry/cold season, in oligotrophic TZ and eutrophic SL, the GPP exceeded the R rates (9.3 ± 3.0 and 27.5 ± 10.5 mg C m⁻³ h⁻¹, respectively), and in mesotrophic SJ, the R rate was higher than the GPP one (−6.6 ± 1.6 mg C m⁻³ h⁻¹).

The SIcalcite values in the water column of the lakes were ~10.0 (TZ), ~9.9 (SJ), and ~10.4 (SL). This suggests that the lakes were supersaturated with respect to calcium carbonate (CaCO₃) and therefore there may be precipitation of it.

3.2. Dissolved CO₂ and DIC

The entire water column of the lakes was persistently supersaturated with CO₂ concerning the atmospheric equilibrium during both seasons (150–600% approximately,

Table 3. Mean values (X) and standard deviation (SD) of the DIC, dissolved CO₂ (C_CO2), as well as the ∇C_CO2 and %, Sat CO₂ in the water column of the lakes in the study. (R/W: rainy/warm season, D/C: dry/cold season, oligotrophic Tziscao (TZ), mesotrophic San José (SJ), and eutrophic San Lorenzo (SL), n: number of observations).

Lake	Season	n	Z	DIC (mmol L−1)		CCO2 (µmol L−1)		∇CCO2	%Sat CO2
			(m)	X	SD	X	SD	µmol L−1 m−1
TZ	R/W	12	0	0.84	0.05	13.9	0.1	0.5	127.4
			34	0.91	0.04	14.8	2.2		135.7
			40	0.99	0.02	42.4	6.2		388.6
	D/C	12	0	0.77	0.02	17.4	1.6	−0.1	159.3
			34	0.77	0.04	16.4	0.6		150.7
			40	0.75	0.04	14.4	0.9		132.4
SJ	R/W	12	0	1.26	0.08	35.2	0.8	−0.2	323.1
			5	1.29	0	39.9	6.6		366.5
			10	1.27	0.01	33.1	1.4		303.9
	D/C	12	0	1.07	0.11	18.8	0.7	0.1	172.2
			5	1.07	0.06	22.2	2.8		204
			10	1.1	0	19.8	1.9		181.4
SL	R/W	12	0	1.11	0.04	29.2	4.3	6.4	267.6
			5	1.06	0.08	46.3	12.9		424.8
			15	1.37	0.33	122.4	65.3		1122.8
	D/C	12	0	0.95	0.08	18.5	7.8	0.4	169.9
			5	1.01	0.05	33.6	6.1		308.3
			15	1	0.03	27.3	13.8		250.7

). C_CO2 and DIC followed a consistent pattern with significant differences between lakes (Kruskal–Wallis: X² = 27.5, p < 0.001 for the C_CO2 and X² = 40.9, p < 0.001 for the DIC). No daily variation was observed in any lake in the DIC or the C_CO2 (Kruskal–Wallis: X² = 0.01, p = 0.0 for the DIC and X² = 0.6, p = 0.41 for the C_CO2).

In the oligotrophic TZ, the DIC concentration in the Z_MIX ranged from 0.63 to 0.95 mmol L⁻¹ with an average of 0.88 ± 0.1 mmol L⁻¹ in the rainy/warm season and 0.76 ± 0.1 mmol L⁻¹ in the dry/cold one (Kruskal–Wallis: X² = 10.5, p < 0.001). The C_CO2 values ranged from 13.0 µmol L⁻¹ to 20.4 µmol L⁻¹ averaging 15.4 ± 2.3 µmol L⁻¹ and without seasonal variation (Kruskal–Wallis: X² = 2.6, p = 0.1). TZ presented C_CO2 positive gradients along the water column (∇C_CO2 = 0.5 µmol L⁻¹ m⁻¹) during the rainy/warm season.

The highest average concentration of DIC was recorded in the SJ Z_MIX (1.16 ± 0.12 mmol L⁻¹). During the rainy/warm season, the DIC concentration ranged from 1.2 to 1.3 mmol L⁻1 and was about 20% higher compared to the dry/cold season (0.95–1.15 mmol L⁻¹) (Kruskal–Wallis: X² = 13.7, p < 0.001). The C_CO2 values ranged from 32.0 to 57.9 µmol L⁻¹ (37.7 ± 8.13 µmol L⁻¹) over the rainy/warm season and 18.1–8.6 µmol L⁻¹ (20.3 ± 2.9 µmol L⁻¹) over the dry/cold one (Kruskal–Wallis: X² = 13.7, p < 0.001). SJ presented the highest C_CO2 in the metalimnion with a negative gradient during the rainy/warm season (~∇C_CO2 = −0.2 µmol L⁻¹ m⁻¹).

In the eutrophic SL Z_MIX, the DIC concentration ranged from 0.9 to 1.14 mmol L⁻¹ with an average of 1.1 ± 0.1 mmol L⁻¹ in the rainy/warm season and 1.0 ± 0.1 mmol L⁻¹ in the dry/cold one (Kruskal–Wallis: X² = 6.9, p < 0.001). The highest mean C_CO2 was recorded in SL (31.9 ± 12.3 mmol L⁻¹) without significant seasonal variation (Kruskal–Wallis: X² = 1.9, p = 0.17). SL presented a heterogeneous vertical profile, with a tendency to increase C_CO2 in the hypolimnion (∇C_CO2 = 6.4 µmol L⁻¹ m⁻¹ in the rainy/warm season and ∇C_CO2 = 0.4 µmol L⁻¹ m⁻¹ in the cold/dry season).

We identified new predictors of C_CO2 by using MRM. On the one hand, in the oligotrophic TZ, the most robust MRM for C_CO2 in the water Z_MIX included pH, N-NO₃, DIC, TSS, and K₂₅ as independent variables and explained 96% of the C_CO2 variance. The variable with the greatest weight within the TZ model was the pH (

Table 4. Multiple regression model for the C_CO2 and the environmental variables. SE = standard error, df = degrees of freedom. (R/W: rainy/warm season, D/C: dry/cold season).

TZ (Oligotrophic)
Regression Summary for the Dependent Variable: CCO2 Adjusted R²= 0.9613, F(5,18) = 89.4, p < 0.0001 Std. Error of Estimate: 2.37
	β	Std. Err. of β	p-Level	Multiple R2	R2 Change	F-to Enter	p-Level
Intercept	224.12	104.69	0.046
pH	−16.43	5.25	0.006	0.8697	0.8697	145.56	<0.001
N-NO3	3.61	0.88	0.001	0.9194	0.0497	109.39	<0.001
DIC	32.88	8.61	0.001	0.9375	0.0181	99.918	<0.001
TSS	−2.05	0.7	0.009	0.9578	0.0203	107.93	<0.001
K25	−0.42	0.34	0.22	0.9613	0.0035	89.399	<0.001
SJ (Mesotrophic)
Regression Summary for the Dependent Variable: CCO2 Adjusted R² = 0.7628 F(4,19) = 15.27, p < 0.001 Std. Error of Estimate: 5.21
	β	Std. Err. of β	p-Level	Multiple R2	R2 Change	F-to Enter	p-Level
Intercept	209.06	50.85	0.001
pH	−4.34	3.05	0.17	0.6481	0.6481	40.51	<0.001
TSS	1.65	0.62	0.015	0.6879	0.0398	23.14	<0.001
K25	−0.42	0.186	0.035	0.7315	0.0436	18.16	<0.001
N-NO3	−1.82	1.146	0.13	0.7628	0.0313	15.27	<0.001
SL (Eutrophic)
Regression Summary for the Dependent Variable: CCO2 Adjusted R² = 0.9471 F(5,18) = 59.26, p < 0.001 Std. Error of Estimate: 11.16
	β	Std. Err. of β	p-Level	Multiple R2	R2 Change	F-to Enter	p-Level
Intercept	−18.89	59.73	0.75
Chl-a	−0.61	0.68	0.39	0.6760	0.676	45.91	<0.001
TSS	3.68	0.737	0.001	0.8116	0.1356	45.23	<0.001
DO	−3.4	1.019	0.004	0.8776	0.066	47.79	<0.001
N-NO3	8.67	2.275	0.001	0.9377	0.0601	71.509	<0.001
DIC	37.98	21.82	0.1	0.9427	0.005	59.26	<0.001
N-NH4	11.9	9.98	0.25	0.9471	0.0044	50.77	<0.001

). In the mesotrophic SJ, the MRM for C_CO2 included pH, TSS, K₂₅, and N-NO₃ as independent variables and explained 76% of the C_CO2 variance. Moreover, the variable with the greatest weight within the SJ model was the pH (Table 4). On the other hand, in the eutrophic SL, the most important predictors in the MRM were Chl-a, TSS, DO, and N-NO₃ which explained 94% of the C_CO2 variance in the Z_MIX (Table 4).

3.3. CO₂ Evasion Rates

CO₂ efflux was recorded in all the lakes with significant variation along the trophic gradient (Kruskal–Wallis: X² = 19.2, p < 0.001). The mean F_CO2 was about 200–600% higher in the eutrophic SL than in the mesotrophic SJ and the oligotrophic TZ (Figure 3A).

In TZ, the F_CO2 values ranged from 0.14 to 1.43 g C m⁻² d⁻¹ during the rainy/warm season and from 0.25 to 1.0 g C m⁻² d⁻¹ over the dry/cold season (Kruskal–Wallis: X² = 4.9, p < 0.05). TZ registered the lowest F_CO2 values of the rainy/warm season. No significant diel variation was observed in TZ (p = 0.07).

During the rainy/warm season, F_CO2 values in SJ ranged from 2.3 to 7.4 g C m⁻² d⁻¹, while during the dry/cold season, they ranged from 1.0 to 2.3 g C m⁻² d⁻¹ (Kruskal–Wallis: X² = 19.86, p < 0.05). The dry/cold season’s highest F_CO2 values were recorded in SJ without significant daily variation (p = 0.96).

The F_CO2 values in SL ranged from 3.5 to 16.6 g C m⁻² d⁻¹ during the rainy/warm season while from 0.03 to 0.5 g C m⁻² d⁻¹ over the dry/cold season (Kruskal–Wallis: X² = 15.0, p < 0.05). SL registered the highest F_CO2 values for the rainy/warm season but the lowest over the dry/cold season. No significant diel variation was observed in SL (p = 0.88).

There was a positive and moderate correlation between the mean C_CO2 in the Z_MIX and F_CO2 (p < 0.6001; Figure 3B). The corrected linear model between the C_CO2 and F_CO2 (F_CO2 = −3.34 + 0.24 C_CO2) explained 50% of the variation in F_CO2 and there was no evidence of autocorrelation in the model (DW = 1.97, p = 0.46).

The three lakes of Montebello were net emitters of CO₂ into the atmosphere in both seasons. When upscaling emissions to the surface area of lakes and considering the lake size disparity (

Table 5. Mean values (X) and standard deviation (SD) of total CO₂ emission rates (F_CO2) and upscaling to the total lake area. (R/W: rainy/warm season, D/C: dry/cold season, TZ: oligotrophic Tziscao, SJ: mesotrophic San José, SL: eutrophic San Lorenzo, n = number of observations).

Lake	Season	n	FCO2 (g CO2 m−2 d−1)		Area	FCO2 (t CO2 d−1)
			X	SD	(m2)	X	SD
TZ	R/W	9	0.88	0.46	3.06 × 10−6	2.68	0.001
	D/C	16	0.59	0.02	3.06 × 10−6	1.81	4.49 × 10−5
SJ	R/W	12	4.63	0.88	6.06 × 10−5	2.8	0.002
	D/C	16	1.42	0.18	6.06 × 10−5	0.86	0.0001
SL	R/W	10	9.13	1.24	1.81 × 10−6	16.54	0.02
	D/C	11	0.29	0.02	1.81 × 10−6	0.53	9.37 × 10−6

), it was estimated that the oligotrophic TZ emitted 2.3 ± 0.6 t CO₂ d⁻¹, the mesotrophic SJ 1.8 ± 1.4 t CO₂ d⁻¹, and the eutrophic SL 8.5 ± 11.3 t CO₂ d⁻¹.

4. Discussion

Our results suggest that the seasonal dynamics of C_CO2 and F_CO2 in karst lakes resulted from a complex interplay of biotic and abiotic variables exposed to strong tropical seasonality. Although the NPP rate and biomass of primary producers (Chl-a concentration) displayed a slight, almost negligible seasonality, the increase in temperature during the rainy/warm season raised C_CO2 and F_CO2 values in all lakes.

Nutrient availability is a key factor in primary production in lakes and the seasonal dynamics of CO₂ [3,8]. Because the thermal stratification depleted the nutrient availability in the Z_MIX of the oligotrophic TZ and mesotrophic SJ, heterotrophic metabolism also determined the amount of CO₂ dissolved in the water column, as usually observed in other lakes [60]. Despite the fact that the increase in temperature stimulates the release of nutrients from the sediments through mineralization, and the increase in runoff results in increased inputs of particulate organic matter and nutrients from the watershed [9], increased nutrient concentration in the lake water column is quickly depleted by the uptake and by the storage in the hypolimnion. The inverse relationship between DO and C_CO2 probably exhibits the importance of heterotrophic respiration in the CO₂ balances of these karst lakes.

NPP rates were about 46% higher during the rainy/warm season compared to the dry/cold season, except for the oligotrophic TZ, where NPP increased in the cold/dry season but independently of phytoplanktonic biomass, which hardly changed in both seasons. According to Chl-a concentrations and NPP rates, the oligotrophic TZ was a heterotrophic system (GPP < R) during the rainy/warm season, whereas it was autotrophic (GPP> R) during the cold/dry; on the contrary, the mesotrophic SJ was autotrophic during the rainy/warm season and heterotrophic over the dry/cold one. The eutrophic lake SL acted as an autotrophic system in both seasons, most likely associated with continuous nutrient availability.

4.1. Variation of C_CO2 in the Water Column: Seasonality and Main Drivers

As expected, the three studied lakes were supersaturated in CO₂ during both seasons. The leading causes for CO₂ supersaturation are usually related to high respiration rates [61], high DIC inputs from the surface or groundwater, high carbonate dissolution [24] high methane oxidation [23], and to a lesser extent with DOC photochemical degradation [62].

There is considerable evidence of the relationship between lake eutrophication and C_CO2 [8,9,13]. According to our results, eutrophication raised C_CO2 consistently throughout the year. The models that best explain the C_CO2 variation in the oligotrophic TZ and the mesotrophic SJ included pH, TSS, DIC, N-NO₃, and K₂₅ which probably evidence the role of the groundwater inflow in the CO₂ balance. In contrast, the variables with more influence within the eutrophic SL model were Chl-a, TSS, DO, and N-NO₃, which probably represent the influence of biological activity on CO₂ dynamics.

Thermal stratification strongly controlled C_CO2 vertical patterns among the lakes [63,64]. Over the stratification, the hypolimnetic anoxia mainly due to organic matter oxidation in the sediments should stimulate the production of CH₄ [65], but probably most CH₄ could have been consumed under the anaerobic oxidation of methane (AOM; [23]) and released as CO₂ to the water column. The strong relationship found between dissolved CO₂ and CH₄ measured in these lakes (C_CO2 = 0.44(C_CH4) + 31.3, R² = 0.88, p < 0.05; unpublished data of Vargas-Sánchez et al.) would suggest around 44% of CO₂ in the water column is released through AOM.

Therefore, over the rainy/warm season, the supersaturation in CO₂ in eutrophic SL should mainly be a consequence of the high rates of anaerobic respiration coupled with AOM. In oligotrophic TZ and mesotrophic SJ, supersaturation should be related more to the groundwater inflow, which is rich in DIC as in most karst landscapes [66] and by the lower NPP. In contrast, during the cold season, while the eutrophic SL exhibits the same pattern but with less intensity because of the temperature, the oligo TZ and mesotrophic SJ remove part of the C_CO2 through primary production. However, the absence of a statistical relationship between Chl-a or NPP and C_CO2 reflects the importance of other abiotic processes, such as groundwater inflow [67] or carbonate precipitation (SIcalcite > 0). However, this approximation does not allow us to know in detail how much carbon leaves the system due to precipitation.

4.2. CO₂ Evasion from Tropical Karst Lakes: Eutrophication Effects

The supersaturation of C_CO2 in the surface water layer (modulated by the physicochemical variables, DIC allochthonous inputs, and metabolic activity) determines the F_CO2 rates [1,3]. F_CO2 exhibited notable variation among the lakes with mean values about 10 times greater in SL compared to TZ over the rainy/warm season and 5 times higher in SJ compared with SL in the dry/cold season. The F_CO2 depended on the C_CO2 in the water column (50% of explained variance). This relationship evidences the importance of the biogeochemical processes that control the CO₂ balance versus those that determine the transfer between water and the atmosphere in this lake ecosystem, supporting the simplified model proposed by Cole et al. [68].

Our results suggest that eutrophication is the primary driver of F_CO2 in the LMNP tropical karst lakes. The eutrophic SL emitted more CO₂ than the oligotrophic TZ and mesotrophic SJ during the rainy/warm season. On the one hand, our results are opposite to Balmer and Downing [20] and Macklin et al. [67], who argued that eutrophic lakes are typically undersaturated in CO₂ and, due to their high primary production, are atmospheric CO₂ sinks [21]. On the other hand, our results are consistent with Huttunen et al. [17], Zhou and Beck [18], and Zhou et al. [19] who affirmed that the increase in OM mineralization leads to higher CO₂ evasion rates. Eutrophication in karst lake SL promotes CO₂ outgassing (9.13 ± 1.4 g C m⁻² d⁻¹), acting as a source of CO₂ to the atmosphere. Under these conditions, we hypothesize that CO₂ uptake cannot offset the methanogenesis and the release of CO₂ through AOM, reversing the sink capacity of these lake ecosystems, along with other additional sources such as groundwater DIC-rich inflow entering the lake.

In contrast, during the dry/cold season, although F_CO2 decreases significantly, oligotrophic TZ and mesotrophic SJ emit more than the eutrophic SL, probably because methanogenesis and AOM are temperature-dependent processes [23,69]. Therefore, the expected temperature rise associated with climate change could further exacerbate CO₂ emissions from eutrophicated karst lakes. Thus, during the cold/dry season, the supersaturation and subsequent degassing of CO₂ is largely controlled by DIC inputs from the surrounding basin [3], in this case, from groundwater.

F_CO2 measured in these tropical karst lakes was in the range of other tropical, temperate, and boreal lakes, demonstrating that the variability between lake ecosystems is not as extensive as we might initially think (

Table 6. Comparison between the F_CO2 rates (g CO₂ m⁻² d⁻¹) obtained in this study and data previously reported for other water bodies with different trophic statuses. O = oligotrophic, M = mesotrophic, E = eutrophic, NM = not mentioned. (R/W: rainy/warm season, D/C: dry/cold season, TZ: Tziscao, SJ: San José, SL: San Lorenzo).

Season	Lake	FCO2		Trophic	Reference
		X	SD	Status
R/W	TZ	0.9	0.5	O	This study
	SJ	4.6	0.9	M	This study
	SL	9.1	1.2	E	This study
Warm	Poza Churince	1.2	1	O	[70]
	Los Hundidos	5.6	2.7	O
	Poza Manantial	15.3	10.2	O
	Poza Azul	6.5	2.7	O
	Poza Becerra	22.9	5.8	O
	Taihu	3.3	1.2	E	[71]
Stratification	Alchichica	0.1	0.4	O	[72]
	Venasjön	0.2	0.8	E	[73]
D/C	TZ	0.6	0	O	This study
	SJ	1.4	0.2	M	This study
	SL	0.3	0	E	This study
Cold	Los Hundidos	3.2	1.6	O	[70]
	Poza Azul	7.7	3.7	O
	Poza Becerra	25.6	7.2	O
	Taihu	0.2	0.8	E	[71]
Mixing	Alchichica	0.1	0.4	O	[72]
	Sau	0.3		E	[74]
	Venasjön	3.3	1.2	E	[73]
NM	Ljusvanttentjärn	0.4	0.1	O
	Parsen	1.1	0.5	M
	Badger	2.2		E	[75]

). Nonetheless, it should be noted that, on the one hand, the mean F_CO2 in the oligotrophic TZ was generally lower than those in other oligotrophic lakes. Contrarily, the mean F_CO2 in mesotrophic SJ and eutrophic SL were higher than in other meso- and eutrophic ones (Table 6).

5. Conclusions

This study is the first to report the C_CO2 and the F_CO2 into the atmosphere of the tropical karst lakes of the LMNP in Mexico. To our knowledge, it is also among the first research documenting eutrophication’s impacts on CO₂ dynamics in tropical karst aquatic ecosystems. Whereas hydrodynamics (stratification, mixing) can drive C_CO2 and F_CO2 in oligo and mesotrophic lakes through the control of NPP, DIC inflows by groundwater discharges exert the main control over the dynamics of CO₂. However, in eutrophic karst lakes, the heterotrophic oxidation of the large amount of organic matter accumulated in the water column and sediments depletes dissolved oxygen from a high fraction of the water column. It promotes anaerobic metabolism using methanogenesis and AOM. Therefore, the high amount of CO₂ entering the lake by groundwater inputs and AOM cannot be removed by phytoplankton uptake and is then largely emitted into the atmosphere. Since cultural eutrophication is due to anthropogenic activities, the excess emissions of greenhouse gases produced globally due to this cause must be considered within the global emissions of anthropogenic origin since we are not counting an important source of greenhouse gases in predicting effects on future climate.