Study area

The Grijalva river headwaters are located in the highlands of Guatemala and the river runs through Chiapas and Tabasco in Mexico, before joining the Usumacinta river at Tres Brazos to flow into the Gulf of Mexico (Fig. 1). For this study, we considered the portion of the Grijalva basin in Mexican territory that covers 51 943.95 km² (^{INEGI 2010}). The Grijalva basin is part of the Hydrologic Region No. 30 and is divided into three subregions: Lower (LG), Middle (MG), and Upper Grijalva (UG) (^{SEMARNAT 2010}), and 50 sub-basins (12 in UG, 12 in MG, and 26 in LG) (^{INEGI 2010}). Main land use and cover in the basin are grassland (29.92 %), agriculture (23.96 %), forest (21.78 %), moist forest (16.19 %), marsh (3.74 %), and urban (1.16 %) (^{INEGI 2016}; Table I).

Fig. 1 The Grijalva basin in Mexico: (a) location of the Grijalva basin, (b) the four hydropower dams into the Grijalva basin (green circles) and (c) study monitoring stations, main rivers and subregions. 

The basin has a population of 5 382 520 inhabitants and an average population density of 104 inhabitants/km² according to the 2020 national population census (^{INEGI 2021}). The largest cities are Tuxtla Gutiérrez (578 830 inhabitants) and San Cristóbal de las Casas (183 509 inhabitants) in Chiapas, and Villahermosa (340 060 inhabitants) in Tabasco (^{INEGI 2021}). Municipal wastewater is a major source of pollution in the Grijalva basin (^{García-García 2005}, ^{INE 2007}, ^{CONAGUA 2019a}), increasing levels of nutrients, fecal coliforms, total suspended solids, biochemical, and chemical oxygen demands (^{SNIARN-SEMARNAT 2016}, ^{CONAGUA 2018a}).

General approach

We examined nutrient concentrations among subregions, stream size, and land uses to determine the need to develop criteria for different areas of the basin as recommended by the ^{USEPA (2000)}. We applied two methods based on frequency distribution to estimate potential nutrient criteria for TP and TN to the whole basin and its subregions, using water quality monitoring data from the Comisión Nacional del Agua (National Water Commission; CONAGUA). We estimated the 75th percentile of the frequency distribution of the reference streams, and percentiles from 5th to 25th of the general stream population. We applied both methods for each subregion and for the whole basin, considering two stream size categories and assuming that environmental conditions were not similar among them. All analyses were performed in RStudio v. 1.3.1093.

Water quality data

Water quality data (physical and chemical parameters, and nutrient concentrations) from 137 monitoring stations in rivers of the Grijalva basin were obtained from CONAGUA. Data were available from January 2012 to December 2018 and had a monthly frequency at some stations, and bimonthly or less at others. Most of the stations included data from periods of low and high flows. For each monitoring station, the subregion, sub-basin, and stream order (1-7) were identified from a hydrographic map (^{INEGI 2010}; Table II). For the analysis, we did not include monitoring stations in lentic systems or near the coast with evidence of saltwater intrusion.

Database contained alphanumeric data for some water quality records such as < LD (lower than detected) and < CMC (lower than minimum measurable), for which we used values equal to 50 % of the reported detection limit (^{Suplee et al. 2007}, ^{Hites 2019}). Records reported as ND (non-determined) and NE (non-estimated) were removed because they represented less than 2 % of whole database (^{Hites 2019}). We identified and assessed suspected outliers by monitoring station to avoid erroneous criterion from atypical data. We did not remove outliers, except in cases in which concentrations were atypically higher (^{USEPA 2010}, ^{Helsel et al. 2020}). After data examination, the database had 125 monitoring stations with 3814 records for TP and 3918 records TN.

Nutrient comparisons among stream size and subregions

We grouped streams in two size categories: low-order (1, 2, and 3) and high-order streams (4, 5, 6, and 7). We compared TP and TN concentrations among subregions and between size categories by using ADONIS (permutation ANOVA), a non-parametric permutation test in the vegan package of R (^{Oksanen 2015}). When interactions were significant, we applied pairwise comparisons using permutation F test to compare the effects of stream size categories into each subregion. We used ADONIS, because data did not meet the normality assumption of ANOVA (^{Oksanen 2015}).

Land use data

We determined land use and land cover at the reach scale from ^{INEGI (2016)}. We used the following land use categories: agriculture, forest, mangrove, grassland, marsh, and urban. Dominant land use at the reach scale for each sampling station was determined by establishing a 2 km buffer around the monitoring station. To explore the relationship between land use and nutrient concentrations, we compared concentrations among the most dominant land use categories (agriculture, forest, grassland, marsh, and urban) at the reach scale and among subregions by applying ADONIS. When the interaction was significant, we conducted pairwise comparisons using permutation F test to compare treatment effects.

Reference streams approach

We identified potential reference sites by subregion through two methods: the best professional judgment (BPJ) and the presence of forest cover as dominant land use. BPJ is used when minimally impacted sites are difficult to find (^{Cunha et al. 2011}), and forested sites usually have optimal conditions to be considered as reference sites (^{USEPA 2000}). By applying BPJ, we selected the less impacted streams based on reach land use and water quality (physical and chemical parameters, nutrients concentrations, biochemical oxygen demand, and faecal coliforms) (^{USEPA 2000}, ²⁰¹⁰, ^{Cunha et al. 2011}). We also selected streams with forest cover at the reach scale to examine its potential use as reference streams. Reference streams for the whole basin were estimated pooling the reference streams selected for each subregion. We calculated the 75th percentile of the frequency distribution of both reference streams (BPJ and forested reference streams) for TP and TN, for each subregion, and for the whole basin.

General stream population approach

We used the whole population of streams and rivers for calculating the concentrations at the 5th and 25th percentiles for TP and TN, for each stream size category, subregion, and for the whole basin. Furthermore, we applied the trisection method (1/6th or 16.7th percentile) also using the whole dataset. This method is recommended when working in impacted regions, where pristine or reference sites are difficult to achieve (^{Cunha et al. 2011}, ^{Huo et al. 2018}).

Nutrient concentrations

TP concentrations ranged from 0.01 to 12.18 mg/L. Higher concentrations were observed in the MG subregion (Alto Grijalva, Cintalapa, De la Venta, Santo Domingo, Suchiapa, and Tuxtla Gutiérrez sub-basins; p < 0.005). There were significant differences in TP concentrations between low-order and high-order streams in LG and MG subregions (p = 0.001; Table III). Greater concentrations were detected in low-order streams in the MG subregion. Low-order streams showed significant differences among subregions, with higher concentrations in the MG subregion than in the LG and UG subregions (p = 0.001). High-order streams also revealed higher concentrations in the MG subregion and lower concentrations in the LG subregion (p = 0.001).

TN concentrations varied from 0.02 to 71.05 mg/L. Higher concentrations were found in low-order streams and in the MG subregion (Alto Grijalva and Tuxtla Gutiérrez sub-basins; p < 0.005). For TN concentrations, there were also significant differences between low-order and high-order streams in the LG and MG subregions (p = 0.001; Table III). For low-order streams, concentrations were higher in the MG subregion than in the LG and UG subregions (p = 0.001). For high-order streams, concentrations were higher in the MG subregion and lower in the LG subregion (p = 0.001).

Land uses

The most predominant land use in the Grijalva basin was grassland (29.92 %). Forest and agriculture were dominant in the MG (29 and 24.20 %, respectively) and UG subregions (36.09 and 30.86 %, respectively), while grassland was in the LG subregion (45.51 %; Table I). In the LG subregion, forested reaches (0.08 ± 0.15 mg/L) had lower TP concentrations than the other land uses (p < 0.012; Table IV); grassland (0.16 ± 0.26 mg/L) and marsh (0.16 ± 0.24 mg/L) showed lower concentrations than agriculture (0.21 ± 0.41 mg/L) and urban (0.25 ± 0.46 mg/L; p < 0.02); and agriculture had lower concentrations than urban (p < 0.036). In the MG subregion, the TP highest concentration was in urban areas (0.95 ± 1.68 mg/L; p < 0.016). In the UG subregion, the TP concentrations were higher in grassland (0.52 ± 0.91 mg/L; p < 0.04).

In the LG subregion, TN concentrations in forest (1.10 ± 0.67 mg/L), grassland (1.13 ± 0.80 mg/L) and marsh (1.15 ± 0.78 mg/L) reaches were lower than in the other land uses (p < 0.04), while in agriculture (1.53 ± 2.34 mg/L) they were higher than in urban (1.42 ± 1.28 mg/L; p < 0.024). In the MG subregion, urban concentrations (6.04 ± 10.82 mg/L) were significantly higher than in other land uses (p = 0.004). In the UG subregion, the TN highest concentration was in grassland (1.79 ± 1.92 mg/L; p = 0.004).

Reference streams

From the 125 monitoring stations in the Grijalva basin, we selected nine reference streams based on BPJ (two in UG, two in MG, and five in LG subregion) and five forested reference streams (one in UG, one in MG, and three in the LG subregion). Only two forested reference streams were included in the BPJ group for the LG subregion. Forested reference streams for the LG subregion were found in Chiapas.

TP and TN concentrations from BPJ reference streams were lower than concentrations obtained for each land use category at the reach scale (agriculture, forest, grassland, marsh, and urban; Table IV). The highest TP concentration was in the MG subregion for both BPJ (0.18 ± 0.31 mg/L) and forested reference streams (0.30 ± 0.68 mg/L; Table V). The highest TN concentration (1.48 ± 0.97 mg/L) was found in the MG subregion in forested reference streams.

TP concentrations at the 75th percentile were higher in BPJ reference streams than in forested reference streams in the LG and UG subregions (0.13 and 0.18 mg/L, respectively; Table VI), while greater values occurred in forested reference streams in the MG subregion (0.25 mg/L). TN concentrations were higher in forested reference streams than in BPJ sites in each subregion. TP and TN concentrations at the 75th percentile for both BPJ (0.13-0.19 and 1.16-1.31 mg/L for TP and TN, respectively) and forested streams (0.09-0.25 and 1.35-1.58 mg/L for TP and TN, respectively; Table VI), were higher than median concentrations of the general stream population (0.08-0.18 and 0.90-1.44 mg/L for TP and TN, respectively; Table III).

General stream population

TP concentrations at 5th and 16.7th percentiles were similar between low-order and high-order streams, while concentrations at the 25th percentile were higher in low-order streams in the MG subregion (Table VII). In the combined stream size category, greater concentrations were observed in the MG subregion and lower concentrations in the UG subregion. At the basin scale, 107 monitoring stations (85.6 %) had values below the 25th percentile (Fig. 2). When analyzed by subregions, 63 stations (91.3 %) in the LG subregion had concentrations below this percentile while six stations (two in Carrizal, one in Cunduacán, one in Cuxcuchapa, one in Pichucalco, and one in Samaria sub-basins) always exceeded (Fig. 3). In the MG subregion, 36 stations (90 %) had data below the 25th percentile while concentrations in four stations (two in Alto Grijalva and two in Tuxtla Gutiérrez sub-basins) always exceeded (Fig. 3). In the UG subregion, all 16 stations showed concentrations below the 25th percentile. When we analyzed concentrations at the 5th percentile at the basin scale, only 52 monitoring stations (41.6 %) had data below this percentile, while by subregions 12 stations in UG, 21 in MG, and two in LG had concentrations below the 5th percentile.

Fig. 2 Nutrient concentrations (mg/L) for the Grijalva basin using 25th percentile approach. Green circles represent monitoring sites that show nutrient concentration values in the 2012-2018 period, equal or lower than the limit concentration for (a) total phosphorus (TP) and (b) total nitrogen (TN). Red circles represent sites with higher values than the limit concentration. 

Similar TN concentrations between stream size categories were obtained for the 5th, 16.7th, and 25th percentiles at LG subregion. Higher values were obtained in low-order streams at the MG subregion (Table VII). When the two stream size categories were combined, greater concentrations were in the MG and UG subregions at the 16.7th and 25th percentiles, while higher concentrations at the 5th percentile were in the MG subregion. At the basin scale, 113 monitoring stations (90.4 %) exhibited values below the 25th percentile (Fig. 2). In the analysis by subregions, 63 stations (91.3 %) in the LG subregion had concentrations below this percentile while concentrations in six stations (one in Almandro, one in Carrizal, one in Chacté one in Cunduacán, one in Cuxcuchapa, and one in Samaria sub-basins) always exceeded (Fig. 3). In the MG subregion, 34 stations (85 %) had data below the 25th percentile while all concentrations at six stations (two located in Alto Grijalva, one in Suchiapa, and three in Tuxtla Gutiérrez sub-basins) were above this percentile (Fig. 3). In the UG subregion, all 16 stations showed values below this percentile. When we analyzed concentrations at 5th percentile at the basin scale, there were 76 monitoring stations (61 %) that had values below this percentile, while by subregions, 14 stations in UG, 23 stations in MG, and 38 in LG had concentrations below this percentile.

Fig. 3 Nutrient concentrations (mg/L) for each subregion of the Grijalva basin using the 25th percentile approach. Green circles represent monitoring sites that show nutrient concentration values in the 2012-2018 period, equal or lower than the limit concentration for (a) total phosphorus (TP) and (b) total nitrogen (TN). Red circles represent sites with higher values than the limit concentration. 

Nutrient concentrations in the Grijalva basin

Nutrient concentrations in the Grijalva basin showed differences between stream size, land use, and subregions. Concentrations were greater for TP and TN in low-order streams, supporting other research findings that have documented a decline in nutrient concentration as stream size increases (^{Alexander et al. 2000}, ^{Smith et al. 2003}, ^{Binkley et al. 2004}, ^{McDowell et al. 2017}). This pattern can be related to greater nutrient inputs from groundwater in small streams, particularly to those located in agricultural areas and to the cumulative effect of nutrient removal processes along the fluvial network (^{Smith et al. 2003}, ^{Kang et al. 2008}). In addition, tropical forested streams can show high nitrate concentrations because of high nitrification rates in soils (^{Neill et al. 1997}, ^{Espinoza-Toledo et al. 2021}) that can also explain higher TN concentrations in low-order streams of the Grijalva basin.

In forested reaches, nutrient concentrations were lower than in other land uses, as expected due to lower inputs from anthropogenic activities (^{Poor and McDonnell 2007}, ^{Cunha et al. 2011}, ^{de Mello et al. 2020}). Forest vegetation improves water quality through filtration, absorption, and interception of runoff nutrients (^{Li et al. 2016}, ^{Duffy et al. 2020}, ^{de Mello et al. 2020}). Agriculture and grassland areas showed greater concentrations than forest areas, likely explained by the application of fertilizers and inputs of livestock waste (^{Aryal et al. 2012}, ^{Jones et al. 2019}).

Higher nutrient concentrations were found in urban areas, probably related to domestic and industrial wastewater inputs. Detergents containing phosphates are still used in Mexico, increasing phosphorous (P) levels in wastewater discharges (^{SNIARN-SEMARNAT 2013}). In addition, in Chiapas and Tabasco only 30.4 and 36.3 %, respectively, of domestic wastewater are treated essentially by primary treatment because secondary treatment is scarce (^{CONAGUA 2018b}, ^2019a), indicating that streams and rivers receive a large proportion of untreated wastewater (^{INE 2007}). Nutrients levels in untreated wastewater can be as high as 4-12 mg/L for TP and 20-70 mg/L for TN, and concentrations in effluents from secondary wastewater treatments plants (WWTP) using activated sludge can be 4-10 mg/L for TP and 15-35 mg/L for TN (^{Carey and Migliaccio 2009}). Low proportion of treated wastewater or inefficient treatment processes can increase nutrient discharge into streams and rivers (^{INE 2007}). In addition, the expansion of urban areas can increase riparian deforestation, decrease infiltration rates, and modify stream and river channels, affecting nutrient dynamics (^{Carey and Migliaccio 2009}, ^{Carpenter et al. 2011}, ^{McMillan et al. 2014}).

The MG subregion showed greater mean TP concentrations, likely related to the presence of large cities such as Tuxtla Gutiérrez and San Cristóbal de Las Casas (^{INEGI 2021}). High TP concentrations in MG (Table III) were similar to concentrations in untreated wastewater (^{Rivera-Vázquez et al. 2007}, ^{Carey and Migliaccio 2009}). The median concentrations in UG and LG subregions (Table III) were in the range of those estimated by ^{Smith et al. (2003)} for US streams and rivers and similar to values reported for rivers draining urban and agricultural catchments in Mexico (^{Jujnovsky et al. 2010}, ^{Arellano-Aguilar et al. 2017}, ^{Fregoso-López et al. 2020}). We expected greater TP concentrations at the MG subregion because of urban wastewater inputs, but median TP concentrations at the LG subregion were unexpectedly lower, despite the greater population density and the presence of a large city (Villahermosa). The location of hydroelectric dams upstream (Fig. 1) and the greater number of WWTP in LG can be related to this result. Dams can retain sediments and nutrients, particularly P (^{van Cappellen and Maavara 2016}); thus, sites located downstream from the dams can be affected by this retention. The LG subregion has more WWTP than MG (70 in LG and seven in MG) (^{CONAGUA 2018b}) that can contribute to lower P inputs into the rivers. ^{Lázaro-Vázquez et al. (2018)} detected an increase in nitrate concentration between 2000 and 2016 in LG but not in TP, suggesting that some removal processes, probably by the dams and the wastewater treatment, prevent an increase in P concentrations. In the UG subregion, higher TP concentrations were found in grassland areas, probably related to fertilizer use, soil compaction, and livestock waste (^{Capece et al. 2007}, ^{Aryal et al. 2012}).

TN showed very high concentrations at the MG subregion in urban reaches, similar to values reported for untreated wastewater (^{Andersen et al. 2004}, ^{Carey and Migliaccio 2009}). This subregion also had intensive agricultural areas, which can provide N inputs from fertilizer application (^{Blann et al. 2009}, ^{Cunha et al. 2011}, ^{de Mello et al. 2020}). In the LG subregion, high TN concentrations were found in agricultural areas, supporting the association between this land use and N inputs (^{Kang et al. 2008}, ^{Billy et al. 2013}). In the UG subregion, in addition to agricultural areas, cattle production probably explains greater TN concentrations due to increased inputs of N from livestock waste (^{Agouridis et al. 2005}, ^{Jones et al. 2019}).

Frequency distribution methods for nutrient criteria

The use of reference streams for establishing nutrient criteria assumes that this group of rivers represents ecosystems with less human disturbance in the basin (^{USEPA 2000}). Average concentrations obtained from BPJ and forested reference streams were lower (Table V) than in agriculture and urban areas (Table IV) and in the general stream population (Table III). However, concentrations at the 75th percentile of the reference streams (Table VI) were higher than values estimated by the 5th and 25th percentiles of the general stream population (Table VII). ^{USEPA (2000)} suggests that concentrations at the 75th percentile of the reference streams and at the 25th percentile of the general stream population can be similar, but several studies have found greater concentrations at the 75th percentile than at the 25th percentile (^{Suplee et al. 2007}, ^{Herlihy and Sifneos 2008}, ^{Evans-White et al. 2013}; Table VIII). For US ecoregions, ^{Suplee et al. (2007)} report that concentrations at the 75th percentile corresponded to percentiles from 4th to 97th of the general stream population, indicating that for some ecoregions, percentiles are not similar. For the Grijalva basin, nutrient concentrations at the 75th percentile of the reference streams corresponded to percentiles from 59th to 64th of the general stream population, suggesting that both distributions are relative similar and probably disturbed streams were included in the reference group (^{Smith et al. 2003}). Nutrient concentrations at the 75th percentile of the reference streams were higher than reported in some nutrient criteria studies, but were also in the range of others (^{Smith et al. 2003}, ^{Suplee et al. 2007}, ^{Hsieh et al. 2016}; Table VIII). Nutrient concentrations at the 75th percentile of forested reference streams were also higher than those reported by ^{Espinoza-Toledo et al. (2021)} for forested streams in LG (0.033 and 0.303 mg/L for TP and TN, respectively).

Human pressures in the Grijalva basin-the location of the monitoring sites-and the low number of reference streams could influence the high concentrations obtained in the selected reference streams. In basins with a wide range of human disturbances such as the Grijalva, reference sites selected based on least-disturbed conditions will not represent minimally disturbed conditions (^{Stoddard et al. 2006}) and rather reflect the consequences of human activities. In this study, reference streams probably included impacted streams, especially in the MG subregion, due to the unavailability of monitoring sites located in rivers with low human impacts. This potentially explains higher concentrations at the 75th percentile of the reference streams than at the 25th of the general population (^{Evans-White et al. 2013}). In addition, estimates based on a small number of reference sites could be affected by natural variability in topography, geology, vegetation, soils, and atmospheric deposition (^{Smith et al. 2003}, ^{Suplee et al. 2007}), and do not represent realistic values of the reference conditions (^{Herlihy and Sifneos 2008}). Although it will be difficult to find undisturbed large rivers in the Grijalva, the identification of small-forested catchments to locate additional monitoring stations will benefit the collection of long-term water quality data.

An alternative for the 75th percentile of reference streams could be the use of the percentiles 5th to 25th of the general population data (^{USEPA 2000}, ^{Rohm et al. 2002}, ^{Suplee et al. 2007}, ^{Herlihy and Sifneos 2008}, ^{Hsieh et al. 2016}). Concentrations at the 25th percentile for the Grijalva basin were 0.03-0.08 mg/L for TP and 0.62-0.83 mg/L for TN, which exceeded values estimated in some studies but were in the range of others (Table VIII). Concentrations at the 5th and 16.7th were similar to values at the 25th percentile of several studies. However, these comparisons must consider that previous studies on nutrient criteria correspond to temperate and sub-tropical systems, and background nutrient concentrations can vary among regions. Factors such as greater N concentrations in tropical than in temperate forested systems (^{Brookshire et al. 2012}) and P as a limiting nutrient in many streams of the region (^{Capps and Flecker 2013}) can influence nutrient criteria values for the Grijalva.

The use of the general stream population could be an alternative for the 75th percentile but it has some disadvantages. The general nutrient data is affected by the degree of disturbance of the study region, and criteria drawn from these data can fluctuate over time, depending on the increase or decrease of nutrient inputs (^{Suplee et al. 2007}, ^{Herlihy and Sifneos 2008}). In areas under high human pressures such as the Grijalva basin, the general stream population will reflect the eutrophication degree, and concentrations at the 25th percentile probably will not be enough to protect aquatic ecosystems (^{Suplee et al. 2007}). On the other hand, the application of concentrations based on the 5th percentile can be very restrictive and difficult to achieve in disturbed areas (^{Hsieh et al. 2016}, ^{Huo et al. 2018}), but some authors suggest that this percentile could represent better background conditions (^{Clune et al. 2020}). Thus, the use of concentrations at the 16.7th percentile could protect the aquatic environments without being too restrictive.

Based on our results, low-order and high-order streams showed similar concentrations at each percentile of the general stream population, suggesting that a single criterion could be used for all stream sizes. Differences among subregions were observed in all percentiles and could be related to differences in natural and human processes among subregions (^{Smith et al. 2003}, ^{Suplee et al. 2007}, ^{Herlihy and Sifneos 2008}). Although it is recommended to establish nutrient criteria for regions with similar characteristics (^{USEPA 2000}), the use of a single criteria for the whole Grijalva basin will facilitate its implementation and the protection of aquatic ecosystems. The application of greater values for upstream (in the UG and MG subregions) than for downstream areas (the LG subregion; Table VII) can compromise nutrient goals in LG. Thus, considering the upstream-downstream connection and human influences in the basin, the application of lower nutrient limits at the upstream subregions based on estimates for the whole basin, probably will benefit aquatic ecosystems.

The use of percentiles can represent an initial guide to develop nutrient criteria in the Grijalva basin but needs to be validated with other methods (^{Herlihy and Sifneos 2008}, ^{Clune et al. 2020}). For example, relationships between nutrients levels and stream biota (algae, invertebrate, fishes) can help to identify thresholds where small increases in nutrient levels cause a large change in the biotic component (^{Robertson et al. 2006}, ^{Smith and Tran 2010}). Models based on basin characteristics or the relationship between land use and nutrient can help to establish reference conditions in areas where undisturbed streams are not available (^{Smith et al. 2003}, ^{Dodds and Oakes 2004}).

Strategies for nutrient control

Implementing and achieving nutrient criteria involve the application of mitigation measures to reduce nutrient discharges into streams and rivers. In Mexican regulations, maximum limit concentrations (monthly averages of 15 mg/L for TP and 25 mg/L for TN) in wastewater discharges into rivers (^{SEMARNAT 2021}) are higher than those established in other countries for TP and TN: 0.005 and 0.1 mg/L, respectively, in Japan (^{MOE 2015}); 1-2 and 10-15 mg/L, respectively, in the EU (^{CNW 2018}); 0.1 and 2-6 mg/L, respectively, in the USA (^{Litke 1999}, ^{USEPA 2013}); 0.5 and 15 mg/L, respectively, in China (^{Li et al. 2012}). Nutrient limits in Mexico are also higher than the TP limit concentration (0.03 mg/L) estimated to avoid eutrophication (^{Dodds 2002}). Values found in this study for the 5th, 16.7th, and 25th percentiles are also lower than current limits set by Mexican regulations, but the 16.7th percentile is not too restrictive and could help to protect aquatic ecosystems in the basin.

In other countries, strategies such as improving WWTP, implementation of maximum daily discharges, optimal manure and fertilizer use, and conservation and restoration of riparian zones has been successfully applied (^{Carey et al. 2013}, ^{Jarvie et al. 2013}, ^{Heatherly 2014}). In the Grijalva basin, sub-basins that always exceeded the 25th percentile were associated with urban areas (Alto Grijalva and Tuxtla Gutiérrez in MG, and Cuxcuchapa, Carrizal and Samaria in LG). Thus, improvement in WWTP is required to reduce nutrient discharges from urban areas. Wastewater treatment in Tabasco and Chiapas is below the national average of 64 % (^{CONAGUA 2019a}), despite an increase of WWTP in Chiapas reported in 2019 ^{(CONAGUA 2019b}). Increasing the number of WWTP and enhancing the treatment process (more secondary and tertiary where possible) (^{Kronvang et al. 2008}, ^{Neal et al. 2010}, ^{Räsänen et al. 2017}) will reduce the amount of wastewater that is directly discharged into the Grijalva streams and rivers. Other options include the implementation of constructed wetlands (^{Dell’Osbel et al. 2020}, ^{Parde et al. 2021}). In addition, reducing the use of detergents containing phosphates will decrease P inputs in urban wastewater (^{Carey et al. 2013}). Decreased nutrient loading from non-point sources can be achieved by the conservation and reforestation of riparian zones, wetland restoration, crop rotation, and setting maximum limits for the application of chemical fertilizers, among other measures (^{Kronvang et al. 2008}, ^{Riemann et al. 2016}, ^{Tromboni and Dodds 2017}).