Spatial modelling of temporal dynamics in stream fish communities under anthropogenic change

Understanding temporal changes in aquatic communities is essential to address the freshwater biodiversity crisis. In particular, it is important to understand the patterns and drivers of spatial variation in local community dynamics, generalizing temporal trends from discrete locations to entire landscapes that are the main focus of management. Here, we present a framework for producing spatially continuous views of community dynamics, focusing on stream fish affected by hydropower development.


| INTRODUC TI ON
In freshwater ecosystems, biological communities are changing worldwide due to multiple anthropogenic pressures (Albert et al., 2020;Reid et al., 2019). These changes involve for instance species extinctions, defaunation and taxonomic homogenization, which tend to occur faster at smaller spatial scales, but then scale up to entire watersheds, regions and even continents (Magalhães et al., 2007;Matthews & Marsh-Matthews, 2016, 2017Villéger et al., 2011;Zbinden, 2020). Therefore, much effort has been devoted to understanding where, how and why local freshwater biological communities change over time, usually through studies conducted for extended periods at a number of discrete locations (e.g. Baranov et al., 2020;Bêche et al., 2009;Erős et al., 2020;Jourdan et al., 2018;Matthews & Marsh-Matthews, 2016). However, relatively few studies have investigated the patterns and drivers of spatial variation in the temporal dynamics of local community, though this would be important to generalize patterns from particular sites to landscapes and regions that are the main focus of management (Erős & Lowe, 2019;Hugueny et al., 2010;Schlosser, 1991;Wiley et al., 1997). Therefore, it is essential to provide conservation and water agencies with spatially continuous views of community dynamics, thereby contributing to assess anthropogenic impacts and to prioritize management action Fausch et al., 2002).
Changes in local biological communities may result from extinctions or colonizations affecting species richness and composition, and from variations in species abundances (Grossman et al., 1990;Magalhães et al., 2007;Matthews & Marsh-Matthews, 2017).
Moreover, changes may follow distinct temporal patterns, being for instance gradual or saltatory, or reflecting variations around loose equilibria, shifts between alternative stable states or gradual directional transitions away from initial community structures (Collins, 2000;DeAngelis et al., 1985;Matthews & Marsh-Matthews, 2017). These changes can be quantified using relatively simple metrics, such as Kendall's coefficient of concordance to estimate constancy in species rank abundances or the coefficient of variation to estimate variability in species abundances (Grossman et al., 1990). However, these measures do not reveal patterns in temporal change, which have often been inferred by visual examination of trajectories in a chosen space of community resemblance (Magalhães et al., 2007;Matthews et al., 2013). Moreover, the temporal convergence/divergence of trajectories at pairs of sites can be used to quantify whether communities are varying in synchrony or converging/diverging from each other (De Cáceres et al., 2019). All these metrics can be used to model spatial patterns in community dynamics.
To generalize a community dynamics metric obtained at discrete locations to a spatially continuous surface, it is necessary to find variables that (a) are correlated with variation in that metric and (b) can be easily mapped at the landscape scale. This may be difficult when community dynamics reflect mainly idiosyncratic variations in local conditions (Erős & Lowe, 2019;Matthews & Marsh-Matthews, 2017), which are hard to extrapolate at larger spatial scales. However, generalization is possible when community dynamics are affected by large-scale gradients, such as for instance the upstream-downstream gradients in rivers or gradients related to sources of human disturbance Gorman & Karr, 1978;Milardi et al., 2019;Schlosser, 1987). Moreover, local dynamics can be influenced by spatial connectivity associated for instance with the topology of stream networks, anthropogenic barriers and habitat fragmentation Crabot et al., 2020;Erős & Lowe, 2019;Hugueny et al., 2010), as it affects meta-community mass effects mediated by dispersal (Heino et al., 2015;Tonkin et al., 2018), as well as the spread of invasive species (Filipe et al., 2017;Gavioli et al., 2019;Milardi et al., 2019;Mota-Ferreira & Beja, 2020). Therefore, spatial modelling of community dynamics requires establishing relations with environmental variables predicting variation in dynamics metrics, and accounting for spatial variables reflecting the effects of connectivity. In the case of rivers, recently developed geostatistical models provide a convenient framework to undertake such modelling exercise, as they account for the complex topology of spatial relations in dendritic networks (Peterson et al., 2013), integrating (Euclidean) spatial dependencies that occur overland, as well as (hydrologic) spatial dependencies along the river network and the effects of flow connection .
Here, we combine community trajectory analysis (De Cáceres et al., 2019) and geostatistical modelling (Peterson et al., 2013) to understand and map spatial patterns of community temporal dynamics in dendritic stream networks. We focused on stream fish communities (sensu Matthews & Marsh-Matthews, 2017) in a watershed where a hydroelectric development was built and started operating during the study (Jackson, 2011;Santos et al., 2017). We expected communities to be more unstable and eventually undergoing directional changes in lotic reaches close to hydroelectric reservoirs, mainly due to the spread of exotic species (Santos et al., 2017). To test this idea, we used data from stream fish monitoring carried out at 30 sites, encompassing from the construction phase (2012-2014), through the filling of the reservoir (2014-15), to the operation phase K E Y W O R D S community dynamics, community trajectory analysis, dendritic networks, exotic species, geostatistics, hydropower generation, stream ecology , and we (a) described community variation in terms of species composition, richness and abundances; (b) quantified dynamics in terms of the velocity and directionality of community change; (c) modelled community trajectory metrics in relation to large-scale ecological gradients (e.g. stream order, elevation), interannual variation in local environmental conditions (e.g. water flow and depth) and the prevalence of exotic species; (d) investigated spatial patterns in community change in relation to environmental and spatial factors; and (e) developed predictive geostatistical models accounting for large-scale ecological gradients and spatial dependencies to produce spatially continuous maps of community dynamics.

| Study area
The study was conducted in NE Portugal, in the River Sabor watershed (N41°09′-42°00′, W7°15′-6°15′; Figure 1), encompassing a wide range of variation in elevation (100-1,500 m above sea level), annual precipitation (443-1,163 mm) and mean annual temperature (6.9-15.6°C). Climate is Mediterranean, with precipitation largely concentrated in October-March and virtually none in hot summer months (June-August). Flow regime is highly seasonal, with most headwater streams drying out or being reduced to pools in summer, while the main watercourse and the largest tributaries are permanent. Two hydroelectric dams (Feiticeiro: 181 ha; Baixo Sabor: 2 820 ha) located near the mouth of the River Sabor started to be built in 2009, with the main reservoir filling in autumn/winter of 2014/2015 (Jackson, 2011;Santos et al., 2017). The Sabor watershed and its fish communities are more thoroughly described by Ferreira et al. (2016) and Santos et al. (2017).

| Fish sampling
As part of a preliminary fish survey (Ferreira et al., 2016) (Ferreira et al., 2016). A relatively small sampling reach was chosen because (a) we were interested in investigating how local community dynamics varied across the watershed, (b) communities are more dynamic at finer spatial grains (Zbinden, 2020) and (c) previous studies demonstrated this reach length to be adequate for capturing responses of Mediterranean stream fish communities to environmental fluctuations (Magalhães et al., 2007). Sites were sampled annually from 2012 to 2019 in June-July (Table S1), when reduced water flows favoured sampling efficiency, but before the peak summer drought when harsh conditions might cause high fish mortality. Each reach was electrofished by the same operator (MMF) using procedures detailed in Ferreira et al. (2016), with consistent effort and methods at each site over the years. Reaches were electrofished for 15-25 min, with longer surveys in wider and deeper streams to enhance detectability of all species in the local fish communities (Ferreira et al., 2016). Fish were identified to species level, measured for total length and returned alive to the stream. Sampling was conducted under licence from the Instituto da Conservação da Natureza e Florestas, which required individuals of exotic species to be euthanized. and annual precipitation [Prec]), which have already been used in predictive species distribution modelling in the Sabor watershed (Ferreira et al., 2016;Filipe et al., 2017;Mota-Ferreira & Beja, 2020;Quaglietta et al., 2018). These variables were extracted using the CCM 2.1 database (Vogt et al., 2007), following procedures detailed in Table S2. We only used precipitation extracted from worldclim (Hijmans et al., 2005), because climatic variables tend to be highly in- water velocity [Vel_cv], which were estimated in the field following procedures detailed in Table S2. We used CV rather than mean values, because we wanted to evaluate how community dynamics were affected by interannual variation in local conditions, and because mean values were correlated with landscape variables such as stream order. To account for the possibility of biological invasions increasing community variability (Erős et al., 2020), we computed the proportion of exotic fish species at each site [Exot] (Table S2).

| Environmental and spatial variables
Spatial data necessary to account for spatial autocorrelation (see below) were obtained in a GIS using the stream network extracted from CCM2.1 and the layer of sampling sites. Estimates included the Euclidean and hydrologic distances (total and downstream hydrologic distances) between each pair of sites . To deal with confluences in tail-up models (see below), we also estimated catchment areas to weight the relative influence of the branching upstream segments (e.g. . Spatial estimates were made using the package "riverdist" (Tyers, 2017) in R and the Spatial Tools for the Analysis of River Systems (STARS) toolbox version 2.0.0  for ArcGiS 10.2 (ESRI, 2016). were very similar across analyses, so we present only the results based on fish abundances, without size structure. We used counts of each species at each site and year instead of species densities, because the later can be driven by changes in fish numbers, habitat area or both (Magalhães et al., 2007). In all analysis, we excluded fish < 5 cm, often corresponding to young of the year (yoy), because they were poorly sampled and their recruitment to the fishing gear might be strongly influenced by the timing of sampling in relation to fish spawning. trajectories. The mean velocity measured how fast the community changed during the study period, while directionality was used to assess whether the community changed over time following a directional pattern or otherwise showed cyclic or random patterns. The pairwise dissimilarity in trajectories was used to understand whether spatial patterns in community temporal variation were associated with spatial environmental patterns and spatial autocorrelation.

| Community trajectory analysis and modelling
To quantify the environmental factors driving temporal community change, we first used generalized linear models (GLMs) with Gaussian errors and identity link, to relate the mean velocity and directionality of trajectories to environmental variables. Pairwise scatterplots were visualized to check for potential outliers and influential points. We found that the single site of order 6 was a potential influential point, and so, it was combined with order 5 sites in a single category. We then screened the pairwise relationships between dependent and predictor variables, considering both linear and nonlinear relations using orthogonal polynomials of second degree. More complex relations (i.e., higher order polynomials) were not considered because of relatively small sample sizes. In subsequent multivariate model building, we considered for each predictor either the linear or polynomial terms that provided the best fit to the data, judged considering the adjusted R-squared and the Akaike information criteria corrected for small sample sizes (AICc). We then screened all combinations of predictors for each dependent variable and, in each case, we retained as best model the combination of predictors minimizing AICc (Murtaugh, 2009).
We also modelled variation in pairwise distances between community trajectories as a function of environmental and spatial distances between sites. First, we computed environmental pairwise distances considering all standardized environmental variables, which were summarized using the symmetric multidimensional scaling implemented in "smacof" (De Leeuw & Mair, 2009). Then, we computed the Mantel correlograms to assess the scale of spatial dependencies in community trajectories and environmental conditions (Legendre et al., 2015), using "vegan" (Oksanen et al., 2012).
Finally, we computed a multiple linear regression on distance matrices (MRM; Lichstein, 2007) with "ecodist" (Goslee & Urban, 2007) and 100,000 permutations, using combinations of environmental and spatial distances between sites, and retaining the model with the largest R 2 . MRM was used, despite the pitfalls of Mantel-based approaches, because these are still considered adequate to analyse dissimilarity matrices, and interpretation was made considering potential problems such as inflated type I errors (Legendre et al., 2015).
All analyses were performed using R software (R Core Team, 2019).

| Geostatistical modelling and mapping
We used geostatistical modelling to relate variables describing community dynamics to both environmental and spatial predictors, considering the spatial structure of dendritic stream networks Peterson et al., 2013;Ver Hoef et al., 2006). These geostatistical models are similar to conventional linear mixed models, with specification in random errors of spatial dependencies as functions of either straight-line distances (Euclidean model) between sites, hydrologic distances between sites connected by the water flow (tail-up model) or hydrologic distances irrespective of water flow connection (tail-down model). The fixed component corresponded to the best linear models (GLMs) developed in previous analysis for the mean velocity and directionality.
The random component was specified considering the full autocovariance structure, which provides the greatest flexibility for representing multiple types of autocorrelation simultaneously . To select the best autocovariance function for each spatial component, we tested all combinations of functions for the models including the three spatial components and selected the one minimizing AICc.
To map spatial variation in community dynamics, we projected the mean velocity and directionality of community change predicted from the geostatistical models on the stream network of the entire Sabor watershed. First, we divided the stream network into segments of a maximum length of 1,000 metres using ArcGIS desktop (ESRI, 2016), and we extracted the value of environmental variables from the centroid of each segment. We then predicted the values of the metrics in each segment using universal kriging within the "SSN" package .

| Changes in species richness and fish catches
The annual number of species recorded in the Sabor watershed remained essentially constant at 10-12 species, but the number of individuals captured (i.e., fish catches) varied widely (Table 1, Figure 2).
Variability over time, as measured by the coefficient of variation, was highest (CV > 1) for rare species ( non-significant for all the other (−.577 < r < .579, p > .10).

| Community trajectories
The fish community trajectories represented in the PCoA biplot indicated major variations across sites, albeit without obvious temporal patterns (Figure 3). There was no visual evidence for directional changes, with communities deviating and later returning to previous states. This was supported by the small directionality values (0.33 ± 0.04, which varied little (0.25-0.43) across the watershed. In univariate analysis for directionality, the best relations were mostly linear, but there were quadratic relations for elevation, stream order and the CV of water depth (Table S3). However, all relations were weak and statistically non-significant. The best multivariate model (Table S4), included a weak U-shaped effect of elevation, suggesting slightly higher directionality at the lowest and highest elevations, and showed a slightly increase in directionality along with slope (Table 2, Figure 4a).
The velocity of community change at each site varied between 0.29 and 0.89 (0.54 ± 0.15). In univariate analysis, the best relations were always linear, except for the quadratic relation with stream order (Table S3). Some univariate effects were statistically significant, with mean velocity declining linearly with elevation, and showing quadratic relations with stream order, with faster changes in 3rd orders and smaller in 2nd-order, 4th-order and particularly 5th/6th-order streams, and with the proportion of exotics, with increases up to about 0.8 and levelling off or slightly declining thereafter. The best multivariate model accounted for 61% of variation, showing faster changes at lower elevation irrespective of stream order, while at any given elevation the velocity of change increased from 5th to 3rd orders, while declining again slightly in 2nd-order streams (Table 2, Figure 4b; Table S5).
The Mantel correlograms showed that pairwise distances between community trajectories were significantly related to Euclidean distances up to about 15 km and to hydrologic distances up to about 30 km ( Figure S1). Environmental and spatial distances were also related to each other, albeit weakly ( Figure S2). The best MRM model (F = 137.018, p < .001) accounted for ~39% of variation in community trajectories between sites, underlining significant effects of environmental (coefficient = 1.6 × 10 −6 ; p < .001) and hydrologic distances (coefficient = 0.051, p < .001).

| Spatial variation in community dynamics
The geostatistical model explained about 40% of variation in community directionality (Table 3) Table 1) small distances. The spatial projection of model predictions produced a map showing that directionality was always low, with minor spatial variation across the watershed (Figure 5a).
Regarding the mean velocity of community change, the geostatistical model explained close to 30% of its variation across the watershed (Table 3)

| D ISCUSS I ON
Our study shows that combining robust descriptors of community change with state-of-the-art geostatistical modelling contributes to understanding and predicting where, how and why the dynamics of local communities vary across stream networks. We found that local fish communities varied widely over the years, but there was no evidence for directional changes, pointing out a state of loose equilibrium (sensu Matthews et al., 2013) across the watershed. However, there was much spatial variation in the velocity of community changes, which were strongly influenced by environmental gradients associated with elevation and stream order.
The spatial patterns of local community dynamics appeared to be affected also by stream network topology, given the strong influence of hydrologic spatial dependencies. Mapping of community dynamics highlighted faster changes in lowland streams affected by hydroelectric development and exotic species. Overall, our framework helps to generalize community dynamics from discrete locations to entire watersheds, providing spatial information needed for freshwater ecosystem assessment and management Fausch et al., 2002).

| Temporal drivers of local community change
As in other Mediterranean-type streams (Bêche et al., 2009;Magalhães et al., 2007), fish communities in the Sabor watershed were highly dynamic, with temporal changes involving mainly fluctuations in species abundances, while much less variation was found in species composition and richness. The high dynamism observed was probably affected by the scale at which the study was conducted, as communities tend to show far more marked changes at the local than at the watershed or regional scales (Magalhães et al., 2007;Zbinden, 2020).
However, community dynamics were probably also driven by strong environmental fluctuations during the study period, particularly the occurrence of extreme droughts that strongly affect stream fish survival and recruitment (Lennox et al., 2019;Magalhães et al., 2003Magalhães et al., , 2007Matthews & Marsh-Matthews, 2017). fish communities take time to recover after environmental extremes (Bêche et al., 2009;Matthews et al., 2013;Resh et al., 2013), which may have contributed further to the variability observed. The spread of exotic species may also have affected community dynamics, due to temporal changes in their own prevalence and abundance, but also due to eventual negative effects on native species (Erős et al., 2020;Gavioli et al., 2019;Milardi et al., 2019;Zanden et al., 2015). This is supported by faster changes in sites with higher proportion of exotic species, and by the high coefficients of variation in the abundance of exotic species when compared to most native species. Finally, random sampling variation may have contributed to the fluctuations observed, though this was minimized through preliminary testing and optimization of sampling methodologies (Ferreira et al., 2016), and by having the same operator applying a constant sampling effort at every site in every year. Therefore, we are confident that the temporal patterns observed are unlikely to result from methodological artefacts. Overall, the local community dynamics recorded in our study seem to be comparable to that reported elsewhere for Mediterranean-type stream fish (Bêche et al., 2009;Magalhães et al., 2007), but also other aquatic organisms (Bêche et al., 2009;Crabot et al., 2020), suggesting that the patterns observed may apply to other study systems.

| Patterns and environmental drivers of spatial variation in community dynamics
Despite the temporal changes observed, there was no evidence for directional community dynamics at any sampling site, suggesting that local fish communities were in a loose equilibrium (Matthews & Marsh-Matthews, 2016;Matthews et al., 2013). This view was supported by consistently small directionality values estimated TA B L E 2 Summary results of the best AICc models explaining variation in mean velocity of community change, directionality and pairwise distances between community trajectories Note: For the intercept and each variable in each model, we provide the regression coefficient, the standard error of the coefficient estimate and the corresponding t-and p-values.

F I G U R E 4
Trend lines (± standard errors) describing the relations inferred from models (Table 2) relating the directionality of community change to elevation for three levels of maximum slope (percentiles 10%, 50% and 90%) (a) and the mean velocity of community change to elevation for each stream order (b)    (Magalhães et al., 2007), and are in line with the view that communities often exist in a dynamic temporal equilibrium (Collins, 2000;DeAngelis et al., 1985). However, the patterns may be considered unexpected given the major anthropogenic changes affecting the Sabor watershed, including the building of a large hydroelectric infrastructure and the increasing spread of exotic species, including fish (Santos et al., 2017) and crayfish (Filipe et al., 2017).
This may be because the study period was too short to capture community trends, as most exotic species were already present in the watershed at the start of our study (Ferreira et al., 2016), and probably not enough time elapsed for the new dams causing major shifts in species composition in nearby lotic areas. This is supported by observations regarding three exotic species potentially spreading from reservoirs (Ribeiro & Veríssimo, 2014;Vinyoles et al., 2007), one of which has been steadily increasing in abundance since 2015 (A. alburnus), while two others only started to be detected towards the end of the study period (M. salmoides and R. rutilus). Directional changes may thus happen in the future, most likely associated with increasing prevalence and abundance of exotic species Milardi et al., 2019).
The velocity of local community changes varied widely across the watershed. This spatial variation was strongly related to structural landscape features such as stream order and elevation, while factors associated with interannual variability in local environmental conditions did not show measurable influences. The effect of stream order might be expected, as it reflects strong longitudinal gradients along rivers in for instance water discharge, and habitat size and heterogeneity (Hughes et al., 2011), which are strong drivers of ecological processes (Vannote et al., 1980) and the distribution of organisms (Harrel et al., 1967;Paller, 1994;Platts, 1979), including in the study area (Ferreira et al., 2016;Filipe et al., 2017;Mota-Ferreira & Beja, 2020;Quaglietta et al., 2018). We found that changes were fast in 2nd-order streams, still a little faster in 3rd-order streams, and then velocity declined in larger order streams. This is in line with studies suggesting higher temporal changes in fish communities in headwaters than further downstream (Schlosser, 1987), though other studies suggest that longitudinal gradients in stream fish community dynamics may vary across watersheds depending on local environmental conditions (Matthews & Marsh-Matthews, 2017). The later was supported by our study, because the observed joint effect of stream order and elevation implied that there were communities in lower order streams at high elevation that varied slower than those in higher order streams at low elevation. Reasons for these patterns are not completely clear, but it is noteworthy that lower orders at high elevation correspond to cold-water mountain streams with species-poor communities dominated by S. trutta, which may be relatively stable over time. In contrast, lower orders at low elevation generally correspond to warm water streams with richer communities dominated by cyprinids and exotics, which during the dry summer months are often reduced to a series of disconnected pools, and thus where fish communities may vary widely from year to year in association with droughts and floods (Bêche et al., 2009;Magalhães et al., 2007). Other possibility is that elevation acted as a surrogate for increasing human disturbance in the lowlands driving higher variability in fish communities, which may be mediated by the increasing prevalence and abundance of exotic species (Erős et al., 2020;Gavioli et al., 2019;Gorman & Karr, 1978;Milardi et al., 2019). This idea is supported by the positive relation observed between the proportion of exotic fish and community variability, and by the inverse relation between the prevalence of exotic crayfish and elevation also found in the watershed (Filipe et al., 2017). However, in multivariate models the effect of elevation was retained but not that of exotic species, possibly because the former may capture the effect of the later, and account in addition for unmeasured ecological processes driving community dynamics. Overall, the velocity of community change appeared to be mainly associated with large-scale proxies, possibly reflecting spatial gradients in more local ecological processes such as biological invasions, which would require further clarification.

| Spatial dependencies in community temporal dynamics
The patterns of community temporal change were also related to spatial dependencies, with more similar community trajectories in sites closer to each other, either overland (Euclidean) or along the waterlines (hydrologic). The effects of Euclidean distances were significant up to about 15 km and possibly reflected similarity between sites in environmental conditions influencing community trajectories. For instance, sites close to each other are likely to be more similar than those farther apart in environmental conditions driven for instance by elevation, which was related to the velocity and, to a lesser degree, the directionality of community change. Hydrologic spatial dependencies were significant up to about 30 km and may be a consequence of similarities between sites associated for instance to stream order and elevation, and to unmeasured spatially structured environmental factors (Legendre & Legendre, 2012). In addition, however, hydrologic spatial dependencies were probably also influenced by mass effects (Heino et al., 2015), with fish dispersal among neighbouring locations homogenizing species composition and synchronizing population fluctuations (Erős & Lowe, 2019;Hugueny et al., 2010;Tonkin et al., 2018). For instance, fish dispersal from larger streams to headwaters may contribute to reduce community fluctuations in the later (Matthews & Marsh-Matthews, 2017). Also, dispersal of exotic fish species across the watershed may contribute to biotic homogenization and to similarities in community fluctuations in sites nearby Milardi et al., 2019).
The geostatistical models further supported the importance of spatial effects and allowed a finer examination of their contribution to community dynamics. In the case of directionality, there were only marked Euclidean effects, with a long range, suggesting influences driven by large-scale spatial gradients overland. This effect should be interpreted with care, given the low values of directionality and its low variability across the watershed. Regarding the velocity of community change, spatial dependencies also explained a large share of variability, with a strong contribution of the tail-down model and a minor contribution of the tail-up model. This is in line with the idea that the tail-down model captures spatial dependencies associated with organisms that can move actively both up and downstream, while the tail-up model mainly reflects spatial dependencies resulting from the passive drift of materials or organisms downstream (Peterson et al., 2013). This, together with the observation that tail-down effects were strong even after accounting for the variable reflecting spatially structured variation in environmental conditions along the stream network (i.e. stream order), further suggest that spatial variation in the velocity of community change was affected by movement of individuals along the waterlines.
Overall, our results support the idea that fish community dynamics is strongly affected by spatial dependencies and the topology of the stream network (Erős & Lowe, 2019;Hugueny et al., 2010;Tonkin et al., 2018).

| Mapping community temporal dynamics to guide management
Mapping of community dynamics highlighted areas across the watershed where larger changes seem to be occurring, some of which may be associated to anthropogenic pressures. In the case of directionality, mapping showed little variation across the watershed, suggesting that at least until now the construction and operation of the Baixo Sabor Hydroelectric Infrastructure did not disrupt the loose equilibrium of fish communities, as there was no evidence for streams closer to the reservoirs showing more directional changes than streams elsewhere in the watershed. It should be noted, however, that our study only encompassed four years after the filling of the larger reservoir, and so, it cannot be ruled out that directional changes will become apparent in the long term. The velocity of community change varied across the watershed, with some evidence for faster changes occurring in streams draining into the reservoirs, and in the Sabor River immediately upstream of the reservoir. This suggests that the presence of large reservoirs may be increasing fish community instability in surrounding lotic environments, through for instance the spread of exotic species (Santos et al., 2017). Notwithstanding, the highest velocity of community change was found in a small watershed (Vilariça) that does not drain into the Baixo Sabor reservoirs. This watershed is affected by a number of anthropogenic pressures, draining into another large dam downstream of Baixo Sabor, flowing through an area of intensive agriculture and being subject to habitat management interventions (Boavida et al., 2018), all of which may have contributed to fast community changes. Overall, the spatially continuous mapping of temporal community dynamics provided a visual representation of the type and spatial extent of anthropogenic impacts on stream fish communities, which would have been more difficult to perceive otherwise.

| CON CLUS IONS
Freshwater biological communities are rapidly changing worldwide due to direct and indirect anthropogenic pressures, making it critical to understand how, where and why such changes are occurring (Albert et al., 2020;Reid et al., 2019). Our study combining community trajectory analysis (De Cáceres et al., 2019) and geostatistical modelling (Peterson et al., 2013) contributes to address these issues, by offering a relatively simple and flexible framework to spatially generalize data on community dynamics collected at discrete sampling locations. Using this framework, we were able to show that local dynamics were affected by larger scale processes operating within the stream network, including both environmental gradients and spatial processes mediated by network topology (Erős & Lowe, 2019). Moreover, we produced maps that helped visualizing community changes across the stream network and that highlighted the effects of a new hydroelectric development in nearby lotic systems. We suggest that our framework may be widely useful to freshwater ecologists aiming to understand spatial variation in local community dynamics under anthropogenic change, while providing a tool for managers to make spatially continuous predictions of community temporal dynamics that can be used in bioassessment and mitigation of anthropogenic impacts on freshwater ecosystems Fausch et al., 2002).