Inference of Fish Community Assembly in Intermittent Rivers Using Joint Species Distribution Models and Trophic Guilds ()
1. Introduction
Non-random patterns of species co-occurrence often reflect both environmental filtering [1] [2] [3] and biotic interactions [4] [5] . However, there is less agreement on the relative importance of these processes across multiple spatial scales [6] [7] . In heterogeneous environments, the processes may synergistically affect species co-occurrence patterns [3] . For example, plant species may negatively associate under favourable abiotic environmental conditions, reflecting competitive exclusion, and positively associate under less suitable conditions when fitness is more strongly influenced by tolerance limits to abiotic factors [8] [9] [10] . In animal communities, predation can result in reduced co-occurrence of predators and prey when prey populations become depleted during periods where and in places where predators are abundant, or when prey is able to avoid areas with predators [11] . Positive species associations may occur when predators are attracted to areas where prey is abundant [12] [13] .
Fish community assembly in rivers varies seasonally with the magnitude, duration, and timing of flood pulses. Biotic interactions tend to be stronger during low-water periods when fishes are concentrated within shrinking aquatic habitats [14] [15] [16] . For example, top-down control of basal resources by benthivorous fish is much stronger during the annual low-water period in piedmont and lowland rivers in Venezuela [17] [18] . Reference [19] observed increasingly non-random patterns of species co-occurrence when species density was higher in a tropical lowland river. The mechanisms that drive species co-occurrence remain poorly documented for rivers, especially within semiarid regions. Rivers, including intermittent rivers in semiarid regions, support important ecosystem services, and often support rich biodiversity [20] [21] . Intermittent rivers experience extremes in discharge, ranging from flooding to complete desiccation. During dry periods, aquatic habitat is reduced and often disconnected, resulting in a mosaic of habitat types with wide variations in physicochemical characteristics and resources [22] . Isolated habitats within a given reach may differ in water residence time as a function of channel geomorphology and bed material [23] .
Evaluation of non-random patterns of community structure is generally based on comparisons between empirical data with distributions generated by null models (e.g. [24] [25] [26] [27] ). To infer mechanisms of community assembly, some ecologists use Joint Species Distribution Models (JSDMs) with species abundance and environmental data [28] [29] [30] . A particularly appealing feature of this approach is the capability to separate environmental filtering from biotic interactions within one model fit [3] [31] .
When a JSDM is fitted without environmental covariates, the model resembles unconstrained ordination [31] . If environmental variables are included in the model, the resulting residual species correlations may represent species co-occurrences without the direct influence of species responses to measured environmental covariates [32] . This means that in a constrained JSDM, the influence of environmental covariates on species occurrence is removed/controlled for [5] . However, residual species correlations may be either a result of biotic interactions or species responses to other covariates that were not included in the model [5] [31] [32] . Therefore, it is not possible to determine the ecological mechanisms that drive residual species associations from these models alone [31] . Interspecific associations depend on several limiting factors and trade-offs, such as traits of prey and predators that influence encounters and capture success and the manner in which environmental variation influences habitat and resources [33] . For this reason, ecologists often analyse species co-occurrence in conjunction with additional species information, such as habitat preferences, migratory patterns, phylogenies, and guild membership [34] [35] . This approach to analyzing species co-occurrence facilitates inferences about the mechanisms driving non-random patterns of co-occurrence [36] .
Designation of guild membership in conjunction with JSDMs can facilitate the inference of probable drivers of co-occurrence patterns in heterogeneous environments. For example, a negative correlation between two species within the same trophic guild could suggest avoidance of interspecific competition [4] [33] . Similarly, a negative correlation between a predator and potential prey species could suggest either predator depletion of local prey (e.g. within an isolated aquatic habitat) or prey avoidance of predators (e.g. within a refugium when prey are capable of dispersal and habitat selection) [11] . Conversely, positive co-occurrence of predatory species and potential prey may indicate predators select areas with abundant prey. Further, the use of guild membership may help with inferences pertaining to those interspecific interactions that are due to unmeasured environmental covariates. For example, a negative co-occurrence involving non-predatory species from different trophic guilds could arise from species differences in tolerance to abiotic conditions. Similarly, inter-guild positive co-occurrence involving non-predatory species could be a result of species aggregating within suitable habitats [4] . Therefore, using trophic guilds in conjunction with JSDMs may help to disentangle the residual correlations and tease apart putative biotic interactions from those interactions that may be influenced by unknown environmental factors.
In this study, we used JSDMs to infer the relative roles of environmental filtering and biotic interactions in the assembly of fish communities from intermittent rivers in the lower Okavango Delta. Previous studies in these rivers suggested that the magnitude of the annual flood pulse has a pervasive influence on the structure of local fish assemblages [37] [38] . However, those studies did not investigate the influence of water level on interspecific associations. Here, we assessed patterns of covariation in species responses to habitat type and physicochemical water parameters. We measured and compared physicochemical water parameters at the mesohabitat scale during multiple stages of the annual flood pulse in the lower Delta. To infer mechanisms influencing species co-occurrence, we used JSDMS to factor out correlations with abiotic environmental variables and then analyzed residual correlations to interpret patterns of species co-occurrence in relation to fish trophic guilds.
We hypothesized that positive residual correlations between species would mainly occur within the same ecological guild, indicating that species tend to aggregate in their preferred mesohabitats. Negative species co-occurrence was expected between predators and their potential prey, an indication of either prey depletion or prey avoidance within isolated habitats.
2. Materials and Methods
2.1. Study Area
The Okavango Delta is located in the semi-arid northwestern region of Botswana (Figure 1). It is characterized by low annual rainfall, which averages around 490 mm [39] . After entering Botswana at Mohembo, the Okavango River forms a relatively narrow perennial channel in the panhandle, which extends for approximately 100 km. In the lower panhandle, the channel branches form an alluvial fan with an extensive network of channels and swamps that cover approximately 13,000 km2 during high floods [40] . Most of the outflow (2%) from the delta into the intermittent rivers occurs through the Boro channel, which forms a confluence with the Thamalakane River at Matlapaneng, eventually draining into Lake Ngami and the Makgadikgadi salt pans [40] . Water temperature and conductivity increase along the fluvial gradient from Mohembo to the lower delta, whereas pH, dissolved oxygen, turbidity and total suspended solids decrease [41] .
![]()
Figure 1. Map showing locations of the study sites in the lower Okavango Delta.
2.2. Fish Sampling and Environmental Descriptors
Fishes were sampled from four sites in the lower Okavango Delta: Borojunction (BJN; 19˚55'01.0''S; 023˚30'52.0''E) a site near the Botswana Defence Force Camp (BDF; 19˚52'28.1''S; 023˚26'45.0''E), Dikgathong (DKN; 20˚08'27.9''S; 023˚22'38.5''E) and Chanoga (CHN; 20˚10'02.4''S; 023˚39'33.6''E), using a multifilament gillnet comprised of eleven panels that were 2.4 m high and 10 m long with different mesh sizes. Mesh sizes were 12 mm, 16 mm, 22 mm, 28 mm, 35 mm, 45 mm, 57 mm, 73 mm, 93 mm, 108 mm and 150 mm. Fish were sampled during the wet (August 2017 to November 2017) and dry (December 2017 to February 2018) seasons. During each survey, the gill net was set once at each site for approximately 12 h between 1800 h and 0600 h the following day. Fish sampling resulted in a total of four gill net samples from each site (two from the low-water season and two from the high-water season) which produced a total of 16 gill-net samples for analysis. Water persisted throughout the sampling period at CHN and DKN, and these sites were classified as permanent habitat. BDF and BJN dried out completely by the end of the dry season, and they were classified as ephemeral habitat.
To facilitate analysis of species response to abiotic environmental conditions at the survey sites, we measured several physicochemical water parameters using a handheld water quality meter at approximately 0600 h. Measured physicochemical water parameters were: dissolved oxygen (measured to the nearest 0.01 mg/l), pH (measured to the nearest 0.01), conductivity (measured to the nearest 0.1 µS/cm) and temperature (measured to the nearest 0.1˚C). We obtained water-level data (measured to the nearest 0.01 m) for each site from the Department of Water Affairs in Maun, Botswana.
2.3. Variation in Abiotic Environmental Factors
First, we summarized the major abiotic environmental gradients associated with the survey sites and dates using a Principal Coordinate Analysis (PCoA). Water physicochemical variables included in the analysis were depth, temperature, Dissolved Oxygen (DO), pH, and conductivity. Prior to inclusion in analyses, values were standardized by taking the Z-score to achieve the same range of variation for each variable. We used Euclidian distance to characterize variation in environmental variables between samples and to generate the distance matrix used as input for PCoA. To identify parameters with the strongest contribution on the abiotic environmental gradient, we correlated each parameter with the scores of the first and second axes of the PCoA. We retained physicochemical parameters with correlation coefficients > 0.6 for further analyses.
Because we were interested in understanding the association of physicochemical water parameters with hydrological seasons and habitat types, we used Redundancy Analysis (RDA). Using a subset of the most influential parameters from the PCoA (DO, depth, and temperature) as response variables, RDA was conducted with hydrological season and habitat type serving as explanatory variables. We tested for differences among physicochemical parameters in relation to hydrological seasons and habitat types using permutation tests after 9999 randomizations.
2.4. Joint Species Distribution Modeling (JSDM)
We collected a total of 8289 fish specimens representing 30 species during the study. Seven species: [Hippopotamyrus ansorgii (3 individuals), Pharyngochromis acuticeps (2 individuals), Parauchenoglanis ngamensis (2 individuals), Enteromius bifrenatus (2 individuals), Enteromius spp. (2 individuals), Clarias theodorae (1 individual) and Enteromius barnardi (1 individual)] had sample sizes < 5 and therefore were removed from the analysis. This resulted in 8276 specimens and 23 species included in the dataset for analysis.
Four non-redundant environmental variables were selected for use as covariates in the analysis. These variables included water physicochemical variables (depth, temperature, DO), and a categorical variable describing habitat type (i.e. ephemeral vs. permanent habitats) that should be associated with fish responses to drought conditions. Hydrological season was not included as a covariate for two reasons. First, hydrological season is redundant with water depth, a metric that provides greater resolution for capturing temporal and spatial variation. Second, hydrologic season represents the temporal component that we address in a different manner (see below).
We used JSDMs to investigate species interactions and individual species responses to environmental covariates. We applied correlated response models that regressed separate Bayesian generalized linear models for each species against environmental covariates and latent variables at the same time. Latent variables are random variables that define the underlying gradients in the data, and they are estimated during model fitting [32] . Therefore, the model simultaneously estimates relationships between multiple species and environmental covariates while approximating species responses to latent variables. This makes it possible to partition the observed correlation into that which is due to species responses to known environmental covariates and that which is unexplained and may be interpreted as indicative of species interactions [28] [31] . The interpretation of species responses to latent variables is not always straightforward, because species may respond to unknown environmental covariates or those that were not measured during the study [31] [42] [43] . To account for the influence of missing environmental covariates in the estimation of species co-occurrence, we included multiple random effects at the level of a sample, assuming that the heterogeneity driving variation in residual correlation at each sampling site was environmental. This ensured that environmental heterogeneity did not influence latent variables and the estimation of residual correlations [44] . To account for temporal fluctuations in data, we considered sampling date as an ordinal variable, and included it in the model as a random effect. This prevents potential dependency of samples taken repeatedly from the same site and ensures that the estimated underlying gradients of species abundances are true estimates of interspecies associations. We also explored potential site-level effects by including the site identity in the model as another random effect. However, the inclusion of this random effect did not change the underlying species responses and it remarkably decreased the explanatory power of the model. For this reason, site identity was excluded from the final model.
2.5. Species Responses to Abiotic Environmental Factors
The JSDM was fitted to species abundance in R [45] using the package boral, which employs Monte Carlo Markov Chain methods for estimating model parameters [28] . Species responses to covariates were modeled assuming Poisson responses with a log link function to control for positively skewed variation in the distributions of species abundances [28] . Simultaneous interspecies associations were approximated using 2 latent variables, which seemed to be more parsimonious than including additional latent variables. The linear model assumptions of normality and homoscedasticity were inspected visually with frequency plots of the model residuals and scatterplots of the residuals versus model predictions. Following [28] , we also employed Dunn-Smythe residual plots to assess additional evidence of violation of model assumptions. All covariates were centered and standardized prior to model fitting, so that coefficients could be compared with each other in terms of relative importance. All model parameters were estimated assuming non-informative normal and uniform priors as the default settings of boral package [28] . We fitted the models by running three chains of 500,000 iterations with the first 10,000 discarded as burn-in. The remaining iterations were thinned by a factor of 30. Model convergence was inspected visually with posterior trace plots of the model parameters. We assumed convergence of MCMC chains if the different chains yielded identical results, if there was no apparent autocorrelation between the sequential interactions in a chain, and if the chains reached a stationary distribution [46] . Species response to each explanatory variable was assessed by inspecting the magnitude and direction of posterior median coefficients for each species and their 95% credible intervals. Variables were considered significant when their corresponding 95% credible intervals did not encompass zero.
2.6. Inference of Interspecific Interactions
To infer potential species interactions, we derived a partial correlation matrix by inverting the residual species correlations obtained from the fitted model. Partial correlations were used because they allow for measuring the influence of one species after controlling for the effects of all others and, therefore, they are more informative about causal interspecies correlations than raw correlations [43] . We tested species partial correlations by inspecting the magnitude and direction of posterior median associations between each species pair and their 95% credible intervals. Species partial correlations were assumed to be significant when their corresponding 95% credible intervals did not encompass zero.
We categorized potential prerequisites for positive and negative species interactions according to [4] and [33] (Table 1). We assigned species to trophic guilds (Table 2) following [47] . Trophic guilds were: 1) Invertivores (feed mainly on insects and other small organisms, such as small snails, diatoms, detritus and mollusks), 2) Predators (feed mainly on insects, crustaceans and fish), 3) Predator/omnivores (feed on fish, insects, shrimps, snails, plant seeds, and fruit), 4) Herbivores (feed mainly on water-plants and algae), and 5) Detritivores (feed mainly on detritus). Types of biotic interactions were inferred by assessing intra and inter-guild pairs [35] . Gillnets were set in deep areas within the main river channel where juveniles are less likely to occur, and the relatively uniform size distributions of captured fish for most species indicated that nearly all were adults (Table 2).
2.7. Variance Partitioning
The variance explained by the model was partitioned as a means to estimate the proportional contribution of environmental filtering, potential biotic interactions, and random processes explaining variation in species abundances. Contribution of environmental filtering was estimated by the proportion of variance attributed to components related to the environmental covariates. The contribution of potential biotic interactions was derived from the variation associated with species responses to latent variables. Variation partitioning was performed separately for each species in the model, and the partitioning for the whole model was obtained by averaging species partitions. The details of this variation partitioning are given in [5] .
3. Results
3.1. Variation and Relationships among Abiotic Environmental Factors
Depth, DO, and temperature were the most important water physicochemical parameters showing high loadings on the PCoA axes (correlation coefficients > 0.6). As expected, variation in water physicochemistry was associated with hydrological season (F = 7.87; DF = 1; P = 0.001) and habitat type (F = 4.06; DF = 1;
![]()
Table 1. Types of interactions, direction, and the inferred assembly mechanisms. Inferences were based on Stachwicz (2001) and Kneitel and Chase (2004).
![]()
Table 2. Sample size (n), mean length (lower [LCI] and upper [UCI] confidence intervals), habitat preference, and feeding guilds of species that were analyzed in this study. Biological and ecological descriptions were obtained from Skelton (2001). Feeding guilds were categorized according to Taylor et al. (2017).
P = 0.02). Permanent habitats during the high-water period were distinct from ephemeral habitats during the low-water period along the first RDA axis. DO and depth were positively associated with high water and permanent habitat, and temperature was positively associated with ephemeral habitat and low water (Figure 2).
3.2. Species Responses to Abiotic Environmental Factors
There were no significant associations among species and habitat types (Figure 3(a)). Six species (Coptodon rendalli, Enteromius paludinosus, Enteromius poechi, Hepsetus cuvieri, Oreochromis macrochir, and Sargocromis carlottae) were significantly associated with shallow water, whereas three species (Brycinus. lateralis, Petrocephalus okavangoensis, and Schilbe intermedius) were significantly associated with deep water (Figure 3(b)). A relatively large proportion of species (57%) (B. lateralis, Clarias gariepinus, C. rendalli, E. poechi, H. cuvieri, Marcusenius altisambesi, O. macrochir, Serranochromis angusticeps, S. carlottae, Sargochromis codringtonii, S. intermedius, Serranochromis macrocephalus, and Synodontis sp.) were significantly associated with warm water, and no species showed significant association with cool water (Figure 3(c)). Only one species (C. gariepinus) was associated with oxygen-poor water. Similarly, one species (E. paludinosus) was significantly associated with oxygen-rich water (Figure 3(d)).
3.3. Inference of Interspecific Interactions
The model was complex, and our sample size was insufficient to achieve statistical power needed to detect significant partial correlations. For this reason, inference of interspecific interactions was based on the analysis of the residual correlation matrix, which resulted in three significant negative residual species correlations (P. okavangoensis vs. C. rendalli, P. okavangoensis vs. H. cuvieri,
![]()
Figure 2. Redundancy analysis showing the relationship between environmental variables, habitat type, and seasonal variation in hydrology.
![]()
Figure 3. Estimated posterior median coefficients for each species correlation with habitat type, depth, temperature and dissolved oxygen with 95% credible intervals. Negative significant species correlations are shown in red and significant positive correlations are shown in blue. (a) Habitat; (b) Depth; (c) Temperature; (d) Dissolved oxygen.
and P. okavangoensis vs. O. macrochir). All positive residual correlations were not significant.
3.4. Variance Partitioning
More than half of the variation in community composition (51%) was attributed to components related to environmental covariates (Figure 4). A relatively small proportion (4%) of the variation was associated with species responses to latent variables. Remaining variation was due to random effects.
4. Discussion
We infer that fish community assembly in the lower Okavango Delta was primarily driven by environmental filtering and random processes with a lesser contribution from biotic interactions. The relatively large contribution of random processes may be attributed to ecological drift [38] [48] . Analysis of species co-occurrence patterns from JSDMs in conjunction with fish trophic guilds suggests that biotic interactions are influenced by abiotic environmental conditions, with a relatively small contribution attributed to predation. Our findings are consistent with the idea that stream fish community assembly, at the intra-annual scale, is dominated by environmental filtering and predator-prey interactions with a weaker influence from interspecific competition [49] [50] .
Two small bodied species (B. lateralis and P. okavangoensis) and a relatively large omnivorous predator (S. intermedius) were significantly associated with deep water. One piscivore (H. cuvieri) and five potential prey species (C. rendalli, E. paludinosus, E. poechi, O. macrochir, and S. carlottae) were significantly associated with shallow water. Enteromius paludinosus also was associated with oxygen-rich habitats, suggesting that dissolved oxygen concentration could have determined the distribution this species within shallow-water habitats. However, many species in both shallow (5/6) and deep water (2/3) also were associated
![]()
Figure 4. Variation partitioning for the model fitted to explain the variation in species abundance in response to habitat type, depth, temperature and dissolved oxygen. Random effects of sampling date were included as residual variation.
with higher temperature, suggesting a strong role of temperature in determining species distribution within both deep and shallow water habitats in the lower Okavango Delta. Our results are consistent with literature reporting the strong influence of temperature on fish spatial distributions at temperate latitudes [51] [52] [53] [54] . Reference [55] proposed similar responses of fish species to temperature and food resources. Warmer aquatic habitats in the lower Okavango delta may provide not only more primary and secondary production, but also better conditions for growth. The significant association of C. gariepinus with hypoxic conditions could be related to its ability to endure harsh conditions via accessory aerial respiration [56] .
Analysis of patterns of species co-occurrence based on species residual correlations after accounting for environmental factors resulted in three significant negative correlations. One negative correlation involved a predator and a potential prey species, which is consistent with the idea that predation influences species spatial segregation [11] [27] . In tropical rivers, one or a few species of predatory fishes can exert a disproportionately large top-down effect on prey populations [57] , by direct mortality and/or indirectly via avoidance of mortality risk [58] [59] [60] . When flows diminish, habitat conditions in intermittent rivers of the lower Okavango Delta become dominated by shallow marginal areas supporting growth of aquatic macrophytes that provide refuge from predators. Increase in habitat structural complexity may result in spatio-temporal variations in water chemistry and nutrient concentrations that may drive species spatial segregation [61] [62] [63] . Interestingly, the predator H. cuvieri occupied such areas and had a negative correlation with P. okavangoensis, possibly inducing this small mormyrid to seek refuge elsewhere in deeper waters. Adult size classes of H. cuvieri were reported to feed mainly on cichlids and mormyrids in the Zambezi River floodplain [64] . Our results suggest that predation could drive negative patterns of species co-occurrence, either by predators locally depleting prey or by prey avoiding predators [57] [58] [65] [66] . This decoupling of predator-prey interactions may be facilitated by habitat complexity, allowing coexistence and increasing ecosystem stability [67] .
Two out of three of the negative species correlations involved non-predatory fishes from different trophic guilds (insectivorous P. okavangoensis vs. herbivorous C. rendalli, and P. okavangoensis vs. detritivorous O. macrochir). Negative correlations among non-predatory species from different trophic guilds imply that spatial segregation was influenced by species differential response to abiotic conditions. This inference is buttressed by the analysis revealing that P. okavangoensis was associated with deep-water habitats, whereas C. rendalli and O. macrochir were associated with shallow-water habitats. Water depth was found to be the principal factor facilitating co-existence in the anuran pond communities [48] . In that case, the inferred mechanism was interspecific variation in breeding site selection. Similarly, habitat selection appeared to be more important than biotic interactions in structuring fish communities in the Macau River, Brazil [34] , and functional redundancy increased as the number of co-occurring species increased at all spatial scales in the Patos Lagoon, Brazil [36] . Our results support the general view that distribution patterns of coexisting species in wetlands are largely derived from species responses to spatial variation in abiotic conditions, i.e. environmental filtering [68] [69] , especially in regions that experience alternating wet and dry conditions [70] .
5. Conclusion
By relating residual species correlations from JSDMs to fish trophic guilds, we were able to infer mechanisms influencing community assembly. However, it was important to note that we derived residual species correlations from latent variables, and these correlations could have been influenced by additional abiotic factors that we did not consider in this study. Nonetheless, our results demonstrate that analysis of residual species correlations from JSDMs, in conjunction with fish trophic guilds, may help ecologists to disentangle residual correlations and tease apart putative biotic interactions from the interactions driven by unmeasured environmental factors in local fish assemblages of intermittent rivers.
6. Recommendations
Environmental management should emphasize the maintenance of a relatively natural environmental flow regime to sustain suitable habitats that can support the persistence of ecological mechanisms that structure fish communities in the lower Okavango Delta. Moreover, investigation of food web ecology using field surveys, dietary analysis, field experiments, and ecological modeling is necessary to provide more insights into how biotic and abiotic factors influence species distributions and co-occurrence in this system.
Acknowledgements
This study was partially supported by the Botswana International University of Science and Technology, through a graduate fellowship awarded to Thethela Bokhutlo. Field sampling was achieved through a director’s permit no. 00000654A obtained from the Department of Wildlife and National Parks.