Drivers of Benthic Macroinvertebrate Assemblages in Equatorial Alpine Rivers of the Rwenzoris (Uganda)

: The Sub-Saharan alpine freshwater biodiversity is currently impacted by human settlements, climate change, agriculture, and mining activities. Because of the limited biodiversity studies in the region, a better understanding is needed of the important environmental variables a ﬀ ecting macroinvertebrate assemblages. In this paper, macroinvertebrate diversity responses to 18 environmental variables were studied at 30 sites along unique Rwenzori rivers at the equator in Uganda. We hypothesized that anthropogenic disturbance and local environmental variables a ﬀ ect macroinvertebrate diversity, irrespective of altitudinal gradients. Based on altitude and climate, the sites were subdivided into three altitude groups consisting of 10 sites each: upstream (US) 1400–1600 m.a.s.l.; midstream (MS) 1091–1399 m.a.s.l., and downstream (DS) 900–1090 m.a.s.l. A total of 44 macroinvertebrate families and 1623 individuals were identiﬁed. The macroinvertebrate diversity patterns were inﬂuenced by temperature, altitude, and latitude. Regression analysis revealed that temperature and nickel, were negative predictors of taxa richness. Nickel, which is released by mining activity, is detrimental to aquatic communities in Sub-Saharan alpine ecosystems. Signiﬁcant longitudinal variation in macroinvertebrate diversity was observed between the sites, which were also a ﬀ ected by mineral and temperature gradients. Our study highlights the need for long-term monitoring in this region to detect and reduce the threats to river biodiversity from anthropogenic activity.


Introduction
The global biodiversity is increasingly under threat due to both man-made and natural stressors [1][2][3]. The Sub-Saharan African biomes (grasslands and forests) are greatly affected by biodiversity losses as a result of human development activities (i.e., mining, agriculture, and industry) and climatological changes (i.e., extreme precipitation events and increasing temperature) [4,5]. Sub-Saharan Africa contains many ecoregions, which are considered to be among the richest in global biodiversity [6][7][8]. One of these biodiverse Sub-ecoregions, the western branch of the East African Rift System (EARS), encompasses the central African Albertine Rift region, where the Rwenzori mountains are situated and which are part of the East African plateau mountain ranges [9][10][11]. The Rwenzori region has a national park that is a UNESCO world heritage site and is earmarked as one of the Global 200 freshwater ecoregions [12,13]. Rwenzori rivers and streams are part of the Lake Edward basin [14], which includes rivers of Nyawamba, Nyamugasani, Mubuku, and Mpanga that are mainly fed by both glacial meltwater and rainfall, and are vital headwaters for the White Nile [15]. These rivers and glacial contribution to stream flow. Due to the influence of lakes, tributaries, and groundwater upwelling, the glacier influence on macroinvertebrate biodiversity reduces with increasing downstream distance from the ice margin [44]. As groundwater contributions increase, range conditions become more stable, and environmental stress is decreased [45]. In alpine systems [45], along with other mountainous regions [46], groundwater has been identified as a moderator of water temperature. The macroinvertebrate community diversity is expected to increase as glaciers retreat, due to upstream taxa migration as a result of increased water temperature caused by reduced glacial influence. Likewise, this shapes the macroinvertebrate distribution and drives the community structure of these alpine ecosystems. Previous longitudinal studies on glacier-fed streams in the Ecuadorian Andes showed that stream temperature and channel stability explained most of the variability in macroinvertebrate taxa composition and richness [41,42,47,48]. Channel stability systematically increased downstream while temperature did not, and the researchers further observed that macroinvertebrate diversity increased with distance from the glacier.
Rivers and streams are also subject to disturbance which can affect taxa richness. The intermediate disturbance hypothesis stipulates that taxa richness is greatest in communities with disturbances that are moderate in both their frequency and intensity [44,49]. Following this hypothesis, both taxon richness [50] and assemblage evenness [51,52] are expected to peak at intermediate levels of disturbance (or stability). The macroinvertebrate families present in a community depend on the physical characteristics, disturbance levels, and biotic factors of the ecosystem. Global studies have further revealed that physicochemical variables (temperature and channel stability) and geographical variables (altitude and longitude) influence macroinvertebrate community patterns in both tropical and temperate aquatic systems [53][54][55][56][57][58]. These studies revealed that temperature influences macroinvertebrate diversity along an altitudinal gradient. There have only been a few studies undertaken on the diversity of macroinvertebrates, with respect to both geographical and anthropogenic disturbance gradients in the study region and Uganda [59][60][61][62]. This study aims at filling this knowledge gap about the Rwenzori alpine rivers, with regard to the importance of both environmental and geographical variables on the macroinvertebrate community. We aimed to study drivers of macroinvertebrate diversity along sites at different altitude ranges with varying levels of anthropogenic disturbance.
To this aim, we sampled macroinvertebrates along rivers and streams on both the leeward and windward sections of the Rwenzoris. Our main goals were: (i) to determine the variables driving the overall macroinvertebrate community diversity of Rwenzori rivers, (ii) to assess the structural variability in faunal univariate diversity metrics between different altitude classes, and (iii) to determine the importance of the site-specific environmental variables on the macroinvertebrate community composition. We hypothesized that: (a) anthropogenic disturbance and local environmental variables drive community diversity irrespective of the geographical altitudinal and longitudinal gradients, and (b) the differences in macroinvertebrate assemblages between altitude groups are minor and show similar relationships with environmental variables.

Study Area
The study rivers and streams are pristine headwaters of the Albert (White) Nile and mainly originate from the glaciers of the Rwenzori Mountain ranges (00 • 23 09 N 29 • 52 18 E). They traverse the western rift valley lowlands and subsequently drain into the wetlands of Lake Edward and Lake George, the latter being a protected Ramsar site of international importance [63]. These mountain ranges are located in southwestern Uganda, in Bundibugyo, Kasese, and Kabarole districts along the border of Uganda and the Democratic Republic of Congo. The study area bridges the equator in the western branch of the East African Rift System, and is bounded by Rwenzori mountains to the north, Lake Edward to the south, and Lake George to the southeast ( Figure 1). The altitude ranges between 900-5109 m above sea level (m.a.s.l.). The main land use activities include tourism, urban settlement, commercial agriculture, mining, and industry, which are conducted between the lowest plateau altitude level of 900 m.a.s.l. to 2000 m.a.s.l. The Rwenzori region is characterized by four climatic types along different altitude ranges (i.e., <1200 m.a.s.l.: Tropical savannah; 1200-2000 m.a.s.l.: equatorial; >2000 m.a.s.l.: temperate) with a snow line at 4500 m.a.s.l. [64]. It has well-defined rainy seasons, with an annual bimodal rainfall regime ranging from 750 mm to 1500 mm occurring from March to June and October to December. However, some highland areas in the region experience heavy rainfall of 1500 mm to 2000 mm all year round [64,65]. The dry seasons are characterized by sporadic rainfall which lasts for a few hours [11].

Site Selection
The 30 sites selected were situated between 900-1600 m.a.s.l. along accessible sections of 20 separate rivers and two Lake George wetland sites. The study area map was created using ArcGIS version 10.0 (ESRI Company, Redlands, CA, USA). The site choice was moderated by accessibility and was distributed across the Rwenzori region (leeward and windward sections) to span the entire altitude range where most anthropogenic activity is conducted (900-2000 m.a.s.l.), in the tropical savannah and equatorial climatic zones. To capture the full range of climatic zone variation using an altitudinal gradient, 10 sites were sampled in the tropical savannah zone (hereafter "downstream sites"), with an additional 10 sites sampled in the transition from tropical savannah to equatorial zones (hereafter "midstream sites"), and lastly 10 sites were sampled from the equatorial zone (hereafter "upstream sites") ( Figure 1).
Additionally, with ground truthing (i.e., direct observation in the field) selected sites had varying anthropogenic disturbance levels. Visible point pollution sources included discharges from municipal wastewater treatment plants, mine leachate from copper mine tailings, and waste from industries, and non-point pollution sources such as runoff from dense urban settlements, agricultural farmlands, and mine drainage. An estimate of the anthropogenic impact based on the land use types (i.e., urban area, mining, roads/paths, residential areas, agriculture, and water abstraction), within 100 m of each site, was visually assessed and scaled from 0 to 3 [66,67]. The four ordinal categories were high (3), moderate (2), low (1), and none (0), with higher values indicating high disturbance. The total disturbance was computed and scored on a scale of 0-10 and based on the scores the sites were then categorized into three levels of disturbance (low, medium, and high). These ranged from least disturbed sites (i.e., non-discernible disturbance) with a score between 0-3 to moderately disturbed (i.e., noticeable impact) with a score of 4-6 and lastly, highly disturbed sites (i.e., visible widespread significant disturbance) with a score of 7-10.
Study sites were both environmentally and geographically diverse and subject to different climatic conditions characterized by spatial and temporal variations in weather. The differences in macroinvertebrate assemblages between the sites could thus be attributed to differences in physicochemical and/or geographical variables. We sampled sites as synchronously as possible between January to February 2017 during the dry season when they were most accessible, and ecologically and hydrologically stable.

Collection of Environmental Data
We measured 15 physicochemical and three geographical variables along each of the 30 sites (see Table S1, in Supplementary Materials). The altitude and exact coordinates were determined using a handheld Global Positioning System (GPS) device (Garmin Legend; Garmin Ltd., Olathe, KS, USA). The physicochemical variables of pH, temperature (Temp), specific conductivity (SpCond), and dissolved oxygen (DO) were measured in the field using a Hydrolab multi-probe meter (H2O Hydrolab multi-probe; Hydrolab Company, Austin, Texas, USA). The current velocity at each site was measured using a handheld Höntzsch probe (HFA-model; Höntzsch, Waiblingen, Germany).
The water samples for chemical analysis were collected using nitric acid pre-washed Wheaton amber glass bottles (1 L). At each site, the sample bottles were washed three times prior to water collection. The glass bottles were carefully filled completely without creating bubbles (zero headspace), which was achieved by filling the bottle until there was a positive meniscus prior to capping and sealing. The bottled water samples were cooled on ice immediately after collection, and then transported in a Coleman cool box (Therapak; Coleman, USA) to the laboratory. The samples were immediately stored at 4 • C until analysis which was carried out within 24 h after collection, according to the World Health Organization (WHO) guidelines [68] and HACH (HACH, 2015) recommendations. These water samples were then analyzed in the laboratory using a HACH Lange DR3500 lab spectrophotometer (HACH Company, Loveland, CO, USA) and matching chemistry test kits to measure, (with their respective minimum to maximum ranges); copper (0.01-1.00 mg L −1 ), nickel (0.05-1.00 mg L −1 ), ammonium-nitrogen (0.015-2.0 mg L −1 NH 4 -N), nitrate-nitrogen (0.23-13.5 mg L −1 NO 3 -N), nitrite-nitrogen (0.005-0.100 mg L −1 NO 2 -N), total nitrogen (1-16 mgL −1 TN), (total and ortho) phosphate (0.01-0.50 and 0.05-1.50 mg L −1 PO 4 ), chemical oxygen demand (5-60 mg L −1 COD) and biochemical oxygen demand, 5 days (1.0-12.0 mg L −1 BOD 5 ).

Aquatic Macroinvertebrate Sampling and Identification
At each site, macroinvertebrates were sampled using standardized multihabitat kick sampling with a dip net (mesh size of 0.5 mm) along all accessible aquatic habitats within a stretch of approximately 30 m [69][70][71]. This included three replicate dip-net samples and visual collections from riverbed substrate (stones, sand, or mud), macrophytes (floating, submerged, or emergent), the immersed roots of overhanging trees, and debris [72]. The sampling effort was evenly distributed across all accessible aquatic habitats. One designated individual sampled each site for 5 min [70,71,73,74]. The three replicate samples per site were combined into one larger sample after collection. A total of 30 macroinvertebrate samples were collected during the daytime, rinsed, sorted, and preserved in 70% ethanol. In the laboratory, the samples were rinsed through a sieve and sorted without magnification. Subsampling was not applied. Using a stereo Olympus SZX10 microscope (Olympus Corporation, Shinjuku, Tokyo, Japan) with LED-light ring at a magnification of 7-45 ×, Rwenzori fauna were identified and counted in the laboratory to family level using the keys of Gerber and Gabriel [75].
Family level identification was considered to be the most suitable for future replicability and uniformity by water stakeholders in the region with limited resources and taxonomic expertise.

Data Analysis
All the analyses were performed using R software version 3.4.0 [76]. The patterns in diversity and environmental drivers of benthic macroinvertebrate assemblages were analyzed using non-parametric multivariate methods. Dufrêne's method was used to identify dominant families from each altitude group, whereby: where Yi is the dominant family, N is the total number of individuals in all sites, N i is the number of individuals of family i in all sites, and f i is the frequency of family i, which is the relative frequency of occurrence of family i in the sites of each group. We carried out indicator value (IndVal) analysis with the three altitude site group combinations to ascertain the association between taxa and the altitude groups using the indicspecies package [78]. In order to identify the group associations for each family, we set the alpha significance to 1 to display the group to which each taxon is associated with. Furthermore, in order to detect the significant indicator taxa for the altitude group demarcations, the alpha significance level was set to 0.5. The indicator value technique integrates approximates of abundance with the frequency of occurrence within a group [79]. A large value indicates that a taxon is abundant at all sites within a group and not present in other groups. The indicator value has two components A and B. Component A is the specificity or positive predictive value as an indicator of the site group while B is known as the fidelity or sensitivity of species as an indicator of the target site group. Component A is the conditional probability that the sampled site belongs to the target site group, in light of the fact that the species has been found. Component B is the conditional probability of finding the species in sites belonging to the altitudinal site group. A Monte Carlo randomization test with 1000 permutations was carried out to determine the index's statistical significance [80,81].
Four frequently used univariate diversity metrics, including taxa richness, total abundance, Pielou's evenness, and Simpson index values, were computed for each of the sites. The differences in diversity indices between the three groups were assessed using the non-parametric Kruskal-Wallis test [82] and Dunn's post hoc multiple comparison tests were run to assess significant differences in the diversity metric values among the altitude groups (p < 0.05).
Principal component analysis (PCA) was executed using the FactoMineR package [83], to determine the primary gradient of variables (predictors) in the data of study sites. We applied Kaiser-Guttmann's criterion and retained the factors that had an eigen value greater than 1. The first six factors were used in the regression analysis (see Table S2, in Supplementary Materials). Using both glmnet [84] and plotmo packages, separate regression models using the least absolute shrinkage and selection operator (LASSO) penalty algorithm [85,86], were run for both total taxa abundance and taxa richness with the 18 predictors from first six dimensions for the three altitude groups. As a shrinkage and variable selection technique for linear regression models, the aim of using the LASSO algorithm was to obtain the subset of predictors that minimizes prediction error for the two quantitative response variables (i.e., taxa abundance and taxa richness). The LASSO algorithm did this by imposing a constraint on the model parameters that causes regression coefficients for some variables to shrink towards zero. Predictors with regression coefficients equal to zero after the shrinkage process, were excluded from the model. Significant predictors with non-zero regression coefficients were retained, as they were strongly associated with the two response variables. Where necessary, prior to analysis, environmental variables were standardized to zero mean and unit variance, to remove disparate measurement scales of the variables.

Macroinvertebrate Community Diversity
A total of 1623 individuals were collected (see Table S3, in Supplementary Materials). We identified 44 macroinvertebrate families, of which Caenidae were the most common family with the taxon recorded at 50% of the sites. The most abundant taxa, constituting 67% of the total individuals identified, were: Simuliidae (26%), Baetidae (14%), Chironomidae (14%), and Caenidae (13%). The midstream (MS) sites had the highest total abundance (793 individuals), with downstream (DS) and upstream (US) sites having lower abundance scores (573 and 257 individuals, respectively). The US sites had the highest mean taxa richness scores (7.5 ± 2.5 SD) with DS and MS sites having lower mean taxa richness scores (5.4 ± 2.3 SD and 5.1 ± 1.9 SD, respectively) ( Figure 2).    Taxa richness ranged from 2-10 taxa for DS sites, 1-7 taxa for MS sites and 5-13 taxa for US sites. The total number of families identified in the three altitude groups were 31 (US), 26 (MS), and 27 (DS). Only 14 families (32% of total families) were common among the three groups. Caenidae and Chironomidae were the most common families, while 17 families (39% of total families) were group-specific. The highest number of group-specific families was detected in US sites (9 out of 44 families, 20% of total families). The lowest number of group-specific families was detected in the MS group (4 out of 44 families, 9% of total families), which included Ecnomidae, Gerridae, Lumbricidae, and Perlidae. Taxa richness, the Simpson index, and the Pielou's evenness diversity index values significantly differed between the three altitude groups and increased along an altitude gradient from DS sites to US sites. The post hoc Dunn's test identified significant differences between the four diversity metrics across the altitude groups (

Indicator Taxa Analysis
Multilevel pattern analysis selected 37 out of 44 families for indicator taxa site group association. As a result, 22 families were associated with only one of the three altitude groups and 15 families were associated with a combination of two altitude groups (Table 1). However, 7 families (Corduliidae, Heptageniidae, Libellulidae, Lymnaeidae, Naucoridae, Nepidae, and Teloganodidae) were not included in the analysis, as these families occurred at sites belonging to all three altitude groups.

291
Multilevel pattern analysis selected 37 out of 44 families for indicator taxa site group association.

292
As a result, 22 families were associated with only one of the three altitude groups and 15 families 293 were associated with a combination of two altitude groups (Table 1). However, 7 families 294 (Corduliidae, Heptageniidae, Libellulidae, Lymnaeidae, Naucoridae, Nepidae, and Teloganodidae)

295
were not included in the analysis, as these families occurred at sites belonging to all three altitude

303
The relationships between all the diversity metrics, geographical and physicochemical variables

322
The first three principle components explained approximately 55% of the variation (see Table   323 S2, in Supplementary Materials). The physicochemical variables nickel, copper, total nitrogen, 324 temperature, orthophosphate, COD, and pH were positively correlated (r >0.5) while the 325 geographical variables latitude, altitude, and longitude were negatively correlated on the first 326 principal component dimension (PC 1) (see Table S2, in Supplementary Materials).

327
The first PCA axis (29.1%) ( Figure 5) was mainly associated with the variables nickel, copper,

Principle Component Analysis (PCA) Ordination
The first three principle components explained approximately 55% of the variation (see Table S2, in Supplementary Materials). The physicochemical variables nickel, copper, total nitrogen, temperature, orthophosphate, COD, and pH were positively correlated (r > 0.5) while the geographical variables latitude, altitude, and longitude were negatively correlated on the first principal component dimension (PC 1) (see Table S2, in Supplementary Materials).
The first PCA axis (29.1%) ( Figure 5) was mainly associated with the variables nickel, copper, nitrite, orthophosphate, and temperature. The second PCA axis (16.1%) ( Figure 5) was mainly associated with two geographical variables (longitude and latitude) and three physicochemical variables (total phosphorus, BOD, and conductivity). Overall, the PCA biplot revealed temperature, nickel and geographical (latitude and longitude) gradients along the sites. The main variables that contributed most to the total variance in the first six principal components were latitude, longitude, total nitrogen, nickel, and conductivity (see Figure S1, in Supplementary Materials).
Water 2020, 12, x FOR PEER REVIEW 10 of 21 contributed most to the total variance in the first six principal components were latitude, longitude, 333 total nitrogen, nickel, and conductivity (see Figure S1, in Supplementary Materials).

Selection of Variables (Predictors) for Macroinvertebrate Community Diversity
The LASSO regression analyses applied on the 18 variables, selected five site-specific environmental variables namely temperature, nickel, conductivity, longitude, and total phosphorus as the strongest predictors for macroinvertebrate community taxa richness. For the total taxa richness in all sites, the two significant (p < 0.05) non-zero coefficient predictors were temperature (−0.269) and nickel (−1.733), while the significant coefficient predictors for taxa richness in the three altitude groups are displayed in Table 2. The taxa abundance had three site-specific environmental predictors in the DS and MS groups. The two predictors for the DS sites were nitrate nitrogen (45) and latitude (2255.65); while temperature (−7.851) was the single predictor of taxa abundance in the MS group.

Influence of Geographical Variables on Macroinvertebrate Diversity
In this study, we assessed a selection of physicochemical and geographical variables influencing Rwenzori macroinvertebrate communities at sites that had varying anthropogenic stress levels. Our study findings revealed temperature, latitude, and altitude to be key drivers of macroinvertebrate assemblages and diversity patterns. Our results were consistent with previous studies in the Ecuadorian Andes [42], Brazilian highlands [87], Tibetan highlands [88], and European mountains [89] which also found altitude, latitude, and temperature to influence macroinvertebrate communities and diversity patterns.
We found spatial structuring to be a key driver of macroinvertebrate community assemblages, as highlighted by the ordination analysis. Sites were structured along geographical gradients of altitude and latitude. Macroinvertebrate communities exhibited varying levels of spatial variation among the sites, possibly demonstrating the macroinvertebrate community responses to different proportions of environmental variation in the study area [90,91]. Our results indicated that latitude contributed most to the variability of the sites, and this was further confirmed by the least absolute shrinkage and selection operator (LASSO)regression that found latitude as a positive predictor for taxa abundance in DS sites. Our findings were supported by a prior study in Artic-alpine boreal headwater streams, that showed latitude as the best predictor of variation in macroinvertebrate assemblages [92,93]. The Artic-alpine boreal study also revealed that the distribution of aquatic ectothermic macroinvertebrates is greatly controlled by stream water temperature and that the latitudinal patterns observed in macroinvertebrate assemblages are additionally underscored by geographical gradients of altitude and ecoregion in stream environmental conditions. The LASSO regression revealed that longitude influenced taxa richness in US sites [48].
The environmentally heterogeneous sites were structured and separated along a latitudinal gradient in the Rwenzori region. The pattern of decreasing macroinvertebrate taxa richness along a decreasing latitudinal gradient (0.66 to 0.01 decimal units) was linked to the variation in the Rwenzori climatic conditions. All the DS are characterized by a tropical savannah climate and are located at lower latitudes close to the equator in the leeward section of the Rwenzori mountains and had the lowest taxa richness score. Latitude and topography also affected the temperature in the Rwenzori sites, as regions that were close to the equator and had low altitudes received the most solar radiation. The seasonal movements of the Inter Tropical Convergence Zone influences rainfall by shifting the wind patterns [64,94,95]. As we sampled sites located in both the leeward and windward sides of the Rwenzori mountain ranges, there was spatial variation in macroinvertebrate taxa composition and environmental characteristics among the different site groups.
Zonal studies in the Andes [96] linked to both altitude and latitude, showed that temperature mainly influences the taxa richness and structure of macroinvertebrate communities and further revealed that taxa richness was linearly related to water temperature, with macroinvertebrate diversity increasing with rising temperature as was observed with taxa abundance that increased from 257 individuals in the low temperature US group to 793 individuals in MS group. Increased water temperatures have been associated with changes in macroinvertebrate community structure, generalist shifts along the river continuum, and reduced resilience of community states [97]. In our study, altitude was negatively correlated with temperature and this was in agreement with previous studies in both glacier and non-glacier fed rivers, which indicated that altitude has a negative relationship with stream water temperature and the dissolved oxygen content [98][99][100]. The atmospheric pressure and oxygen partial pressure decrease with increasing altitude eventually leading to a decrease in the dissolved oxygen content [101]. Numerous studies have considered the influence of dissolved oxygen availability on the behavior, growth, energetics, and overall survival of aquatic macroinvertebrates. They revealed a dramatic variation in hypoxia tolerance among taxa and thus demonstrating a close association between dissolved oxygen availability and the presence or absence of particular macroinvertebrate taxa within a community [102]. For example, Simuliidae tend to be present in fast flowing well oxygenated waters [103], and our results were consistent with this general habitat association of the family as Simuliidae were the most abundant taxon. Altitude influenced the macroinvertebrate assemblages as it was significantly positively correlated with the diversity indexes of taxa richness, Pielou's evenness, and Simpson index. Taxa richness significantly differed between the three site groups, while the taxa abundance significantly differed between the DS and MS site groups. Out of the three altitude groups, US sites had the highest taxa richness score and this contradicted earlier studies in high altitude streams in northern Ecuador [101,104,105] and the Himalayas [106]. Their findings indicated that the decrease in taxa richness with increasing altitude was due to a decrease in the dissolved oxygen saturation and they further concluded that low oxygen availability could limit aquatic biodiversity at high altitudes. The differing findings from our study, could be attributed to the different ecoregions in Himalayan and Ecuadorian montane regions. Furthermore, we studied sites at altitudes less than 2000 m.a.s.l. where the climatic and environmental conditions are less harsh.

Indicator Families and Macroinvertebrate Diversity
Indicator value analysis ascertained the altitude site group preferences of the macroinvertebrate families. The analysis revealed that Corixidae was the only family that was significantly associated with DS sites which were mainly characterized by low velocity streams. Our findings are supported by previous studies, that indicated low velocity streams are the preferred habitat for Corixidae in Africa and North America [107,108]. Dytiscidae, another taxon found in our study, was significantly associated with US sites which are characterized by high velocity streams and rivers. An earlier study in the Iberian peninsula supported our finding and revealed that certain species of Dytiscidae are confined to mountainous regions and high velocity streams [109]. Dytiscidae are the most diverse water beetle family and are highly adapted for aquatic life [110], by having enlarged prothoracic and elytra edges, an ideal shape for good swimming in high water velocity environments [111][112][113].

Physicochemical Variables Influence on Macroinvertebrate Diversity
Regression analysis indicated that the two site-specific local environmental variables, temperature and nickel, were key in structuring the overall taxa richness. Nickel was identified as a key variable shaping the macroinvertebrate communities along a mineral gradient, as some of the MS and DS sites were located along the copper belt of western Uganda [114]. The Kilembe mines area has several old copper piles that leach into the water system [23,114,115]. These mines have deposits of copper-nickel ore and the mine drainage water contains copper and nickel ions that leach into the water system and this negatively influences the physiological function of macroinvertebrate taxa and are toxic at high concentrations. The concentrations of copper and nickel were highest at sites from DS and MS groups and this was due to the stream water velocity as well as the weathering action of the copper mine tailings causing increased concentrations of copper and nickel at these sites in comparison to the US sites which recorded extremely low nickel concentrations. The negative impact of heavy metals on macroinvertebrate assemblages was supported by studies performed in high altitude streams revealing the sensitivity of macroinvertebrates to metal concentrations [116][117][118]. Acid mine drainage reduces taxa diversity in aquatic ecosystems through both acute and chronic mechanisms, including the toxicity of metals in stream water and sediment [119][120][121][122]. In the Danube delta, an overload of copper, nickel, and zinc has been shown to have inhibitory effects on the development of molluscs and gastropods [123]. The metallic compounds were found to disrupt the oxygen level, byssus (attachment fibers) formation, reproductive processes, and the overall mollusc development. The lethal effects of these metals in aquatic organisms is related to the inhibition of enzymes involved in cellular respiration [123].
Temperature is of prime importance to the living conditions of aquatic fauna with the thermal tolerance limits playing a vital role in reducing macroinvertebrate community diversity, along with various ecological variables of topography and climate, which may thus explain the temperature dependent patterns of taxa richness in our study [124][125][126]. Water temperature was a strong predictor of macroinvertebrate taxa richness. This was in agreement with earlier studies in high-altitude Ugandan rainforest streams [59] revealing that forested habitats that had low stream temperature, were dominated by pollution sensitive macroinvertebrate families such as Ephemeroptera, Plecoptera, Trichoptera, and Odonata (EPTO). These EPTO families require habitats with higher dissolved oxygen, low levels of turbidity, a neutral pH, and cool water temperature, therefore making them more susceptible to physicochemical changes caused by water pollution. The sites that were characterized by both high stream temperature and human impact were dominated by tolerant taxa such as Hemiptera, Coleoptera, and Diptera.
Previous studies in glacier-fed rivers and streams have revealed that glacial meltwater influences environmental variables of conductivity and temperature [32,127,128]. In our study, the positive correlation between conductivity and temperature can also be attributed to higher glacial contributions to streamflow resulting in lower conductivity levels. The influence of conductivity on macroinvertebrate diversity is in agreement with our previous study in the Rwenzori region [11] as the LASSO regression analysis that revealed conductivity to be a key predictor for taxa richness in the US and DS site groups. Water temperature is lowered with heightened glacial contributions, and varying amounts of glacial melt inputs received by different sites thus cause conductivity and temperature to differ among the sites. The study revealed positive correlations between pH, copper, and temperature, which is supported by other previous studies of copper in water which showed that temperature and pH were the key variables influencing copper concentrations in surface water [129]. Furthermore, this pattern is seemingly driven by the area catchment geology as revealed by the positive correlations between pH, copper, and temperature as previous studies revealed that the catchment geology, anthropogenic activities, and geographic location influence the variability of copper in surface waters [130,131]. In our study, the DS group sites that had high concentrations of copper were also characterized by high anthropogenic disturbance levels and were located downstream of the Kilembe copper mines, however, we did not include catchment geology in our study. Future research should include catchment geology as our study area is located in a mining region with heightened weathering of the Kilembe mine tailings and rocks which contribute to metal concentrations in the surface water [17,114].
Family abundance was highest in MS sites, as the majority of the sites were mainly characterized by moderate and low levels of anthropogenic disturbance as explained by the intermediate disturbance hypothesis. However, the marked decrease in taxa abundance at a few DS sites contrasted with previous studies on alpine aquatic systems, which revealed that taxa abundance increased with distance from the source [132]. A reduction in the abundance could be attributed to the toxicity exerted by total phosphorus [133][134][135], nitrate nitrogen, and nitrite nitrogen on macroinvertebrates [136]. Our observed study pattern is also related to the medium to high levels of anthropogenic disturbance at the majority of the DS sites. This has been demonstrated by earlier studies in Gabonese streams that revealed an inverse relationship between anthropogenic disturbance and arthropod abundance [137]. The low family abundance in US sites was in agreement with earlier tropical alpine studies [138,139] that showed that macroinvertebrate abundance decreased with altitude. However, the low taxa abundance can also be attributed to the high and moderate levels of anthropogenic disturbance from mining, urbanization, and water abstraction at these sites, which were similar to the activities in our study. The majority of US sites are currently facing the anthropogenic pressures of human activities, agriculture, and water abstraction [11], and previous studies have demonstrated linkages to anthropogenic pressures in other tropical high altitude streams contributing to a decreasing macroinvertebrate abundance [99]. In the tropics, human settlements, agriculture, and other anthropogenic activities can extend to relatively high altitudes because of the benign climate, as is the case in our study region whereby the main human activities take place between 900 to 2000 m.a.s.l. In turn, these tropical high-altitude streams can be affected by organic pollution, agrochemicals, sedimentation, and mining wastes. This will eventually have a negative effect on the aquatic ecosystem diversity as sensitive taxa will drastically decrease due to the aforementioned threats.
Macroinvertebrate community diversity is presumed to increase as glaciers retreat, due to upstream taxa migration as a result of increased water temperature and higher channel stability caused by reducing glacial influence [32,34,127,140]. More research is needed to better understand the macroinvertebrate community responses to declining glacial meltwater contributions in the Rwenzori aquatic ecosystems across a wider spatial and temporal scale. In future studies, incorporating the glaciality index (computed as a function of glacier area and distance from the glacier), in addition to the mineral and temperature gradients, is necessary to clarify the relative impact of the factors on macroinvertebrate biodiversity.

Future Research and Implementation of Insights in Management
Based on our sample size using both categorical and continuous analyses, it was possible to link the macroinvertebrate community diversity to the environmental drivers for the three altitude groups thereby adding to the existing knowledge of the Rwenzori region. Our findings indicate that family-level monitoring of key indicator taxa along with measurements of temperature and copper-nickel concentrations can be effectively used to assess river health in these unique Rwenzori aquatic ecosystems. The categorical analysis provided clear indications as to which variables were important, reduced complexity of the analysis and interpretation of results. Based on our study, biomonitoring and conservation programs would benefit from spatial stratification according to different climatic zones and altitude ranges where land use activities are conducted. This would guarantee that (a) the monitoring and conservation programmes of the Rwenzori aquatic systems are conducted in an ecologically systematic and meaningful way, and that (b) predictive models for river bioassessment are more robust and accurate than those generated from larger and more heterogeneous areas. This could enable cost effective regular monitoring by water resources stakeholders to mitigate water pollution and sustainably manage the aquatic ecosystems in the region. As macroinvertebrate assemblages are strongly driven by both temperature and mineral gradients, these unique equatorial alpine ecosystems can be expected to shift in taxa composition and distributions, as they respond profoundly to climate change induced glacial retreat and other anthropogenic activities [141][142][143].

Conclusions
In our study, we investigated the environmental drivers of macroinvertebrate diversity in Rwenzori region. Both geographical and physicochemical variables influenced the overall macroinvertebrate diversity in the Rwenzori region. In spite of the contribution of anthropogenic and natural stressors, the outcomes of this field study highlighted the significance of these stressors in determining macroinvertebrate diversity in the Rwenzori region. Family richness responded to physicochemical variables tested in the study, revealing temperature and nickel as important predictors of macroinvertebrate community diversity. Our findings at family level indicate characteristic patterns of structural variability in faunal univariate diversity metrics between the US, MS, and DS altitude groups, and thus definitive results have been determined in the context of drivers of macroinvertebrate diversity. We determined that in studied sites with similar climate and altitude, the macroinvertebrate community responds differently depending on site-specific environmental characteristics. The study findings indicate that in future biomonitoring programmes, indicator taxa, site-specific environmental variables are important and should be taken into account given the intricacy of generalizing categorical responses, along the altitude site groups. Thus, site-specific environmental variables are integral when carrying out conservation and restoration programmes in these alpine rivers. On a practical level, there are important protocols that can be deduced from this study, which can be immediately applied to assess river health and key drivers of macroinvertebrate community in sub-Saharan equatorial alpine rivers.