Water Quality Analysis in a Subtropical River with an Adapted Biomonitoring Working Party (BMWP) Index

: Subtropical rivers in developing countries often lack adequate monitoring, which makes it difﬁcult to comprehensively determine their water quality when faced with different anthropic impacts. There are no proper protocols in the regulations to incorporate indicators and adapt them to different biogeographic regions, limiting the potential success of conservation and restoration of river ecosystems. This study proposes implementing macroinvertebrates as bioindicators of water quality in river ecosystems, and modifying the calibration of the widely used Biomonitoring Working Party (BMWP) index for its adaptation in a subtropical river. The Duero River, Mexico, was used as an example in this study. Data were explored with multivariate statistics, and the water quality and habitat values were averaged to obtain the families’ bioindication values and the index categories. The BMWP adequately described a deterioration gradient from the origin to the river mouth (from fair to extremely polluted), with some intermediate recovery points related to the presence of springs. Its performance was compared with other biological indices and exhibited a positive relationship with all of them. In addition, how BMWP changed over time was analyzed by examining previous samples, and highlighted increased river deterioration over time. A calibrated BMWP will allow for long-term monitoring at a low cost.


Introduction
One of the most commonly used biological indices to monitor the quality of lotic aquatic ecosystems is the biological monitoring working party (BMWP) [1][2][3][4][5]. The BMWP is a biotic approach because it includes taxonomic groups, considering their sensitivity or tolerance to pollution, and both aspects are incorporated into an index [3]. The index describes and analyzes the macroinvertebrate community at the taxonomic family level [6]. The characteristics that make macroinvertebrates good bioindicators are their wide distribution, limited mobility, numerical abundance, sensitivity and response to distinct types of environmental conditions and stressors, and ease of finding, quantifying and standardizing them [7]. The application of the BMWP and other indices that use bioindicator organisms provides complementary information on biotic and abiotic conditions of lotic ecosystems in addition to traditional monitoring techniques (physical, chemical, and bacteriological variables). The development of such indices, arises from the need to systematically reflect changes in riparian and fluvial ecosystems and to express the environmental and habitat factors in an integral way, with the expectation of long-term management [8].
The BMWP assigns a bioindication value for each family and is the minimum perceived value of the tolerance of macroinvertebrates to organic contamination. It ranges from 1 (very tolerant) to 10 (very sensitive). Consequently, the total index score for a sample is defined as the sum of the minimum perceived value of the tolerance of all the families present [9]. In general, if the sum is approximate or greater than 100, the rivers and streams are deemed healthy, but if the sum is less than 10 they are considered highly polluted [6]. The BMWP has been used in various regions worldwide, and has been adapted, mainly due to the absence of some taxa used in the original version and the presence of others not initially included [10,11]. In addition, the modifications include the combination of families due to taxonomic difficulties and changes in some families' values, which is sometimes related to their frequency of occurrence and the degree of saprobity [1,12]. A more recent calibration and validation of the index has been proposed to relate it more closely to a river's specific environmental characteristics [13].
Despite these advantages, there are also some disadvantages. The extrapolation of the BMWP in large regions or countries results in poor generalizations because the sensitivity of some organisms changes over space and time [14]. There is subjectivity in the allocating tolerance intervals because fine taxonomic levels can bias the index, since the macroinvertebrate genera and species within the families have different sensitivities to distinct types of environmental degradation [9]. Therefore, a calibration process has been implemented, which was also suggested in the original proposal to establish suitable bioindication values [4,15]. Accordingly, once an adequate monitoring scheme is established and the protocol for obtaining and reviewing samples in the whole system is developed (e.g., Cornejo et al. [16]), with an adequately adapted index proposal (e.g., Ruiz-Picos et al. [17]), then the results will better reflect water quality condition of the aquatic ecosystem. Particularly, it is important to implement frequent monitoring and measure the BMWP regularly to satisfactorily describe the ecosystem conditions.
In this context, different countries, mainly in North America, Europe, Asia, and Australasia, have adopted the use of macroinvertebrates in their environmental regulations and have long-term monitoring programs [18][19][20][21]. More recently, some developing countries have also assumed this approach [16,22,23], but several nations still lack official biomonitoring programs. This deficiency is not necessarily related to the absence of biological community studies, but to the inadequate standardization and sporadic monitoring depending on the river typology, habitat type, and time of year [24].
Although a protocol has been developed to sample macroinvertebrates in Mexico and includes the application of the BMWP to be adopted in the official normative about ecological flow (NMX-AA-159-SCFI-2012, [25]), it concludes that the bioindication values and the water quality ranges should be adapted according to the sampled system. This adaptation should consider the different biogeographic, physiographic, environmental, and climatic characteristics as well as the anthropogenic impacts that affect aquatic ecosystems in the country. Some isolated studies in tropical systems of Mexico have already implemented the BMWP using the values and ranges from other places [26,27] or implemented the calibration and validation process according to environmental characteristics [17].
We focused on the subtropical Duero River because it is located within a different climate region than the previous studies and belongs to the so-called Mexican Transition Zone, which joins the biogeographical Nearctic and Neotropical regions [28]. This biotic crossroad is characterized by the presence of several hydric resources with high biodiversity and endemicity and is also related to ancient prehispanic human intervention. Additionally, previous studies on the Duero River have analyzed water quality and entomofauna during different time periods since the 1980s [29][30][31]. In this river, environmental conditions have shown contrasting results, with some upstream localities and sites related to springs showing good water quality (e.g., control sites); in contrast, other localities experience different impacts with high levels of organic and bacteriological contamination, and are associated with urban development and agricultural and livestock production [32]. Moreover, another study included the addition of other bioindicators, such as fish, to describe the system's biotic integrity in specific decades [33]. In this context, the current study's main aim was to develop a calibrated BMWP according to the river environmental conditions and compare the results in different periods and with other biological indices. We hypothesized that (1) the ecosystem presents a general spatial gradient from least polluted at the origin to severely polluted at the mouth, related to human impacts, and contains intermediate recoveries by occasional tributaries from the springs, (2) the category of the BMWP scores does not necessarily coincide with the category of the index of water quality and the index of habitat quality because some habitat modifications (not included in the index of water quality) and pollution (not included in the index of habitat quality) affect the macroinvertebrates distribution and abundance directly, and (3) the river has experienced a degradation in water quality over time.

Study Area
The Duero River is located on the Central Plateau of Mexico and it is an affluent of the Lerma River, which finally enters Chapala Lake, the largest natural lake in the country. The river has a catchment area of 2531 km 2 , a length of 85 km, and a total annual flow of 457.8 hm 3 [34]. The source of the river is at 1860 m.a.s.l. and has an altitude change of 340 m [35]. The catchment of the Duero River is divided into four sub basins ( Figure 1). The upper sub basin (Cañada) consists of volcanic rocks (andesites, basalts), and it is the primary recharge zone, while in the middle and lower zones there are three sub basins containing valleys (Guadalupe, Zamora, and Ciénega) and where vertisol soils predominate, supporting agricultural land use (143,333 ha) [31]. Although there are 52 springs, and theoretically there is enough water available to meet the needs of different users, this is not the case due to inadequate management and distribution; consequently, the Duero River has a deficit of 44.7 hm 3 s −1 [34].

Sampling and Data Analysis
We selected 14 sites, some of which have good water and habitat conditions, whereas others were affected by anthropogenic impacts associated with principal human activities and urban development ( Figure 1). Most of the sites are located in the main river course, but three sites are near springs in tributaries, La Toma (site 3), the national park Camécuaro Lake (site 6), and Verduzco (site 9). Water and macroinvertebrates were sampled under different flow conditions: low velocity in the dry season (May), high in rainy season, and intermediate in transition season (November) in 2019.
We used the multihabitat protocol proposed by Cornejo et al. [16] to sample macroinvertebrates. In a 100 m stretch of the river or tributary, we sampled the main aquatic habitats (e.g., bed substrate, aquatic and riparian vegetation, areas of different water characteristics). The number of samples per habitat depended on the percentage coverture. Samples covered an area of approximately 1 m 2 , and the organisms were captured with two nets; in fast-flowing riffles and pools, we used a kick net and for floating, submerged, and riverine vegetation a type-D net (both with a mesh size of 500 µm). We preserved the macroinvertebrates collected with 70% alcohol. Organisms were identified to the family level according to different taxonomical keys [37][38][39].

Sampling and Data Analysis
We selected 14 sites, some of which have good water and habitat conditions, whereas others were affected by anthropogenic impacts associated with principal human activities and urban development ( Figure 1). Most of the sites are located in the main river course, but three sites are near springs in tributaries, La Toma (site 3), the national park Camécuaro Lake (site 6), and Verduzco (site 9). Water and macroinvertebrates were sampled under different flow conditions: low velocity in the dry season (May), high in rainy season, and intermediate in transition season (November) in 2019.
We used the multihabitat protocol proposed by Cornejo et al. [16] to sample macroinvertebrates. In a 100 m stretch of the river or tributary, we sampled the main aquatic habitats (e.g., bed substrate, aquatic and riparian vegetation, areas of different water characteristics). The number of samples per habitat depended on the percentage coverture. Samples covered an area of approximately 1 m 2 , and the organisms were captured with two Some of the water parameters were used to determine the water quality index developed by the National Sanitation Foundation (NSFWQI) [40]. The habitat quality index (HQI) related to the physical structure of the river and its surroundings was to describe each sampling event following the criteria proposed by Barbour et al. [41]. The biological quality of the river was evaluated with the use of macroinvertebrates according to three indices: the percentage of Ephemeroptera (minus the Baetidae family), Plecoptera, Trichoptera individuals index (EPT-B) [42], the Hill's numbers that describe biological diversity [43], and the BMWP adapted to Mexican rivers [13]. The EPT index was modified because Masese and Raburu [42] proposed the subtraction of the Baetidae, Caenidae, and Hydropsychidae families; however, we only subtracted Baetidae because this family was the only one that showed low bioindication values in different tropical areas [1,12,16,17]. Hill numbers describe the effective number of species based on a sampling framework with rarefaction and extrapolation methods and derives species richness (q = 0), Shannon's entropy index exponential (q = 1), and the inverse of Simpson's concentration index (q = 2) [43].

Index Calibration and Comparison
Although the BMWP has been used in different rivers in Mexico, Ruiz-Picos et al. [13] emphasize the necessity to adapt the index, not only to the macroinvertebrate families present in the region, but also to the specific environmental characteristics of the ecosystem studied. They propose a complete protocol to calibrate the BMWP, which includes four main steps. First, the physicochemical quality index (Pcq) is calculated. From the average physicochemical data set (seasons), the variables that better describe the sites are chosen using statistical methods and according to the parameter's range and the way it describes impacts. The water parameters selected are standardized and combined on a scale of 0 to 10, from poor to better water quality, and sites are arranged within this scale. Second, the bioindication value is obtained. The macroinvertebrate abundance is structured into five classes from zero (no individuals) to 5 (>100 individuals). Every family abundance class is assigned to the different Pcq intervals according to their presence or absence in the sites. The bioindication value of every family is then derived from the 5th percentile of each abundance class. Third, these bioindication values are replaced at every site and in every study period for each family found, and they are summed to evaluate the BMWP index. Fourth, the general scale of the index in the whole system is established from the median and tenth percentile values of the sites with better water quality. Above the median is excellent quality, in the tenth percentile is good quality, and below this value, the amount is divided by four to get the rest of the categories: regular, polluted, very polluted, and extremely polluted ( Figure 2).  [13]. In red and italics, the modification implemented in the present study.
We followed this protocol to obtain BMWP-calibrated values for the Duero River, but we also implemented two modifications. In the first step, we average the Pcq calculated with the HQI values to combine the information from both indices. The HQI was rescaled from 0 to 10 to have compatible values with the Pcq because the original scale was from 0 to 200 [41]. We assumed this initiative because the physical and chemical parameters of the water do not always entirely reflect the characteristics of the sites (i.e., good water quality but bad habitat conditions). Consequently, this index reflected not only the organic pollution but also the habitat quality in every site. In the fourth step, from the multivariate analysis of the differences and trends in the environmental and community variables, the season and sites were selected to establish the general BMWP scale for the river. We used the results from the calibration to compare different sites at distinct time periods (1984 and 1986 [29]; and 2013 [30]). Moreover, we related the BMWP with the physical and chemical parameters to identify which of them better explain the index variations, and with other biotic indices, such as the EPT-B, and Hill's numbers [43], to identify the similar or dissimilar response to pollution [6,44].  [13]. In red and italics, the modification implemented in the present study.
We followed this protocol to obtain BMWP-calibrated values for the Duero River, but we also implemented two modifications. In the first step, we average the Pcq calculated with the HQI values to combine the information from both indices. The HQI was rescaled from 0 to 10 to have compatible values with the Pcq because the original scale was from 0 to 200 [41]. We assumed this initiative because the physical and chemical parameters of the water do not always entirely reflect the characteristics of the sites (i.e., good water quality but bad habitat conditions). Consequently, this index reflected not only the organic pollution but also the habitat quality in every site. In the fourth step, from the multivariate analysis of the differences and trends in the environmental and community variables, the season and sites were selected to establish the general BMWP scale for the river. We used the results from the calibration to compare different sites at distinct time periods (1984 and 1986 [29]; and 2013 [30]). Moreover, we related the BMWP with the physical and chemical parameters to identify which of them better explain the index variations, and with other biotic indices, such as the EPT-B, and Hill's numbers [43], to identify the similar or dissimilar response to pollution [6,44].

Statistical Analyses
From the total samples (N = 42), two methods were used to select the environmental variables that best described the sites and to obtain the Pcq index. First, Spearman's rho correlation was used to identify redundancy and collinearity among the variables (≥0.9 and p < 0.0001) [45]. The selection of one of the variables with similar contributions was made, on the one hand, according to their magnitudes and, on the other, to reduce the influence of specific water characteristics in the subsequent analyses (several variables measuring similar aspects, like ion composition). Second, we performed a principal component analysis (PCA) to reduce the number of variables and to identify those variables that have a higher contribution to the variance in the dataset (a threshold to the chosen principal components > 75% and a variable loading > 0.7). From 15 variables, 5 were excluded: salinity, total dissolved solids, chloride, depth, and percentage saturation of DO.
We used a distance-based permutational multivariate analysis of variance (PER-MANOVA) [46] in a two-way factorial design to test differences among sub basins and seasons. The information was log-transformed to reduce the effect of distinct units and magnitudes, and then we used the Euclidian dissimilarity coefficient in the environmental matrix; in the community dataset we used the Bray-Curtis dissimilarity due to the high frequency of zeros [47]. Prior to the PERMANOVA and because the data had an unbalanced design (different number of sampled sites per sub basin), we ran a method to test the homogeneity of multivariate dispersions (PERMDISP), since heterogeneity could heavily influence the PERMANOVA results [48]. In both datasets, the PERMDISP showed no significant difference: environment (sub basins F = 1.11 and p = 0.35; Seasons F = 1.57 and p = 0.21) and community (sub basins F = 2.58 and p = 0.07; Seasons F = 0.02 and p = 0.99).
To find the predominant trends in the dataset and to identify the spatial ordination of the sites in different seasons in relation to the macroinvertebrate families, we used nonparametric multidimensional scaling (NMDS) ordination with Bray-Curtis dissimilarity. We calculated three-dimensional solutions from 250 random starts of real data, with up to 1000 iterations, to evaluate stability, and obtained a final stress of 0.14. In the plot, we incorporated the environmental variables; although they were not part of the NMDS analysis, they were integrated as vectors, scaled by their correlation with the axes, and their significance was assessed using permutations [49].

Results
We found 62 families at the different sites and in different seasons (Table A1) Zamora differed from the other zones but Cañada, which is located in the upper portion of the basin (with Guadalupe: F = 3.49, p = 0.006; with Ciénega: F = 3.75, p = 0.006) and Cañada only with Ciénega (F = 2.5, p = 0.018; Figure 1).
The NMDS ordination showed an indirect gradient, and the first axis separated in the left to right direction, with sites reflecting better (higher oxygen and transparency, and the presence of more sensitive families) to worse conditions (higher Turbidity, BOD 5 , nitrogen nutrients, and mostly tolerant families; Figure 3). In this context, according to the amount of information, we duplicated the plot in the same figure, emphasizing this distinction.
To the left, we mainly found that most sites in the different seasons related to the origin of the river in the Cañada sub basin (Kunio, Tacuro, La Toma, but El Rastro and Etúcuaro) or to the presence of springs (Verduzco and Camécuaro). On the right side of the plot, the sites were mainly located near the river mouth in the Ciénega sub basin (Estanzuela, Capulín, San Cristóbal, Cumuato, but Ibarra). The environmental variables which most correlated to the axes and that showed significance using permutations were temperature (r 2 = 0.17, p = 0.04), conductivity (r 2 = 0.19, p = 0.001), dissolved oxygen (r 2 = 0.34, p = 0.001), ammonium (r 2 = 0.34, p = 0.001), nitrate (r 2 = 0.14, p = 0.001), turbidity (r 2 = 0.47, p = 0.001), and BOD 5 (r 2 = 0.18, p = 0.02).  The water parameters showed a general increasing trend in turbidity, BOD5, nitrogen nutrients, conductivity, and temperature, but dissolved oxygen and transparency, which were also related to the indirect gradient in the first axis. In addition, fluctuations in the physicochemical parameters along the river are common because of the presence of springs in the different sub basins, aspect that is evident by low values of turbidity, BOD5, and nitrogen nutrients (Table 1). Of particular note are the high values of fecal coliform bacteria present at most sites (average 9488-7,334,062 MPN/100 mL), which indicates the discharge of wastewater without treatment, and even very high values were found in one of the springs (Camécuaro from 19,000 to 5,333,333 MPN/100 mL).  The water parameters showed a general increasing trend in turbidity, BOD 5, nitrogen nutrients, conductivity, and temperature, but dissolved oxygen and transparency, which were also related to the indirect gradient in the first axis. In addition, fluctuations in the physicochemical parameters along the river are common because of the presence of springs in the different sub basins, aspect that is evident by low values of turbidity, BOD 5 , and nitrogen nutrients (Table 1). Of particular note are the high values of fecal coliform bacteria present at most sites (average 9488-7,334,062 MPN/100 mL), which indicates the discharge of wastewater without treatment, and even very high values were found in one of the springs (Camécuaro from 19,000 to 5,333,333 MPN/100 mL). To calibrate the BMWP for the Duero River, each family's bioindication scores were obtained, representing the minimum tolerance value to organic pollution ( Table 2) [17]. Despite the scale ranging from 1 to 10, no families were classified as a '1 because no site was highly polluted with few macroinvertebrate families, or as a '10 due to a lack of sites without apparent alteration. Additionally, the general quality categories for the index were defined. Temporally, we decided to establish this general BMWP scale of the river according to the measurements from the dry season (May) because it does not differ from the other seasons in the physicochemical and macroinvertebrate variables according to the PERMANOVA results, and it had a continuous moderate flow and stable conditions, unlike the rainy or intermediate seasons. Spatially, we selected four sites with better water quality, from which the median and tenth percentile were obtained. These sites were located at different parts of the river: two represent small creeks from springs (La Toma and Verduzco), and the other two were located in the main river with little contamination (Kunio and Antes Central de Abastos). Above the median, the water quality was deemed excellent, and between the tenth percentile and the median, the water quality was deemed good. Then, the tenth percentile is divided by four to obtain the rest of the water quality categories: regular, polluted, very polluted, and extremely polluted (Table 3) [52]. Table 2. Bioindication values were modified and adapted to the different families found in the Duero river as part of the BMWP calibration.

Families Score
Calamoceratidae, Corydalidae, Leptoceridae, Naucoridae, Nepidae, and Petaluridae 9 Anthomyiidae, Leptophlebidae, Perlidae, Polycentropodidae, Rhyacophilidae, and Sciomyzidae 8  The water, habitat, and biological quality described by the different indices during the dry season are presented in Table 4. No sites showed good conditions, and most were regular, but two were extremely polluted. Regarding the water quality variables, the BMWP had a positive correlation with dissolved oxygen (Spearman = 0.87, p = 0.00001) and a negative correlation with temperature (Spearman = −0.6, p = 0.02), conductivity (Spearman = −0.57, p = 0.03), and ammonium (Spearman = −0.64, p = 0.01). In the correlation of the BMWP with the biological indices, we found positive correlation with the EPT (Spearman = 0.67, p = 0.008) and the Hill numbers of Shannon's entropy index (Spearman = 0.88, p = 0.00001) and the inverse of Simpson's index (Spearman = 0.88, p = 0.00001). When we compared the historical BMWP values to try to identify changes over time, the analysis did not show a statistically significant difference among the years (χ 2 = 7.54, p = 0.11). However, when plotting the values at the different sites, sites at the origin of the river tended to decrease their water quality in recent years ( Figure 4).
(Spearman = 0.67, p = 0.008) and the Hill numbers of Shannon's entropy index (Spearman = 0.88, p = 0.00001) and the inverse of Simpson's index (Spearman = 0.88, p = 0.00001). When we compared the historical BMWP values to try to identify changes over time, the analysis did not show a statistically significant difference among the years (χ 2 = 7.54, p = 0.11). However, when plotting the values at the different sites, sites at the origin of the river tended to decrease their water quality in recent years (Figure 4).

Discussion
The calibrated BMWP was defined according to the community structure (abundance intervals of macroinvertebrates) associated with the river's environmental condi-

Discussion
The calibrated BMWP was defined according to the community structure (abundance intervals of macroinvertebrates) associated with the river's environmental conditions. We characterized the general physicochemistry and nutrient contents related to organic loads as defined in other studies and protocols [53,54]. However, the Pcq index showed higher values at sites with an evidently transformed habitat (e.g., river rectification, humanmade ponds, and small dams in the riverbed), such as Etúcuaro and Antes de la Central (sites 5 and 7, respectively), which placed them in a higher quality interval than expected. In contrast, we obtained small values of Pcq at sites with good habitat quality, such as the national park Camécuaro (site 6). Consequently, the combined Pcq and HQI value, better reflected the condition of the river at the different sites.
Adapting the biological index from the original version is a regular procedure because it establishes a better description of the conditions found in the ecosystem analyzed [55,56]. However, the way through which BMWP is calibrated is an important process, reducing subjectivity and taking into account the way in which macroinvertebrate families respond to environmental conditions in every system [17]. In this context, we used and modified the calibration procedure by Ruiz-Picos et al. [13] to propose a quantitative protocol. We based the analysis spatially on a low but representative number of sites, describing each sub basin's hydrology, water quality, habitat structure, and human impacts. In each site, we implemented a rigorous multi-habitat sampling to characterize all habitats in different stretches of the river. Temporally, the study covered the most contrasting climatological conditions present in the Duero river to describe better the water and habitat quality and the macroinvertebrates community. In addition, this represents the adapted BMWP in an important biogeographical zone (transition zone) within a highly productive and crowded basin in the central part of Mexico (Lerma River basin).
As expected from our first hypothesis, the ecosystem showed a general decrease down the stream in the water quality gradient, with some intermediate recoveries in the environmental variables and the macroinvertebrate community ( Figure 3). The BMWP clearly revealed this pattern, since in more transparent, oxygen-rich, and fast current sites, we found higher scores and index categories (regular; Table 4). In these sites we found macroinvertebrate families with higher bioindication values such as Ephemeroptera (Leptophlebiidae, Tricorythidae), Plecoptera (Perlidae), and Trichoptera (Glossosomatidae, Lepidostomatidae). These families are considered the most sensitive to pollution, and most have also been described in other subtropical and tropical rivers with equivalent bioindication values [13,16,57]. In contrast, the water quality decreased mainly in the final section of the river, where we found sites with high turbidity, pollution (high nutrients, DBO 5 , and fecal coliforms), low flow, and with some of the lowest BMWP categories (very polluted and extremely polluted; Table 4). In this area, although some of the macroinvertebrates found were present in all sites along the river, their dominance was higher in this section of the river, showing low bioindication values in the Duero River and other rivers, such as in the case for Diptera (Chironomidae), Oligochaeta (Lumbriculidae), and Isopoda (Asellidae) [54,58]. Other organisms also had low bioindication values at the river mouth, but they are from families with low frequency and abundance values, such as Hirudinea (Erpobdellidae, Glossiphonidae) and Diptera (Culicidae).
Two aspects are important to mention. First, some macroinvertebrates showed a different bioindication value compared with other studies [13,16,57]. For instance, the Ephemeroptera families, Heptageniidae and Leptophlebiidae, had lower and higher values in the Duero River (three and eight, respectively). However, the latter family value was similar to another study where they calculated the bioindication score according to the degrees of saprobity [12]. These results are part of the adaptation process of the BMWP and because of the calibration. For instance, Heptageniidae, which had low frequency and abundance was absent in most sites with better water quality but site 8 (Después de la Central). However, it had higher abundances at site 14 (Ibarra), the site with the lowest Pcq-HQI combined index. Second, regarding the intermediate recoveries, the springs improved the water quality at adjacent sites. For instance, the national park Camécuaro (site 6) has the highest discharge value in the watershed (Q = 1.75 m 3 s −1 , [34]), and it is related to the better water quality in Antes and Después Central de Abastos (sites 7 and 8, respectively).
For our second hypothesis, the results showed mixed trends among indices. For instance, Camécuaro showed contrasting results between the water and habitat quality indices (high values) and the EPT and BMWP (low value; Table 4). Low values could be related to habitat modifications for recreational use, such as the removal of submerged macrophytes to promote swimming and the introduction of non-native fish species to promote sport fishing; these aspects are also linked to the disappearance of other organisms at the site, such as the native fish Skiffia multipunctata [59]. Additionally, natural aspects in different zones, such as the higher water flow and the presence of roots and shade by ahuehuetes (cypress), restricts the growth of submerged vegetation and diminishes habitat diversity. Another site with different values was Ibarra (site 14) at the mouth of the river. This area of the river has been modified to be a small reservoir and the flow is controlled by a gate adjacent to the union with the Lerma River. The higher BMWP value at this site, in contrast to the other indices, could be related to continuous flow (the gates were open in all seasons) and the presence of water hyacinths that contribute to overall biotic diversity because they offer different habitats for macroinvertebrates [60].
In the third hypothesis, the lack of statistically significant seasonal or annual differences in water quality may result from a temporally stable gradient of water quality from mild pollution in the upper basin to severe pollution in the river mouth. However, on one hand, historically the Duero River had a macroinvertebrate composition of 72 families and only 62 are found currently. On the other hand, when comparing each site, there was a marked difference in the first sub basin related to habitat modification. Some springs were channelized to provide more water to nearby towns, or isolated with infrastructure to avoid degradation because several of these places are used as recreation areas, or are used to provide water for agriculture. According to the BMWP categories obtained in the present study, the main difference is that Kuinio, Tacuro, and La Toma (sites 1 to 3) showed that the average water quality was excellent and good in 1984 and 1986, respectively, and only regular in 2019. The higher values in the middle of the river, related to Antes and Después de la Central (sites 7 and 8), could be affected more by the sampling plan than by the conspicuous restoration of the habitat.
Unfortunately, this river has been used since pre-Hispanic times, and there are no major initiatives for its restoration, so there are no very high-quality sites. However, those that were used of better quality as controls allowed to make the calibration and the BMWP adequately reflect the Duero River conditions. In future studies, it is important to validate the BMWP by including additional sites and samples collected from multiple years to prove that the model fits well in the Duero River. A critical aspect of Duero River conservation is the maintenance of water flow to promote the river's biotic integrity because, in some years, different parts of the river dry up and lose continuity [33].