DOI: 10.5061/dryad.brv15dvkb
title: Data from: Community size and disturbance history jointly explain the interplay between stochastic and deterministic community variation in benthic invertebrates

description: Community ecology has so far struggled to integrate both deterministic and
 stochastic processes into a global model of community variation. To
 address this issue, we aimed to characterise the ecological conditions
 that determine the relative importance of deterministic and stochastic
 community variation in benthic macroinvertebrates. We sampled
 macroinvertebrate and algal periphyton communities, and microhabitat
 conditions at monthly intervals over a 5 year period in two sampling sites
 in the regulated Durance River and one site in the Asse River (a natural
 tributary of the Durance). The relative importance of stochastic and
 deterministic processes was estimated over time as functional and
 taxonomic spatial β-deviation (i.e. β-diversityobs – mean β-diversitynull
 / SD β-diversitynull) for each sampling date and site. We additionally
 quantified deterministic variation as the trait-environment relationship
 and predicted a positive correlation with β-deviation metrics. We found
 evidence of overdispersal in both functional and taxonomic β-diversities
 compared to null expectations. Spatial β-deviation was highly variable
 over temporal scales and was jointly explained by invertebrate community
 size, disturbance (i.e. flooding), periphyton diversity and to a lesser
 extent, microhabitat diversity. The key factor that explained β-deviation
 was recent river discharge (&lt; 90 days), which had a strong negative
 effect on community size but was also directly associated with higher
 β-deviation, likely related to recolonization processes post-flood.
 Lastly, we found that when β-deviation was high, the hydraulic
 trait-environment relationship was stronger. This indicates that spatial
 heterogeneity in hydraulic conditions was the main driver of deterministic
 variation in invertebrate community structure. This study demonstrates
 that the importance of stochastic processes can vary greatly over short
 (~months) and long (~years) temporal windows with significant consequences
 for trait-environment relationships. Our study provides an example of how
 stochastic community variation can be modelled to better interpret
 deterministic patterns of community structure in natural communities.
description: <strong>Hydrological conditions</strong>
 We obtained hourly hydrological data (Q
 m<sup>3</sup>s<sup>-1</sup>) for the period of
 December 2012 to August 2018. Hydrological data for the two Durance sites
 was obtained from the Escale hydroelectrical dam (28.3km & 34.7km
 upstream of the Dur-1 and Dur-2 sites, respectively). We also obtained
 hydrological data for the Bléone River (906 km² drainage basin) as both
 Durance sites were situated downstream from its confluence with the
 Durance (Figure 1). The Bléone River is also regulated by a dam but as
 discharge measurements were not available. We obtained simulated discharge
 values inferred from reservoir water levels for the sample period from the
 French electrical company (EDF). These simulated discharge values were
 then added to the Durance hydrograph. For the Asse site, we obtained
 hydrological data from the Clue Chambrières limnology station (14.7km
 upstream from the Asse study site). We also added the input of the Asse to
 the Dur-2 hydrograph, as Dur-2 is located downstream from the confluence
 between the Asse and Durance rivers. Lastly, we standardised river
 discharge (Q m<sup>3</sup>s<sup>-1</sup>) to
 account for differences in the size of catchment basins between our
 sampling sites. River discharge was converted to mm/h and then
 standardised by catchment area using the following formula: Q(l/s)*T(s) /
 (S km² *1000*1000). Where Q is discharge (in l/s), T is the number of
 seconds in the desired time unit (i.e. 3600 seconds in an hour) and S is
 the surface area of the catchment basin (Dur-1 = 8131km², Dur-2 = 9185km²
 and Asse = 538km²). <strong>Macroinvertebrate and
 periphyton sampling protocol</strong> We
 performed monthly sampling beginning in late winter (~February) and ending
 in early summer (~July) from 2014 to 2018 (Table 1). We selected this
 period in order to capture the main growth period of invertebrate and
 periphyton communities in our study sites (Cazaubon & Giudicelli,
 1999). Sampling was designed to be representative of
 the various hydrologic ambiances found in riffle habitats. To this end, we
 defined three hydraulic ambiances: 30-50 cm s<sup>-1</sup>,
 70-90 cm s<sup>-1 </sup>and &gt; 100 cm
 s<sup>-1</sup>, wherein we obtained three macroinvertebrate
 and periphyton samples at each sampling date (i.e. three samples per
 ambiance = nine samples total). Samples were distributed haphazardly
 within each ambiance with a minimum of three meters between two samples.
 The length of the sampling reach varied according to the availability of
 the three ambiances, ranging from 50m to ~1000m in length.
 Macroinvertebrate samples were obtained using a Surber sampler (0.10 m²)
 with a 0.25 mm mesh. Samples were immediately stored in 70% ethanol,
 awaiting morphological identification. In the laboratory,
 macroinvertebrate samples were rinsed and further sieved (0.63 mm mesh)
 before being counted and identified under binocular microscope. Species-
 or genus-level Identification was obtained for the majority of specimens
 with the exception of Diptera (which were identified to the family-,
 sub-family- or tribe-levels) and other groups like Oligochaeta and
 Hydrachnidia which could not be reliably identified. Nine periphyton
 samples were obtained in parallel to each of the nine macroinvertebrate
 samples in the same sample area. A periphyton sample was constituted of
 three medium-sized stones found beside the Surber sample. The biofilm on
 each stone was scrapped from a standardised 25 cm² area (drawn on by
 pencil), using a scalpel and/or nylon brush and transferred to a container
 with a solution of formaldehyde 2-5%. The resulting solution from each of
 the three stones (total area sampled = 75 cm²) was then pooled into a
 single container. Note that the area sampled changed after 2014 (when 100
 cm<sup>2</sup> area was sampled on each stone) to account for
 the large amount of material that was collected. In the laboratory, each
 sample was brought to a fixed volume of liquid (50 mL) and put into
 suspension to homogenize the samples. Counting and identification of the
 periphyton was performed using a microscope and expressed as the number
 (of cells) per cm<sup>2</sup>.
 <strong>Micro-habitat description</strong>
 We described microhabitat conditions that are pertinent to the
 ecology of benthic organisms (granulometry, clogging, depth, velocity,
 etc.) at the scale of each Surber sample (i.e. 0.10 m²). The granulometry
 of the riverbed was quantified using semi-quantitative size-classes
 described by Malavoi & Souchon, 1989. The substrate size-class
 scale ranged from 0 (silt: 0.0039 – 0.0625mm diameter) to 10 (bedrock:
 &gt; 1024mm diameter; Supporting information). We used substrate
 size-classes to estimate the surface and underlayer substrate richness,
 largest substrate size and dominant substrate size. Substrate clogging was
 evaluated visually using a semi-quantitative scale ranging from 1 to 5
 (Supporting information) designed to describe the level of substrate
 embeddedness and clogging of the interstitial space by silt and algae
 (Archambaud et al., 2005). Current velocity was measured at regular
 intervals beginning at 3cm from the river bottom and at approximately
 0.20, 0.40 and 0.80 water depth, using an electromagnetic current meter
 (OTT MFT Pro). 
description: # Data from: Community size and disturbance history jointly explain the
 interplay between stochastic and deterministic community variation in
 benthic invertebrates
 [https://doi.org/10.5061/dryad.brv15dvkb](https://doi.org/10.5061/dryad.brv15dvkb) ## Description of the data and file structure ## Files and variables #### File: FunctionalTraits.xlsx **Description:** Functional trait database adapted from Tachet et al. 2010 ##### Variables **FunTraits excel sheet:** *  **Final_assignment**: Final taxonomic assignment of macroinvertebrate taxa. * **FinalFun**: Final taxonomic assignment for functional traits. This is a tradeoff between the resolution of the Tachet et al. functional trait database and the final taxonomic assignment (Final_assignment). When final assignment is more precise than the functional trait database, taxa were homogenised to the next best taxonomic level. For example, *Baetis* spp. were homogenised to the genus level. * **TaxLevel**: The taxonomic resolution aquired for FinalFun * **NbGene_1 - Vit_4:** Fuzzy-coded functional traits for macroinvertebrates. See the trait info sheet for additional information. **Trait info excel sheet:** * **Trait name**:  Fuzzy-coded functional traits for macroinvertebrates.Corresponds with names from the FunTraits sheet. * **Trait_type:** Indicates the type of functional trait, note F refers to fuzzy-coded traits * **Fuzzy_name**: Indicates the fuzzy trait, e.g. Nb.gen, number of generations * **Extra_info:** Provides additional information for each fuzzy trait #### File: CommunityData.xlsx **Description:** Invertebrate and periphyton community sampling ##### Variables * **Sample_ID**: Sample ID can be broken down into: site name =&gt; Hydrological ambiance (1,2 or 3), Sample (A, B or C), Sampling campaign (C0 to C6; see Table 1) ==&gt; Sampling year * **Site**: Sampling site * **Campaign**: Sampling campaign, refers to the order of sampling in a given year * **Year**: Sampling year * **Date:** Sampling date * The remaining columns provide abundance data for invertebrate (individuals per 0.10m²) and periphyton (cells per 0.01m²) taxa. #### File: Microhabitat.xlsx **Description:** Description of microhabitat conditions, samples correspond to invertebrate and periphyton sampling ##### Variables * **Sample_ID**: Sample ID can be broken down into: site name =&gt; Hydrological ambiance (1,2 or 3), Sample (A, B or C), Sampling campaign (C0 to C6; see Table 1) ==&gt; Sampling year * **Site**: Sampling site * **Campaign**: Sampling campaign, refers to the order of sampling in a given year * **Year**: Sampling year * **Date:** Sampling date * **Ambiance**: Hydrological ambiance of the sample. Hydrological ambiances are broken down into: (1) 30-50 cm s-1, (2) 70-90 cm s-1 and (3) &gt; 100 cm s-1 * **Largest_substrate**: Largest substrate size-class in sampling area. Refer to the README panel in the excel file for additional information. * **Smallest_substrate**: Largest substrate size-class in sampling area. Refer to the README panel in the excel file for additional information. * **Dominant_substrate**: Dominant (most abundant) substrate size-class in sampling area. Refer to the README panel in the excel file for additional information. * **Substrat_accessory1-6**: Presence of accessory substrate size-class in sampling area in order of abundance. Up to 6 accessory substrates were recorded if present. Refer to the README panel in the excel file for additional information. * **Substrate_richness**: The total number unique substrate size-classes in the sampling area. Refer to the README panel in the excel file for additional information. * **Largest_underlayer**: Largest underlayer substrate size-class in sampling area. Refer to the README panel in the excel file for additional information. * **Smallest_underlayer**: Largest underlayer substrate size-class in sampling area. Refer to the README panel in the excel file for additional information. * **Dominant_underlayer**: Dominant (most abundant) substrate size-class in sampling area. Refer to the README panel in the excel file for additional information. * **Undercover_accessory1-3**: Presence of accessory underlayer substrate size-class in sampling area in order of abundance. Up to 3 underlayer accessory substrates were recorded if present. Refer to the README panel in the excel file for additional information. * **Underlayer_richness**: The total number unique substrate size-classes in the sampling area. Refer to the README panel in the excel file for additional information. * **Clogging**: The degree (0-5) of interstitial substrate clogging by fine sediment and/or algae. Refer to the README panel in the excel file for additional information. * **Depth**: River depth (cm) * **Water velocity (cm s-1) and associated depth measurements:** | v3cm.s | \~ | Water velocity at 3cm from riverbed | | :------ | :- | :----------------------------------- | | h\_0.2h | \~ | Height (from river bed) at 20% depth | | v0.2h | \~ | Water velocity at 20% depth | | h\_0.4h | \~ | Height (from river bed) at 40% depth | | v0.4h | \~ | Water velocity at 40% depth | | h\_0.8h | \~ | Height (from river bed) at 80% depth | | v0.8h | \~ | Water velocity at 80% depth | #### File: Hydrology.xlsx **Description:**  | Hourly discharge values (Q m3/s) for Asse, Dur-1 and Dur-2 sites | | :--------------------------------------------------------------- | ##### Variables * README is included in the xlsx file ## Code/software All statistical analyses and data manipulation were performed in R v4.20 (Core Development Team, 2023) running in RStudio (RStudio, 2011). ### Spatio-temporal variation in macroinvertebrate and periphyton communities We performed Principal Components Analysis (PCA) to characterise spatiotemporal variation in macroinvertebrate and periphyton communities. PCAs were based on total invertebrate and periphyton community density (inv/0.10m² & cells/cm² for invertebrates and periphyton, respectively), taxonomic richness (based on the highest possible taxonomic identification), Shannon and Simpson diversity indices, and taxonomic evenness. Lastly, as benthic communities have a tendency to be dominated by a few hyper abundant taxa (e.g. Battle et al., 2007; Minshall et al., 1992), we calculated the percentage contribution of the most abundant taxon to the total sample density (hereafter dominance). Functional diversity indices were also calculated for the macroinvertebrate community. Functional diversity indices were based on fuzzy-coded traits from Tachet et al., 2010, and included ecological (microhabitat, temperature and water velocity preferences, and feeding strategy) and biological (respiration, diet, locomotion, dispersion strategy, lifestyle, reproduction, number of generations per year) characteristics (Supporting information). Differences in taxonomic assignment between the invertebrate and functional datasets were resolved by summarising groups to the highest shared taxonomic-level (genus, tribe or family; taxa identified to order-level or higher were excluded from functional analyses) and averaging fuzzy trait values (Supporting information). Functionally identical taxa were regrouped into functional entities using the function *sp.to.fe *from the *mFD *package (Magneville et al., 2022). Fuzzy-traits were then transformed into two continuous variables (axes 1 & 2) using Principal Coordinates Analysis (PCoA) based on Gower distances (i.e. functional dissimilarity between taxa). We then calculated functional richness, evenness and divergence using *fd_fric*, *fd_feve* and *fd_fdiv *functions, respectively from the *fundiversity* package (Grenié & Gruson, 2023). The *fundiversity* package calculates functional metrics according to Villéger et al., 2008, wherein functional richness corresponds to the volume of the convex hull of multidimensional functional space in each sample. Functional evenness corresponds to the minimum spanning tree which links all the species in the multidimensional functional space. Functional divergence is quantified by the extent to which taxa diverge in their distance from the centre of functional space, and is weighted by taxa abundance. ### Calculating β-deviation The objective of the β-deviation metric was to quantify the degree of non-random invertebrate community structure at each sampling site and date. We calculated β-deviation as the extent to which observed spatial dissimilarity (β-diversity) within macroinvertebrate communities deviated from null expectations (i.e. based on random community reshuffling). β-deviation was estimated for both taxonomic and functional (see section 3.5.1 for trait information) β-diversity indices, separately, with the aim of providing complementary information on the processes that drive community assembly. β-diversities were calculated separately for each sampling date and site (i.e. among the nine Surber samples per campaign) as pairwise Hill numbers (q = 1) based on a Jaccard-type complement. Taxonomic β-diversities were calculated using the *pair_dis *function from the *hilldiv *package* *(Alberdi & Gilbert, 2019)* *and functional β-diversities were calculated using the *beta_fd_q* function (tau = mean, package *mFD; *Chao et al., 2019; Magneville et al., 2022) based on Gower distances between functional entities (function *funct.dist, *package *mFD;* (Magneville et al., 2022). The choice of null model is particularly important as it imposes constraints on the null distribution generated by the model, thus modifying the ecological significance of β-deviation (Gotelli & Ulrich, 2012; Harvey et al., 1983; Vellend et al., 2014). We therefore selected the most widely used null model based on the reasoning of Kraft et al., 2011 which randomly assigns individuals to sites while preserving the overall species abundance distributions and the total number of individuals in each site (e.g. Shi et al., 2021; Siqueira et al., 2020; Tucker et al., 2016). Abundance-based null model approaches like the Kraft et al. null model have been shown to be more robust than β-deviation estimates based on presence/absence data (Tucker et al., 2016). It should be mentioned that the Kraft et al. null model has been criticised, notably, for constraining species abundance distributions. Qian et al., 2013 has pointed out that species abundance distributions actually arise largely due to deterministic processes (e.g. Alonso et al., 2008). However, as the abundance distribution of invertebrate communities varies greatly over seasonal timeframes, we justify our use of the Kraft et al. model by our need to explicitly control for such temporal variation in species abundances. β-deviation was calculated as the standardised effect size of a two-sided test comparing the simulated and observed β-diversity values using the *oecosimu *function from the *vegan* package (n = 1000, method = *r2dtable*; Oksanen et al., 2008; Patefield, 1981). There are three possible outcomes for this test: (i) the mean β-diversity value falls within 95%CI of simulated values, indicating that dissimilarity could have arisen solely from stochastic processes, (ii) the mean β-diversity value exceeds the 95%CI of simulated values, indicating overdispersal of community variation and (iii) the mean β-diversity is significantly lower than the 95%CI of simulated values, indicating clustered community variation. Positive (i.e. outcome ii) and negative (outcome iii) β-deviation values are expected to arise from a dominant role of heterogenous and homogenous selection processes, respectively. It should be noted that non-significant β-deviation values (outcome i) can also arise due to opposing deterministic forces of community assembly. For example, an equal contribution of heterogenous and homogenous selection would be expected to result in non-significant β-deviation. In the case of a non-significant β-deviation value (i.e. *p* &gt; 0.05), the value was set to zero. This two-sided test assumes that the null distribution is normally distributed. An assumption that is often invalidated in the case of null model generated metrics (Bernard‐Verdier et al., 2012), including the present study. We therefore rank transformed our β-deviation metrics into 10 distinct classes in order to reduce the magnitude of differences among standardised effect sizes driven by non-normal null distributions. ### The trait-environment relationship We quantified the trait-environment relationship as an additional indicator of deterministic community structure (Cottenie, 2005), and tested its relationship with β-deviation. The trait-environment relationship was quantified as the goodness-of-fit between invertebrate functional β-diversities and environmental β-diversities. Invertebrate and environmental β-diversity were both determined within the context of each sampling date at each site. Invertebrate assemblages were compared using functional β-diversities (see section 3.5.2). Environmental β-diversity was separated into two abiotic components including (i) substrate and underlayer characteristics (dominant substrate, size-class richness, smallest and largest substrate size-class, and clogging class), and (ii) hydraulic characteristics (water velocity mean, max, min, coefficient of variation and at 3cm from the riverbed), hydraulic ambiance (30-50 cm s-1, 70-90 cm s-1 or &gt; 100 cm s-1) and water depth. To calculate abiotic environmental β-diversity we first scaled all microhabitat variables to zero mean. Microhabitat samples were then plotted using Principal Components Analysis (function* dudi.pca*, *package ade4; *Dray & Dufour, 2007) and β-diversity was determined as the Euclidean distance between points. We also quantified (iii) periphyton β-diversity as a biotic dimension of environmental β-diversity. To do so, periphyton community data (log + 1 transformed) were compared using pairwise Bray-Curtis distances. The trait-environment relationship was quantified as the Pearson’s r correlation between functional invertebrate β-diversities and environmental β-diversities (i.e. β-diversityInv. ~ β-diversityenv) for each sampling date at each site and environmental β-diversity component (i.e. substrate, hydraulic and periphyton), separately. This resulted in three trait-environment relationships related to: substrate, hydraulic and periphyton β-diversities. ### The relationship between ecological conditions, β-deviation and the trait-environment relationship   To characterise the relationship between ecological conditions, β-deviation and the trait-environment relationship we fitted causal models. We considered four distinct categories of ecological variables: disturbance, phenology, abiotic conditions and biotic conditions (periphyton and invertebrate communities; see Supporting information for *a priori *model structure). However, due to incomplete periphyton data (Table 1) we fitted two causal models (model *i*) one using all data types (n = 56) and (model *ii*) a second excluding periphyton data (n = 80). Ecological variables were all zero-mean standardised. Prior to model selection we verified the normality of taxonomic and functional β-deviation metrics, which were both skewed towards zero. We therefore log + 1 transformed both β-deviation metrics prior to rank transformation (see section 3.5.2) to meet the assumption of normality. The first step was to select the best performing ecological predictors of β-deviation. Variable performance was evaluated for each ecological compartment according to RandomForest analysis (mtry = 10, ntree = 5000, 100 runs, package *randomForest; *Liaw & Wiener, 2002), based on mean square error divided by its standard deviation. The two best performing disturbance, microhabitat, periphyton and macroinvertebrate variables were selected for causal analysis. Causal analysis was performed using the R package *piecewiseSEM* (Lefcheck, 2016). The global model was composed of multiple mixed effect models, all of which included a year and site random effect. The final model was obtained by removing all non-significant causal paths. We also tested the significance of dependencies between covariates (function *dSep*), and included significant dependencies in the final model. We tested the relationship between β-deviation and the trait-environment relationship (i.e. invertebrate β-diversity ~ environmental β-diversity corelation statistic) using the same mixed model structure as indicated above. The β-deviation ~ periphyton trait-environment relationship was tested in model *i *(n = 56), while the β-deviation ~ hydraulic and substrate trait-environment relationships were tested in model *ii *(n = 80). Lastly, to detect site-specific effects, we visualised and tested the relationship between the trait-environment relationships and β-deviation (Pearson’s r correlation and associated *p-*value) by site and type. ## Access information Data was derived from the following sources: * Functional traits: Tachet, H., Richoux, P., Bournaud, M., & Usseglio-Polatera, P. (2010). *Invertébrés d’eau douce: Systématique, biologie, écologie* (Vol. 15). CNRS éditions Paris