Diversity and distribution of benthic invertebrates dwelling rivers of the Kruger National Park, South Africa

water flowed downstream through the park, suggesting a beneficial effect of protected river reaches on benthic invertebrate diversity. However, for the Crocodile River, which makes up the southern border of the park, this trend was less conspicuous, suggesting that this river may experience the greatest threats. More generally, benthic invertebrate communities were driven by the concentrations of phosphates, sulphates, ammonium and organic matter and by substrate characteristics. Conservation implications: Meiobenthic organisms are very abundant in KNP rivers and react to environmental gradients; thus, they should be more considered for bio-monitoring or conservation of comprehensive assemblages of animals. Interestingly, protected reaches tended to show a reduced dominance of the invasive T. granifera and a higher diversity of benthic invertebrates.


Introduction
Rivers represent only a minor part of Earth's water resources, but they play a key role for the connectivity of ecosystems and the maintenance of biological diversity at a global scale (Dudgeon et al. 2006). However, most river ecosystems have been degraded by human activities throughout history (e.g. fisheries, commercial routes, irrigation, production of energy, receptacles for wastewater and introduction of invasive species), resulting in a current situation where rivers are one of the most threatened ecosystems worldwide (Malmqvist & Rundle 2002). As an example, South African rivers are intensively exploited to ensure water security for urban areas, agriculture and industry, where water use often exceeds natural water availability (O'Keeffe 1989). Those rivers further serve as outlets for urban wastewater, excess nutrients and xenobiotics from agriculture, as well as metals draining from mine exploitation or industrial activities (e.g. Jackson et al. 2013;Rimayi et al. 2018;Riddell et al. 2019). Nel et al. (2007) analysed the health status of South African rivers and found that 84% of river reaches (and especially those located in the largest, permanently flowing river systems) were degraded to the point they could be considered as highly threatened. Nature conservation regions including relatively large portions of river catchments can help mitigate pollution inputs and restore water quality (Nel et al. 2007). Although some xenobiotics like organochlorine pesticides Meiobenthos (or meiofauna) are microscopic invertebrates that inhabit biofilms and interstitial spaces in rivers. They are diverse and extremely abundant, and they perform essential ecological functions by linking microbial production to higher trophic levels (e.g. macrobenthic invertebrates and fishes). However, meiobenthic communities remain poorly studied in Africa. Here, we sampled meio-and macrobenthic invertebrate communities associated with biofilms and sediments across an upstream-downstream gradient along the Olifants, Sabie and Crocodile rivers flowing through the Kruger National Park (KNP). We expected to link differences in community structure to environmental gradients as those rivers show different degrees of anthropogenic stress as they enter the park. Both meio-and macrobenthic communities differed across rivers and also structured along an upstream-downstream gradient. The upstream sites, which were the closest to the park borders, consistently showed a lower diversity in all three rivers. There, the invasive snail Tarebia granifera strongly dominated (making up 73% -87% of the macrobenthos), crowding hard substrates, while concomitantly the abundances of biofilm-dwelling meiobenthos like nematodes and rotifers were substantially reduced. Nevertheless, the diversity and evenness of communities then tended to increase as water flowed downstream through the park, suggesting a beneficial effect of protected river reaches on benthic invertebrate diversity. However, for the Crocodile River, which makes up the southern border of the park, this trend was less conspicuous, suggesting that this river may experience the greatest threats. More generally, benthic invertebrate communities were driven by the concentrations of phosphates, sulphates, ammonium and organic matter and by substrate characteristics.
can still drift over long distances, their bio-accumulation potential in aquatic organisms of conservation areas is correlated with their proximity to pollution sources (Wolmarans et al. 2021). This means that river reaches further downstream from pollution sources should show a better ecological status.
To reduce threats and develop proactive management practices, one requires an in-depth understanding of the response of natural communities and ecosystem processes to water quality. However, in the case of South African rivers, management of water resources is often decided in a context of incomplete knowledge of the response of the different biological communities present in the field, leading to uncertainty in decision-making (see, e.g., Roux, Kleynhans & Thirion 1999a;Roux et al. 1999b). For example, management decisions concerning the conservation actions applicable in the large protected area of the Kruger National Park (KNP) have been traditionally based on the assessment of abiotic rather than biotic factors (Solomon et al. 1999). This is based on the assumption that resultant biotic patterns are likely to be correlated with abiotic components. However, to measure, protect and restore natural river integrity in the KNP, it would be more relevant to include multiple biological indicators of water quality (Rogers & Biggs 1999). Pollutants not only alter the chemistry of the water column but also persist in sediments (Gerber et al. 2015) and within tissues of organisms (Gerber et al. 2016;Seymore 2014;Wolmarans et al. 2021). The consequences of altered environmental conditions have been found to translate into conspicuous modifications of the structure of bacteria, diatom, fish and insect communities in KNP water (Farrell et al. 2019;Rasifudi et al. 2018;Riddell et al. 2019;Shikwambana et al. 2021). Another cause of stress to the native aquatic communities of KNP is the accidental introductions followed by problematic blooms of invasive pathogenic micro-organisms, fishes, crayfishes and molluscs species (e.g. Jones et al. 2017;Macdonald 1988;Petersen et al. 2017).
While the response of vertebrates to abiotic and biotic threats is conspicuous, the response of benthic invertebrates (living hidden in biofilm matrices and aquatic sediments) is mostly overlooked when it comes to assess ecosystem health status. Among those benthic invertebrates, one may further make a distinction between the macrobenthos, such as aquatic insects, snails and leeches that are large enough to be visible and usually are retained on 500 μm meshes, and the meiobenthos, such as nematodes, rotifers, copepods and tardigrades that are so tiny that they are invisible and usually pass through 500 μm-meshes but are retained on 20 μm meshes (Ptatscheck, Gehner & Traunspurger 2020a). Meiobenthic invertebrates are little studied although they show complex behaviours and extraordinary physiologies that allow them to colonise most, if not all, benthic habitats (Brüchner-Hüttemann, Ptatscheck & Traunspurger 2020;Rebecchi, Boschetti & Nelson 2020) where they can reach remarkable abundances (between 10 5 -10 6 ind. m -2 on any submerged substrate) (Majdi, Schmid-Araya & Traunspurger 2020;Traunspurger, Wilden & Majdi 2020). Meiobenthos play an important role in riverine food webs as intermediaries: grazing on microbes (e.g. bacteria, protozoans and micro-algae) and serving as prey for a variety of macro-invertebrates and fish juveniles (Majdi & Traunspurger 2015;Ptatscheck et al. 2020b;Schmid-Araya et al. 2002). Some studies have highlighted the potential of meiobenthic invertebrates to reflect aquatic ecosystem health in South Africa: for instance, Gyedu-Abadio et al. (1999) observed that the density and diversity of nematode assemblages were affected by the concentration of metals in a subtidal portion of the Swartkops estuary. Further experimentation in the laboratory showed that nematode genera like Axonolaimus, Theristus and Paramonhystera were tolerant to metal pollution (Gyedu-Ababio & Baird 2006). In Europe, a nematode 'species at risk' index has been developed to monitor river sediment pollution (Höss et al. 2011). This index adds value to routinely used macrobenthicbased indices because nematodes are ubiquitous and experience the effects of pollutants during their wholelife cycle within the sediment (Brüchner-Hüttemann et al. 2021). Whether such meiobenthos-based indices could be tested, developed and added to the bio-indication toolbox outside Europe is a question of general interest for ecologists and stakeholders, but lacks support from an insufficient collection of field data so far.
The specific objective of this study is to provide fresh insights into the diversity and distribution of benthic invertebrates (meio-and macrobenthos) in three large rivers (Olifants, Sabie and Crocodile) flowing through an emblematic protected area of austral Africa, the KNP. We sampled the main benthic habitats in those rivers (i.e. biofilms associated with hard substrates and surface sediment) to examine the structure of resident communities. The three rivers were expected to show a gradient of environmental stress as they enter the park: Olifants and Crocodile rivers experience different types of pressures (mostly related to mining practices for Olifants and to agricultural practices for Crocodile), while Sabie is quite pristine in comparison. Thus, we expected that the structure of benthic invertebrate communities could reflect gradients of pollution experienced by the three rivers. Finally, an upstream-downstream gradient was also expected to be manifest through a beneficial effect of preserved areas of the park on the taxonomic richness of benthic assemblages.

Study sites
The landscape in the southern part of the KNP consists of plains showing a gentle eastward slope and being drained by three major river systems (Crocodile, Sabie and Olifants rivers; Figure 1). The Crocodile (catchment area: 10 455 km 2 ) and Sabie rivers (6252 km 2 ) start off in the Drakensberg region. The Olifants River has the largest catchment area (54 434 km 2 ) and starts in the Highveld region (Muller & Villet 2004;Venter & Bristow 1986). It generally shows the highest annual run-off followed by the Crocodile and Sabie rivers (Muller & Villet 2004) but, at the time of sampling, all rivers had similar, relatively low discharges (3.82 m 3 s -1 , 3.87 m 3 s -1 and 4.11 m 3 s -1 in Olifants, Crocodile and Sabie rivers, respectively). This was the result of a steady decrease in flow over 1 month prior to sampling (Figure 1-A1).
The three rivers drain residential and agricultural areas and then flow eastwards through well-preserved ecozones in the park, except for the Crocodile River that delineates the southern border of the park, one of its bank draining the large agricultural zone of the Lowveld (Figure 1). The Crocodile River is considered highly threatened, draining untreated industrial and household wastewater from large cities such as Johannesburg and Tshwane, its water being increasingly used for irrigation and receiving agricultural effluents. The Olifants River drains the industrial fertiliser complex of Phalaborwa as well as phosphate rock mining facilities (foskorite and pyroxenite) before entering the park and is considered one of the most polluted rivers in South Africa (Gerber et al. 2015 and references cited therein). In contrast, the Sabie River is recognised as a remarkable hotspot of aquatic biodiversity and one of the most pristine river systems of South Africa (Riddell et al. 2019).
Three sites were selected in each river system, with the help of South African National Park officials and KNP game rangers. Each site was chosen according to its accessibility, the availability of riverine microhabitats (hard and soft substrates) and its coherence regarding our aim to examine the effects of an environmental gradient through the park. As far as possible, we selected river reaches that showed the most similar hydro-morphologies and without riparian canopy cover (e.g. Figure 2-A1). Upstream sites were located close (< 500 m) to the park's border, and to ensure substantial spatial coverage, each sampling site was located 10 km -60 km apart from each other ( Figure 1).

Sampling
Samples were collected over a 4-day sampling campaign in May 2019 during austral autumn, a period of low flow conditions (Figure 1-A1). To homogenise sampling, at each site, a 25-m reach was selected along which five 1 m 2 shallow (< 30 cm water depth) plots located approximately 5 m apart were set (e.g. Figure 3-A1). Each plot comprised both hard substrates (in the form of large bedrock or cobble stones) and soft substrates (mostly in the form of fine, sandy sediment). To avoid trampling disturbance, plots were always sampled from downstream to upstream (Figure 3-A1). Samples were placed at 4°C in the dark upon collection and further preserved with fixatives or frozen until further laboratory analyses could take place.
At each plot, 500 mL of water was collected for nutrient analyses; and water temperature, pH, conductivity, dissolved inorganic salts, total dissolved solids and turbidity were measured in situ with a portable probe (HI9813-5, Hanna Instruments Inc., Bedfordview, South Africa). Sediments were sampled to measure total organic carbon by pushing one 2-cm diameter Perspex corer (Figure 1-A1) through 5 cm of the upper sediment layer at one location chosen randomly within each plot. Sediments were sampled for meiobenthos by pushing two 2-cm diameter Perspex corers through 5 cm of the upper sediment layer at two different locations within each plot. The two sediment cores collected were merged in the same tube and the resulting sample was fixed in a final concentration of 4% buffered formaldehyde. A calibrated brush-sampler (Figure 3-A1; Peters et al. 2005) was used for underwater sampling of 3.14 cm 2 areas of biofilm (and its associated meiobenthos) growing on hard substrates such as rocks and cobbles. This process was repeated three times at different random locations within each plot, using a 20-μm-aperture nylon mesh sieve each time to concentrate the collected subsamples. The contents of the sieve were then poured into a tube and the resulting total sample (representing an area of 9.42 cm 2 ) was fixed in 4% buffered formaldehyde.  After meiobenthos sampling, a 'kick-sampling' procedure was used to obtain a semi-quantitative sample representative of the macrobenthic community dwelling in each plot.
Stones, macrophytes and superficial sediments were thoroughly agitated by hand and foot all over the plot for 3 min. The resulting suspended organisms were collected in a    sturdy hand net (500 μm meshes) held downstream of the plot. The content of the net was poured in a jar and preserved in 70% ethanol.

Sample processing
Water chemistry variables were measured in the laboratory using a Merck Spectroquant Pharo 300 Spectrophotometer and appropriate test kits following standard protocols (Merck KGaA 2014). Variables were measured on defrosted, unfiltered water samples and included Cl -(further referred to as chloride) (Merck protocol #14897), P-PO 4 3- . Subsamples of frozen sediments were defrosted, dried in an oven at 60°C for 96 h, weighted, burned to ashes at 600°C for 6 h, then weighted once again to measure the percentage of total organic content (TOC) as a ratio between ash-free dry weight and dry weight of the sediment.
To extract meiobenthos from sediment samples, a densitycentrifugation procedure was used in the laboratory involving the flotation of organic particles on Ludox HS-40 (specific gravity set at 1.14), following the method of Pfannkuche and Thiel (1988). The supernatant, containing the meiobenthic invertebrates, was poured through 20-μm meshes, preserved in 4% formaldehyde and stained with a few drops of Rose Bengal. The pellet (i.e. inorganic sediment particles from which the meiobenthos had been extracted) was further passed through a series of stacked sieves (1000 μm, 500 μm, 250 μm, 125 μm, 63 μm and 32-μm-aperture meshes). Each sediment size fraction retained on the respective sieve (further coined F1000, F500, F250, F125, F63 and F32) was dried in an oven at 60°C for 96 h and weighted to estimate the sediment particle size distribution of each sample.
Meiobenthic invertebrates were counted and identified under a 40-80x magnification stereomicroscope (Nikon SMZ1500). Morpho-taxonomic features were used to classify the meiobenthic individuals into coarse taxonomic groups: nematodes, rotifers, gastrotrichs, copepods and their Nauplii larvae, ostracods, oligochaetes, tardigrades, mites, molluscs (mostly veliger and pediveliger larval stages) and larval stages of plecopters, ephemeropters and dipters (detailing Chironomidae and Ceratopogonidae families). The weight of sediment from which meiobenthos had been extracted was used to express abundances per gram sediment dry weight. Meiobenthic organisms dwelling on hard substrates were not extracted using the Ludox procedure and they were directly counted and identified as described above, except that abundances were expressed per area of biofilm.
Individuals of macrobenthic invertebrates were counted in each sample under a 1-5x magnification dissection microscope (Nikon C-LEDS) and identified to lowest possible taxonomic level based on morphological characteristics and regional identification key for aquatic macro-invertebrates.
Total abundances were expressed as number of individuals per 1 m 2 plot. However, relative abundances were used for community structure analyses to avoid potential misinterpretations because of the semi-quantitative nature of the 'kick-sampling' procedure used to sample (Table 1-A1).
Effects of river identity (Olifants, Crocodile or Sabie) and reach location (upstream, midstream or downstream) on abiotic parameters, faunal abundances and macrobenthic diversity indices were tested using a two-way analysis of variance (ANOVA) method followed by a post-hoc Tukey's HSD test for pairwise comparisons. The univariate data were checked for normality and homoscedasticity using Shapiro-Wilk test and Levene's test, respectively. If data were not normally distributed and homoscedastic, the non-parametric Kruskal-Wallis rank sum test followed by a multiple comparison test was used instead. For macrobenthos data, effects of river identity and reach were tested on Shannon's diversity index (H') and Pielou's evenness (J').
Multivariate effects of river and reach on the benthic community structure were tested using permutational multivariate analysis of variance using distance matrices (PERMANOVA: adonis function with 9999 permutations).
The assumption of multivariate homogeneity of variances was first checked using the PERMDISP2 procedure (betadisper function), which is a multivariate analogous to Levene's test. The function rankindex was also used before tests to rank the performance of different dissimilarity indices. In our case, Bray-Curtis distance was the most relevant dissimilarity index.
Non-metric multidimensional scaling (NMDS) based on Bray-Curtis distances (metaMDS function) was used to ordinate samples and species scores. To highlight the most important environmental factors describing the variability in community structure of macrobenthos and meiobenthos, envfit function (that fits vectors of continuous environmental variables) was used. The significance of environmental variables (vectors) fitted onto the ordination was assessed using a goodness-of-fit (R 2 ) statistics and empirical p-values based on a permutation test (N = 999 permutations). The direction of vector showed the direction of the gradient, and the length of the arrow was proportional to the correlation between the variable score and the ordination space. Only significant vectors were further displayed on NMDS plot using 'p.max=0.05' as a filter argument when plotting vectors on the ordination space.
The response of individual taxon abundances to river identity and reach location was examined using multi-level pattern analysis after De Cáceres and Legendre (2009), which is an extended method of indicator species analysis. Multi-level pattern analysis (multipatt function) is a permutation test examining which taxa are significantly responsible for differences among a group of samples. Multi-level pattern analysis was chosen over other species indicator methods because it is less sensitive to the weight of over-dominant species, and it is not based on the relative contribution to differences (De Cáceres & Legendre 2009).

Ethical considerations
This article followed all ethical standards for research without direct contact with human or animal subjects.

Water and sediment environment
Overall, water pH was slightly alkaline but did not differ significantly across rivers or reaches (Table 1). However, conductivity was significantly higher in the Olifants River (on average 663 μS cm -1 ) in comparison to Crocodile (450 μS cm -1 ) and Sabie (104 μS cm -1 ) rivers (Kruskal-Wallis test: χ 2 = 38.7, p < 0.001). The concentrations of sulphate and chloride were likely the dominant contributors to variance in conductivity values as strong positive correlations were observed between conductivity and sulphate and chloride concentrations (Figure 4-A1). Turbidity was significantly lower in Sabie in comparison to Crocodile and Olifants rivers (ANOVA, F 2,36 = 31, p < 0.001), and in the latter two more turbid rivers, turbidity was found to be significantly higher at upstream sites (ANOVA, F 2,36 = 9.8, p < 0.001). Nitrate and ammonium concentrations did not show significant associations with river or with reach, although they were found to be positively correlated with phosphate concentration (Figure 4-A1). Phosphate was significantly more concentrated in the Crocodile River (ANOVA, F 2,36 = 37.1, p < 0.001) and interestingly its concentration tended to increase in downstream areas in all rivers (ANOVA, F 2,36 = 187.7, p < 0.001).
The sediment of the Olifants River showed a significantly greater proportion of organic matter in comparison to the other rivers (Kruskal-Wallis test: χ 2 = 8.2, p = 0.015). In terms of granulometry, relatively coarse sediment dominated (fractions > 500 μm represented on average 63.5% of the sediment grain size distribution). However, it is worth noting that the downstream site on Olifants River showed unusually low proportions of the largest grain size categories (Table 1). The proportion of coarser sediment was negatively correlated with the proportion of finer grain size categories (Figure 4-A1). The coarsest sediment fraction (> 1000 μm) showed a significant river*reach interaction (ANOVA, F 4,36 = 6.2, p < 0.001), essentially confirming that Sabie's riverbed showed coarser sediment as we moved downstream.

Abundance and diversity of benthic invertebrates
Macrobenthos was at least three orders of magnitude less abundant than the meiobenthos: on average, 130 macrobenthic individuals per m 2 against 332 000 ind. m -2 for biofilm-dwelling meiobenthos (Table 2). Macrobenthos abundance appeared significantly lower in the Sabie River (ANOVA, F 2,36 = 8.4, p < 0.001). The macrobenthic community consisted of 82 taxa from 53 families or coarser taxonomic groups. Fifty-seven taxa could be further resolved to genus or species level (Table 1-A1). The taxon dominating the macrobenthic community was Tarebia granifera, making up 36.95% of the macrobenthic community on average. However, T. granifera made up 73% -87% of individuals in upstream reaches (Figure 2a). As a result, reach position had a significant effect on the relative abundance of T. granifera (ANOVA, F 2,36 = 101, p < 0.001). The larval stage of the ephemeropteran genus Caenis sp. was the second most common macrobenthic taxon, representing 16.8% of the assemblage on average and being especially abundant in mid-and downstream reaches of the Sabie and Crocodile  (Table 2); however, taxon richness was significantly lower in upstream reaches in comparison to mid-or downstream reaches (Table 2, ANOVA, F 2,36 = 12.7, p < 0.001). Shannon's diversity index (H') was also lower in upstream reaches (Table 2) but was further found to be significantly lower in the Crocodile River as well (ANOVA, F 2,36 = 11, p < 0.001). Pielou's equitability index followed a similar pattern as H' with communities generally more equitable in mid-and downstream reaches (ANOVA, F 2,36 > 11.2, p < 0.001), except for the downstream reach of the Crocodile River, which was uneven and strongly dominated (87%) by Caenis sp. (Table 2, Figure 2a).
Meiobenthos abundance in sediment did not show any significant trend (Table 2). However, biofilm-dwelling meiobenthic organisms were significantly less abundant in upstream reaches (ANOVA, F 2,36 = 9.3, p < 0.001). Nematodes and rotifers were the numerically dominant representatives of the meiobenthos in all rivers and reaches (Figure 2b and Figure 2c), but while nematodes made up on average 21% of the meiobenthic assemblage both in the sediment and in biofilms, rotifers were relatively more abundant in biofilms (45.7% of individuals) than in sediments (31.9%; Table 2).

Structural responses of the macrobenthic community to environmental gradients
The structure of the macrobenthic community was significantly affected by river identity (PERMANOVA, F 2,42 = 29, R 2 = 0.12, p = 0.002). Multivariate homoscedasticity was met without data transformation (PERMDISP2, p = 0.72); however, the effects of reach location on macrobenthic community structure could not be tested as multivariate homoscedasticity of variances was not met, even after reasonable data transformation techniques. This was caused by a largely heterogeneous distribution of distances to centroids among reaches: samples from upstream reaches had a very skewed community structure because of the large dominance of T. granifera (Figure 2a), resulting in a very low dispersion of sample scores around group centroids (N = 15, average distance to centroid: 0.15 for upstream samples), while for mid-and downstream sections, average distances of samples to centroids were 0.51 and 0.54, respectively.
The multi-level pattern analysis identified 14 taxa significantly associated with one river group. Those groups are highlighted with stars on the NMDS biplot ( Figure 3). Namely, only macrobenthic nematodes were found to be significantly associated with the Crocodile River. The ephemeropteran larvae Tricorythus sp. and Elassoneuria sp., the trichopteran larvae Cheumatopsyche sp. and Macrostemum sp., the dipteran larvae of Simulium sp., the dragonfly larvae of Pseudagrion sp. and coleopteran adults of Berosus sp. were significantly associated with the Olifants River. The ephemeropteran larvae of Machadorythus maculatus, the hemipteran Laccocoris sp., the dipteran larvae of subfamily Chironominae, the dragonfly larvae of Lestinogomphus sp. and Notogomphus sp., and the coleopteran larvae of subfamily Larainae were significantly associated with the Sabie River. The clam Corbicula fluminalis was associated with both Crocodile and Olifants rivers. Dragonfly larvae Coenagrionidae and Paragomphus sp., and ephemeropteran larvae of the Leptophlebiidae family were associated with both Crocodile and Sabie rivers (Figure 3).
When performed to examine the associations with a specific river reach, the multi-level pattern analysis identified two taxa, namely, the larvae of Ceratopogonidae and copepods being associated with both mid-and downstream reaches ( Table 2).
The structure of the meiobenthic community dwelling biofilms growing on hard substrates was significantly affected by river identity (PERMANOVA, F 2,36 = 4.9, R 2 = 0.15, p < 0.001), river reach (PERMANOVA, F 2,36 = 5.1, R 2 = 0.16, p < 0.001) and by the river*reach interaction (PERMANOVA, F 4,36 = 1.8, R 2 = 0.11, p = 0.023). Multivariate homoscedasticity was achieved without data transformation (PERMDISP2 river , p = 0.79; PERMDISP2 reach , p = 0.65). This may have been caused by the relatively coherent pattern of lower abundances in upstream reaches coupled with relatively higher abundances of chironomids in the Sabie River (Figure 2b and Table 2). Concentration of ammonium in the water was found to be the strongest driver of biofilm-dwelling meiobenthic community structure (R 2 = 0.17, p =0.02; Figure 5), followed by TDS and conductivity (R 2 = 0.145, p = 0.035 and R 2 = 0.146, p = 0.032, respectively). The multi-level pattern analysis identified three taxa significantly associated with the Sabie River, namely, oligochaetes, tardigrades and gastrotrichs; NMDS, non-metric multidimensional scaling. while one taxon (Chironomids) was significantly associated with both Olifants and Sabie rivers ( Figure 5). When performed to examine associations with a specific river reach, the multi-level pattern analysis identified four taxa being associated with both mid-and downstream reaches, namely, ostracods and larval stages of Ceratopogonidae, Chironomidae and Ephemeropterans.

Discussion
Here we studied for the first time the entire assemblage of benthic invertebrates dwelling hard and soft substrates of three major rivers flowing through the KNP, South Africa. Overall, we observed a great diversity of benthic invertebrates of various sizes and morphologies belonging to a total of seven major animal groups: Arthropoda, Mollusca, Nematoda, Rotifera, Gastrotricha, Tardigrada and Annelida. We found the highest diversity of families in the Sabie River (42 families), followed by the Olifants (36 families) and Crocodile River (25 families). These results generally confirmed a better biodiversity status for the Sabie River as was already observed in previous studies focusing on macrobenthic invertebrates dwelling rivers of the KNP. For example, Muller and Villet (2004) identified 58, 51 and 46 families in the Sabie, Crocodile and Olifants rivers, respectively. The Crocodile and Olifants rivers are known to experience substantial anthropogenic pressures in comparison to the Sabie River (Riddell et al. 2019). We found different structures for both meio-and macrobenthic communities across the rivers, confirming our hypothesis

TDS Cond
NMDS, non-metric multidimensional scaling. that benthic communities would react to environmental gradients across rivers. Crocodile River's water had more phosphates and showed the largest alterations of benthic communities, suggesting that this system faced acute stress. Olifants River's water showed more sulphates and a higher conductivity and turbidity, and the sediment texture was finer housing higher amounts of organic material. But as Olifants ran through the KNP, diversity and evenness tended to improve. Sabie showed lower concentrations of nutrients and the most diverse assemblages.
However, we found that upstream sites located at the park border were consistently less diverse, supporting that the richness of benthic assemblages benefits from the course of the river through protected areas of the park. The Crocodile River, making the southern border of the park, did not show such improvement presumably because it receives effluents from intensive agriculture plots all along its southern catchment outside the KNP (Van der Laan, Van Antwerpen & Bristow 2012). One of the most conspicuous effects of the upstream-downstream gradient was the alteration of community structure because of the strong dominance of T. granifera at all upstream sites. This snail species originates from Asia and was first reported from South Africa in 1999 in northern KwaZulu-Natal (Appleton & Nadasan 2002). It has further rapidly spread into the Mpumalanga province, KNP and Swaziland and will probably continue its expansion to northern sub-tropical lowlands in Zimbabwe and Mozambique (Appleton, Forbes & Demetriades 2009). Fifteen years ago, Wolmarans and De kock (2006) sampled T. granifera from the Sabie and Crocodile rivers, but did not find it in Olifants. In our study we found substantial numbers of T. granifera in all river sites, including all river sites along the Olifants. T. granifera can attain high densities in invaded areas, thus impacting other benthic fauna and flora (Appleton et al. 2009;Jones et al. 2017). In invaded estuarine reaches, it has been shown that T. granifera mostly exploits algal biofilms and detritus (Miranda, Perissinotto & Appleton 2011). Thus, their grazing activity on hard substrates might particularly affect resource availability for biofilm-dwelling meiobenthos, as was evidenced empirically in the case of other aquatic snails dampening biofilm-associated assemblages (Peters, Hillebrand & Traunspurger 2007).
Here, we found a significant reduction of biofilm-dwelling meiobenthos in upstream reaches, reduction that could have been caused by a more severe biofilm grazing activity by T. granifera (we observed that ostracods and small instars of insect larvae were rare in the locations strongly dominated by T. granifera). However, stronger evidences are needed to show direct causality so we recommend experimental evaluation of the direct top-down or indirect bottom-up effect of the invasive T. granifera on local biofilmdwelling meiobenthos. Interestingly, our results showed that T. granifera did not over-dominate other downstream reaches in the KNP. A possible rationale is that the more diverse communities in mid-and downstream reaches could have reduced the invasion success of T. granifera through increased competition, predation, or parasitism (Miranda et al. 2016). This would also need further experimental evaluation. Nevertheless, our results contribute to support that the success of invasive species is all the more increased when indigenous fauna is impoverished because of stress like habitat deterioration (Van der Velde et al. 2002).
The positive association of sulphates with the Olifants River is not surprising as this river has historically suffered from salt enrichment because of sulphates that enter the river from effluent discharge which is one of the main contributors to the high salinity measurements in this river (Goetsch & Palmer 1997;Riddell et al. 2019). Furthermore, we found a positive association of phosphates in the Crocodile River as this system receives a great deal of organic effluents through fertilisers from agricultural activities taking place along its catchment outside the KNP ( Van der Laan et al. 2012). Similar results were obtained by Riddell et al. (2019) who attributed impacts of orthophosphates on rivers in KNP to agricultural activities taking place immediately adjacent to the park. In our study, increased sulphate and phosphate levels, along with other organic and inorganic variables emerged as significant drivers of the meio-and macrobenthic community structure. Biota such as nematodes, the water scavenger beetles Berosus spp., the Hydropsychid caddisfly Cheumatopsyche spp., Caenis spp. mayflies, dipteran Simulidae and Chironominae are adapted to tolerate a wide range of water quality parameters, such as organic enrichment and high salinity (Erasmus et al. 2021;Malherbe, Wepener & Van Vuren 2010), which is likely why they were associated with the Olifants and Crocodile rivers.
Habitat likely also played a role in the distribution of macrobenthos in the rivers as most of the biota specifically associated with either the Olifants or the Sabie river also have particular habitat and water flow preferences. For instance, the mayflies Elassoneuria sp. and Tricorythus sp. and the caddisflies Macrostemum sp. were significantly associated with the Olifants River presumably because of their habitat preference for rocky surfaces in fast-flowing water (Gerber & Gabriel 2002), a habitat that was more represented at sampling sites in this river. Similarly, the Sabie River habitats consisted of dense reed vegetation (Figure 2-A1), and most of the macrobenthos that were significantly associated with the Sabie River have affinities with such dense vegetation patches, particularly Laccocoris, or show affinities with the sandy to muddy substrate in slower flowing water at the edges of streams, particularly the gomphids Lestinogomphus sp. and Notogomphus sp. (Gerber & Gabriel 2002).
Although they are omnipresent in limnetic habitats, meiobenthic communities are poorly considered in ecology and conservation studies ). This knowledge gap is even more preoccupying for tropical and subtropical regions because studies about freshwater meiobenthos ecology and taxonomy have traditionally focused on temperate regions from the northern hemisphere (e.g. Fontaneto et al. 2012;Traunspurger et al. 2020;Zullini 2014). Here we found an abundant meiobenthic community dominated by nematodes and rotifers, and the taxonomic composition and abundance values reported in the present http://www.koedoe.co.za Open Access study are quite comparable to meiobenthic communities found in streams and rivers in Europe (e.g. Brüchner-Hüttemann et al. 2020;Majdi et al. 2012;Majdi, Threis & Traunspurger 2017). Some meiobenthic groups have been previously investigated in the KNP: for example, 30 years ago, Botha and Heyns (1991, 1992, 1993) identified a total of 33 nematode species (including species new to science) in sediment samples collected from Crocodile and Olifants rivers. Here we also observed an abundant nematode community, particularly in fine sediments of the Olifants River. In South African estuaries, Pillay andPerissinotto (2009) andNozais, Perissinotto andTita (2005) have observed that the abundance and community structure of the meiobenthos respond to river inflow and drought conditions (mostly through a decrease in physiologically sensitive taxa under harsher conditions). Here we found that biofilm-and sediment-dwelling meiobenthos exhibited different responses: nutrient and water quality proxies affected the biofilm-dwelling organisms, while sedimentdwelling meiofauna was also shaped by sediment texture and the availability of benthic organic material. Overall, the observed patterns are in agreement with previous studies showing that freshwater meiobenthic communities can become relevant bio-indicators of sediment pollution, water eutrophication or sediment texture (e.g. Haegerbaeumer et al. 2017;Schenk et al. 2020;Traunspurger et al. 2020). However, further studies should detail the response of dominant meiobenthic phyla (e.g. nematodes or rotifers) to the genus or species level in order to identify sensitive and tolerant species and develop subtler bioindication tools.

Conclusion
The KNP is known worldwide for its unique and emblematic terrestrial megafauna but less so for its remarkable freshwater fauna of benthic invertebrates. In this study we observed that the diversity of those invertebrates increased as we moved away from park borders, highlighting a potential beneficial effect of protected river reaches. Furthermore, in protected reaches, the assemblage was more even, less affected by the domination of the invasive snail T. granifera. The effects of this species on native assemblages should be further assessed, our results suggesting that T. granifera could affect the structure of biofilm-dwelling communities. Other drivers (like the concentration of nutrients in water and sediment granulometry) also showed significant effects on the distribution of benthic invertebrates. We recommend that the interplay between biotic and abiotic drivers should be further studied for a more comprehensive management and conservation of aquatic resources. Finally, our results stress that the Crocodile River showed the poorest and most unbalanced communities all the way. This suggests that further environment-damaging project on the Crocodile River (such as current applications for largescale coal-mining on the southern boundary of the park) could have catastrophic impacts on an already stressed river system.