Aquatic Macroinvertebrate Assemblages in Rivers Under the Inuence of Mining Activities

Mining is one of the main pollution issues worldwide, causing the greatest disturbances to the environment. Industrial and artisanal mining activities are widespread in Mexico, being a major global producer of various metals. This study aimed to assess the ecological impairments resulting from mining activities using the aquatic macroinvertebrates assemblages (MA). A multiple co-inertia analysis (MCOA) was applied to determine the relationships between environmental factors, habitat quality, heavy metals, and aquatic macroinvertebrates in two rivers of the Central Plateau, Mexico. The results revealed three contrasting environmental conditions and different MA. High concentrations of heavy metals, nutrients, and salinity limit the presence of various families of seemingly sensitive macroinvertebrates, these factors were identiﬁed as the drivers of structural changes in the MA, showing that not only mining activities, but also agriculture, and villages in the basin, exert negative effects to the macroinvertebrate communities. Diversity indices showed that the lowest diversity matched with both, the most polluted and the most saline rivers. The rivers studied displayed a high alkalinity and hardness, which can lead to the formation of metal precipitates and thus acting as a protection to aquatic biota. Aquatic biomonitoring in rivers, impacted by mining and other human activities, is critical for detecting the effect of metals and other pollutants to improve management and conservation strategies. This study supports the design of cost-effective and accurate water quality biomonitoring protocols in developing countries. nutrients act as drivers of the community structure. contrasting and taxa richness) between Extoraz Escanela a set tolerance of some taxa to the presence of heavy metals. found high concentrations of heavy Cd Hg) sites some cases, life


Introduction
Mining is one of the main pollution issues worldwide, causing the most significant disturbances to the environment, involving air, terrestrial, and aquatic ecosystems [1][2][3] . Particularly in Latin America, among the main stressors of aquatic ecosystems are mining activities, which lead to metal pollution 4,5 . Additional significant sources of pollution in natural water bodies include agriculture (which drives deforestation, causing changes in land cover and land use) and the massive use of agrochemicals 6 . Furthermore, industry and human settlements discharge wastewater containing complex mixtures of pollutants 7 . All these stressors cause degradation of water quality of rivers, provoking adverse effects on aquatic life -particularly in tropical regions where human populations overgrow - 8 . Mexico is the second most populated country in Latin America 9 and one of the world s most biodiverse countries 10 . However, Mexico faces several serious pollution challenges from anthropogenic activities such as mining and agriculture, which affect biodiversity and ecosystems 11,12 . The country is the largest producer of silver in the world and a major global producer of gold, copper, and zinc, amongst other minerals 13 . Industrial and artisanal mining activities are widespread in the country, including opencast mining and mercury extraction, the latter associated with the use of mercury for gold amalgamation 14 . Artisanal mercury extraction is a prevalent activity in developing countries 15,16 . Efforts to protect and preserve biological diversity are usually focused on natural protected areas (NPAs), which are the cornerstone of conservation strategies worldwide 17 . The federal government has decreed 182 NPAs to date to preserve the biological diversity of Mexico. NPAs in México harbor a high diversity of various taxonomic groups 18 , most of which are not systematically monitored today; this is especially true of freshwater species 19 . Nowadays, even though NPAs are subject to a conservation protocol, these areas face several challenges in Mexico, such as mining, deforestation, and agriculture, which cause adverse effects 20,21 . Recent studies by Armendáriz-Villegas et al. 11 in Mexican NPAs, including the Sierra Gorda Biosphere Reserve (SGBR), have identified industrial and artisanal mining activities (exploration and exploitation) as important disturbance factors in Mexican NPAs, including the Sierra Gorda Biosphere Reserve (SGBR). There are currently more than 140 mines in the region of Sierra Gorda (in Central México), 60% of which occur on the SGBR polygon (most of them as artisanal mining) 22 ; however, data of impacts from mining affecting aquatic life in NPAs are currently unknown. Neither the detection of metals in water nor the impact of metals on aquatic ecosystems is systematically monitored in Mexico. Furthermore, biomonitoring is not part of the government programs in Mexico for environmental protection to date. Macroinvertebrates are recognized as the most suitable organisms for biomonitoring as they have low mobility, are in contact with sediments and the water column being exposed to pollutants in both compartments, and display a wide range of tolerance to contaminants (including highly sensitive species to very tolerant to polluted conditions) 23,24 . Thus environmental disturbances can modify the structure of macroinvertebrate assemblages in response to several stressors. The modifications include changes in species and abundance in impacted areas with a predominance of tolerant species, while sensitive species occur only in the most negligible impact or unimpaired conditions 25 . Nevertheless, the composition of aquatic macroinvertebrate communities in Mexican NPAs has been little studied, and scarce information is available on the habitat requirements of the aquatic macroinvertebrates or their particular responses to different stressors. Some experimental studies about the impact of heavy metals on the distribution and diversity of macroinvertebrates have been carried out [26][27][28][29] and, in some cases, community structure was included 27,30,31 . However, additional information is needed to fully understand the response of macroinvertebrates to heavy metals and other minerals from the geological strata to assess the disturbances caused by human impacts on aquatic wildlife, especially in an NPA affected by mining activities 32 . This study aims to identify how anthropogenic pollutants, particularly those associated with mining (metals), agriculture (nutrients), and villages (wastewater), affect aquatic macroinvertebrate communities in rivers with a high mineral content running across an NPA in Central Mexico. We first assessed the biological diversity of macroinvertebrates inhabiting the rivers in the study area and then identified the environmental drivers that determine the composition and abundance of aquatic macroinvertebrates. As the study area is affected by mining activities, we emphasized the evaluation of heavy metals and taxa identification that might indicate conditions of inorganic pollution. In addition, we examined the impact of organic pollution on aquatic macroinvertebrates.

Aquatic Macroinvertebrate Communities
A total of 77 000 specimens in 93 families, 23 orders, and six phyla were identified. Taxonomic richness ( 0 D) was lowest (37 ± 4 families) in the Extoraz River (sites PB, EP RQ, BC) and highest (50 ± 1.7 taxa) in the Escanela-Jalpan and Santa Maria rivers ( Fig. 1 and Supplementary Table S1). There were statistically significant differences between the Extoraz and Escanela-Jalpan rivers (Tukey HSD; p ≤ 0.05). The exponential of Shannon's diversity index ( 1 D) showed the lowest values in the Extoraz (except for site PB, where diversity was high) and Escanela rivers (sites EN, ES, and particularly site AH, where the lowest value was recorded). In contrast, the highest values were recorded in the Jalpan and Santa Maria rivers (sites PI, JL, PA, AY, and SM) ( Fig. 1 and Supplementary Table S1). The dominance index (the inverse of the Gini-Simpson's 2 D) was lowest in the Escanela River and the Conca River, where dominance is high. The Extoraz River showed medium values of 2 D (in contrast to its low 0 D values), denoting a moderate dominance of a few taxa (Diptera). High 2 D values were observed at site SM (Santa Maria River), where the abundance of the different taxa showed a high evenness (see Fig. 1 and Supplementary  Table S1). The Shannon's diversity index (H') showed minor differences between sites (2.25 ± 0.42). The lowest diversity values (1.63 ± 0.5) were recorded in the Escanela River (sites ES, EN, and, particularly, site AH) ( Fig. 1 and Supplementary  Table S1) and the highest (2.67 ± 0.1) in the Santa Maria River (sites SM, AT).

Environmental Factors
A few environmental factors showed a high degree of collinearity (Spearman's rho ≥ 0.80, Supplementary Table S2). The variables total coliforms, color, electrical conductivity, sulfates, total suspended solids, turbidity, and iron were excluded from further analyses because they correlated with fecal coliforms, ammonium, salinity, and aluminum, which were kept for subsequent analysis. The decision to keep these variables was based on the assumption that these factors directly affect benthic invertebrates than the others. As shown by their hardness, alkalinity, chloride, and salinity values (Table 1), some streams (Extoraz and Concá) are highly mineralized. A high concentration of fecal coliforms and high oxygen demand (BOD 5 ) revealed domestic wastewater pollution in the Escanela-Jalpan river. Other sites (e.g., sites on Extoraz river, VC and CN) are impacted by nutrients from various sources. Overall, all study sites were well-oxygenated (Table 1). High concentrations of heavy metals were found in the rivers of SGBR. The order of predominance of heavy metals in the rivers was observed as follows: Extoraz basin Al > Fe > Zn > Cd > As > Hg > Mn > Sb > Co > Cr > Cu, Escanela-Jalpan Zn > Fe > Al > As > Cd > Hg > Sb > Cu > Cr > Co > Mn, Ayutla basin As > Fe > Al > Zn > Sb > Hg > Cd > Co > Cr > Cu > Mn, Concá basin Al > Fe > Hg > Zn > As > Cd > Sb > Co > Cr > Cu > Mn, and Santa María basin Al > Fe > Zn > Cd > As > Sb > Hg > Co > Cr > Cu > Mn.

Ordination of environmental factors
The biplot of the co-inertia analysis showed a clear ordination of study sites ( Fig. 2b and Fig. 2k); three groups of sites can be clearly differentiated: 1) sites ES, EN, AH, and PI at the Escanela-Jalpan River, characterized by high BOD 5 and fecal coliforms due to the input of domestic wastewater, in addition to high concentrations of Zn, pH, dissolved oxygen, and the best habitat quality; 2) sites PB, RQ, EP, and BC (at the Extoraz River), and VC and CN (Conca and Santa María River), characterized by high hardness, alkalinity, salinity, total nitrogen, ammonia, nitrates, nitrites, and orthophosphates, indicative of high mineralization and nutrient enrichment, in addition to high concentrations of Al, Cd; and 3) sites AY, SM, AT, PA, and JL (at the Ayutla River and the Santa Maria River, and one site at the Jalpan River), characterized by high concentrations of TP, Hg, As, water discharge, and temperature (air and water). Particularly, Jalpan (JL) River was related to the higher concentrations of Cr, Sb, and Cu ( Fig. 2b and Fig. 2k).

Identification of environmental drivers of macroinvertebrate community structure
The MCOA analysis revealed remarkable differences between the rivers studied due to their co-structure of richness and abundance of the macroinvertebrate community and the environmental variables. The first two PCA axes accounted for 37.33% and 18.31% of total inertia. Most taxonomic groups (except for Lepidoptera) were highly correlated with diverse environmental factors; RV values ranged from 0.64 to 0.80 ( Fig. 2a and Fig. 2c-k). Trichoptera and Coleoptera showed the highest correlations with some environmental factors (Fig. 2g, i). Study sites were arranged according to biological and environmental features along the F1 and F2 axes (Fig. 2b). The relationships between environmental factors and taxa are shown separately for each taxonomic order of aquatic macroinvertebrates (including its various families), plus a group of miscellaneous taxa including those orders represented by a single family (See Appendix for Family codes). Taxa in the miscellaneous group were highly correlated with heavy metals (i.e., Hg, Cd, Cr, and Al), TP, water discharge, and temperature ( Fig. 2c, k). Highest concentrations of heavy metals such as Cr, Sb, and Cu were correlated with Perlidae (A21), Planariidae (A1), and Hirudinidae (A3); the metalloid As and the heavy metal Hg were correlated with several mollusk families (particularly Thiaridae, A10) (identified with arrows in Fig. 2c, k). High values of DO, BOD 5 , fecal coliforms, and habitat quality were highly correlated with Trombidiformes (A17); high concentrations of several nutrients and heavy metals such as Al and Cd were associated with the occurrence of Lumbriculidae (A5). The Ephemeroptera families, Heptageniidae (B3), and Baetidae (B2), were associated with well-oxygenated water and high BOD 5 , fecal coliforms, and pH. The family Leptophlebiidae (B5) was associated with high concentrations of Cr, Leptohyphidae (B4) was correlated with high concentrations of nutrients, Cd, and Al. The Odonata families showed a high affinity for heavy metals. A high nutrient concentration, Al, and Cd were tolerated by Gomphidae (C4), and Libellulidae (C5) showed tolerance to high BOD 5 and fecal coliform concentrations. Coenagrionidae (C2) and Platystictidae (C6) showed a relationship to Cr and As (Fig. 2e, k). The Hemiptera families (particularly Naucoridae, D10) were related to high Cd, Al and nutrient concentrations (Fig. 2f, k). The occurrence of Trichopterans such as Hydropsychidae (E7), Hydrobiosidae (E6), and Philopotamidae (E11) was correlated with Cr, As, and Hg. High oxygen concentration, good habitat quality, and fecal coliforms were related to Xiphocentronidae (E13) and Policentropodidae (E12) (Fig. 2g, k). Lepidoptera families Crambidae (F1) and Pyralidae (F3) were highly related to rivers with Cr, As, and Hg and high water discharge ( Fig.  2h, k). High heavy metal concentrations (Cr, As, Hg, Cd, and Al) were tolerated by some Coleoptera families, especially Psephenidae (G13), Dytiscidae (G4), and Lampyridae (G11). Well-oxygenated sites with good habitat quality were related to Dryopidae (G3), Hydrophilidae (G10), and most Coleoptera families (Fig. 2i, k). Most dipterans occurred in well-oxygenated sites; however, the Chironomidae family (H2) was related to the highest concentration of fecal coliforms and high BOD 5 . High concentration of Cr was associated with Empididae (H6), while As and Hg with Culicidae (H3) (Fig. 2j, k).

Discussion
Our study aimed to identify the environmental factors and heavy metals that determine the MA community structure in rivers with significant influence by mining activities located in their upper reaches. A set of environmental factors was examined to find their correlations with biological data; to this end, we used the multivariate analysis MCO proposed by Dolédec and Chessel 33 , which helped elucidate the main environmental factors influencing the biological composition structure. The MCOA combines separate ordinations into a single analysis based on the cross-covariance matrix 34 . We first identified and removed the environmental variables that were redundant to improve the analysis. Thus, any potential biases in our results were avoided by selecting a set of variables of similar prediction efficacy 35,36 . This study highlights the high diversity of aquatic macroinvertebrates inhabiting the SGBR (see Appendix). High diversity values were found mainly in rivers located in the central and northwestern parts of the SGBR (Escanela-Jalpan and Santa Maria rivers). These rivers overall show suitable conditions for establishing macroinvertebrates, i.e., good habitat quality and differences in vegetation between the river reaches. However, the wastewater treatment facility (WWTF) located between sites EN and AH on the Escanela river produces a drop in diversity ( 0 D, 1 D, and Shannon entropy H') that can be attributed to unsuitable conditions caused by habitat deterioration and unfavorable water quality affecting the establishment of various 3/15 macroinvertebrates. According to Vinson 37 , factors including habitat fragmentation, changes in water discharge, temperature, and poor water quality are strong environmental disturbances for establishing aquatic macroinvertebrates in rivers. Water discharge is reduced (even interrupted, causing fragmentation) downstream of the WWTF, and, therefore, the community of macroinvertebrates showed a marked reduction in taxa richness 0 D and other diversity measures ( 1 D and Shannon entropy H'). Macroinvertebrate diversity is likely associated with several factors, among them, is the vegetation. The Escanela-Jalpan and Santa Maria basins are dominated by pine forest and tropical rainforest, respectively, both basins reached high diversity measures. In contrast, the Extoraz basin is largely covered by xerophytic vegetation. The Extoraz River showed the lowest diversity ( 0 D and 1 D) and medium dominance ( 2 D) in the entire SGBR. Despite the relatively low diversity in the Extoraz River, we identified 32 different taxa (mostly at the family level); by comparison, Rocha et al. 38 found 25 macroinvertebrate families in similarly arid conditions in Brazil. Studies carried out by Bogan et al. 39 in the southern United States of North America recorded 148 taxa (most of them at species level) across a network of six rivers in dry ecosystems under similar conditions to those in the Extoraz River. Several authors have pointed out the importance of riparian vegetation and mature trees to provide cover and habitat for the establishment of fish and macroinvertebrates 40,41 . Such is the case of the Escanela-Jalpan and Santa Maria Rivers, contrasting with the Extoraz River, where large mature trees are scarce. The diversity numbers ( 1 D and 2 D), showed an increasing trend downstream as far as the Santa Maria River. According to Vannotes's concept 42 and Malmqvist and Hoffsten 43 aquatic diversity tends to be higher in downstream reaches due to the tributaries that flow into the mainstream along its gradient and the heterogeneity of the habitat. We found marked differences between Shannon's diversity index (H') values and Hill numbers. The Shannon's index was not as sensitive to small differences in diversity between and within the rivers studied as the Hill numbers were. As Shannon's index is based on the amount of information provided by the identity of an individual chosen at random, rare or uncommon species are not expected. Jost et al. 44 mentioned that linear diversity calculations (as Shannon's index) are not entirely suitable to determine diversity. For this reason, we calculated Hill numbers to estimate the "effective number of species". The corrected diversity values provided by Hill numbers showed a pattern of increasing diversity downstream, which was most evident in the Concá and Ayutla rivers. In most cases, community diversity is determined by a significant environmental component 45 . In our study, environmental factors such as heavy metals, BOD 5 , oxygen, and nutrients act as drivers of the macroinvertebrate community structure. We assumed that the contrasting diversity (dominance and taxa richness) between the Extoraz and Escanela Rivers is defined by a set of environmental factors and the tolerance of some taxa to the presence of heavy metals. We found high concentrations of heavy metals (particularly Cd and Hg) in the sites studied; in some cases, these concentrations exceeded the aquatic life criteria set by the US EPA 46 . These results reflect the effects of mining activities inside the SGBR, which is the most common human activity in the upper reaches of the rivers of SGBR 47 . According to Robles et al. 48 , cinnabar (HgS) attains high concentrations due to illegal mining. The SGRB has been historically referred to as a well-mineralized place, mainly by cinnabar (HgS), which is the most common ore in the SGBR and the primary mercury-containing mineral. Mercury production in the area has attained a high volume in the past century through illegal mining 49,50 . Then, erosion processes may lead to high mercury concentration in water (maximum values of 0.049 mg L −1 in the present study). As per the criteria for aquatic life set by the US EPA 46 , the maximum mercury concentration in saltwater is 0.0018 mg L −1 . However, concentrations recorded in SGBR water were twice as high. Thus, the SGBR is at imminent risk of heavy metal pollution, particularly by mercury, as Hernández-Silva et al. 49 have reported for maize crops in the SGBR. Cadmium is another metal whose levels in SGBR rivers are twice the limit (0.00025 mg L −1 ) set by the US EPA 46 . Other heavy metals (Al, As, Cr, Cu, and Sb) measured in this study did not exceed the aquatic life thresholds. Nevertheless, their concentration might increase in the future if erosion continues at the current rate. Our analyses revealed that a well-mineralized environment prevails in the rivers of the SGBR. High CaCO 3 concentrations, evidenced by alkalinity and hardness (Table 1), are indicative of a particular geological stratum. Moreover, soil types, including Regosol, Vertisol, Litosol, Rendzina, and Cambisol, were identified by Carrillo-Martínez and Suter-Cargneluti 51 . These findings suggest that long-term erosion of soluble rocks naturally occurs in the area 47 . We found several forms of nitrogen directly linked to current human activities in the area, such as urbanization and agriculture. Mostly rural settlements occur in the SGBR (population density is currently 25 persons per km 2 ), in which sewage treatment is scarce, and only three treatment plants are available for municipal wastewater from the nearest towns. However, wastewater discharges from the larger villages (as in the Escanela-Jalpan River) are starting to affect the nearby rivers. Total nitrogen and total phosphorus reached concentrations higher than 4 mg L −1 and 0.7 mg L −1 , respectively, usually associated with extensive deforestation and agriculture in impacted basins 52,53 . Khatri and Tyagi 54 have discussed the impacts of rural and urban sewage, where farming and runoff from human settlements are the primary factors increasing the input of coliforms and organic compounds into nearby stream water. We identified villages such as Ahuacatlan (AH) and Jalpan (JL), located in the Escanela River, as well as river modifications such as the reservoir built on the course of the Jalpan River, that may cause adverse effects on the ecosystem 52, 55 . The MCOA analysis revealed a significant strong relationship between key environmental factors and macroinvertebrate families pooled in groups in a spatial dimension. The MCOA yielded RV values for each taxonomic group, which were generally high (RV from 0.46 to 0.74). Dalu et al. 56 examined approximately 5 000 specimens in short-term studies and found RV values that rarely reached a maximum of 0.30. We examined more than 70 000 specimens that made it possible to achieve higher correlation values between environmental factors and macroinvertebrate community. Overall, we identified three different conditions for aquatic life in the SGBR: a) water with high concentrations of heavy metals caused by mining activities, b) water with high nutrient concentrations from villages and primary activities carried out in the area, and c) well-oxygenated water with high numbers of coliforms (Fig. 2k). Using the MCOA, we were able to correlate particular macroinvertebrate families with each of the three water quality conditions identified in the SGBR. Several families in the miscellaneous group of aquatic macroinvertebrates were related to heavy metals; the vectors representing most of them appear on the upper part of the PCA biplot (Fig. 2c-i) and show a high affinity for heavy metals. Groups such as Plecoptera, mollusks, and worms showed a close relationship with high concentrations of heavy metals (Cu, Cr, Sb, As, and Hg). In a study in high mountain streams of the Gangqu River, China, Qu et al. 57 found an inverse relationship between heavy metal concentration and diversity of macroinvertebrate community, particularly, affecting those belonging to Ephemeroptera-Plecoptera-Trichoptera taxa (EPT). Other studies have shown an association between aluminum in the water column and Perlidae such as the genus Anacroneuria 26, 58 ; however, we found a different response in the SGBR. The family Perlidae showed sensitivity to aluminum and highnutrient water regardless of the concentration of other heavy metals, including mercury. Other groups such as mollusks, worms, and flatworms were highly correlated with high concentrations of heavy metals. Ankley 59 and Croteau et al. 60 found bioaccumulation of Cd, Cu, Pb, Ni, and Zn from sediments in several mollusks, Perlidae, and worms, confirming that those taxa are tolerant to high concentrations of heavy metals. We found Leptophlebiidae and Psephenidae living in water with Cd, Hg, and Al. High Cu, Cr, Cd, and As concentrations have been related to high abundances of some Ephemeropterans (Leptophlebiidae), and Coleopterans (Psephenidae) 61,62 . Therefore, we might consider these taxa as tolerant to heavy metals. Some Odonata (Coenagrionidae and Platystictidae) were also related to high concentrations of heavy metals. Studies by Corbi et al. 63 and Michailova et al. 28 suggest that some Odonata taxa are tolerant to heavy metals in water. Our results are consistent with those studies, as we found a high correlation between Odonata (Coenagrionidae and Platystictidae) with Cd and Al. We found similar responses in some Hemiptera and Trichoptera. However, Hemipterans are air-breathing, skate on the water surface, and, eventually, leave the aquatic environment 64 , conditions that, taken together, allow them to live in polluted water. Trichoptera belongs to the EPT group, considered as less tolerant to pollution than non-EPT insects 65 . Nonetheless, some Trichoptera families, e.g., Hydropsychidae and Philopotamidae, have been recognized as tolerant to heavy metals in water 29,30,57 . We found several families of Trichoptera, particularly Hydropsychidae, Philopotamidae, and Hydrobiosidae, that showed a high correlation with As and Hg. In contrast, Lepidopterans did not reach a sufficient weight to relate their occurrence to heavy metals pollution, only the family Crambidae showed tolerance to heavy metals. Ephemeroptera such as Baetidae and Heptageniidae showed a high affinity for well-oxygenated water, while Leptophlebiidae displayed high tolerance to heavy metals (As and Hg) (Fig. 2d). Several authors mention that the Ephemeroptera life cycle depends on water temperature in running water 66,67 , and that they prefer cool well-oxygenated water. As Jacobsen et al. 68,69 found in previous studies, this condition is generally true for Ephemeroptera, Plecoptera, and Trichoptera. Overall we observed the same pattern in some Ephemeroptera, but not in all of them. Heptageniidae and Baetidae are highly sensitive to the presence of heavy metals in water, but Baetidae tolerates water with coliforms. Those responses have also been observed in several studies under laboratory conditions; for example, Clements 27 and Courtney and Clements 70 demonstrated that Heptageniidae is highly sensitive to heavy metals such as Zn, Cd, and Cu. Overall, Ephemeroptera is sensitive to heavy metals; however, Baetidae and Leptophlebiidae are the most diverse families of mayflies, and their responses are equally diverse 71 . Baetidae (as the genus Baetis) is a family well-known for tolerating little to highly polluted water 66 . Although Baetids and Heptageniids tolerate human impacts by coliform bacteria, they have the potential to be good indicators of water quality, particularly as indicators of the impact from mining activity in the SGBR. The families Lumbriculidae, Leptohyphidae, Gomphidae, and Naucoridae were closely related to nutrient-enriched water; also, Oligochaeta includes nutrient-tolerant taxa; Ristau et al. 72 demonstrated changes in their density driven by phosphorous concentration peaks. We found the same pattern with Leptohyphidae, which showed tolerance to water containing high nutrient concentrations.
This study reveals high metal concentrations in the studied rivers of the SGBR; however, a high diversity of aquatic macroinvertebrates was also detected. The geologic strata of the reserve (mainly limestone shale and dolomite limestone) contribute to the high hardness and alkalinity in water. Both factors are likely to modify the bioavailability and the toxic effect of metals on macroinvertebrates present in the rivers. A buffer effect is currently influencing due to high concentrations of carbonates (CO 3 − and HCO 3 − ) in the study area, forming soluble or insoluble complexes and/or precipitates 73 . Thus, higher water alkalinity decreases the toxicity of metal ions either by active surface competition for binding sites in tissues 74 or by reducing their concentration through the formation of insoluble precipitates. Furthermore, it is likely that due to ions Ca 2+ and Mg 2+ , in the form of carbonates, may both compete with other divalent metal ions for organisms' binding sites 75 . Carbonate functions as a blocker to the entry of metals into the organisms; hence, the alkalinity and hardness of the SGBR rivers are likely acting synergistically decreasing the availability and toxicity of heavy metals and thus protecting the aquatic biota. However, this effect was less evident in the Extoraz River, which has xerophytic vegetation, and also is the zone with the higher number of 5/15 mines in their upper reaches and higher salinity, all of which could contribute to the lower diversity of this river. The results of this study using the MCOA revealed diverse macroinvertebrate assemblages that can be ranked along a stressor gradient (heavy metals concentration) ranging from families that are highly sensitive to others that are tolerant to high metal concentrations. This allowed us to identify taxa that might be used as indicators of the ecological health of rivers influenced by mining activities. Our study is one of the first efforts in Latin America to characterize the response of macroinvertebrates in support of their use as biomonitoring tools in natural protected areas under mining impact.

Study area
The SGBR is one of the most important NPAs in Mexico, with a surface area of 3 800 km 276 , located in the Central Mexican Plateau, between Neartic and Neotropical biogeographic regions 77 . This NPA is located in a wider mountain range called "Sierra Gorda" with two major rivers: the Santa Maria river in the northeast and the Extoraz river in the south. Both rivers flow into the Panuco River (to the Gulf of Mexico). Data were collected at 15 sites across the SGBR, on the Extoraz and Santa Maria rivers and the tributaries (Conca, Ayutla, and Escanela-Jalpan rivers) (Fig. 3). Two main human settlements in the SGBR (Ahuacatlán and Jalpan towns) are located along the Escanela-Jalpan River; a wastewater treatment facility is located on the periphery of Ahuacatlan; the river then flows into a reservoir before joining the Santa Maria River. Several small towns with some 93,000 inhabitants occur in the SGBR; small parts of the area are used as cropland 78 . The main geological strata in the SGBR is composed mostly of limestone 22 , shale and dolomite. Sierra Gorda is home to one of the most diverse vegetation types among Mexican protected areas. The main vegetation types in the SGBR are pine forest (on the Escanela river), tropical rainforest, oak forest, deciduous tropical forest (on the Santa Maria river), and xerophytic shrubland (on the Extoraz river). The climate ranges from semi-warm to warm and sub-humid to semi-dry with an annual mean temperature mostly above 18 • C. Precipitation averages 313 mm in the dry season and 883 mm in the summer rainy season 79 . Currently, 140 mines are located in Sierra Gorda: 83 inside the SGBR; however, all of them occur in the upper reaches of the five sub-basins that converge into the mainstreams (Extoraz, Escanela, Jalpan, Conca, Ayutla, and Santa María rivers). The influence of mines extends to the entire basin downstream from the mining area where waste occurs ( Fig. 3b and 3c). Most mines provide minerals such as Au (46%), Ag (27%), Hg (8%), Pb (5%), and Sn (2%), and non-mineral materials (Sb, barite, fluorite, phosphorite, marble, and gypsum, 12%). The total of mines located occur as follows: Extoraz sub-basin 46, Escanela Jalpan sub-basin 40, Ayutla sub-basin 28, Concá sub-basin 1, and Santa María sub-basin 2. In most cases, the exploitation is artisanal 22 .

Environmental Factors
The geographic coordinates and elevation (in masl) of each sampling site were recorded with a GPS Garmin device. At each sampling site, environmental variables, water and air temperature ( • C), dissolved oxygen (mg L −1 ), pH, turbidity (NTU), conductivity (mS cm −1 ), salinity (PSU), were recorded (using a Quanta probe); water samples were collected twice in both the dry and the rainy seasons (February 2017 and January 2018 and July and October 2017, respectively). In all cases, seasonal averages were calculated for each sampling site. The concentrations of nutrients, major ions, and other physicochemical factors (Supplementary Table S3), were measured using HACH 80 and APHA 81 procedures. Heavy metal and metalloid concentrations were measured after microwave acid digestion (EPA 3015A) with the microwave oven Anton-Paar Multiwave Go and following the methods recommended by the relevant Mexican regulations 82 , using Inductively Coupled Plasma Optical Emission Spectrometry (ICP OES) (Perkin Elmer Optima 4300DV). Metal concentrations below the detection limit were arbitrarily assigned a value of half the corresponding detection limit 83 . In each reach, habitat quality -a summary description of the variety of habitats and their features -was assessed through habitat surveys, scoring physical features that included stream morphology, substratum, and riparian coverage, using an ordinal categorical scale. The scores were then used to compute a habitat score following the procedures described by Barbour et al. 84 ; habitat scores range from 0 to 20, with higher values indicating more heterogeneous habitats. River discharge was estimated according to Michaud and Wierenga 85 . For each environmental factor, the mean of four measurements was calculated.

Macroinvertebrates
Aquatic macroinvertebrates were sampled at each study site in February (dry season) and July (rainy season) 2017, using multi-habitat methods recommended by the United States Environmental Protection Agency (US-EPA), Barbour et al. 84 , and the AQEM 86 . Four 5-minute subsamples were collected at each site, then pooled into a 20-minute sample for each sampling site. Two subsamples were collected in riffle sections using a kick net, and the other two in the most dominant habitat using a scoop net; both nets had a 500 µm mesh size. Specimens were preserved in 70% alcohol for subsequent sorting and identification in the laboratory. These specimens were identified using specialized literature for North American and Neotropical areas: Merrit and Cummins 87 , Thorp and Covich 88 , Bueno-Soria 89 , Springer et al. 90 , and Hamada et al. 91 . Taxa were primarily determined to family level. Taxa from the two seasons were pooled together into a single list to compute mean values for each study site.

Data Analyses
The diversity of macroinvertebrates at each study site was evaluated using four biodiversity indices. The Hill diversity numbers or "true diversities" and the Shannon entropy index (H') were calculated using the methods described by Jost et al. 92 . In this method, community diversity measures are converted into "effective number of species" and three Hill diversity numbers can be assessed at different orders (q); the parameter q controls the sensitivity to common and rare species. When q = 0, abundance values are raised to the power of 0, rare and abundant species have the same weight; thus, the Hill diversity number of order 0 ( 0 D) represents species richness. When q = 1, abundance values are raised to the power of 1, and the Hill diversity number of order 1 ( 1 D) represents the exponential of Shannon's diversity index. When q = 2, abundance values are raised to the power of 2, thus increasing the weight of dominant species, and the Hill diversity number of order 2 ( 2 D) represents the inverse of the Gini-Simpson's dominance index 92 . A distinctive advantage of this approach is that all diversity numbers can be interpreted as effective numbers of species. Diversity profiles for each river were summarized as heatmaps. Diversity indices were computed with the package iNEXT for R 3.1.0. A one-way ANOVA was used to test for differences in diversity indices between rivers; the data were tested for normality and homoscedasticity before analysis. Data on environmental factors (water and habitat quality) were tested for collinearity using the Spearman's rank correlation coefficient in order to exclude redundant variables from further analysis. The relationship between environmental factors and aquatic macroinvertebrates at the different study sites was examined. First, an environmental data matrix containing the mean values of each of the 25 environmental factors recorded at the 15 study sites was compiled. Macroinvertebrate abundance data were log-transformed, i.e., ln (x+1), to reduce the effect of dominant taxa 19 and compiled in a species-X-sites matrix. The two matrices were then subjected to a multiple co-inertia analysis (MCOA) (999 permutations) following Dolédec and Chessel 33 . We perform a multi-table STATIS analysis to identify the relationship of taxa (orders and a group of miscellaneous families) and environmental factors on a spatial scale, according to Lavit et al. 93 . Finally, to evaluate the significance of the relationships between taxa and environmental factors, we calculated RV (Rho values). These analyses were carried out using the packages ADE4 22 and Vegan 2.0. 10     Results of the MCOA: a) inter-structured correlation circle; b) ordination of study sites on the plane de ned by the rst two axes; c) to j) contribution (vectors) of aquatic macroinvertebrate taxa (for easier visualization, separate biplots are shown for each taxonomic order); and k) contribution (vectors) of environmental factors. See Appendix for taxa codes. Note: The designations employed and the presentation of the material on this map do not imply the expression of any opinion whatsoever on the part of Research Square concerning the legal status of any country, territory, city or area or of its authorities, or concerning the delimitation of its frontiers or boundaries. This map has been provided by the authors. the part of Research Square concerning the legal status of any country, territory, city or area or of its authorities, or concerning the delimitation of its frontiers or boundaries. This map has been provided by the authors.

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download.