Structure of Benthic Macroinvertebrate Communities in the Rivers of Western Himalaya, Nepal

According to River Continuum Concept (RCC), channel morphology, including sediment loads and channel width, river habitat, flow regimes and water quality, differs from the tributary to the downstream river’s mainstem, allowing shifts in faunal composition from dominance of shredders to collectors downstream, respectively. Tributaries are responsible for contributing organic carbons, nutrients and water. However, such knowledge is still limited in the monsoon-dominated river systems of the Himalaya. The study was conducted in the river’s mainstem and tributaries of the Karnali River Basin, which are glacier and spring-fed river systems, respectively, in the western Himalaya, Nepal. A total of 38 river stretches in the river’s mainstem and tributaries were sampled during post-monsoon and pre-monsoon seasons in the years 2018 and 2019. Water quality parameters, such as pH, temperature, electrical conductivity, total dissolved solids, dissolved oxygen, alkalinity and hardness, and the benthic macroinvertebrates were studied. Ten subsamples of benthic macroinvertebrates were collected following the multi-habitat sampling approach at each site. High taxa richness was recorded in tributaries compared to the river’s mainstem while abundance was similar between river types. Non-metric multidimensional scaling (NMDS) formed two distinct groups, reflecting high similarities in benthic macroinvertebrate composition within the tributaries and river’s mainstem rather than between river types. Redundancy analysis (RDA) indicated water temperature and pH as major environmental predictors for benthic macroinvertebrate variability between river types. Therefore, river type-based conservation efforts that account for upstream–downstream linkages of aquatic biota and resources in freshwater ecosystems can ensure the ecological integrity of the whole river basin.


Introduction
Tributaries in a landscape determine the physical and ecological condition of a river's mainstem. Their number and characteristics shape the river's mainstem geometry, substrate types, hydraulics, nutrients, organic matters and water quality [1,2]. Based on the River Continuum Concept (RCC), tributaries receive allochthonous inputs, and thus the structure of the aquatic communities-the functional feeding groups-is different from that of the downstream river's mainstem, in which the energy cycle depends on autochthonous inputs [3].
Tributaries play a crucial role in maintaining macroinvertebrate populations in the downstream river's mainstem by flow-induced drifting along the longitudinal gradient [4,5]. Depending on the river's flow regime, macroinvertebrates can drift from 100 m to potentially over 1 km downstream [6,7]. Katano et al. [8] analyzed the contribution of a tributary in terms of resources and aquatic fauna to the Figure 1. Study area, located in the western Himalaya, Nepal, south Asia (top left image). The area's physiographical zones (bottom left image); the two major national parks (Khaptad NP and Bardia NP) and the sampling sites across the Karnali River's tributaries and mainstem are also depicted (right image). NP stands for National Park.

Physico-Chemical and Hydro-Morphological Variables
At each site, the physico-chemical variables (pH, water temperature, electrical conductivity (EC), total dissolved solids (TDS) and dissolved oxygen saturation (DO, %)) were measured using a HANNA (Model:HI9829) multi-parameter probe. Alkalinity and hardness were measured on site. Flow velocity was recorded only for wadeable rivers (i.e., tributaries) with a cross-sectional approach at 0.6 x water depth (measuring from the surface) by using a Global Flow Probe (Xylem brand) at 1 m intervals across the wetted length of the cross section. Water discharge was calculated by multiplying each flow velocity value with the relevant cross-sectional area and summing up all the values across the wetted perimeter. The proportion of surface flow types (%riffle, %run, %rapid and %pool) and river width was visually estimated at each site prior to macroinvertebrate sampling.

Benthic Macroinvertebrate Sampling and Processing
At each site, 10 sub-samples of benthic macroinvertebrates were collected from a 50-100 m river stretch [24]. Sub-samples were taken from available micro-habitats with respect to their relative habitat coverage of 10%. Micro-habitat less than 10% coverage were not considered for sampling of benthic macroinvertebrates. Benthic macroinvertebrates were collected using a hand net with a 500 µm mesh size and a frame of 25 cm × 25cm, where the stream bottom was disrupted for a minute, which dislodged the benthic organisms, and transported by the water current into the net. The collected sub-samples of each site were made on one composite sample, preserved in 95% ethanol and brought to the laboratory. In total, benthic samples covered an area of 0.625 m 2 in each sampling site.
In the laboratory, each macroinvertebrate sample was rinsed with tap water and the entire sample was sorted. Macroinvertebrates were identified and counted at the genus (mainly Thirty-eight sites were sampled along the river's mainstem and the tributaries of the Karnali River during the pre-monsoon (April/May) and post-monsoon seasons (November) of 2018 and 2019. The sites were distributed in three watersheds, namely Tila, Middle-Karnali and Lower-Karnali across four physiographic zones: (i) Terai, (ii) Siwalik, (iii) Mid-Hill and (iv) High Mountain ( Figure 1). Sampling sites in the tributaries were located in the downstream sections where the channel width gets wider and lack riparian vegetation.

Physico-Chemical and Hydro-Morphological Variables
At each site, the physico-chemical variables (pH, water temperature, electrical conductivity (EC), total dissolved solids (TDS) and dissolved oxygen saturation (DO, %)) were measured using a HANNA (Model:HI9829) multi-parameter probe. Alkalinity and hardness were measured on site. Flow velocity was recorded only for wadeable rivers (i.e., tributaries) with a cross-sectional approach at 0.6 x water depth (measuring from the surface) by using a Global Flow Probe (Xylem brand) at 1 m intervals across the wetted length of the cross section. Water discharge was calculated by multiplying each flow velocity value with the relevant cross-sectional area and summing up all the values across the wetted perimeter. The proportion of surface flow types (%riffle, %run, %rapid and %pool) and river width was visually estimated at each site prior to macroinvertebrate sampling.

Benthic Macroinvertebrate Sampling and Processing
At each site, 10 sub-samples of benthic macroinvertebrates were collected from a 50-100 m river stretch [24]. Sub-samples were taken from available micro-habitats with respect to their relative habitat coverage of 10%. Micro-habitat less than 10% coverage were not considered for sampling of benthic macroinvertebrates. Benthic macroinvertebrates were collected using a hand net with a 500 µm mesh size and a frame of 25 cm × 25 cm, where the stream bottom was disrupted for a minute, which dislodged the benthic organisms, and transported by the water current into the net. The collected sub-samples of each site were made on one composite sample, preserved in 95% ethanol and brought to the laboratory. In total, benthic samples covered an area of 0.625 m 2 in each sampling site.
In the laboratory, each macroinvertebrate sample was rinsed with tap water and the entire sample was sorted. Macroinvertebrates were identified and counted at the genus (mainly Ephemeroptera, Plecoptera, Trichoptera, Mollusca and Oligochaeta, except for Tubificidae), family (Coleoptera, Heteroptera, Odonata, Diptera, Lepidoptera and Megaloptera) and subfamily (Chironomidae and Psephenidae) levels, using available region-specific keys [25][26][27][28]. The identified samples have been preserved in vials containing 95% ethanol and stored at the Aquatic Ecology Centre of Kathmandu University.

Statistical Analyses
A Mann-Whitney U test was carried out to differentiate variation in taxa richness and abundance of macroinvertebrates between the Karnali River's mainstem and its tributaries for each watershed, i.e., Tila, Middle-Karnali and Lower-Karnali.
Non-metric multidimensional scaling (NMDS) was used to group the samples based on abundance of benthic macroinvertebrates across river types (river's mainstem and tributaries) and physiographical zones. NMDS ordination is a robust ordination technique for exploring similarities or dissimilarities in biological data as it does not require any assumptions of multivariate normality and yields good results even when large numbers of data sets have zero values [29]. Benthic macroinvertebrates abundance was log (x + 1)-transformed prior to the NMDS analysis. Rare taxa occurring only in 3 sites were removed from the analysis. Permutational Multivariate Analysis of Variance (PERMANOVA) was carried out with the Adonis function in R/vegan to test whether benthic macroinvertebrate assemblages differed significantly across the physiographical zones and river types. The Bray-Curtis distance was used as a distance measure in the benthic macroinvertebrates abundance data.
A detrended correspondence analysis (DCA) was employed on benthic macroinvertebrates abundance data to examine whether the data followed a linear or unimodal distribution [30]. As the length of the first DCA axis (2.2 SD, longest gradient) was larger than 2.0 but smaller than 3.0 SD, the linear ordination method, redundancy analysis (RDA), was chosen to find the explanatory variables governing the distribution of the benthic macroinvertebrates. Prior to RDA, multicollinearity test was conducted and "TDS" was removed due to high collinearity (r > 0.70) with electrical conductivity in the final RDA plot. All the analyses were carried out in R software [31].

Results
A total of 128 taxa of macroinvertebrates belonging to 84 families and 22 orders were recorded in the study sites of the Karnali River Basin. Taxa richness of benthic macroinvertebrates was higher in the tributaries compared to the river's mainstem for Tila (p < 0.05) and Middle-Karnali (p < 0.01), but not for the Lower-Karnali (p = 0.79) watershed ( Figure 2). Abundance was significantly higher in the tributaries (compared to the river's mainstem) of the Middle-Karnali watershed (p < 0.01), but it was not different for the Tila (p = 0.21) and Lower-Karnali (p = 0.12) (Figure 3). High variability across sites was documented in the river's mainstem, except for the Tila watershed.  Abundance was significantly higher in the tributaries (compared to the river's mainstem) of the Middle-Karnali watershed (p < 0.01), but it was not different for the Tila (p = 0.21) and Lower-Karnali (p = 0.12) (Figure 3). High variability across sites was documented in the river's mainstem, except for the Tila watershed. Abundance was significantly higher in the tributaries (compared to the river's mainstem) of the Middle-Karnali watershed (p < 0.01), but it was not different for the Tila (p = 0.21) and Lower-Karnali (p = 0.12) (Figure 3). High variability across sites was documented in the river's mainstem, except for the Tila watershed.  Taxa richness was higher for all macroinvertebrate groups in tributaries compared to the river's mainstem ( Figure 4). The faunal richness was mainly dominated by Ephemeroptera, Trichoptera and Diptera in the river's mainstem and tributaries of the Tila and Middle-Karnali watersheds (Figure 4a,b), while Mollusca was one of the three dominant groups in the river's mainstem and the tributaries of the Lower-Karnali watershed (Figure 4c).
Geosciences 2020, 10, x FOR PEER REVIEW 6 of 14 Taxa richness was higher for all macroinvertebrate groups in tributaries compared to the river's mainstem ( Figure 4). The faunal richness was mainly dominated by Ephemeroptera, Trichoptera and Diptera in the river's mainstem and tributaries of the Tila and Middle-Karnali watersheds (Figure 4a,b), while Mollusca was one of the three dominant groups in the river's mainstem and the tributaries of the Lower-Karnali watershed (Figure 4c). Composition of benthic macroinvertebrates was dominated primarily by Ephemeroptera and Diptera for the river's mainstem and tributaries of the Tila and Middle-Karnali watersheds, while in the Lower-Karnali, Heteroptera, Diptera and Mollusca were highly abundant ( Figure 5). Taxa richness and abundance were not different between the pre-monsoon and post-monsoon seasons, within and among watersheds.
Geosciences 2020, 10, x FOR PEER REVIEW 8 of 14 tributaries of each watershed. PERMANOVA showed significant differences in macroinvertebrate composition across sites between the tributaries and the river's mainstem for Mountain (F = 6.09, df = 1, R 2 = 0.11, p = 0.001) and Lowland (F = 3.58, df = 1, R 2 = 0.16, p = 0.003). Both the richness (p > 0.05) and abundance (p > 0.05) of the functional feeding groups did not vary across sites in the river's mainstem and tributaries, for each studied watershed: Tila, Middle-Karnali and Lower-Karnali; also in general. For the Mountain zone, the richness of scrapers, predators and collector-filterers significantly varied across sites between the river's mainstem and tributaries but not for collector-gatherers (p = 0.08) and shredders (p = 0.18) (Figure 8a). Similarly, abundance, except for shredders (p = 0.18) and scrapers, were significantly different across sites between the river's mainstem and the tributaries for collector-gatherers, predators and collector-filterers (Figure 8b). Both the richness (p > 0.05) and abundance (p > 0.05) of the functional feeding groups did not vary across sites in the river's mainstem and tributaries, for each studied watershed: Tila, Middle-Karnali and Lower-Karnali; also in general. For the Mountain zone, the richness of scrapers, predators and collector-filterers significantly varied across sites between the river's mainstem and tributaries but not for collector-gatherers (p = 0.08) and shredders (p = 0.18) (Figure 8a). Similarly, abundance, except for shredders (p = 0.18) and scrapers, were significantly different across sites between the river's mainstem and the tributaries for collector-gatherers, predators and collector-filterers (Figure 8b). In the RDA, the first two axes, associated with water temperature, pH, dissolved oxygen, conductivity, alkalinity and hardness, explained only 23.24% of the taxonomic variance; the RDA's 1st axis was found to be highly significant (RDA1, df = 1 F = 7.1531, p = 0.001). The sum of canonical eigenvalues was 15.7. A strong negative correlation existed between the first axis and water temperature (r = −0.49), while oxygen saturation (DO%) was −0.82. Alkalinity (0.37) and hardness (−0.57) were positively correlated with axis 1 (Figure 9). In the RDA, the first two axes, associated with water temperature, pH, dissolved oxygen, conductivity, alkalinity and hardness, explained only 23.24% of the taxonomic variance; the RDA's 1st axis was found to be highly significant (RDA1, df = 1 F = 7.1531, p = 0.001). The sum of canonical eigenvalues was 15.7. A strong negative correlation existed between the first axis and water temperature (r = −0.49), while oxygen saturation (DO%) was −0.82. Alkalinity (0.37) and hardness (−0.57) were positively correlated with axis 1 (Figure 9).
In the RDA, the first two axes, associated with water temperature, pH, dissolved oxygen, conductivity, alkalinity and hardness, explained only 23.24% of the taxonomic variance; the RDA's 1st axis was found to be highly significant (RDA1, df = 1 F = 7.1531, p = 0.001). The sum of canonical eigenvalues was 15.7. A strong negative correlation existed between the first axis and water

Discussion
The tributaries in our study have high habitat heterogeneity and diverse flow types, which resulted in increased taxonomic richness and abundance in the tributaries compared to the river's mainstem. Due to the high water discharge and volume in the river's mainstem, mineral habitats are tightly bounded making colonization of macroinvertebrates unfavorable. It is not surprising that we did not observe seasonal variation in the macroinvertebrate community structure across sites of the tributaries and the river's mainstem in all watersheds, in contrast to the results reported for highly fragmented and polluted rivers of Nepal (see [19,22]). This might be due to the fact that the study sites were in near pristine conditions that did not lead to large variation in community structure between seasons. It is also true that the results were based on limited number of sites along the river's mainstem of only one river, i.e., Karnali River. The study conducted in tributaries of Mahakali and Karnali River Basin (in particular West Seti and Bheri Rivers) by Tachamo Shah [22] indicated that richness and abundance of macroinvertebrates varied between seasons in the rivers of the western Himalaya under flow modification due to water abstraction projects [22].
Mac Nally [32] found lower densities of macroinvertebrates in tributaries compared to the river's mainstem during low flows in the Acheron River, Australia, but densities were similar during high flow; in turn, taxarichness was similar during low flows but increased in the river's mainstem during high flows. Studies showed that discrepancies in benthic macroinvertebrate community composition between tributaries and the river's mainstem are huge with increased distance between the tributaries and mainstem [5].Our data showed variation in richness and abundance between the river's mainstem and the tributaries only for the Middle-Karnali watershed. This could be that the Tila and Lower-Karnali, situated in the High-Mountain and Lowland zones, respectively, had a similar water quality status (in particular temperature, p > 0.05) between the tributaries and the river's mainstem. The water quality parameters alone explained over 20% variation in benthic macroinvertebrate composition for the watersheds (Figure 9). In contrast, for the Middle-Karnali situated in Middle mountain, the water quality parameter value (in particular temperature, p < 0.01) was found to be different between the tributaries and the river's mainstem. The majority of variations in macroinvertebrate composition are governed by other factors that were not included in this study, such as altitude, slope, aspects, channel morphometry, flow regimes, water depth and river channel width; these could all be of high importance in structuring the benthic macroinvertebrate community. Research findings in other parts of the world revealed that slope, angle of intersection, size of watershed and the amount of rainfall shaped macroinvertebrate communities in river mainstems [33][34][35]. In the low-flow season, low-order headwater streams may not support high taxonomic richness and abundances due to river contraction [36], which decreases the amount of available habitat [37] and food resources (distant from riparian vegetation) for most aquatic biota compared to conditions found at higher flows. In Canada, macroinvertebrate abundance was found high in mainstems but no difference in taxonomic richness was observed during the high-flow season [38].These findings indicate that macroinvertebrate diversity and density vary temporally depending on the availability of habitats and food resources.
The River Continuum Concept (RCC) foresees that macroinvertebrate community composition changes gradually from the headwaters to the downstream river's mainstem [3]. For instance, the relative proportions of functional feeding groups change from shredder dominance in headwaters to collector dominance in the lower reaches of large rivers. These findings are broadly applied to our study river system. The abundance of collector-gatherers was higher in tributaries compared to the river's mainstem ( Figure 9); this was also observed for scrapers in contrast to Vannote [3]. This might be due to the fact that our sites in the tributaries were located in the downstream section and devoid of riparian vegetation along the banks, supporting abundant food resources and reduced competition among taxa that might have influenced the results. Moreover, in the river's mainstem, macroinvertebrate abundances were more evenly distributed among different functional feeding groups, which aligned with the results of Heino et al. [39]. In general, the abundance of the macroinvertebrates increased with increasing stream size, and some species indicated restricted distribution only in the tributaries, which is in contrast to the findings of Heino et al. [39]. Taxa belonging to grazing snails were absent from the tributaries and mainly occurred in river reaches of the Lowland.
Stream size seemed to be a major factor influencing richness and abundance of functional feeding groups of macroinvertebrates in our study (Figure 8), concurring with earlier findings from the central Himalayas and other parts of the world [20,40].The changes in community composition along the river channel gradient suggests that conservation of aquatic biodiversity might be achieved by maintaining tributary habitats as drifting of benthic macroinvertebrates occur from the tributaries to the river's mainstem. Even if a river's mainstem is degraded, tributaries can play a vital role for restoration of faunal assemblages in the river's mainstem. Moreover, tributaries act as refugia for the river's mainstem in case of temporary disturbances in the river's mainstem (e.g., [8,41,42]). Similarly, tributaries provide excellent spawning habitats for migratory species [43]. Therefore, protecting headwaters from habitat degradation may preserve several aquatic species of different groups, including endemic species and threatened species of the drainage system. Furthermore, headwaters are the habitats of regionally rare species and rare assemblage types that are not typically found in river mainstems [44][45][46][47]. It has been suggested that headwaters should be regarded as key zones to focus conservation efforts upon in freshwater ecosystems, as they are important contributors to the ecological integrity of whole drainage systems [48]. Since river size is a key environmental driver determining the taxonomic and functional composition of macroinvertebrates, stream order-based river typology should be considered in river conservation plans [20].
Our results provide a basis for future studies on the benthic macroinvertebrate composition and various environmental variables across river mainstems and tributaries. It is evident from this study that tributaries can restructure both the habitats and macroinvertebrate community composition of the downstream river mainstems; however, the rate of change may depend on the relative size of the tributaries and mainstems. The scientific water resources development in tributaries and river mainstems could enhance upstream-downstream linkages of biota and other resources, including substrates that directly or indirectly influence ecosystem services to the people living in the vicinity.
Author Contributions: R.D.T.S., D.N.S., S.S. and D.R. conceived the study; R.D.T.S. and D.N.S. designed the study and conducted field works; R.D.T.S. analyzed the data and wrote the paper; S.S., D.N.S. and D.R. reviewed and edited the manuscript. All authors have read and agreed to the published version of the manuscript.
Funding: This research was funded by USAID-Paani program -Grant no. G-KAT-013.