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²) and Sabie rivers (6252 km²) start off in the Drakensberg region. The Olifants River has the largest catchment area (54 434 km²) 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³ s⁻¹, 3.87 m³ s⁻¹ and 4.11 m³ s⁻¹ 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).

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² 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² 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²) 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.

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₄^3-(phosphate) (#14848), SO₄^2- (sulphate) (#14791), N-NH₄⁺ (ammonium) (#14752), N-NO₃^- (nitrate) (#14773). 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 density-centrifugation 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² 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).

Statistical analyses were performed in R (v 4.0.3.; R Development Core Team 2020) using packages ‘vegan’ and ‘indicspecies’ (De Cáceres & Legendre 2009).

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²) 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).

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

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⁻¹) in comparison to Crocodile (450 µS cm⁻¹) and Sabie (104 µS cm⁻¹) rivers (Kruskal–Wallis test: χ² = 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: χ² = 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.

Macrobenthos was at least three orders of magnitude less abundant than the meiobenthos: on average, 130 macrobenthic individuals per m² against 332 000 ind. m⁻² 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 rivers (Figure 2a, Table 1-A1). Other macrobenthic taxa collectively represented 46.2% of the assemblage on average (Figure 2a), with only 14 taxa making up individually more than 1% of the assemblage (Tricorythus sp., Corbicula fluminalis, Paragomphus sp., Acanthiops sp., Leptophlebiidae, Cheumatopsyche sp., Notogomphus sp., Hydroglyphus sp., Laccocoris sp., Pseudagrion sp., Simulium sp., Chironominae Anisops sp., Elminae larvae; see Table 1-A1). On average, 9.2 different macrobenthic taxa were identified from each sample (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).

The structure of the macrobenthic community was significantly affected by river identity (PERMANOVA, F_2,42 = 29, R² = 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.

Phosphate was found to be the first driver of macrobenthic community structure (R² = 0.29, p < 0.001), followed by variables associated with sediment granulometry: F125 (R² = 0.24, p = 0.002), F1000 (R² = 0.21, p = 0.009), F250 (R² = 0.20, p = 0.01) and F63 (R² = 0.15, p = 0.034). Sulphate was also a significant driver (R² = 0.14, p = 0.036) although its direction was opposed to that of phosphate (Figure 3).

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).

The structure of the meiobenthic community dwelling sediments was significantly affected by river identity (PERMANOVA, F_2,36 = 2.3, R² = 0.08, p = 0.003), multivariate homoscedasticity being met without data transformation (PERMDISP2_river, p = 0.21; PERMDISP2_reach, p = 0.4). It was significantly affected by the interaction river*reach (PERMANOVA, F_4,36 = 2.3, R² = 0.17, p < 0.001). This was caused by the high nematode and rotifer abundances found in the downstream reach of Olifants River (Figure 2b). Sediment TOC was found to be the first driver of sediment-dwelling meiobenthic community structure (R² = 0.38, p < 0.001; Figure 4), followed by variables associated with fine sediment granulometry: F63 (R² = 0.34, p < 0.001), F32 (R² = 0.28, p < 0.001). Water quality variables such as total dissolved solids (TDS) and conductivity (Cond) were also significant drivers (R² = 0.18, p = 0.012, and R² = 0.17, p = 0.017, respectively). The multi-level pattern analysis identified only one taxon: nematodes, as being significantly associated with the Olifants and Sabie rivers (Figures 2b, Figure 4 and Table 2). 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² = 0.15, p < 0.001), river reach (PERMANOVA, F_2,36 = 5.1, R² = 0.16, p < 0.001) and by the river*reach interaction (PERMANOVA, F_4,36 = 1.8, R² = 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² = 0.17, p =0.02; Figure 5), followed by TDS and conductivity (R² = 0.145, p = 0.035 and R² = 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; 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.