Assessment of Water Quality and Heavy Metal Environmental Risk on the Peri-Urban Karst Tropical Lake La Sabana, Yucatán Peninsula

Abstract

In the Yucatán Peninsula, anthropogenic activities such as urbanization and the final disposal of solids and wastewater critically impact aquatic systems. Here, we evaluated the anthropogenic-related environmental alteration of Lake La Sabana, located in the northern limits of one of the main cities of the Mexican Caribbean. We evaluated lake water quality, using physical, chemical, and microbiological indicators, and heavy metals in surficial sediment and fish tissue to evaluate the potential environmental risk. Multivariate analyses revealed that environmental conditions in La Sabana are spatially and temporally heterogeneous. Medium to bad water quality was determined within different basins by the National Sanitation Foundation water quality index, related to the degree of anthropogenic influence at each zone. The central-south zones displayed critical microbiological values largely exceeding national standards. Heavy metals in sediment and fish tissue such as Zn and Hg were relatively low, but Hg concentrations threaten the ecological environment. Incipient wastewater treatment and its final disposal in La Sabana are mainly responsible for the changes in the trophic status and availability of nutrients, which in turn may have promoted changes in the biological structure and aquatic plant invasions. Lake La Sabana can be considered a model of the potential and sequential effects of anthropogenic alterations in the oligotrophic karst tropical aquatic systems in the Yucatan Peninsula.

1. Introduction

The Yucatán Peninsula in southern Mexico is a continental karstic platform composed of highly permeable limestone deposits [1,2,3]. It is a unique hydrological region characterized by the most prominent subterranean aquifer in the world (extending over 165,000 km²), more than 8000 cenotes (freshwater dolines), and little availability of surficial waters with the practical absence of rivers and fewer than 10 lakes > 10 km² [4,5]; most of the inland aquatic environments interact to a certain degree with the surrounding marine environments such as the Gulf of Mexico and the Caribbean Sea [6,7]. Freshwater biodiversity is highly endemic in the Yucatán Peninsula, with more than 40 hypogean endemic species [8], while in epigean environments, groups such as Ostracoda (crustacea), display about 70% of endemism to different extents, from single lakes to the entire region [9].

One of the most important environmental and ecological problems in the Yucatán Peninsula is the management of wastewaters and the pollution aquatic environments [10,11]. This region is inhabited by about five million people, of which 52% are concentrated into five main cities: Mérida, Cancun, Playa del Carmen, Campeche, and Chetumal [12]. In these cities, for example, only between 35 and 30% of the population use a wastewater treatment system, whereas the remaining 70% of the population in cities and all the population in towns discharge wastewater into the ground in septic tanks or similar systems [11,13]. Given the geomorphology and high karstification of this region, the final disposal sites of treated or untreated water are the subterranean aquifer, lakes, wetlands, or the sea, facilitating the dispersal of pollution. In subterranean waters, water quality monitoring has been intensive, and it has demonstrated the presence of contaminants such as pharmaceuticals (acetaminophen and ibuprofen) [13], personal care products (triclosan and synthetic musks) [14], polycyclic aromatic hydrocarbons [15], herbicides (chlorophenoxy), heavy metals [16], pesticides (organochlorine) [17], bacteria (coliforms) [18], and viruses. In surficial waters, little effort has been made to continuously monitor the effects of anthropogenic activities on the water quality and ecological interactions. Surficial waters such as lakes and aguadas (permanent and ephemeral freshwater ponds) are key for several ecological processes, such as migration, reproduction, and as feeding grounds, and have also been fundamental for human population development, since prehispanic times. In the Yucatán Peninsula, lakes and aguadas are mainly distributed in the southern region, being flooded depressions dependent on precipitation or fed by the groundwater flow from the aquifer [1,2].

La Sabana is a freshwater lake, with the highest degree of anthropogenic alteration in the southern Yucatán Peninsula. This lake belongs to the Bacalar hydrological system, characterized by oligotrophic, well-oxygenated, and alkaline waters with high concentrations of calcium (Ca²⁺), bicarbonates (HCO3⁻), magnesium (Mg²⁺), sulfates (SO₄²⁻), and chloride [19]. La Sabana is a peri-urban lake located along the northwest limits of Chetumal city. Chetumal is one of the main cities in the Mexican Caribbean with a population growth rate of about 2.9% in 2010 [20], which was higher than the Mexican national average of 1.4%. In 2020, Chetumal city had a population of about 233,648 inhabitants [21]. Since 1999, about 120 lt s⁻¹ of treated water has been discharged into Lake La Sabana, which corresponds to the wastewater produced by ~40% of the city’s population [22]. During the last 10 years, alterations of the natural conditions, related with increase urbanization, has occurred in La Sabana, including the modification of littorals, blocking of natural surficial flows, the removal of littoral native vegetation and changes in the trophic status caused by nutrient increase. Currently, how wastewater pollution andenvironmental alterations in La Sabana have modified water quality and hydrological dynamic, and what has been the responses of aquatic biota is practically unknown; thus, governmental actions toward the restoration and conservation of the lake are void. Understanding the environmental and ecological responses to pollution and the rapid nutrient increase in La Sabana is relevant for modelling the potential effects of contamination in the surficial aquatic environments at a regional scale in the Yucatán Peninsula.

Environmental monitoring is the key to infer changes in the natural conditions in aquatic ecosystems. The multiparameter monitoring from sources such as the water column to infer water quality, and from sediments to detect contaminants such as heavy metals, is the fundamental to capture broad aspects of the environment, and to infer more clearly the potential effects on nature and society [23].

In this study, we comprehensively evaluated the anthropogenic-induced environmental alteration in Lake La Sabana by spatially analyzing (1) the water quality using physical, chemical, and microbiological indicators; (2) heavy metals in sediments as evidence of leeching, and (3) the presence of heavy metals in fish tissue as a consequence of contaminant transfer to food webs. We estimated the current environmental risk of heavy metals in sediment and fish tissue and documented how the increase in nutrient availability has promoted sequential alterations in the lake ecosystem structure and functionality. A bathymetric map of the Lake La Sabana is provided.

1.1. Study Site

Lake La Sabana is a freshwater body located along a geological fault zone in southern Quintana Roo, Mexican federal state [24]. The lake has an elongated form with approximately 6.5 km length (considering the associated wetland) and less than 600 m width (Figure 1). La Sabana is currently divided into four sub-basins, given the construction of three dams used to access from one side of the lake to the other, each basin having variable conditions regarding water color, transparency, and biota [25,26]. La Sabana belongs to the Bacalar hydrological system and, it is assumed, has a subterranean connection with it, as a surficial connection is lacking. The main basin (northern basin) is about 3 km long and <600 m wide; the southern basins correspond to the wetlands and are more sensitive to dry seasons and can dry out sporadically. The climate in the region is tropical (A-climate “tropical/megathermal climate”) and warm subhumid (Aw) according to the Köppen classification [27]. The mean annual temperature is 26 °C and is related to the descending limb of a Hadley cell, centered at 20° N [28]. Precipitation is the most fluctuating variable in the region, and it depends on the seasonal migration of the Intertropical Convergence Zone (ITCZ). The rainy season occurs in summer, associated with the northern position of the ITCZ in the American continent, whereas the dry season occurs during winter–spring. The mean annual precipitation in the Yucatán Peninsula is 1600 mm but ranging monthly from 300 to 1800 mm [29,30].

1.2. Water Physical, Chemical and Microbiological Parameters, and Quality Assessment

Lake La Sabana was divided into four zones: north, center, south, and wetland (Figure 1). At each zone, three sampling sites were established, except for the wetland, in which six sampling sites were considered because of its size. Seven environmental parameters were quantified at each site: temperature, conductivity, pH, and dissolved oxygen were measured in situ with a Hach HQ40D multiparametric probe; transparency was measured with a turbidity tube; 5-day biochemical oxygen demand (BOD₅)) was measured with a portable dissolved oxygen and BOD meter HI98193 from Hana instruments. Total suspended solids (TSS) were quantified gravimetrically using a glass fiber filter of 2 µm. For nutrient content determination, 1.5 L of water was collected in plastic bottles, previously washed with distillated water. Collected samples were preserved in ice in situ and then transported to the Laboratorio de Química analítica from El Colegio de la Frontera Sur, Chetumal unit. Nitrites, nitrates, orthophosphates, and sulfates were measured with a Shimadzu UV-1700 spectrophotometer. Samples were measured in triplicate and for all runs two blank samples were used. Calibration curves reached an r² = 0.990–0.995. All variables were measured, near surface at 0.5 m depth, during February and August 2021, coinciding with the dry and rainy season, respectively.

For the microbiological analysis, three surficial water samples (at 0.5 m depth) were collected at each lake zone during May (coinciding with the dry season), and July 2021 (coinciding with the rainy season). We applied the most probable number (MPN) technique to estimate total and fecal coliforms in La Sabana, as threshold values of the Mexican national regulation NOM-210-SSA1–2014 are expresses in MPN units. This technique is used to estimate microbial populations in water and to determine the probability of containing disease-producing organisms that make the water unsafe for consumption or recreational activities [31,32]. We used sterilized Whirl-Pak bags and collected about 100–150 mL of sub-surficial waters (0.5 m depth). Samples were labeled and stored in ice and then transported to the Laboratory of Microbiology of food and water at the Laboratorio Estatal de Salud Pública del Estado de Quintana Roo for analysis. Presumptive tests for both total and fecal coliforms were first conducted to detect coliform presence in samples. We used the method proposed by Peeler et al. [33], which is based on lake water inoculation in a series of lactose-rich media at different concentrations. Given that all presumptive tests for both total and fecal coliforms were positive, we conducted confirmatory tests. For total coliforms, from each tube with presence of gas, a loopful of the sample was streaked onto the selective medium brilliant green bile lactose broth and incubated at 35 +/− 0.5 °C for 24–48 hrs. Confirmatory tests for fecal coliforms were conducted by inoculating subsamples from positive tubes into the selective medium EC broth (Escherichia coli Broth). The determination of viable numbers of bacteria in the samples was based on the Bacteriological Analytical Manual (BAM) Appendix 2: Most Probable Number from Serial Dilutions of the U.S. Food and Drug Administration (FDA) [34]. For comparison purposes, we used the threshold values of the Mexican Standard NOM-210-SSA1–2014 for total coliforms, which established guidelines for water consumption and recreational activities. Fecal coliform reference values were obtained from the Laboratorio Estatal de Salud Pública, Laboratorio de Patógenos. Reference values for both total and fecal coliforms are expressed in MPN 100 mL⁻¹.

The National Sanitation Foundation water quality index (NSF-WQI) was applied to obtain relative values of quality in La Sabana. Nine environmental parameters were used: pH, dissolved oxygen, turbidity (derived from transparency tube values, expressed in nephelometric units), fecal coliforms, BOD₅, orthophosphates, nitrates, temperature, and TSS. The NSF-WQI weighted average of each parameter and its relative importance was summed to obtain values between 0 (zero) representing very bad quality and 100 representing excellent quality [35].

1.3. Heavy Metals in Sediment, Fish Tissue and Ecological Risk Assessment

A total of eight sampling sites were selected for heavy metal determination in sediment. Two sampling sites were taken in the north (N1, N3), center (C1, C3), south (S1, S3), and wetland (W1, W2). Samples were collected during August–September 2020 with the aid of an Ekman dredge. Only the uppermost 3 cm of each grab were used for analysis. For heavy metal determination in fish tissue, a total of 15 individuals of the non-native fish Oreochromis niloticus were collected in La Sabana during November 2020 and January 2021 using an individual cast net. After the capture, individuals were kept in ice and immediately transported to the Laboratorio de Zoología of the Instituto Tecnológico de Chetumal and stored at −5 °C. The total lengths and weights of the fishes were recorded. For heavy metal determination, fish liver was used because it effectively accumulates metals that cannot be metabolized or assimilated by the individuals [36]. Heavy metal concentration in the liver can be comparable or even higher than in other tissues such as muscle or gills [37]. For the fish of La Sabana, the whole liver was removed and stored in Eppendorf tubes with 70% ethanol until analysis.

For both sediment and fish liver, Pb, Cd, Zn, and Hg were measured, as they are the most frequent contaminants in sediment [38,39] and biological groups [40,41] in the region. The atomic absorption spectrophotometry technique was used. For sediment digestion, standard methods 3051A [42] and 3015A [43] were used. Fish liver acid digestion was performed as described in Buenfil-Rojas et al. [40]. Graphite furnace atomic absorption spectrophotometry (GFAAS) (Avanta PM–105 GF3000, GBC) was used for Pb, Cd and Zn and hydride generation atomic absorption spectrophotometry (HGAAS) (Avanta PM–HG 106 3000, GBC) for Hg. Quality control for heavy metals in sediment was performed as in Tun-Canto et al. [39]. The recuperation percentage was between 90 and 110% for all elements. The accuracy of the analysis of the heavy metals in the fish tissue was determined using certified reference material TORT-2 (lobster 72 107 hepatopancreas) from the National Research Council of Canada (NRC-CNRC) and sample blanks. For reference material, Hg concentration was 0.27 ± 0.04 μg g⁻¹ dry weight and our results were on average 0.25 ± 0.009 μg g⁻¹ dry weight, with a recuperation percentage of 95 ± 3%.

Two approaches were used to evaluate the relative ecological risk of heavy metal in sediment: the total metal content measure, based on the Canadian Sediment Quality Guidelines (CSQG) for the protection of aquatic life developed by the Canadian Council of Ministers of the Environment, that establishes threshold values for ecotoxicological effects on aquatic biota [44]; and the total content index that uses the ecological risk factor (Er). This index is calculated using the following equations:

$$\mathrm{E}r={\mathrm{T}}_{r }\times {\mathrm{C}}_f$$

(1)

$$\mathrm{C}f=C_i/{\mathrm{C}}_n^i$$

(2)

where T_r is the toxic response factor values for each different metal, and C_f is the contamination factor, determined by the content of target metal in the sediment (C_i), divided between the background value of metal in the study area or consensus reference values for a given substance (${\mathrm{C}}_n^i$) [45]. For our samplings, given the lack of heavy metals reference values in La Sabana and the surrounding environments, the consensus preindustrial reference values were taken: Pb = 70, Cd = 1, Zn = 175, and Hg = 0.25 [45]. Toxic response factors are as follows: Pb = 5, Cd = 30, Zn = 1, and Hg = 40 [45,46]. RI consists of five levels of environmental risk: RI < 40 suggests low risk; 40 ≥ RI < 80, moderate risk; 80 ≥ RI < 160 considerable risk; 160 ≥ RI < 320 high ecological risk; and RI > 320, very high ecological risk [43,47].

The relative contamination of heavy metal in the fish liver of O. niloticus from La Sabana was assessed by comparison with national standards and with the metal pollution index (Pi), which is a mathematical method to assess and monitor contamination of an individual metal or a set of heavy metals in aquatic ecosystems [48]. The Mexican Standard NOM-242-SSA1-2009, which determines sanitary specifications for fishing fresh products for human consumption was used for comparison. The Pi was calculated using the following equation [49]:

$$\mathrm{Pi}=\mathrm{C}n/\mathrm{GB}$$

(3)

where Cn is the content of metal in sediment and GB is the geochemical background value. For this study, the geochemical background value of Hg was taken from the NOM-242-SSA1-2009, in fish samples of 0.5 mg kg⁻¹ wet weight. Pi values classify four pollution levels, Pi < 0.2 no significant pollution; 0.2 ≥ Pi < 0.6 minor pollution; 0.6 ≥ Pi < 1 moderate pollution; and Pi > 1, severe pollution. The Pearson correlation test was used to identify the relation of fish length and weight with heavy metal content in liver.

1.4. Spatial and Temporal Environmental Variability between Lake Zones

For the statistical analysis, the physical and microbiological data with different measured units were normalized. In order to determine if the samples could be assigned to statistically significant groups regardless of the known factor levels, a hierarchical cluster analysis (CLUSTER) coupled with a similarity profile test (SIMPROF) with 999 permutations and a significance level of 5% was performed based on Euclidean distances. To explore the multivariate composition of the samples, a non-metric multidimensional scaling (nMDS) was used. Significant groups identified by cluster and SIMPROF were overlaid as ellipses in the nMDS plot.

A test of homogeneity of dispersions (PERMDISP) and a two-way permutational multivariate analysis of variance (PERMANOVA) were used to detect the significant differences in the dispersion and location, respectively, among the four levels of the lake zone factor (north, center, south, wetland) and the two levels of the season factor (dry, rainy). These analyses were performed with 999 permutations based on the Euclidean distances.

The relation between the studied samples and the environmental variables was evaluated by means of the ordination obtained in a principal components analysis (PCA). All calculations were performed in the statistical software PRIMER v7.0.13 [50].

Bathymetry of Lake La Sabana

A rigid-bottom inflatable boat of about 2.8 m long and 1.2 m beam, equipped with a 6.5 hp outboard motor, was used as platform to determine the bathymetry of La Sabana. Given the elongated morphology of the lake, we established seven main longitudinal transects and three transversal transects in each of the basins of the lake. Additional transects were performed in the lake margins. All transects were equidistant, with about 60 m of distance between each other and a total of 726 depth measurements were taken from La Sabana. The resulting mesh had irregularly distributed sampling points. Latitude, longitude, and depth of sampling points were recorded with a Garmin (Lenexa, Kansas) eco-sounder model Echomap Chirp 42CV. GPS had an accuracy of 1–3 m, and the transducer had a precision of 95%. A polygon of the coastline was delineated with satellite imagery, using the google satellite base map and then converted to a vectorial file in the open source QGIS software 3.22.1. The XYZ data (longitude, latitude, and depth) were visually inspected in an excel work sheet. Then, using the Surfer software (Golden Software Inc., Golden, CO, USA), we applied the Kriging interpolation algorithm to obtain the bathymetric map and the isobaths to represent water depth.

2. Results

2.1. Coliform Bacteria and NSF Water Quality Index in La Sabana

The spatial and temporal patterns of the distribution and concentration of coliform bacteria was variable in La Sabana. During the dry season, the coliforms were almost homogenously distributed in the central-south part of the lake, where they reached concentrations >16,000 MPN 100 mL⁻¹ and > 1000 MPN 100 mL⁻¹ for total and fecal coliforms (Figure 2a,b), respectively. The north was characterized by drastic reductions in microbial concentrations with an average of ~1000 and ~100 MPN 100 mL⁻¹ for total and fecal coliforms (Figure 2a,b), respectively, whereas in the wetlands, the values were variable with generally low values of ~1100 and ~120 MPN 100 mL⁻¹. During the rainy season, the interpolated maps showed that coliforms were almost restricted to the center of the lake (Figure 2c,d). The total coliform concentration significantly increased to > 50,000 MPN 100 mL⁻¹ and the fecal coliform concentration to >18,000 MPN 100 mL⁻¹ compared with the dry season. Lower bacterial values were observed in the wetlands where values remained <2500 and <260 MPN 100 mL⁻¹ for total and fecal coliforms, respectively (Figure 2c,d). Compared with reference values of the Mexican Standard NOM-210-SSA1-2014, which indicates a threshold of <16,000 MPN 100 mL⁻¹ for total and <450 MPN 100 mL⁻¹ for fecal coliforms, for recreational use of waters, La Sabana largely exceeded the threshold in almost all sections during the two climatic seasons. The wetlands were the exception as coliform values remained below the threshold, during both dry and rainy seasons (Figure 2a–d).

The NSF WQI was relatively homogeneous spatially and temporally, with values fluctuating from 41 to 61 (Figure 2e,f and Table S1). The north showed values between 48 and 56, suggesting that water quality fluctuates between bad to medium conditions throughout the year. The center displayed mostly medium water quality conditions (NSF WQI= 53 average) during the dry season but dropped to bad conditions in the rainy season (NSF WQI = 48). The south displayed mostly bad conditions throughout the climatic seasons with NSF WQI values ranging from 46 to 48. The wetlands demonstrated the higher values of the index ranging from 57 to 61, and thus suggesting medium water quality conditions throughout the year.

2.2. Heavy Metals in Surficial Sediment and Ecological Risk Index

The heavy metal analysis showed two available elements in the surficial sediment of Lake La Sabana, Zn and Hg (

Table 1. Heavy metal concentrations and basic statistical metrics in surficial sediment in La Sabana and associated ecological risk factor.

Site	Heavy Metals in Surficial Sediment (in mg kg−1)				Ecological Risk Factor (Er)
	Pb	Cd	Zn	Hg	Zn	Hg
N1	<LOD	<LOD	73.2	0.25	0.41	41.29
N3	<LOD	<LOD	52.1	0.27	0.29	43.88
C1	<LOD	<LOD	590.1	0.44	3.3	70.99
C3	<LOD	<LOD	40.4	0.12	0.23	19.55
S1	<LOD	<LOD	41.3	0.14	0.23	22.94
S3	<LOD	<LOD	664.5	0.30	3.79	49.24
W1	<LOD	<LOD	104.2	0.34	0.59	55.37
W2	<LOD	<LOD	42.3	0.31	0.49	49.6
Mean	-	-	202.1	0.27
St dev	-	-	264	0.10
Cv (%)	-	-	130.6	37.9

), both detected in 100% of the samples. The Cd and Pb were below the limit of detection (LOD = 0.125 and 0.12 mg g⁻¹, respectively) in all samples (Table 1). Therefore, Cd and Pb were excluded from further analyses. The concentration of Zn was relatively higher than Hg and displayed a higher variability between the sites. Table 1 shows Zn and Hg concentrations (in mg kg⁻¹ wet weight) per site, and the basic statistics of total content in La Sabana. The CSQG revealed that Zn in most sites was below the threshold effect levels (TELs, Zn = 124 mg kg⁻¹), except for the sites C1 and S3, for which concentrations were above the probable effect level (PEL = 271 mg kg⁻¹), suggesting that adverse effects for aquatic biota may occur in those sites. For Hg, the concentrations in all sites were between TELs (0.13 mg kg⁻¹) and PEL (0.7 mg kg⁻¹), and therefore falling within the range where adverse effects may occasionally occur [51]. The Er which is based on a score index, highlights that Zn concentrations were a low environmental risk, whereas for Hg, moderate environmental risk was detected.

2.3. Heavy Metal in Oreochromis niloticus and Metal Pollution Index

The total length of O. niloticus ranged from 10 to 26 cm, while the weight ranged from 170 to 300 gr (

Table 2. Size measurements and Hg concentration in Oreochromis niloticus from Lake La Sabana.

Weight (gr)	Total Length (cm)	Hg Concentration (mg kg−1)	Metal Pollution Index
300	26	0.0085	0.17
200	22	0.019	0.38
225	21	0.0247	0.494
300	25.5	0.0205	0.41
200	22.5	0.1289	2.578
190	20	0.0233	0.466
350	26	0.0254	0.508
200	23	0.0145	0.29
330	25	0.0278	0.556
210	25	0.0227	0.454
200	22	0.0129	0.258
180	19	0.1196	2.392
200	21	0.1362	2.724
190	19	0.0451	0.902
170	18	0.2073	4.146

). Of the four nonessential heavy metals suspected to be available in the lake water and to be accumulated in the fish livers, Zn, Cd and Pb were below LOD (0.065, 0.125 and 0.12 mg g⁻¹) of the analyte, thus suggesting negligible concentrations in the liver. Mercury (Hg), in contrast, was detected in all specimens, with values ranging from 0.008 to 0.2 mg kg⁻¹. The Hg values were mostly below the threshold of 1 mg kg⁻¹ of Hg and 0.5 mg kg⁻¹ of methylmercury, established by Mexican Standard NOM-242-SSA1-2009, which determines sanitary specifications for fresh fishing products for human consumption. The metal pollution index, however, demonstrated that most individuals fell in the moderate pollution category and in four individuals, severe pollution was inferred. The Pearson correlation test revealed negative correlation between fish length and Hg content in liver (r = −0.59, p-value 0.01 at 0.05 significance level) suggesting that smaller individuals have higher Hg concentrations, whereas statistical significance is lacking for fish weight and Hg concentration (r = −0.45, p-value 0.08 at 0.05 significance level).

2.4. Lake La Sabana Bathymetry, Spatial and Temporal Environmental Variability

The bathymetric map revealed that La Sabana is a shallow lake that does not exceed 1.9 m depth during the dry season (Figure 3). The northern basin displays deeper waters, whereas in the wetlands the average depth is 1.1 m (Figure 3). During the rainy season, lake water levels can increase by 40 cm, and then a water outlet is formed in the northernmost section, discharging waters northward toward the Chetumal bay (Caribbean Sea). The Lake area is 1.48 km², the perimeter is 18.81 km and the length through the middle of the lake is 6.43 km.

The SIMPROF test projected on the nMDS ordination plot, with a good level of stress (0.09), clearly showed five significant groups (Figure 4). The first group comprised samples from the center zone of the lake in the rainy season (Figure 4). The samples from wetland formed two well-defined and coherent groups, one corresponding with samples obtained in the rainy season and the other one comprising samples from the dry season (Figure 4). The remaining two groups included samples from the north, south and central zones of the lake, differentiated by the climatic seasons (Figure 4).

PERMDISP showed that in the case of the seasonal factor, the multivariate distances from each group’s centroids were homogeneous (F = 3.024, p perm = 0.131), but in the case of the lake zone factor, the dispersions were heterogeneous (F = 7.637, p perm = 0.001). Two pairwise compared lake zones resulted with significantly different dispersions after the Benjamini–Hochberg correction [52]: central wetland and south wetland (p perm = 0.001 and p perm = 0.002, respectively). Given that multivariate dispersions were not homogeneous and that the design was unbalanced, the pseudo-F ratio as well as the permuted p value of PERMANOVA were computed using Type III sums of squares and permuting residuals under a reduced model [53]. For pairwise PERMANOVA tests, Monte Carlo simulations were also computed, taking into account the limited number of permutations available in some cases.

The results of the PERMANOVA showed that significant differences existed between the levels of the lake zone factor and between the levels of the season factor (p ≤ 0.001 in both cases;

Table 3. Table of the principal results of permutational multivariate analysis of variance (PERMANOVA) test. The total variation is partitioned according to four sources: lake zone, season, lake zone × season—interaction between lake zone and season factors, residuals. Bold values indicate significant results (** p ≤ 0.01, *** p ≤ 0.001).

Source	df	SS	MS	Pseudo-F	p Perm	Unique Perms
Lake zone	3	100.25	33.418	10.834	0.001 ***	998
Season	1	70.617	70.617	22.893	0.001 ***	999
Lake zone × Season	3	23.175	7.725	2.504	0.003 **	999
Residuals	22	67.862	3.085
Total	29	261

; S2); in the interaction between both factors, significant differences were detected as well (p ≤ 0.01; Table 3). In the pairwise tests, under the term lake zone, all comparison with wetlands resulted significant at α = 0.05 level after the Benjamini–Hochberg correction, both with Monte Carlo and permuted values of p (S3), while north-center and north-south tests resulted not significant when taking the Monte Carlo values of p; albeit, they were significant regarding the permuted probabilities (Table S2). As expected from the PERMANOVA main test, the only pairwise test under the term season (dry–rainy) was significant (p ≤ 0.001) with both Monte Carlo and permuted values of p (S3). Under the interaction term lake zone × season, all the comparisons for pairs of levels of the lake zone factor, both within dry and rainy levels of the season, resulted significant against wetland at α = 0.05 after Benjamini–Hochberg with Monte Carlo values of p (Table S2); due to the reduced number of unique possible permutations in only one of those comparisons (north-wetland within dry level) the test was also significant using the permuted p value (Table 3).

Under the interaction term lake zone × season for pairs of levels of the season factor, although the number of unique permutations available was very restricted, in three out of four cases, taking the Monte Carlo p values instead of the permuted p values after the Benjamini–Hochberg correction, all the comparisons between dry and rainy seasons within any zone level, except one (within south level) were found to be significant at α = 0.05 (Table 3).

The first two principal components obtained by the principal components analysis (PCA) captured 66.9% of the total variation (PC1: 35.6%, PC2: 31.2%). The variables total coliforms, BOD₅, transparency, and fecal coliforms were the most correlated to PC1 and contributed to segregating the samples under the lake zone factor; samples of the wetland zone appeared to be associated with higher values of transparency, while the remaining zones, particularly in the rainy season and, remarkably, samples from the center zone, showed a strong negative association with transparency and strong positive association with total coliforms, BOD5, and fecal coliform values (Figure 5). The most correlated variables to PC2 were temperature, pH, conductivity, and TSS, and they strongly contributed to segregate the samples under the two levels of the season factor; samples obtained during the dry season appeared to be associated to high pH, conductivity and TSS values in contrast to samples obtained in the rainy season, which instead showed higher temperatures (Figure 5).

3. Discussion

3.1. Anthropogenic Alterations in La Sabana: Water Quality, Heavy Metals and Environmental Risk

La Sabana is a shallow lake that does not exceed, in the main basin, 1.9 m during the dry season and 2.3 m in the rainy season. Shallow conditions are typical of the Bacalar hydrological system as surrounding lakes such as lakes Milagros and Bacalar do not exceed 5 and 20 m, respectively [54,55]. Southern basins of La Sabana corresponding to wetlands (Figure 5) are more sensitive to flooding and desiccation [6], whereas in the northern basin, the increase in water levels forms a temporal water outlet that discharges toward the Chetumal bay (Caribbean Sea). The overall shallow conditions in La Sabana suggest that it is sensitive and vulnerable to anthropogenic activities, particularly to contaminant loads. The NSF-WQI did show bad to medium water quality conditions within La Sabana (Figure 2e,f). The center and south showed predominantly bad-quality conditions through most of the year (Table S1). The PCA revealed that total coliforms, BOD₅, low transparency, and fecal coliforms are the most relevant variables of the central south of La Sabana (Figure 5). Fecal and total coliform values were substantially high in this same zone during the dry and rainy seasons, largely exceeding the threshold of Mexican standard NOM-210-SSA1–2014 of <16,000 MPN 100 mL⁻¹ for total and <450 MPN 100 mL⁻¹ for fecal coliforms, which indicates the maximum values for waters for recreational use (Table S1). High values in fecal coliforms are usually associated with the presence of other enteric pathogens [56] and therefore, direct or indirect interaction with the water of the central south of La Sabana must be made with caution as it may cause waterborne diseases such as gastrointestinal diseases, ear, eyes and wound infections, typhoid fever, or dysentery [57]. The area of continuously bad water quality and maximum amounts of bacterial content coincides with the discharge area of the wastewater treatment plant “El Centenario”. Contrastingly, in the wetlands corresponding to the southern basins, where the influence of “El Centenario” is limited, the NSF-WQI shows better water quality (medium quality) and the bacterial concentration drastically dropped in both climatic periods, compared with the central south of the lake. The PCA illustrated that wetlands are associated with waters of higher transparency and lower levels of variables related to contamination, such as fecal coliforms and BOD₅ [58]. The wetlands are the section of the lake where the anthropogenic influence is most limited. Extensive zones surrounding the system still maintain native vegetation and settlements have not been developed yet. Man-made dams contribute to isolating the wetlands from water exchange with the polluted zones of the lake, favoring better water quality.

Wastewater treatment plants discharging to lakes and even coastal areas are known to represent source of contaminants and organic matter [59,60]. In Latin America, one of the most important problems for wastewater management is the lack of infrastructure for wastewater treatment plants; consequently, only about 20% of wastewater is properly treated, and the remaining 80% is partially or not treated at all in urban centers [61,62]. In Chetumal city, “El Centenario” processes about 120 lt s⁻¹ of wastewater and provide services to about 40% of the population of the city, but the treatment capacity is exceeded, and high percentage of wastewaters is only partially treated. Nitrogen and phosphorous are the main constituents of wastewater and promoters of eutrophication in lakes [63,64]. In La Sabana, the nitrogen available in the form of nitrites+ nitrates (NN, Table S2) revealed relatively high concentrations (> 1.6 mg L⁻¹), particularly in the south, and phosphorus in the form of orthophosphates (0.19–0.23 mg L⁻¹, Table S1) is exceeding the limits of international standards such as USEPA (0.05 mg L⁻¹) [65], suggesting that the water conditions in La Sabana are above the maximum acceptable (0.1 mg L⁻¹) and susceptible to rapid eutrophication. Currently, La Sabana is categorized as within eutrophic to hypereutrophic based on the N and P concentrations [25]. Bad water quality and over-enrichment (eutrophication) of La Sabana may, therefore, be a consequence of the combined effects of the insufficient capacity of “El Centenario” to process wastewater, resulting in its release of partially treated water to the lake, and the high percentage of population without access to water treatment services.

Heavy metal concentrations were relatively low in both sediment and fish tissue (Table 2), with almost all values not exceeding the national norms, such as NOM-242-SSA1–2009, and international standard thresholds, such as the CSQG. Most values of the risk index (RI) and the pollution index (PI) of sediment and fish similarly indicate low environmental risk. The RI, however, demonstrated moderate ecological risk for Hg in sediment, and the metal pollution index showed values of severe pollution with Hg in 26% of individual fish. The presence of Hg in aquatic ecosystems is a global concern as it can be toxic even at low concentrations [66,67]. Hg is easily incorporated into trophic chains and can also be potentially toxic to humans through fish consumption [68,69,70]. In La Sabana, fishing is an important activity, with products being locally traded. The consumption of fish from La Sabana must be carefully evaluated, as the results obtained in our study were conducted in the liver, an organ not typically consumed, and that may display relatively higher values of Hg than muscle [71,72]. Concentrations of Hg in muscle and liver, however, can be equaled or even reverted given the overall Hg availability in the lake system [73,74]. For instance, in our study, putative juvenile fish (with smaller length) displayed higher values of Hg, which is indicative that the Hg contamination may be extended to other wildlife and lower levels of the trophic chain, as the juvenile stages of O. niloticus are more specialized in zooplankton, insects and detritus, compared with adults characterized by more herbivorous feeding habits [75,76].

The incorporation of heavy metals in La Sabana was a priori considered to be related to leachate from the open-air landfill of Chetumal city located less than 5 km from the northern basin. In our study, low values or the absence of heavy metals typically generated by landfills, such as Cd, Zn, Pb, Cr [77,78,79], were detected, suggesting that the leaching influence of the landfill was minor. This could be explained by the hydrological dynamic of the zone, as the landfill is located in the drain area of the lake during high water levels. Therefore, in the case of availability, the leachate may be moved toward the Chetumal Bay (Caribbean Sea). The source of the Hg incorporation in La Sabana is relatively difficult to determine, as activities that increase its incidence [80,81], such as metallurgical processes, chlor-alkali industries, battery production and agriculture, are absent in the surrounding areas. Fossil-fuel burning and airborne emissions from activities in landfills such as the compaction of solid waste and burning [82,83,84] may be the most important sources, but it remains to be proved.

The hydrological connectivity in the region, including the groundwater flow, the surficial aquatic environments and their interactions with the Caribbean Sea, suggests that the anthropogenic impact and contaminants of La Sabana may not only be focalized in the core area, but largely extended in the region, affecting a broad range of ecosystems. For instance, both coliform bacteria [85] and heavy metals in sediment and aquatic species [41] have been recorded at a regional scale in the Chetumal bay (Caribbean Sea).

3.2. Ecosystem Effects of Wastewater Discharge and the Nutrient Increase in La Sabana

Lake La Sabana is relevant for hydrological balance in the region and particularly for Chetumal city. Although its waters are not used for the water supply, most meteoric water received by the city is discharged to both La Sabana and Chetumal Bay and, therefore, these systems control flooding. The growth of the city caused the lake to be partitioned into four sub-basins, leading to changes in the dynamic and functionality of the lake. The nMDS and SIMPROF tests revealed the presence of five well-defined groups, corresponding either to basins, or climatic periods. Wetlands were the most different regions in comparison with any other zone of the lake. Climatic variables such as temperature, pH, conductivity, and TSS were more variable between climatic periods, coinciding with what has been observed for the region by different authors [86,87]. The relative isolation of the basins thus focalizes contaminants and drives a different biological structure. The north, center and south of the lake are the zones more drastically altered by anthropogenic activities; partially treated wastewater loads in these zones can be tracked back to at least 10 years ago, during the population growth of the city. The first effect of such anthropogenic influence was the change in the trophic status of the lake. Natural mesotrophic conditions in La Sabana changed in a relatively short period of time (<10 years) to hypereutrophic conditions [25]. In the hypereutrophic basins, the zooplankton community exhibited population changes, likely as a response to increased nutrients. Currently, non-native species dominate the northern basin, as a typical response from highly perturbed aquatic environments, with groups such as Asplanchnidae (Rotifera), Cyclopidae (Copepoda) and Cypridopsis (Ostracoda); whereas in the wetlands, regional distributed groups are more common, e.g., Alicenula, Chlamydotheca (Ostracoda) [25]. The La Sabana fisheries, which were previously based on native species, have changed to only O. niloticus, an introduced and broadly environment-tolerant species, suggesting an additional effect on the fish assemblage.

The coverage of P. stratiotes on La Sabana is an additional and sequential consequence of the over enrichment of lake water. Pistia stratiotes is an invasive species in tropical and subtropical regions of the world that colonizes water bodies rich in N and P [88,89] and with prevailing eutrophic conditions [90,91,92]. Due to its geographic location and climate, the Yucatán Peninsula is susceptible to invasion by aquatic plants, such as P. stratiotes. The region was, however, almost free of this plant species given the overall oligotrophic conditions of most water bodies [93,94,95]. In La Sabana, the colonization of P. stratiotes from 2017 to 2019 ranged from absent to more than 14 hc [22], revealing the high risk to the region of aquatic plant invasion if aquatic bodies become enriched with organic matter originated from anthropogenic activities.

Wastewater production, treatment, and final disposal is a serious concern in the Yucatán Peninsula, because a large percentage of wastewater in the population center and production areas does not undergo any treatment [11,13,14,15,16]. Consequently, pollution easily enters aquatic environments and spreads beyond its source of origin. Given the high connectivity of the aquatic environments on the Peninsula, pollution is currently widespread, affecting all aquatic ecosystems including the marine environment [13,96]. The regional tendency of an increase in nutrient availability, therefore, makes it necessary to understand how aquatic environments respond to pollution and anthropogenic alterations. Peri-urban lakes such as La Sabana are key to understanding such responses, because these are the most drastically anthropogenically altered aquatic environments in the region with accelerated processes of eutrophication [6,25]. In the southern Yucatán Peninsula, the water quality and environmental risk assessment in La Sabana illustrate that wastewater discharge without proper treatment can critically affect water quality in surficial environments and represent a risk for the human population. Moreover, in La Sabana, it becomes evident that increased nutrients and changes in trophic status introduce profound changes in the ecosystem functionality. In La Sabana, previous studies detected the elimination of native fauna by generalist species and the invasion of aquatic macrophytes. At a regional scale, the current status of water quality, nutrient availability, heavy metals and ecosystem alterations in La Sabana must be considered a model of what can be expected in the future in the Yucatán Peninsula, if waste loads to aquatic environments continue or are increased. The environmental restoration of la Sabana is fundamental to maintain the regional ecosystem services and the wellness of society in Chetumal city. It must be conducted from a holistic perspective, addressing the root causes of pollution and with the participation of government, academics and society. An effective restoration program in La Sabana may also benefit aquatic environments with similar environmental conditions.

4. Conclusions

Lake La Sabana is a fundamental system for the ecological, dynamic and hydrological balance of the region, particularly benefiting the management of meteoric water in Chetumal city. The diagnosis of Lake La Sabana highlighted that within its four basins, the environmental conditions are variable, both spatially and temporally, with water quality going from medium in the wetlands, to bad in the northern basin, in accordance with the National Sanitation Foundation water quality index. Heavy metals were relatively low in sediment and fish tissue, except for mercury which showed a high potential environmental risk determined by the ecological indices. The wastewater treatment plant “El Centenario” was recognized as the major source of contamination on La Sabana, because of partially treated wastewater loading to the lake, thus contributing to the increase in pathogenic bacteria in the surficial waters and altering the balance of N/P related to eutrophication. Lake La Sabana must be considered a model of the environmental and ecological consequences of improper wastewater management and anthropogenic alteration in the aquatic environments of the Yucatán Peninsula. Plans for restoration are necessary to reestablish ecosystem processes and maintain the wellness of society in Chetumal city.