Impact of Changes of Land Use on Water Quality , from Tropical Forest to Anthropogenic Occupation : A Multivariate Approach

Worldwide, it is acknowledged that changes of land use influence water quality; however, in tropical forests, the relationship between land use and water quality is still poorly understood. This study assessed spatial and seasonal variations in water quality, and the relationship between water quality and changes of land use in the Bobos-Nautla River, whose upper course runs across a patch of a tropical cloud forest. Spatial and seasonal variations in water quality and land use were assessed with multivariate tools. A cluster analysis, as well as a Principal Component Analysis (PCA-3D), identified three groups of sites: (1) an upper portion, which showed the best water quality and the broadest natural vegetation coverage; (2) a middle course, with high nitrogen and phosphorus concentrations associated with extensive agricultural uses; and (3) a lower course, characterized by the highest levels of total and fecal coliforms, as well as ammonia nitrogen, associated with the highest percentage of urbanization and human settlements. Our findings demonstrate the impact of changes of land use on water quality of rivers running through cloud forests in tropical zones, which are currently endangered ecosystems.


Introduction
It is recognized that climate change is a process that evolves slowly, reflecting a set of modifications that are evident in the long term.However, changes of land use have short-term, frequently drastic effects, leading to considerable impacts on environmental and hydrological processes (infiltration, groundwater recharge, base flow, runoff, water quality) [1].Understanding the relationship between land use and freshwater quality is key for effective water management, both currently and in the future [2,3].
The quality of freshwater reflects the combined effects of multiple natural processes and the anthropogenic changes of land use [4,5].Runoff from catchment areas into water bodies is the main source of nutrients and pollutants [6]; consequently, the relationship between water quality in lotic environments and changes of land use has raised a growing interest within the scientific community, as landscape patterns regulate physical and biogeochemical processes in the basin [7].Agriculture and urbanization are the main sources of nutrients and xenobiotics that degrade water quality in water bodies [8].The conversion of native vegetation to agriculture and human settlements has significantly increased the human well-being at the expense of the degradation of many ecosystem services and biodiversity [9,10].Particularly, the annual deforestation rate has increased dramatically in tropical zones [11,12], and although the deforestation rate showed a 50% decrease in 2015 relative to 1990 worldwide, in Mexico the loss of forest areas still prevails [12].According to the Secretariat of Environment and Natural Resources (SEMARNAT), 29% of the original natural vegetation coverage in the Mexican territory had been lost by 2011 [13], with mountain cloud forests among the ecosystems most vulnerable to changes of land use [14].In developing countries, the combination of a vigorous population growth rate, industrialization, and expanding agricultural areas has exerted and ever-increasing pressure on water [15] and forest resources, impacting ecological processes.In Mexico, cloud forests are seriously threatened by deforestation and habitat fragmentation [16].For this reason, environmental policies and management plans for the conservation of water and forest resources should be developed.
Cloud forests cover a mere 1% of the Mexican territory, ranking second among the ecosystems with the highest species richness in Mesoamerica [17].This forest is characterized by several arboreal strata with a high abundance of ferns and epiphytes, as well as a high biodiversity associated with its unique climatic conditions: high cloudiness, fog at the canopy level, and high relative humidity [18,19].However, models on the effects of changes of land use on climate change have predicted an estimated loss of up to 70% of cloud forest area by 2080 [20].Changes of land use in ecosystems as vulnerable as cloud forests affect biodiversity and ecosystem services, including the hydrological cycle [21].The influence of changes of land use on water quality of rivers running through these forests is currently poorly understood, particularly in areas with a mosaic of land uses along the river course.Some research [22,23] has demonstrated that cloud forest vegetation can function as a natural filter against pollutants dumped in watersheds; hence, vegetation cover, diversity, and structure are key drivers of water quality.Although changes of land use associated with loss of water quality have been reported for several countries in various works [22][23][24], the information on this issue in Mexico is scarce.First-order streams with different land uses including cloud forest, grassland, and coffee plantations, which were related to eight ions in water (Ca 2+ , CO 3 2− , NO 3 − , Cl − , Na + , Mg 2+ , and K + ) and total suspended solids, and to the structure of vegetation, were analyzed [25], indicating that the ions studied increased in grassland and coffee plantations but were lower in forest cloud forest, and a loss of vegetation diversity was detected.In this sense, there is a large number of tools for the detection of the effects of climate change and land uses on water quality [1,3].Among them, modeling has proven to be a reliable tool for predicting the long-term effects of climate change on water quality.Some models are based primarily on spatial analysis using Geographic Information Systems (GIS) through satellite imagery and digital elevation that incorporate various watershed attributes in addition to land use areas and water quality parameters (slope, precipitation, air temperature, topography and soil parameters, river flow, elevation, and anthropogenic factors influencing land use, among others) [1,3,[26][27][28]; other models include time series, linear models, multiple linear regression, and constrained least squares models [26,29].Models are also used to explore the behavior of pollutant levels in stream water under different scenarios.The current heavy deforestation and high rate of land use change in tropical forests cause short-term effects in the water quality of rivers flowing across them; therefore, studies should be conducted to identify the main factors affecting water quality in these areas.
Water quality is controlled by multiple factors, including the connectivity and interaction of the river with the watershed.Multivariate statistical methods (cluster analysis, CA; principal component analysis, PCA) are key methods to evaluate the correlations among the various water quality parameters and land uses to identify the key drivers affecting water quality and influencing the aquatic ecosystem.CA and PCA have been extensively used to characterize and evaluate freshwater quality, revealing temporal and spatial variations caused by natural and anthropogenic factors [30][31][32][33].CA is an unsupervised pattern-recognition technique that reveals the intrinsic structure or underlying behavior of a data set without making a priori assumptions about the data, in order to sort the components of a system into categories or clusters according to their similarity [34].In the clustering procedure, study sites (components) are sorted in such a way that similar components are clustered into the same class sharing a particular set of properties [35].Thus, cluster analysis allows interpreting the data and identifying patterns [34].For its part, PCA is a powerful tool for pattern recognition and detection of correlations (or co-variances) among a set of variables or samples.It provides information on the most significant parameters associated with spatial and temporal variations; it also describes whole data sets by excluding the least significant parameters with minimum loss of the original information [36].It is extremely useful for the analysis of data involving a large number of variables.It has been extensively used as an unbiased method that reveals associations between samples and variables [37].Furthermore, the potential capabilities of GIS combined with multivariate statistical analyses for assessing water quality have been explored [38].Thus, the combined use of GIS with CA and PCA facilitates the interpretation of complex databases to gain a deeper insight on water quality in watershed studies, being key in the development of appropriate strategies for effective management of water resources [34].Although several studies in temperate latitudes use CA/PCA to assess the relationship between changes in land use and their impact on water quality, as is the case of [34], this is the first study that uses both tools (CA/PCA) in combination with GIS to explore the effect of changes of land use in a tropical cloud forest on water quality of a stream running through it in Mexico.
This study aims to investigate the temporal and spatial variations in water quality as related to land uses along the Bobos River, Veracruz, a hydrological system that flows through an important cloud forest area in the Gulf of Mexico slope.In order to analyze land uses and water quality indicators, we used CA and PCA to detect patterns (clusters of study sites) and to identify the key parameters (land uses and water quality indicators) associated with seasonal and spatial patterns.

Study Area
The Bobos River is located in the northeastern portion of the Gulf of Mexico slope, within the Nautla River basin; a section of the main river channel flows across the natural protected area named "Filobobos River and its Surroundings" (Río Filobobos y su Entorno; Figure 1).The basin is used for multiple activities, including mining in the upper course, the production of banana, citrus, and coffee in the middle course, and fisheries, humid-soil agriculture, induced pastures, and human settlements in the lowlands.The region is also used for tourism activities such as rafting.The largest town, Martínez de la Torre (>250,000 inhabitants), is located in the lowlands, a region that also includes less densely populated towns (<250,000 inhabitants) such as Altotonga, Jalacingo, and Tlapacoyan [39].

Field Work
This study involved 12 study sites located along the Bobos and Nautla rivers.In the upper course, the sites selected were Pimiento (PI) and Huapala (HU).Along the Alseseca river tributary, the sites selected were Jalacingo (JL), Tomata (TM), and Encanto (EN), located close to Tlapacoyan and within the natural protected area.Last, the sites Tezcapa (TZ), Manantial (MA), Puente Filo (FI), Palmillas (PA), Rojo Gómez (RG), Martínez de la Torre (MZ), and Paso Largo (PL) are located along the mainstream of Bobos-Nautla rivers, with RG, MZ, and PL being close to the largest town in the region (Martinez de la Torre).The study involved six monitoring events that encompassed the dry, wet, and cold seasons, the latter influenced by northerly winds (cold fronts): August (rainy) and December (cold) 2013, February (cold), June, and September (rainy) 2014, and September (rainy) 2015.In addition, the altitudinal profile of each stream in the basin was delineated, locating study sites according to altitude and distance traveled, in kilometers, from the highest to the lowest study site; a topographic map was included.
The following variables were recorded at each site and study season: dissolved oxygen (mg/L), water temperature (°C), conductivity (μs/cm), turbidity (NTU), and pH (using a Quanta ® probe, Hydrolab, Austin, TX, USA).In addition, water samples (100 mL and 500 mL) were collected from each site for microbiological and chemical testing in the laboratory.

Field Work
This study involved 12 study sites located along the Bobos and Nautla rivers.In the upper course, the sites selected were Pimiento (PI) and Huapala (HU).Along the Alseseca river tributary, the sites selected were Jalacingo (JL), Tomata (TM), and Encanto (EN), located close to Tlapacoyan and within the natural protected area.Last, the sites Tezcapa (TZ), Manantial (MA), Puente Filo (FI), Palmillas (PA), Rojo Gómez (RG), Martínez de la Torre (MZ), and Paso Largo (PL) are located along the mainstream of Bobos-Nautla rivers, with RG, MZ, and PL being close to the largest town in the region (Martinez de la Torre).The study involved six monitoring events that encompassed the dry, wet, and cold seasons, the latter influenced by northerly winds (cold fronts): August (rainy) and December (cold) 2013, February (cold), June, and September (rainy) 2014, and September (rainy) 2015.In addition, the altitudinal profile of each stream in the basin was delineated, locating study sites according to altitude and distance traveled, in kilometers, from the highest to the lowest study site; a topographic map was included.
The following variables were recorded at each site and study season: dissolved oxygen (mg/L), water temperature ( • C), conductivity (µs/cm), turbidity (NTU), and pH (using a Quanta ® probe, Hydrolab, Austin, TX, USA).In addition, water samples (100 mL and 500 mL) were collected from each site for microbiological and chemical testing in the laboratory.

Water Quality
Water samples were tested for nitrates (mg/L), nitrites (mg/L), ammonia nitrogen (mg/L), total nitrogen (mg/L), orthophosphates (mg/L), total phosphorus (mg/L), color (Pt-Co U), and hardness, using the spectrophotometric techniques of the Hach ® DR 3900 instrument (Hach, Loveland, CO, USA); on the other hand, chlorides (mg/L), BOD 5 (mg/L), alkalinity (mg/L CaCO 3 ), and total and fecal coliforms (MPN/100 mL) were determined using American Public Health Association (APHA) techniques [40].The water quality index (WQI) was calculated from the multiplicative weighted index proposed by Dinius [41], ranging from 0 to 100, where zero corresponds to extremely poor water quality and 100 to excellent water quality, considering five different potential uses of water (public water supply, fish and shellfish, agricultural, industrial, and recreational uses).The WQI considers the percentage of oxygen saturation (%DO), BOD 5 , NO 3 concentration, water hardness, conductivity, color, total and fecal coliforms, alkalinity, chlorides, pH, and atmospheric and water temperature, using the following equation: where WQI = Water Quality Index (0 to 100); I i = subindex of the ith parameter (0 to 100); W i = weighting value of the parameter (0 to 1); n = number of parameters.

Land Uses
For the characterization of study sites based on land uses, an interactive map was developed based on the geographic layers of land use at 1:250,000 scale provided by INEGI (National Institute of Statistical and Geography), using the QGIS software (Las Palmas version) (Open Source Geospatial Foundation, Chicago, IL, USA).For each study site, "buffers" (i.e., zones influenced by land use) were set measuring 2 km in length upstream of the monitoring site by 500 m at both sides (2 km × 500 m).Additionally, land use was assessed watershed upstream from each study site.In this case, we calculated the percentage of each land use by means of the Simulador de Flujos de Agua de Cuencas Hidrográficas (SIATL) (3.2) [42]; these were subsequently grouped into major land-use categories.

Statistical Analysis
A dissimilarity analysis was conducted between study sites through Euclidean distances, with the dendrogram drawn according to the Ward's hierarchical clustering method using averages of physicochemical parameters for each study site to identify sections along the river and to define clusters for comparison.
Homoscedasticity testing was performed with cluster data; afterwards, analysis of variance (ANOVA) (α = 0.95) and multiple comparison tests between clusters were performed using the Student-Newman-Keuls method.
The physicochemical characterization of study sites and their relationship with land uses were performed through a principal component analysis with three dimensions (PCA-3D) using the Spearman's correlation coefficient.To this end, a matrix was built with study sites as rows and environmental parameters and land uses as columns.This analysis considered averages of physicochemical parameters for sampling seasons and percentage of each land use.The data set was standardized using the formula X = log (x + 1); for percentages, the formula was 2π × arcsin √ p, where p is the relative value of each individual land use.All statistical analyses were run using the software XLStat-2015 (Addinsoft, Paris, France).

Clustering of Study Sites According to Their Physicochemical Properties
The dendrogram of the dissimilarity analysis of study sites (Figure 2) revealed three clusters.The first (C1) grouped together sites PI, HU, MA, TZ, and FI and was defined as the upper course.C2 clustered sites RG and PA (the middle course of the Bobos River) and sites JL, TM, and EN (from the Alseseca River tributary); finally, a third cluster (C3) included sites PL and MTZ (the lower course).

Longitudinal Profile of the River
The longitudinal profile of the river, including topographic relief and slope of each river section, as defined by the dissimilarity analysis, led to the identification of a marked change in gradient (Figure 3A,B); as expected, C1 showed the steepest slope (4.76%), C2 attained an intermediate slope (3.18%), and C3 had the lowest slope (0.44%).

Longitudinal Profile of the River
The longitudinal profile of the river, including topographic relief and slope of each river section, as defined by the dissimilarity analysis, led to the identification of a marked change in gradient (Figure 3A

Longitudinal Profile of the River
The longitudinal profile of the river, including topographic relief and slope of each river section, as defined by the dissimilarity analysis, led to the identification of a marked change in gradient (Figure 3A,B); as expected, C1 showed the steepest slope (4.76%), C2 attained an intermediate slope (3.18%), and C3 had the lowest slope (0.44%).
The mean WQI attained the highest value in the upper course (78.05 ± 0.97), which was significantly different from the lower course (72.38 ± 1.48; p < 0.05).Furthermore, in accordance with the water quality categories of [41], WQI values for the upper course can be classified as follows: "Acceptable for all water sports; "Acceptable for fish and shellfish"; "Minor purification for crops requiring high-quality water"; "Treatment needed for public water supply."By contrast, WQI values in the lower course correspond to the following categories: "Becoming polluted to practice any water sports"; "Marginal quality for sensitive fish and shellfish"; "No treatment needed for most crops and industries", "Treatment required for public water supply." On the other hand, averages by season (Figure 4) yielded the lowest values in the rainy seasons of August 2013 and June 2014 (71.30 ± 0.77 and 73.82 ± 1.04, respectively), which were significantly different (p < 0.05) versus December 2013 and February 2014, when the highest values were recorded (79.36 ± 0.97 and 84.42 ± 0.93, respectively).According to WQI values, the seasonal patterns show differences only in the case of public water supply, which fell into the category "Minor purification required" in the dry season, but into the category of "Necessary Treatment for Public Water Supply" in the rainy season.The global WQI for the whole study area ranges from mildly polluted for human use to excellent quality for industrial, recreational, and fishing uses.
Water 2018, 10, x FOR PEER REVIEW 7 of 16
The mean WQI attained the highest value in the upper course (78.05 ± 0.97), which was significantly different from the lower course (72.38 ± 1.48; p < 0.05).Furthermore, in accordance with the water quality categories of [41], WQI values for the upper course can be classified as follows: "Acceptable for all water sports; "Acceptable for fish and shellfish"; "Minor purification for crops requiring high-quality water"; "Treatment needed for public water supply."By contrast, WQI values in the lower course correspond to the following categories: "Becoming polluted to practice any water sports"; "Marginal quality for sensitive fish and shellfish"; "No treatment needed for most crops and industries", "Treatment required for public water supply." On the other hand, averages by season (Figure 4) yielded the lowest values in the rainy seasons of August 2013 and June 2014 (71.30 ± 0.77 and 73.82 ± 1.04, respectively), which were significantly different (p < 0.05) versus December 2013 and February 2014, when the highest values were recorded (79.36 ± 0.97 and 84.42 ± 0.93, respectively).According to WQI values, the seasonal patterns show differences only in the case of public water supply, which fell into the category "Minor purification required" in the dry season, but into the category of "Necessary Treatment for Public Water Supply" in the rainy season.The global WQI for the whole study area ranges from mildly polluted for human use to excellent quality for industrial, recreational, and fishing uses.

Land Uses
Nine different land uses were identified in the assessment of land use watershed upstream from each study site (tropical mountain cloud forest, TMCF; coniferous forest, CF; high evergreen forest, HEF; humid soil agriculture, HSA; rainfed agriculture, RA; induced pasture, IP; cultivated pasture, CP; human settlements, HS; urban zones, UZ) (Figure 5).Land-use percentages obtained from this assessment (Figure 5A) show that the tree types of forests (TMCF, CF, and HEF) are represented in all study sites, ranging from 24 to 38% in C3 and C2, and from 41 to 65% in CI.Agriculture uses (RA, IP, CP, and HSA) ranged from 62 to 76% in C3 (particularly HSA) and from 35 to 58% in both C2 and C1.In all clusters, HS and UZ showed values below 1%.On the other hand, the assessment of land use including the buffer zone (2 km × 500 m) (Figures 5B and 6) showed two types of forest (TMCF and HEF); in C1, the percentages ranged from 23 to 100%, although forests were not represented in three study sites of this cluster (TZ, MA, and FI); HEF was observed in only one C2 site (PA: 8%) but was completely absent in C3.HS and UZ ranged from 25 to 42% in C3; in C2, UZ varied from 5 to 20% only in study sites PA and RG; C1 shows no influence of this land use.
Water 2018, 10, x FOR PEER REVIEW 8 of 16

Land Uses
Nine different land uses were identified in the assessment of land use watershed upstream from each study site (tropical mountain cloud forest, TMCF; coniferous forest, CF; high evergreen forest, HEF; humid soil agriculture, HSA; rainfed agriculture, RA; induced pasture, IP; cultivated pasture, CP; human settlements, HS; urban zones, UZ) (Figure 5).Land-use percentages obtained from this assessment (Figure 5A) show that the tree types of forests (TMCF, CF, and HEF) are represented in all study sites, ranging from 24 to 38% in C3 and C2, and from 41 to 65% in CI.Agriculture uses (RA, IP, CP, and HSA) ranged from 62 to 76% in C3 (particularly HSA) and from 35 to 58% in both C2 and C1.In all clusters, HS and UZ showed values below 1%.On the other hand, the assessment of land use including the buffer zone (2 km × 500 m) (Figures 5B and 6) showed two types of forest (TMCF and HEF); in C1, the percentages ranged from 23 to 100%, although forests were not represented in three study sites of this cluster (TZ, MA, and FI); HEF was observed in only one C2 site (PA: 8%) but was completely absent in C3.HS and UZ ranged from 25 to 42% in C3; in C2, UZ varied from 5 to 20% only in study sites PA and RG; C1 shows no influence of this land use.

Multivariate Analysis
The principal component analysis (PCA-3D) accounted for 69.51% of the cumulative variance of the data set in the first three axes (F1: 35.56%,F2: 20.91%, and F3: 13.03%); factor loadings are shown in Table 1.The first component (F1) showed the importance of the altitudinal gradient and the variables associated with changes along the river course: altitude and WQI are inversely proportional to turbidity, temperature, nitrogen (NH3 and TN), coliforms (total and fecal), total suspended solids, alkalinity, chloride, color, and land uses HS and UZ; F2 is related to nutrient enrichment (NO3 and O-PO4) and mineralization (conductivity, SO4, and hardness; land uses related to F2 were cultivated and induced pastures, all of which were directly proportional to F2); finally, F3 is related to DO, BOD, and rainfed agriculture, which directly proportional to F3 (Table 1).The 3-D plot (Figure 7) shows clusters of study sites: Cluster I comprises sites PI, HU, TZ, MA, and FI, characterized by the highest conductivity, alkalinity, hardness, dissolved oxygen, and sulfates.These sites also reached the highest WQI values of the six periods studied; in addition, these showed the highest percentages of natural vegetation (tropical mountain cloud forest and high evergreen forest), as well as rainfed agriculture from TZ to FL, along with clearing of natural vegetation.Cluster II includes sites located along the Alseseca river tributary (JL, TM, and EN) and sites in the middle course (RG and PA), characterized by the highest levels of nitrates, nitrites, orthophosphates, total phosphorus, chlorides, BOD5, and fecal coliforms.These sites were closely related to rainfed agriculture, as well as to induced and cultivated pastures.Finally, Cluster III includes sites located in the lower course (MZ and PL), associated with the highest percentage of human settlements and urban zones, humid soil agriculture, and the highest values of total nitrogen, ammonia nitrogen, pH, total coliforms, turbidity, suspended solids, and color (Table 2).

Multivariate Analysis
The principal component analysis (PCA-3D) accounted for 69.51% of the cumulative variance of the data set in the first three axes (F1: 35.56%,F2: 20.91%, and F3: 13.03%); factor loadings are shown in Table 1.The first component (F1) showed the importance of the altitudinal gradient and the variables associated with changes along the river course: altitude and WQI are inversely proportional to turbidity, temperature, nitrogen (NH 3 and TN), coliforms (total and fecal), total suspended solids, alkalinity, chloride, color, and land uses HS and UZ; F2 is related to nutrient enrichment (NO 3 and O-PO 4 ) and mineralization (conductivity, SO 4 , and hardness; land uses related to F2 were cultivated and induced pastures, all of which were directly proportional to F2); finally, F3 is related to DO, BOD, and rainfed agriculture, which are directly proportional to F3 (Table 1).The 3-D plot (Figure 7) shows clusters of study sites: Cluster I comprises sites PI, HU, TZ, MA, and FI, characterized by the highest conductivity, alkalinity, hardness, dissolved oxygen, and sulfates.These sites also reached the highest WQI values of the six periods studied; in addition, these showed the highest percentages of natural vegetation (tropical mountain cloud forest and high evergreen forest), as well as rainfed agriculture from TZ to FL, along with clearing of natural vegetation.Cluster II includes sites located along the Alseseca river tributary (JL, TM, and EN) and sites in the middle course (RG and PA), characterized by the highest levels of nitrates, nitrites, orthophosphates, total phosphorus, chlorides, BOD 5 , and fecal coliforms.These sites were closely related to rainfed agriculture, as well as to induced and cultivated pastures.Finally, Cluster III includes sites located in the lower course (MZ and PL), associated with the highest percentage of human settlements and urban zones, humid soil agriculture, and the highest values of total nitrogen, ammonia nitrogen, pH, total coliforms, turbidity, suspended solids, and color (Table 2).

Discussion
The water quality index proposed by Dinius [41] has been used by several authors to indicate the effect of anthropogenic impacts and used as an easy-to-interpret tool for assessing water quality [43].The results obtained in this work demonstrate a marked temporal variation, with the dry cold season (February 2014) showing the highest water quality; however, as regards spatial variation, this index is insensitive to variations between sites, seemingly unsuitable for identifying sections with different characteristics along the river.Some authors [44] stress that this analysis should be conducted in parallel with an assessment of the changes of land use, since study sites with different land uses yielded very similar scores.The above derives from the fact that the WQI does not include water-quality parameters that are affected by changes of land use, such as nutrients (orthophosphates, ammonia, nitrites, and nitrates) [44][45][46].However, the combined analysis of the 23 physicochemical variables and all land uses considered in this study and using multivariate methods (CA and PCA) allowed for the identification of the relationship between changes of land use and physicochemical water quality, and revealed spatial variation patterns.
The CA identified three regions along the course of the river, which share local issues requiring particular attention or management.For instance, C1 (upper course) shows effects of the changes of land use (HP, CP, and RA) according to the PCA, although it showed the highest coverage of TMCF vegetation.These changes of land use were most evident in sites TZ, MA, and FL, characterized by nutrient enrichment.C1 sites, despite their location in the best preserved area, should be closely monitored and deforestation should be regulated.Furthermore, the steep slope along the river course in this zone promotes self-purification.Nilsson et al. [47] indicate that the upper course of any stream is key for the conservation of water quality in the rest of the basin, as it plays a key role in maintaining the buffering system and diluting pollutants downstream.
For its part, the PCA evidenced that C2 is closely related to agricultural practices (RA, CP, and IP account for the higher nutrient levels in sites TM, JL, RG, and PA), contributing to high nitrogen (NO 2 and NO 3 ) and phosphorus (O-PO 4 ) inputs.This has been highlighted by various authors [47][48][49][50], who point out that increased nutrient levels in water bodies are directly related to the use of agrochemicals, with higher concentrations during the rainy season as a result of leaching and runoff.This fact evidences the need to implement best management practices to prevent excess nutrient discharges into the river.
Last, C3 (MTZ and PL) was associated with human settlements (UZ, HS, and HSA) and showed wastewater indicators (fecal and total coliforms, ammonium, total N, and chlorides).Due to the low gradient in this section of the river, its self-cleaning capacity decreases sharply, as reported by [51] for the lower course of the Potrero de los Funes river, Argentina.The zone of the river associated with C3 requires the establishment of wastewater treatment facilities before discharge into the river.In all cases, the studies relating land-use changes at different spatio-temporal scales show that forest clearing or replacement of natural vegetation coverage within a basin leads to the deterioration of water quality [48,49,52].
The combination of two multivariate tools was useful in this study.The cluster analysis produced groups based on similarities and differences; for its part, the PCA allowed identifying the variables, producing a new set of variables (principal components) that are consistent with the clusters derived from CA.In this study, the PCA confirmed its validity as a method for monitoring changes in water systems and the key drivers associated with them.
In contrast with the assessment of the effect of land use in the watershed upstream from each study site, the interpretation of the effects of land use on water quality from the analysis of buffer zones involves scale and interpretation implications according to the research hypothesis established a priori.It has been pointed out [53,54] that the significant effects on water quality result from changes of land use in the areas closest to the study sites, since it is in these areas that the watershed is being washed off of materials, which enter river water and cause immediate effects.Our findings reveal that, when the influence of land uses is explored according to the watershed upstream from each study site, C2 and C3 include percentages of forests (TMCF, HEF, and CF) with an altitudinal range lacking an evident relationship with the river course; thus, these may not directly affect water quality.For their part, AH and ZU display a minor relative influence, since both represent a small and constant percentage across all study sites (<1%), particularly in C3, which is the area where human settlements are located in the proximity of the river; consequently, wastewater indicators are evident.On the other hand, our results show that changes at a local scale (buffer land uses) significantly affect water quality, as stressed by [7].In this respect, our results show that study sites located in the upper course (PI and HU), with extensive coverage of natural vegetation (92% and 100%, respectively), show good water quality, as nutrient concentrations were low, while C2 and C3 showed the highest impact by agricultural practices and human settlements.
Gove et al. [55] found a greater influence of upstream vs. lateral inputs as the length of the catchment area increases.The fact that these increases reach a point of diminishing returns indicates that processing within the river eventually dampens the memory of upstream events relative to more recent longitudinal and lateral inputs.Therefore, if the total cumulative upstream land use governs water quality at any point in the basin, water quality measurements would exhibit little variations downstream.In this study, when watershed upstream from each study site are considered, land uses that are distant from the river course become evident, offsetting the effect of land uses adjacent to the river.On the other hand, the use of a 2 km × 500 m buffer zone use reveals a shift of land uses with an immediate effect on the water body.The cluster analysis showed three different groups based on similar water quality characteristics, which are consistent with Figure 5B.In this sense, and considering the objectives of this study, the use of information derived from small-scale buffers seems more appropriate, since these actually provide information on the effects of land uses adjacent to the river course.
The analysis of results discussed in this work proposes a baseline assessment methodology for intertropical latitudes, as environmental monitoring data are scarce in many Latin American countries [25], and Mexico has poor records of historical changes of either land use or water quality.
Based on the above, developing countries in tropical zones should start to identify the direct effects of different land uses using approaches similar to the one described herein, in order to investigate cause-effect responses.According to [55], the incorporation of the effect of the watershed upstream of each study site could be masking/neutralizing the effects of land uses adjacent to the river course.Working with small-scale buffers (e.g., 2 km × 500 m) and building a monitoring network produce baseline data on the effects of land uses in different basins, and produce useful parameters for incorporation into various modeling tools such as SWAT.In turn, this facilitates building models on potential changes of land use and the expected hydrological changes, which should be taken into account by policy makers focused on sustainable water resource management and climate change, as suggested in [56].In turn, these tools predict the effects of changes of land uses; examples include [3], wherein future changes in water quality within a watershed dominated by agriculture are anticipated, [57], wherein various hydrological land-use scenarios in the upper portion of the Mississippi river are modeled to 2091 and the impact of land uses on hydrologic processes is determined, and [1], wherein different scenarios of the effects of climate change on water quality are assessed.In this sense, our results can be used as a small-scale quantitative component for use in future predictions about changes of land use and land cover, as defined in [58], and in the development of management plans and public policies for conservation.

Conclusions
The loss of tropical cloud forest is mainly due to the change of land cover to anthropogenic land uses.Our results revealed that streams located in this basin are highly vulnerable to changes of land use.These changes alter the proper functioning of the tropical cloud forest and the ecosystem services, including water quality.Hence, conservation measures targeting cloud forests and high evergreen forests, along with reforestation programs in the rest of the basin, should be implemented; wastewater treatment plants should also be established in the lower course of the river-the section most affected by direct pollution issues and where the river self-purification capacity is not as efficient as in the upper course.Our results evidence that the natural protected area in the river urgently requires the implementation of reforestation programs to improve the condition of the local ecosystems (streams and forest).

Figure 2 .
Figure 2. Dendrogram of dissimilarity of Euclidean distances between study sites based on average physicochemical values of all study periods.

Figure 3 .C1Figure 2 .
Figure 3. Altitudinal profile and topographic relief of the Bobos-Nautla River.(A) Altitudinal profile of streams in the Bobos-Nautla River basin.The relative slope (%) for each section is shown in brackets according to the clustering from the dissimilarity analysis.(B) Topographic relief of the Bobos-Nautla River watershed.

Figure 2 .
Figure 2. Dendrogram of dissimilarity of Euclidean distances between study sites based on average physicochemical values of all study periods.

Figure 3 .C1Figure 3 .
Figure 3. Altitudinal profile and topographic relief of the Bobos-Nautla River.(A) Altitudinal profile of streams in the Bobos-Nautla River basin.The relative slope (%) for each section is shown in brackets according to the clustering from the dissimilarity analysis.(B) Topographic relief of the Bobos-Nautla River watershed.

Figure 4 .
Figure 4. Mean values of water quality index (WQI) by study site, season, and river portion.a ≠ b (p < 0.05).

Figure 4 .
Figure 4. Mean values of water quality index (WQI) by study site, season, and river portion.a = b (p < 0.05).

Figure 5 .Figure 5 .
Figure 5. Percentage of land use in each study site.(A) Percentage of land use for watershed upstream from each study site.(B) Percentage of land use for the upstream 500 m × 2 km buffer zone.

Figure 6 .
Figure 6.Percentage of land uses from buffer and mean WQI of study sites along the Bobos River.Tributaries and streams are shown in colors as in Figure 1.

Figure 6 .
Figure 6.Percentage of land uses from buffer and mean WQI of study sites along the Bobos River.Tributaries and streams are shown in colors as in Figure 1.

Figure 7 .Table 2 .Figure 7 .
Figure 7. Principal component analysis for study sites, physicochemical variables, WQI, and land uses.(A) Environmental and land-use variables; (B) clusters of study sites.
Water 2018, 10, x FOR PEER REVIEW 11 of 16

Table 2 .
Mean values (± Standard Error (SE)) of components related to N and P, WQI, and predominant land uses (%) for each study site.