Disentangling the effect of climatic and hydrological predictor variables on benthic macroinvertebrate distributions from predictive models

Lotic freshwater macroinvertebrate species distribution models (SDMs) have been shown to improve when hydrological variables are included. However, most studies to date only include data describing climate or stream flow-related surrogates. We assessed the relative influence of climatic and hydrological predictor variables on the modelled distribution of macroinvertebrates, expecting model performance to improve when hydrological variables are included. We calibrated five SDMs using combinations of bioclimatic (bC), hydrological (H) and hydroclimatic (hC) predictor datasets and compared model performance as well as variance partition of all combinations. We investigated the difference in trait composition of communities that responded better to either bC or H configurations. The dataset bC had the most influence in terms of proportional variance, however model performance was increased with the addition of hC or H. Trait composition demonstrated distinct patterns between associated model configurations, where species that prefer intermediate to slow-flowing current conditions in regions further downstream performed better with bC–H. Including hydrological variables in SDMs contributes to improved performance, it is however, species-specific and future studies would benefit from hydrology-related variables to link environmental conditions and diverse communities. Consequently, SDMs that include climatic and hydrological variables could more accurately guide sustainable river ecosystem management.


Introduction
Species distribution models (SDMs) are ecological predictive models, increasingly used to inform and complement large-scale distribution analyses to aid conservation efforts (Araújo et al., 2011;Guisan et al., 2013;Eaton et al., 2018). In river ecosystems, SDMs have been less often applied due, in part, to (1) complex interactions between the numerous driving factors of river systems and (2) insufficient abiotic and biotic data to effectively describe hydro-ecological relationships in streams.
The hydrological regime is said to be the ''master variable'' (Power et al., 1995) of lotic habitats and critical to the ecological stability of river ecosystems . It is highly variable, both spatially and temporally, making it a core driver of the physical structure of river habitats and a regulator of species distribution and abundance (Resh et al., 1988;Poff et al., 1997). For many river species, flow regimes are very important for, at least, part of their life history (Lytle & Poff, 2004). Numerous macroinvertebrate species depend on flow-related cues that directly or indirectly initiate, for example, breeding period (Hancock & Bunn, 1997), development (Gray, 1981) and emergence & metamorphosis (Peckarsky et al., 2000;Lytle, 2002). Therefore, species have evolved to the heterogeneity in rivers caused by the variability in flow regime, e.g. average/low/high flows, intermittent and ephemeral flows. It is essential to understand the influence stream flow has on the distribution of species, so that suitable recommendations can be made to restore or conserve river systems successfully. This is particularly important with the increasing changes in hydrological regime due to climate change, e.g. severity and frequency of floods and droughts (IPCC, 2007;Döll & Zhang, 2010;Kiesel et al., 2019;Gudmundsson et al., 2021).
There are still only a limited number of studies that directly investigate the influence of flow regime on riverine species distribution. Some recent attempts have been made to include, at least, some aspects of flow regime, e.g. high flow days (Kuemmerlen et al., 2015a) as well as aggregated flow statistics, e.g. mean annual flow Pyne & Poff, 2017), which were shown to be of high relevance to macroinvertebrate distribution. With data describing climate and hydrology becoming increasingly available it becomes possible to include both aspects to SDMs. However, the question remains, if including both will affect model results, leading to improved predictive accuracy and/or differing predicted ranges? Considering both climate and hydrology variables in SDMs would possibly allow for a distinct examination of the abiotic drivers of species distribution and facilitate management decisions to develop wise conservation or restoration strategies.
We hypothesized that SDMs including hydrological variables next to commonly used climatic variables would improve SDM performance (H1). To investigate this, we applied SDMs on a community of benthic macroinvertebrates with combinations of three predictor datasets including (1) climate only variables, (2) hydrology only variables and (3) both climate and hydrology-related variables. We evaluated the influence of each predictor dataset on species' distribution by investigating individual variance as well as shared variance explained by each dataset. We also compared differences in the functional traits; current preference and stream zonation preference. By investigating trait composition of, for example, species' that perform better with hydrology variables, we aimed at better interpreting the model's results, from a perspective closely related to river hydrology. We investigated which model configuration performed better by evaluating differences in model performance. Further, we analysed how the choice of predictor datasets influences the predicted species distributions. By comparing differences in model performance, functional traits and explained variance, we expected to determine the individual influence of both climate and hydrology datasets, and to what extent these datasets influence predicted species' distributions. If combinations of datasets that included hydrology outperformed the climate only combinations, we would accept H1. If performance across datasets was similar, H1 would be rejected.

Study area
The study area comprised the Ems (17,934 km 2 ) and Weser (46,306 km 2 ) catchments located in Germany ( Fig. 1). Selecting two catchments fully within Germany ensured consistent data availability. The two catchments are adjacent to each other and both are split into two distinct ecoregions: central plains (lowland) and central highlands (mountain, sensu Illies 1967).
The stream network for the study area was based on a layer in GeoTIFF raster format of a modeled 1 km 2 gridded stream network with a total of 13,749 cells. The network was downloaded from earthenv.org/ streams, which was derived from Hydrosheds (www. hydrosheds.org; 30-arc-s spatial grain; Lehner et al., 2008), which in turn is based on the SRTM dataset (srtm.csi.cgiar.org; Jarvis et al., 2008). This spatial resolution was applied due to the requirement of spatially analogous environmental variables in SDMs, i.e. hydrology, bioclimate and hydroclimate at 1 km 2 (Araújo et al., 2019).

Biological data
Macroinvertebrate species datasets were obtained from German Federal State agencies. These macroinvertebrate samples were collected using a 25 9 25 cm hand net (500 lm mesh size), following the AQEM STAR sampling methodology, in which samples consist of 20 microhabitats, sampled based on their relative cover at the sampling site (AQEM 2002;Haase et al., 2004). To be included in the study, each species had to be identified to species level and to avoid issues with modelling a small sample size, each species had to have at least 20 occurrences within the study area (Stockwell & Townsend Peterson, 2002). A total of 91 species occurrences at 1,258 sites remained from this process, in a presence only format. Sampling frequency occurs in a 3-year cycle, all sites were sampled at least once within the period between 2005 and 2013.

Environmental predictors
The predictor variables for bioclimate, hydroclimate and hydrology used in this study were all openly sourced data and freely available to the user. All datasets are available in raster GeoTIFF format.

Bioclimate
There are 19 bioclimatic variables openly available from worldclim.org (Hijmans et al., 2005;Lehner et al., 2008;Fick & Hijmans, 2017) describing temperature and precipitation. These data are applied frequently in SDMs and other predictive modelling studies. The variables were downloaded at a 1 km (30arc-s) spatial resolution. All 19 variables were masked to the base layer 1 km 2 stream network described above. The bioclimate data were measured at local grid cell scale and represents our measurement of climate.

Hydroclimate
The hydroclimate variables were downloaded from earthenv.org/streams (Domisch et al., 2015a) and are based on, among others, the 19 bioclimatic variables described above. However, it differs in that the bioclimatic information is related to the stream network by accumulating information from the upper subcatchment in every point (i.e. grid-cell) along the stream network. Flow accumulation was used as the mechanism to relate these environmental variables with the stream network (see Domisch et al., 2015a for details). This accumulative feature causes high correlation among many of the variables and with stream flow (Kuemmerlen et al., 2014(Kuemmerlen et al., , 2015a and hence includes an aspect of hydrological information. This dataset therefore contains information describing both climate and hydrology, embedded within the values.

Hydrology
The dataset from Irving et al. (2018), includes 53 of the Indicators of Hydrological Alteration (IHA) metrics that describe the magnitude, frequency, timing, duration and rate of change of flow events (Olden & Poff, 2003). IHA metrics are commonly used in flowecology assessments (e.g. Kakouei et al., 2018) and environmental flow research (Poff et al., 2010;Peres & Cancelliere, 2016) as they comprehensively describe hydrological flow regime. These metrics have only recently been included in river SDMs (but see Irving et al., 2020). As the information contained in the IHA metrics is directly related to the hydrological regime of rivers, it is logical to suggest that this data could affect predicted species' distributions.
In brief, stream flow was extrapolated for the German stream network through a weighted linear regression using accumulated seasonal precipitation from earthenv.org/streams (Domisch et al., 2015a). The daily stream flow (m 3 s -1 ) was then applied as input to calculate the IHA metrics via the R package Eflowstats (www.github.com/USGS-R/EflowStats; Henriksen et al., 2006;Archfield et al., 2014). These data were purposefully created at the same spatiotemporal resolution as the bioclimate and hydroclimate data for explicit use in SDMs. The data have been found to be effective for use in predictive modelling  within the same study area (Irving et al., 2020). All 53 IHA metrics were computed for the time period 1985-2013 to capture enough variability to produce accurate metrics following Kennard et al. (2010) who recommends at least 15 year's worth of data, and informs that there is negligible change in variability over 30 years. In addition, we wanted to ensure the hydrological data covered both the biological (2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013) and bioclimate (1950-2000) data.
These metrics were originally derived from both gauging station daily stream flow data and seasonal precipitation sourced from earthenv.org/streams (Domisch et al., 2015a). It is therefore important to note that some degree of correlation is to be expected with the hydroclimate variables used in this study. Nonetheless, the additional high-resolution stream flow data adequately adds various aspects of flow regime at a high temporal resolution and therefore better represents hydrology.  Table S1). By comparing the explained variance of each predictor set based on variance partitioning analysis, we identified the influence of each predictor set (bC, hC and H) on the spatial distribution of macroinvertebrates. Here, the full model (hC-bC-H) represents the full coverage of environmental predictors used in our study (n = 0) therefore, it is intended for comparative purposes only.

Variable selection process
The variables were selected by adapting the procedure from Irving et al. (2020) using boosted regression trees (BRTs), in a two-step process. First, each predictor set was applied in BRTs separately (hydroclimate n = 19, bioclimate n = 19, hydrology n = 53) for every species within the community (n = 91). BRTs calculate the variable importance of each predictor from the number of times each variable was chosen by the algorithm (Elith et al., 2008). The variable importance was averaged (mean) across all species to find the variable importance for the community. The average variable importance was used to determine the most important individual variables (i.e. 30% of variables with the highest relative importance) from each predictor set (hydroclimate n = 6, bioclimate n = 6, hydrology n = 19). The remaining variables from all predictor sets (n = 31) were then applied collectively into the 2nd run of BRTs with the same criteria as above.
A pair-wise Pearson's correlation analysis was undertaken for each model configuration with the threshold 0.7 (Dormann et al., 2013). The variable importance from the 2nd run of BRTs was used to determine which of the correlated variables were chosen to remain in the analysis. In order to compare the datasets effectively, we retained a minimum of two predictor variables per dataset. Variables chosen for each model are outlined in Tables 1 and 2. As the variables are related, i.e. hydroclimate was derived from the bioclimate dataset, and hydrology was derived, in part, from the hydroclimate dataset, it was likely that a high level of correlation would be observed between datasets (see Table S2 for correlation matrix).
The selected variables (n = 9) were included in our final full model (hC-bC-H). The variables contained in the full model were then distributed according to the remaining model configurations: bC (n = 5), hC-H (n = 7), bC-H (n = 7), hC-bC (n = 7).

Species distribution models
All SDM analyses were undertaken in the R package sdm (Naimi & Araujo, 2016). We applied the four model configurations outlined above to each species within the community in separate SDMs. Each SDM was applied with an ensemble of five algorithms: artificial neural network (ANN), generalized linear model (GLM), flexible discriminant analysis (FDA), BRT, and classification tree analysis (TREE). As the species occurrence data consisted of presences only, we applied 2,000 randomly placed pseudo-absences in geographical space as background absences. Each model was repeated 10 times by bootstrapping. This resulted in 50 models per species, per configuration. For validation, the data were randomly split into training and testing datasets in 70:30 ratios. The true skills statistic (TSS) and the sensitivity values were derived from the validation as a measurement of model performance. Sensitivity is a measure of true positives, i.e. where the model correctly places a presence, and the TSS is derived from both the sensitivity and the specificity (the number of true negatives). The TSS values are reported as weighted mean ± standard error (mean ± SE). In addition, TSS values for individual species were used as a grouping measure for the functional traits analysis (see below).
To predict probability of occurrence for each species and model configuration, an ensemble model was produced by averaging all 50 previously mentioned individual models (5 algorithms 9 10 repetitions) weighted by each model's performance as given by the TSS values. This weighting emphasizes the higher performing models, without excluding lower performing models (Araújo & New, 2007). As output, the ensemble model produced a probability of occurrence map for the entire stream network of the study area. To convert the probability map to binary presence/absence predictions (1, 0), we applied a threshold determined from maximising sensitivity and specificity (Liu et al., 2005(Liu et al., , 2013.

Model performance
All analyses were undertaken in R version 3.5.2 (R Core Team, 2018). Model performance was analyzed using the TSS values calculated through the training and testing validation datasets. Pairwise Wilcoxon tests were applied to the 50 TSS values of each species to test for differences between model configurations. Any values P \ 0.05 were considered significantly different. Each configuration modelled 91 species resulting in a total of 364 Wilcoxon tests. We summarize the outcome as a percentage of significance for each model configuration, i.e. % S = (number of significant models/364) * 100.

Predicted distribution
The predicted distributions were compared according to range size and percentage overlap. Range size was defined as the number of raster cells predicted by the model as a presence within the study area, after converting the predicted probabilities to binary presence/absence predictions by maximizing the sensitivity and specificity produced by the ensemble prediction (Liu et al., 2005(Liu et al., , 2013. Range size was determined by counting the number of species presences predicted by each model configuration. To test for differences in range size between model configurations, pairwise Wilcoxon tests were applied. To compare the predicted distribution of the community, and how they were similar or different in geographical space, pairwise range overlap values were calculated by counting the number of grid cells that contained a same species' presence predicted by each respective model configuration. The proportion of shared grid cells in relation to the predicted range of the pairwise model configurations, was calculated as percentage overlap (mean ± SE) for the community. The predicted distributions of all species per model  (2014), where the proportion of variance explained by the first predictor set can be described as total variance explained minus the proportion of variance explained by the second predictor set in the configuration. The shared variance of both predictor sets can be described as the total variance explained minus the sum of both predictor sets in the configuration. For example, the proportional variance of the hydroclimate data in the hydroclimate and hydrology configuration would be; R 2 hydroclimate = 1 -R 2 hydrology , the proportional variance of hydrology; R 2 hydrology-= 1 -R 2 hydroclimate , and the amount of shared variance: R 2 shared = 1 -(R 2 hydrology ? R 2 hydroclimate ). These proportional variance values were then used as input into the varPart function in modEva package (Barbosa et al., 2016) to calculate the proportional variance partition. This procedure was applied to every species and proportional variance was averaged across all species and reported as community proportional variance.

Analysis of traits
To better understand the influence of the different datasets on the composition of species' communities, we grouped species according to their performance for each configuration based on TSS values from the SDMs. The highest TSS value determined a model configuration group for each species, that is, if 5 species showed higher TSS values for the bC configuration than all other configurations, those species were grouped as ''bC''.
The percentage of trait composition for each group (e.g. bC group, bC-H group) were compared visually. In order to disentangle the influence of each dataset on species traits, we did not include the full model in this analysis. To compare trait composition to predictor variables, we applied pairwise Wilcoxon tests of difference to compare values located at the species occurrences associated with each group.

Variable selection process
The variable selection process resulted in nine predictors: five bioclimate, two hydroclimate and two hydrology variables (Table S3 for BRT coefficients). These were applied in the full model configuration and distributed to each 2-set model (Table 1). Interestingly, variables describing mean temperature of both wettest quarter and of driest quarter from both the hydroclimate and the bioclimate predictor sets were included in the model (see Table 1). It could be expected that the same variable from both climaterelated predictor sets would be highly correlated, however the correlation was negligible: Bioclimate/ Hydroclimate 08 from bC and hC, corr = 0.37, Bioclimate/Hydroclimate 09 from bC and hC, corr = 0.30 (Table S2). Mean values for all chosen predictor variables are outlined in Table 2.

Variance partitioning
We applied variance partitioning to the predicted distributions from the SDMs (presence/absence). The variance partitioning of the 2-set model configurations showed that bioclimate also had the highest explained variance compared with hydroclimate (0.4, Fig. 3a) and hydrology (0.57, Fig. 3b). Hydrology showed the lowest amount of explained variance in all 2-set models, compared with hydroclimate (0.16, Fig. 3c) and bioclimate (0.09, Fig. 3b). Hydroclimate showed lower explained variance than bioclimate (0.12, Fig. 3a) but higher explained variance than hydrology (0.44, Fig. 3c). The shared variance in all 2-set models was relatively low (Fig. 3a-c), the highest being between Bioclimate and Hydroclimate (0.1, Fig. 3a). This demonstrates that the predictor sets have a separate influence on species' distribution.
Model configurations hC-H and bC-H showed a negative shared variance meaning that the two predictor sets explain the variance in different directions, i.e. both a positive and negative relationship. Unexplained variance was lowest in bC-H model configuration (0.36, Fig. 3b), compared with hC-H (0.44, Fig. 3a) and hC-bC (0.38, Fig. 3c).
However, differences were apparent in the geographical location of predicted ranges as shown by the range overlap. The smallest percentage range overlap was between model configurations bC and hC-H (42.3 ± 1.6%, n = 945, Table 4). The largest range overlap was between model configurations bC and hC-bC (73.2 ± 1.3%, n = 1,477, Table 4). The differences in predicted ranges are shown when mapping species richness (Fig. 4) in geographical space. The pairwise Wilcoxon tests on species richness showed a significant difference between hC-H and all other model configurations (P \ 0.001). In addition a significant difference was found between bC and hC-bC-H (P = 0.03) and a near significant difference between bC and bC-H (P = 0.08).

Analysis of traits
We compared the influence of model configurations groups; bC, bC-H and hC-bC on species trait composition. Group hC-H was removed from this analysis as it comprised only three species, hence was not data sufficient to investigate patterns in trait composition. There then remained models configuration groups of 10 species in bC, 44 species in bC-H and 37 species for hC-bC. This analysis derived some noteworthy patterns (Fig. 5). Trait composition for current preference showed rheophil and rheo-limophil (Fig. 5a) comprised species from model group bC (both 30% trait composition) and hC-bC (both 27% trait composition). Species in model group bC-H showed a high trait composition within limo-rheophil category (36%). Trait composition for stream zonation (Fig. 5b) for all model configurations tended to increase from upstream to downstream until a peak was reached, then the values decreased further downstream. Distinct peaks were associated with some models, i.e. bC (20%) and hC-bC (18%) species peaked within lower-trout region, whereas bC-H (20%) peaked in barbel region further downstream. Patterns of trait composition from current preference and stream zonation are comparable in that bC and hC-bC species were associated with higher current velocity and regions further upstream in comparison to bC-H species that were associated with lower current velocity and downstream regions. Trait composition for all models showed a high association for littoral region, i.e. bC (17%), bC-H (26%) and hC-bC (16%). The Wilcoxon tests of difference between species occurrence associated with each model configuration group showed a significant difference between groups bC and bC-H for predictors; Bioclimate 04 (mean ± SE; bC = 628 ± 0.4, bC-H = 619.6 ± 0.3, P = 0.01) and Bioclimate 08 (mean ± SE; bC = 12.3 ± 0.1, bC-H = 14.5 ± 0.05, P = 0.0005).

Discussion
We compared three datasets, combined in five dataset configurations, to evaluate their influence on macroinvertebrate distribution predictions and traits using SDMs. We found that the dataset bC had the most influence on the model in terms of proportional variance, however model performance was increased with the addition of hC or H. After the full model, the configuration bC-H performed best, while hC-H configuration performed least well. The different model configurations predicted similar range sizes, however were not always consistent in geographical space. The composition of species' functional traits, demonstrated distinct patterns between the species with and without the dataset describing hydrology.

Variable selection
One precipitation variable was included in each model configuration, i.e. precipitation seasonality (Bioclimate 15, Tables 1, 2) from the bioclimate dataset. The inclusion of this variable indicates that precipitation- related data are important to include in SDMs as a separate entity to hydrology-related parameters (such as high flow volume), suggesting that local precipitation should not be used as a substitute for hydrologyrelated data (e.g. Domisch et al., 2019). This is supported by the high variance explained individually, as well as the low amount of shared explained variance demonstrated in all model configurations through variance partitioning.

Model performance and variance partitioning
The full model performed best overall. The full model, however, was only used a comparison in our study, therefore from the model configurations containing 1 or 2 datasets the model bC-H performed best. We therefore argue that these two predictor sets (bC and H) complement each other well and that adding hydrology into the model configuration can improve model performance, while keeping the number of predictor variables low. The only model configuration that did not contain bC variables, i.e. hC-H, performed least well overall, confirming the importance of including climate when evaluating species distribution.
Variance partitioning of the individual predictor sets both suggest that bC, i.e. spatially static temperature and precipitation-related variables, have the most influence on the studied macroinvertebrate community, at this scale. Nonetheless, hC and H, i.e. stream and flow-related variables improved model performance, even though they demonstrated a smaller influence in terms of variance partitioning. Our finding is not consistent with current flow-ecology theory that suggests hydrology (H) to be as important as climate (bC) in determining species distribution (Pyne & Poff, 2017).
The scale (13,749 km 2 ) and resolution (1 km 2 grid cells) at which the predictors are applied may partially explain these results (Pearson & Dawson, 2003;Randin et al., 2009;Lenoir et al., 2013;Domisch et al., 2015b;Record et al., 2018) as it may be too coarse to fully capture the influence of all the dimensions of the hydrological regime on macroinvertebrate distribution. The IHA metric MH21; high flow volume is the discharge of the highest flow, represents the more extreme aspects of hydrological regime. High flow volume is likely related to catchment size, i.e. peak flows driven by climate and geology, at large-scales (i.e. catchment scale) (Poff,Fig. 3 Proportional variance partitioning of all four models; a hC-H; hydroclimate and hydrology, b bC-H; bioclimate and hydrology, and c hC-bC; hydroclimate and bioclimate 1997). At small scales, i.e. reach scale, large-scale variables describing the hydrological regime, such as flood and droughts, are able to induce changes in river biota communities by influencing small-scale habitat, e.g. hydraulics, riverbed substrates and stream channel morphology (Allen & Vaughn, 2010;Soranno et al., 2010). However, local-scale hydraulics, e.g. pool/ riffles, and their resultant impact on physical microhabitats, e.g. creation of refugia, influence the distribution of biota, which reduces the influence of largescale drivers (Poff, 1997;Frissell et al., 1986). Of the few freshwater SDM studies that include variables directly describing hydrology, together with climate, most show a high importance of either both climate  Bold is mean percentage overlap in geographical space across all species (n = 91) ± standard error. Italics is absolute averaged number of overlapping predicted presences and hydrology, or hydrology as the most importance variable (Bond et al., 2011;Kuemmerlen et al., 2014;Huang et al., 2016), although a lower influence of hydrology has been documented in stream headwaters (Monbertrand et al., 2019). Hydrology is certainly important at finer spatial scales (Kuemmerlen et al., 2014), and climate becomes more important with increasing scale (Friedrichs-Manthey et al., 2020). This scale-dependency is a recognised challenge in SDM research (Domisch et al., 2015b;Friedrichs-Manthey et al., 2020). Nonetheless, hydrology metric RA3; Fall Rate, describing the rate of negative changes in flow between consecutive days, can drive potential drought stress. This negatively affects different communities in the ecosystem such as local plant communities (The Nature Conservancy, 2009), and initiates cascading effects like resource scarcity. Therefore, even on large-scales, factors impacting local-scale changes may still be influential.
In addition, the disproportionate number of variables could drive the differences in explained variance, that is, five variables were applied from the Bioclimate dataset, whereas two variables were applied from each of Hydroclimate and Hydrology datasets. Although we ensured at least two variables from each dataset for a robust comparison, the remaining variable selection was chosen by the selection process. As the selection process had a high influence in determining the most important variables for this community and study area, we believe this approach is less biased than forcing the exact same number of variables from each dataset.
Despite these challenges, incorporating hydrology at the scale of this study does not hinder predictive ability (Araújo et al., 2019) and resulted in an improvement in SDM performance compared with climate alone. Applying the same predictors at a smaller spatial (e.g. reach scale) with finer resolution (e.g. \ 100 m 2 : Kuemmerlen et al., 2014) may result in hydrology demonstrating a stronger influence (Friedrichs-Manthey et al., 2020), on macroinvertebrate distribution, however further investigation would be needed for confirmation.

Predicted distributions
The model configuration bC-H predicted the largest range size overall. Conversely, hC-bC predicted the smallest range size. However, we only found a significant difference in range size between model configurations hC-H (the model that performed least well) and hC-bC.
Despite similar range size, some model configurations predicted species' distributions in different geographical locations as suggested by the range overlap and differences in species richness between model configurations. The ranges predicted by model configuration hC-H overlap between 42 and 57% with all other models and were significantly different to all other configurations in terms of species richness. The correlation analysis and variance partitioning suggests that the datasets contain distinct environmental information, rendering the differences in predicted range size and location plausible. Here, the SDMs are relating different environmental conditions to the species' known occurrences and predicting suitable habitat accordingly. Nonetheless, as hC-H performed least well, and does not include the most important dataset (bC) it is reasonable to suggest that these predicted ranges do not contain enough relevant information for accurate predictions.
For other model configurations 67-77% range of each model overlapped. We found a significant difference in species richness between our most basic model (bC) and the full model (hC-bC-H). This difference may be driven by the addition of hydrologyrelated variables (i.e. hC and H datasets), or it may be driven by the full model containing the most information, hence higher complexity (Araújo et al., 2019). As the full model, in this study, was evaluated for comparison purposes only, it is challenging to disentangle the reason behind this result. Considering the near significance of the differences between bC and bC-H it is possible that the addition of pure hydrology variables (H) influences species' predicted range, however we were unable to confirm this conclusion with our analysis.

Analysis of traits
The models that contained specific information on hydrology (bC-H) performed best when applied on species that prefer intermediate to slow-flowing current conditions in regions that were further downstream. Species modelled with more information regarding climate (bC, hC-bC) showed preference for high to intermediate flow conditions located further upstream. Model configurations hC-bC and bC-H were comparable in terms of model performance, explained variance and predicted ranges. However, the species that performed better with each configuration, demonstrate distinct differences in their trait composition. These findings agree with Irving et al. (2020) where species that were modelled with additional descriptive hydrology variables prefer relatively slower current conditions and downstream regions.
Variables describing climate (mean temperature of wettest quarter and temperature seasonality) varied significantly for species occurrences grouped in bC and bC-H models, with higher temperature and lower seasonality shown for bC-H. These differences suggest that the variability in climate is strong enough upstream to drive the distribution of certain macroinvertebrate species, however in downstream regions, hydrology plays a much larger role. It is possible that the species located in upstream regions, at least in relative terms, are better adapted to colder temperatures than species located downstream (Monbertrand et al., 2019), and hence may respond more strongly to climate than hydrology.
Flow regime in both the Ems and Weser catchments is generally driven by precipitation (i.e. pluvial) with seasonal high flows driven by snowmelt and precipitation (Koeniger et al., 2009). In the Weser catchment, increasing low flows as well as decreasing low flow duration have been documented over the last few decades (Bormann and Pinter, 2017). Although increased precipitation due to climate change may have influence (Hisdal et al., 2001), increased flows are likely attributed to reservoir management (Bormann & Pinter, 2017). The climate datasets (hC and bC) applied in this study may therefore be unable to capture changes in flow due to anthrogenic influences. In contrast, the hydrology dataset is derived, among other variables, from discharge gauging stations, which measure real-time stream flow (m 3 s -1 ), and therefore account, to some extent, hydrological variability not driven by precipitation, e.g. anthropogenic alteration, urban run-off and groundwater influence. The hydrologic variables (and the anthropogenic disturbance to hydrological regimes) are likely to assist in better depicting the varied influences on the physical habitat of the river biota (Resh et al., 1988;Poff et al., 1997), which shapes the structure and function of river ecosystems. In contrast, the climate datasets (bC and hC) are applied as indirect surrogates for water temperature (Moore et al., 2005;Caissie, 2006) and hydrology (e.g. Domisch et al., 2019).
It is therefore plausible that the species with upstream preferences, may be located in areas with limited influence of reservoir management (i.e. upstream of reservoirs), and hence performed better with climate-related variables. In contrast, species downstream of reservoirs are more influenced by direct measures of hydrology. However, more investigation is needed to confirm this theory.
Role of hydrology in SDMs: implication for practice Our findings have important implications when applying such models to inform conservation efforts: omitting flow regime variables in SDMs may lead to an underrepresentation of macroinvertebrate species that are sensitive to changes in flow. For example, SDMs applied under current, and future climate conditions have predicted range shifts of e.g. macroinvertebrate distribution to higher altitudes (Domisch et al. 2011(Domisch et al. , 2013. Adding complementary factors describing flow regime may result in more accurate predicted range, both for current and future potential distributions, subsequently requiring adjustments in mitigation strategies for conservation efforts. Many rivers in Germany (and globally) are highly altered from their natural regime. Water as a resource is in high demand continuously abstracted, recycled and discharged back into the system, or rivers serve as heavily managed waterways. By focusing on climate or climate-related surrogates, predicting species ranges may not capture the full diversity of functional traits, and hence may ignore important relationships between flow-sensitive species and their environment. This factor is especially important when predicting species' distributions under future conditions, as it can provide insights to the ultimate drivers of range shifts (Monbertrand et al., 2019;Irving et al., 2020), and also should be considered when applying SDMs on other aquatic species (Ihlow et al., 2011;Markovic et al., 2014;Ruiz-Navarro et al., 2016;Rodriguez-Merino et al., 2019).

Concluding remarks
The main findings of our study suggest that by including environmental predictors describing flow regime in SDMs applied on macroinvertebrates can potentially increase model performance, despite a low contribution of hydrology to explained variance. The IHA metrics applied in this study are partially derived from real-time stream flow data from gauging stations, which incorporate the principal factors that control river hydrology. The metrics describe direct influencing factors of river habitat, including some aspects of anthropogenic disturbance. These characteristics are not described by either of the climate-related predictors included in our study. In addition, our analysis of functional traits provides insights into the driving factors of macroinvertebrate distribution. Improvement in model performance is species-specific and likely linked to the relevance of hydrology to particular species.
Our study intentionally focused on climatic and hydrologic variables however, we acknowledge that to provide a holistic assessment of benthic macroinvertebrate distribution, additional elements known to influence river systems need to be considered. It is therefore possible that the inclusion of data describing the surrounding land use (that influences stream flow and nutrient influx), topography (that influences river hydraulics) and connectivity (that influences dispersal and habitat availability) of the study area could increase predictive accuracy and influence range size.
Our study highlights improvements to predictive ability, and provides insights into drivers of community composition through functional traits. We recommend that given its fundamental importance, variables describing flow regime must be considered in SDM studies applied on river biota.
Author contributions All authors contributed to the study conception and design. Material preparation, data collection and analysis were performed by KI. The first draft of the manuscript was written by KI and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript.
Funding Open Access funding enabled and organized by Projekt DEAL. This work has been supported by the German Federal Ministry of Education and Research (BMBF; 01LN1320A and 033W034A).
Data availability Environmental datasets are available open access. Biological data were collected from German Federal State agencies and not publicly available.
Code availability Not applicable.

Declarations
Conflict of interest All authors declare that they have no conflict of interest or financial ties to disclose.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.