DOI: 10.5061/dryad.prr4xgxj5
title: Data from: As time goes by: 20 years of changes in the aquatic macroinvertebrate metacommunity of Mediterranean river networks

description: Aim: to analyse temporal metacommunity dynamics in river networks in
 relation to hydrological conditions and dispersal. Location: 15 river
 reaches from the Llobregat, Besòs and Foix catchments in the North-Eastern
 Iberian Peninsula. Taxon: aquatic macroinvertebrates belonging to 99
 different families. Methods: we sampled aquatic macroinvertebrate
 communities during spring in 20 consecutive years. We built two
 environmental distances (one related with water chemistry and another one
 with river flow regime) and two spatial distances (network distance and
 topographic distance). Then we used Mantel tests (accounting for spatial
 autocorrelation) to relate macroinvertebrate dissimilarity with
 environmental and spatial distances. Additionally, we determined the dry
 and wet years using the Standardised Precipitation Index (SPI) and we
 classified macroinvertebrate families based on their ability to fly and to
 drift. Finally, we ran a linear regression model including the correlation
 value (r) of each Mantel test as response variable and distance type
 (environmental or spatial), SPI, dispersal mode, their pairwise
 interactions and a three-way interaction as predictor variables. Results:
 metacommunity organization varied over time and it was significantly
 affected by precipitation, which can be related to river network
 connectivity. The environmental filters, mainly the flow regime, were
 generally more important than the spatial filters in explaining community
 dissimilarity over the study period. However, this depended on the
 dispersal abilities of the organisms. Network fragmentation due to flow
 intermittence during the dry years significantly reduced the dispersal
 capacity of strong aerial dispersers, leading to spatially structured
 metacommunities. For strong drift dispersers, community dissimilarity
 patterns were generally best explained by environmental filters regardless
 of SPI. Main conclusions: a significant temporal variation in
 metacommunity organization can be expected in highly dynamic systems (e.g.
 Mediterranean rivers) and it might depend on the dispersal modes and
 abilities of the organisms, since they determine the response to changes
 in environmental and landscape filters.
description: We studied 15 river reaches from the Llobregat, Besòs and Foix
 catchments in the North-Eastern Iberian Peninsula. The annual
 precipitation ranged from 500 to 1400 mm, and the average annual
 temperature from 7 to 15 °C. The mean flow ranged from 0 to 12,000 l/s,
 thereby including both perennial (i.e. surface flow is maintained
 throughout the year) and temporary (i.e. surface flow ceases during dry
 periods) river reaches. All river reaches were under low human pressure as
 they fulfilled most of the 20 criteria to meet reference conditions in
 Mediterranean rivers. These criteria include a wide range of human uses
 and disturbances on rivers and streams (e.g., diffuse sources of
 pollution, invasive species, land use intensity, riparian vegetation,
 river geomorphology, habitat conditions and hydrological alterations) and
 some general aspects of naturalness, and they have already been used to
 assess the impact of stressors in Mediterranean rivers.
 At each site and sampling date, we recorded water temperature,
 conductivity, oxygen (concentration and percentage of saturation) and pH
 using a multi-parametric digital probe YSI® Pro Plus. We calculated river
 flow (l/s) measuring the river section (width x depth) and water velocity
 with a digital anemometer Schiltknecht® MiniAir2. Additionally, a water
 sample was collected, filtered through glass fibre filters (GF/F; Whatman,
 Maidstone, UK), transported to the laboratory in ice, and frozen for water
 chemistry analysis. At the laboratory, we analyzed major anions (chloride,
 sulphates, nitrites and nitrates) by high-pressure liquid chromatography
 and estimated soluble reactive phosphorus and ammonium concentrations
 using standard colorimetric methods. We collected
 aquatic macroinvertebrates using a circular hand net of 250 µm mesh size.
 Sampling consisted of an initial 3-min kick sample from all available
 habitats. We examined the initial kick sample in the field and we
 collected successive samples until no additional macroinvertebrate
 families were found. Since samples were originally collected for
 biomonitoring purposes, family-level resolution was used. However,
 families should be good surrogates of species-level assemblage patterns.
 Samples were preserved in 10% formaldehyde or 70% ethanol solution
 following the ECOSTRIMED protocol. We sorted the samples and identified
 macroinvertebrates to family level in the laboratory. The abundance of
 each occurring family was quantified and ranked as follows: 1 for 1 to 3
 individuals, 2 for 4 to 10 individuals, 3 for 11 to 100 individuals and 4
 for more than 100 individuals. This ranking followed the application of
 biological indices for assessing water quality, which was the initial
 purpose of this database. We collected all samples
 during the spring season (mostly in May) in 20 consecutive years
 (1997-2017).  <i>Standardised Precipitation Index
 (SPI)</i> We determined the dry and wet years
 using the Standardised Precipitation Index (SPI). SPI represents the
 standardized deviation from a reference series. Positive and negative SPI
 values indicate precipitation greater or lower than the mean of the
 reference series, respectively. It is usually computed using a three-,
 nine- or 12-month span, but we chose the three-month SPI  (i.e.
 accumulated precipitation of the last three months before each sampling
 date) as it reflects short- to medium-term soil moisture
 characteristics and it has been successfully used as a surrogate of river
 flow in our study region. Therefore, we assumed that lower SPI values (dry
 years) corresponded to a higher probability of dry river reaches, whereas
 higher SPI values (wet years) corresponded to a higher probability of high
 flows (supporting information). We calculated the three-month SPI for each
 year from 1997 to 2017 using available monthly rainfall data measured in
 six weather stations of the Meteorological Service of Catalonia from 1950
 to 2017. According to the SPI values, the study periods
 could be classified as extremely dry (SPI ≤ -2), severely dry (-2 &lt;
 SPI ≤ -1.5), moderately dry (-1.5 &lt; SPI ≤ -1), neutral (-1 &lt;
 SPI &lt; 1), moderately wet (1 ≤ SPI &lt; 1.5), severely wet (1.5
 ≤ SPI &lt; 2), extremely wet (2 ≤ SPI). The SPI
 values are standardised and normalised with the following
 formula: SPI = (P<sub>i</sub>−P<sub>m
 </sub>) ÷ S<sub>p </sub> where
 P<sub>i</sub> is the accumulated precipitation of a defined
 time span (i), P<sub>m</sub> is the mean of rainfall of the
 period analysed, and S<sub>p</sub> is the standard deviation
 of the series of precipitation from the same period.
 <i>Classification of the stream flow
 regime</i> We used the TREHS (Temporary Rivers
 Ecological and Hydrological Status) open access software to classify each
 site as perennial or temporary according to their natural hydrological
 regime . We calculated the coefficients Mf (% of months in a year with
 flow), Mp (% of months in a year with isolated pools) and Md (% of months
 in a year with a dry riverbed), which integrate information about the flow
 regime of the rivers over the years. According to this information 8 river
 reaches were permanent and 7 were temporary.
 <i>Classifying macroinvertebrates into dispersal
 groups</i> We classified macroinvertebrate
 families according to their ability to fly and to drift using information
 from a set of biological traits and expert opinion as described in the
 DISPERSE database. For the ability to fly, we considered traits that
 indicate dispersal mode (passive/active in the aquatic and the aerial
 environment), but also other related traits from the larvae and adult
 stages (e.g., female wing length, adult life span). Female wing length
 included 8 categories from &lt;5 mm to &gt;50 mm, whereas adult
 life span included 4 categories from &lt; 1 week to &gt; 1 year.
 For the ability to drift, we considered propensity to drift, which
 included three categories from rare to frequent intentional drift.
 Finally, considering these traits and expert opinion, we assigned an
 affinity to fly and to drift between 1 and 3 (1 for low affinity to 3 for
 strong affinity) to each family, and created two dispersal groups: strong
 flyers (i.e. strong fly affinity) and strong drifters (i.e. medium or
 strong drift affinity). This affinity accounted for intrafamily trait
 variability (i.e. associated to different genera or species) following the
 fuzzy coding approach. We considered taxa with medium drift affinity as
 strong drifters because only three families had the maximum drift
 affinity. The strong flyers group comprised 23 families (mainly Odonata,
 Trichoptera, Heteroptera and some Coleoptera and Diptera) and the strong
 drifters group comprised 20 families (mainly Ephemeroptera, Trichoptera,
 Plecoptera and some Diptera, plus the Crustacea Gammaridae).
 <i>Calculation of distances</i> We
 built a Euclidean matrix (water chemistry distance matrix) based on
 log-transformed and standardised values (mean=0, SD=1) of river flow,
 oxygen, ammonium, nitrites, nitrates, soluble reactive phosphorus,
 sulphates and chloride to characterise local environmental conditions. We
 also built a hydrological dissimilarity matrix through a Gower index based
 on the hydrological regime (Mf, Mp, Md), river flow of each sample and
 river flow of the previous summer. The latter intended to include
 information on preceding hydrological conditions that can be relevant for
 aquatic macroinvertebrates (i.e. the macroinvertebrate community can be
 different if it comes from a wet of dry year). This matrix expressed the
 differences in local hydrological conditions between sites and it was
 termed “flow regime difference”. The dispersal-based
 distances included two physical distances (geographical and network
 distance) and a landscape resistance distance (topographic distance),
 which were used to approximate different potential pathways used by
 organisms to disperse between sites. The geographical distance considers
 the straight line between sites and was calculated as the Euclidean
 distance between geographical coordinates. The network distance measures
 the minimum distance between pairs of sites by following the river
 network. It was measured using a digitized stream network map (ECRIN
 database EcrRiv.mdb) and the Network Analyst extension in ARCGIS 10.0.
 Since the geographical and the network distance were highly correlated, we
 only considered the latter in our analyses, because it better represents
 potential dispersal distances for aquatic organisms. The topographic
 distance considers dispersal through the landscape taking into account the
 topography and it was calculated using a 200-m DEM and the program
 CIRCUITSCAPE. CIRCUITSCAPE uses electronic circuit theory to calculate the
 resistance of the landscape to movements of individuals between each pair
 of sites based on map pixel values (elevation in our case). Because our
 sites were located in different catchments we were unable to calculate the
 network distance between all the sites. Thus, we assigned the maximum
 network distance value found in our dataset (182.15 km) to pairwise
 comparisons between sites of different catchments.
 <i>Statistical analyses</i> We
 calculated alpha diversity (local number of families) for each site and
 year, and beta (pairwise Bray-curtis dissimilarity distance) and gamma
 (total number of families) diversity for each year as broad descriptors of
 macroinvertebrate variation across years. We used
 traditional Mantel tests for spatial distances (topographic and network
 distance) and Mantel tests corrected by spatial autocorrelation through
 Moran spectral randomization (MSR) for environmental distances (water
 chemistry and flow regime difference) (999 runs for each test) to identify
 the relative effect of environmental and spatial filters on community
 composition in each year. The MSR correction was able to remove the
 spurious spatial dependence from our environmental distances, providing a
 correlation value that reflected the net environmental importance. Each
 Mantel test included the Bray-Curtis dissimilarity based on
 macroinvertebrate abundances as response variable and one environmental
 (water chemistry, hydrological) or spatial (topographic, network) distance
 as predictor. We also ran separated Mantel tests for community
 dissimilarity exclusively including taxa either with strong aerial
 dispersal or strong drift dispersal to assess the role of dispersal mode
 in modulating the importance of environmental and spatial filters. To
 estimate the relative importance of environmental and spatial distances in
 explaining community dissimilarity, and given that Mantel tests only
 allowed testing one distance at a time, we summed the Mantel statistic r
 of environmental (water chemistry + flow regime difference) and spatial
 (topographic + network) components. We replaced negative Mantel statistics
 by zero to avoid artificial distortions resulting from a lack of
 independence between the two studied matrices. Before summing Mantel
 r-values, we tested the collinearity between distances. Each pair of
 distances (water chemistry vs. flow regime, and topographic vs. network)
 was only weakly correlated (Pearson <i>r</i>: -0.20 to 0.25),
 suggesting that they explained independent and complementary changes in
 community dissimilarity. To estimate if standardized
 deviations in precipitation (i.e. SPI) modified the importance of spatial
 and environmental filters, we performed a linear regression model
 including the correlation value of each Mantel test (Mantel test
 statistic) as response variable and SPI (linear and quadratic terms),
 distance type (environmental or spatial) and their pairwise interaction
 (SPI x distance type) as predictor variables. Also, to test if the
 importance of spatial and environmental filters varied with dispersal mode
 and precipitation, we ran a linear mixed-effect regression model including
 the correlation value of each Mantel test as response variable and SPI
 (linear and quadratic terms), distance type (environmental or spatial),
 dispersal mode, and their pairwise interactions (SPI x distance type, SPI
 x dispersal mode, distance type x dispersal mode) and three-way
 interaction (SPI x distance type x dispersal mode) as predictor variables.
 This model included year as a random factor. Before building the models,
 we applied a logit-transformation to environmental and spatial correlation
 values (response variables) to reduce variable distribution skewness and
 improve linearity. For each model we adopted a multi-model inference
 approach to quantify the regression coefficients, statistical support and
 importance of the predictor variables, using the function dredge from the
 MuMIN R package. This function ranks all the models generated from all the
 possible combinations of predictors included in the global model through
 the Akaike’s Information Criterion for small sample sizes (AICc). All the
 generated models were ecologically plausible. FWe retained all the models
 differing in less than 7 AICc units (∆AICc≤7) from the model ranking first
 (hereafter, top model). We also estimated the weighted-averaged response
 of top models (weighted by their Akaike weight) and plotted the fitted
 values to facilitate model output interpretation. All models were
 validated by visually checking their residuals for normality and
 homoscedasticity. For some models, we found a temporal dependence in the
 model residuals. In these cases, we added an autoregressive integrated
 moving average (ARIMA) term to the models to account for the lack of
 temporal independence in the residuals. All analyses were conducted using
 the R version 3.5.2.
description: R code and functions to reproduce the results
 0_functions.R<br> Additional R functions
 R script to perform Mantel Test to explore the importance  of
 environmental and spatial distances in explanining community
 dissimilarity <br> 1_Mantel_tests.R R script to
 perform Multi-model inference to explore the importance  of environmental
 and spatial distances in explanining community dissimilarity <br>
 2_final_models_plots.R