Analysing the contribution of intermittent rivers to beta diversity can improve freshwater conservation in Mediterranean rivers

In Mediterranean climate regions, intermittent rivers (IRs) harbor highly dynamic communities with species and trait composition changing over time and space. Considering simultaneously multiple biodiversity facets and a spatiotemporal perspective is, therefore, key to develop effective conservation strategies for these ecosystems. We studied the spatiotemporal dynamics of aquatic macroinvertebrates in rivers of the western Mediterranean Basin by analysing (1) the taxonomic and functional richness and the local contribution to beta diversity (LCBD; measured considering taxonomic and functional facets) of perennial rivers and IRs over �ve sampling times, and (2) their relation with �ow intermittence, local environmental uniqueness, and the number of anthropogenic impacts. Both analyses were also conducted for the subset of data including only IRs to compare values between their �owing and disconnected pool phases. According to our results, taxonomic and functional richness tended to be higher in perennial rivers than in IRs, while IR sites made the greatest contribution to taxonomic and functional LCBD. When comparing among IRs sites over time, higher values of taxonomic and functional LCBD corresponded mostly to their disconnected pool phase. Flow intermittence, the number of impacts and the environmental uniqueness were signi�cant predictors of taxonomic and functional richness, but only �ow intermittence was an important predictor of taxonomic LCBD. For the IR-only data subset, disconnected pool permanence was the main predictor explaining spatiotemporal patterns. Our results highlight the importance of IRs to biodiversity conservation of Mediterranean climate rivers, especially during the disconnected pool phase, suggesting that these ecosystems cannot be ignored in conservation planning strategies.


INTRODUCTION
Incorporating relevant information on community composition analysis is crucial to develop effective biodiversity and ecosystem conservation strategies (Pereira et al. 2013;Hill et al. 2016).Commonly, conservation efforts have focused on protecting the number of species within a site (taxonomic alpha diversity; α) or region (taxonomic gamma diversity or regional diversity; γ), while efforts to characterize and incorporate the variation in species composition (taxonomic beta diversity; β) are relatively more recent (Koleff et  Functional diversity can provide complementary insights into the impacts of disturbance on ecosystem functioning (Flynn et al. 2011;Gross et al. 2017).For instance, the use of functional diversity allows clarifying the role of each species in ecosystem processes and their resistance and/or resilience capacity to environmental changes, either natural or anthropogenic (Tobias and Monika 2011;Villéger et al. 2013Villéger et al. , 2017;;Soria et al. 2020).As the increase of multiple anthropogenic impacts is threatening the stability of ecosystems, understanding how both taxonomic and functional diversity in uence ecosystem functioning can contribute to better predict the ecological consequences of biodiversity loss (Hooper et  Intermittent rivers (IRs) are non-perennial watercourses that typically shift among owing, disconnected pool and dry phases, and are therefore highly dynamic ecosystems (Gallart et al. 2012(Gallart et al. , 2017)).Flow intermittence exerts a primary control on IRs' ecosystem structure and function over time and, consequently, biodiversity patterns typically follow these changes (Datry et   were also conducted for the subset of data including only IRs to compare values between their owing and disconnected pool phases.Perennial rivers and IRs were sampled along gradients of natural ow intermittence and anthropogenic impacts to investigate their in uence on LCBD.We used the LCBD because it allows identifying which sites and times contribute the most to total β diversity.We expected that TRic, FRic, T-LCBD and F-LCBD would be higher in IRs than in perennial rivers because macroinvertebrate communities in IRs are subjected and adapted to higher spatiotemporal variability (Bonada and Resh 2013; Tornés and Ruhí 2013; Leigh et al. 2019).We also hypothesized that T-LCBD and F-LCBD values would be even higher during IRs' periods of disconnected pools than during their owing phase because these ecosystems harbour unique species adapted to non-owing conditions (Bonada et al. 2006(Bonada et al. , 2020)).Finally, in addition to TRic, FRic, T-LCBD and F-LCBD being related to ow intermittence, we also expected them to be positively correlated with environmental uniqueness (Castro et

Study sites and sampling design
The study was conducted in 20 rivers located in the Mediterranean climate region of the Iberian Peninsula.A 100 m site was de ned in each river.Ten sites were intermittent and ten were perennial.The

Biological dataset
Macroinvertebrates were collected during owing and disconnected pool phases by sampling each site ve times during a six-week interval between April and December 2015: April-May (henceforth t1; spring in the northern hemisphere), June (t2; spring), July-August (t3; summer), September (t4; summer), and December (t5; autumn).Therefore, sampling included owing (t1), drying (t2-t4) and rewetting (t5) periods with their corresponding owing, disconnected pools, and dry riverbed aquatic phases.Because ve sites were dry between one and three occasions, a total of 91 samples were obtained.Our sampling procedure followed the o cial quantitative standardized protocol used by water agencies in Spain (MAGRAMA 2013).Samples were collected using a 250 µm-mesh D-net across all available microhabitats and preserved in 4% formaldehyde.Macroinvertebrates were identi ed to the lowest taxonomic resolution possible, usually genus, but with some Chironomidae and Ceratopogonidae identi ed to subfamily or tribe.Overall, 194 macroinvertebrate taxa were identi ed (Soria et al. 2019(Soria et al. , 2020)).
Resistance and resilience traits related to ow intermittence and anthropogenic impacts were obtained to study the functional facet of biodiversity (see details in during the 30-week study period (hereafter ZF T ) and the number of days in the disconnected pool phase since the last sample was taken (hereafter DP i ) (Soria et al. 2020).As the ZF T hydrological predictor informed on how intermittent an IR was, it was used to test the spatiotemporal patterns of taxonomic and functional richness and LCBD between perennial rivers and IRs.Instead, the DPi was used when analysing IRs' spatiotemporal dataset, as it provided information on the permanence of the disconnected pool phase.
For each sampling site and time, the following water quality parameters were measured in situ and analysed in the laboratory when required: conductivity (µS/cm), pH, temperature (ºC), dissolved oxygen (mg/l), chlorophyll-a (mg/m 3 ), HCO 3 , Ca, TOC, Mg and SO 4 (µg/l).Chlorophyll-a concentration, however, was not used because it was correlated with the number of impacts (Soria et

Statistical analysis
The TRic and FRic of each site and sampling time were calculated.FRic was obtained from Soria et al.

RESULTS
Our results suggest that both TRic and FRic tended to be higher in perennial rivers than in IRs over time (Fig. 1a,b and Table S1), while T-LCBD and F-LCBD showed the opposite pattern (Fig. 1c,d and Table S2), especially in summer.When comparing the contributions of perennial rivers and IRs to overall taxonomic and functional β diversity through the ve sampling times, IR sites were the largest contributors.Speci cally, signi cant site contributions to T-LCBD were found in ve IRs corresponding to their disconnected pool phase (t3-t5; Fig. 2a and Tables S1 and S2).For F-LCBD, signi cant site contributions were observed at one speci c sampling time (t2) of a perennial river, three IRs during its owing phase (t1 and t3) and four IRs during its disconnected pool phase (t2, t4-t5; Fig. 2b and Tables S1 and S2).
When considering perennial rivers and IRs, TRic and FRic showed a signi cant negative correlation with all three predictors: ZF t , environmental uniqueness and the number of impacts (Table 1).Our results also showed a signi cant and positive relationship between T-LCBD and ZF t , but no signi cant results were observed for environmental uniqueness and the number of impacts (Table 1).None of the explanatory variables were correlated with F-LCBD (Table 1).

Table 1 Partial regression coe cients of the explanatory predictors used in linear mixed-effects models.
The response variables were the local taxonomic and functional richness (TRic and FRic, respectively) and the taxonomic and functional local contribution to β diversity (T-LCBD and F-LCBD, respectively).ZF t : total zero-ow days.Signi cant coe cients (P < 0.05) are indicated in bold.For the subset of data including only IRs, TRic showed a signi cant negative correlation with all three predictors, but FRic was only signi cantly and negatively correlated with the number of impacts (Table 1).T-LCBD and F-LCBD were signi cantly correlated with DP i , but no signi cant patterns were observed for the other predictors (Table 2).

DISCUSSION
Overall, our results showed that IRs had higher taxonomic LCBD than perennial rivers, despite their lower local taxonomic and functional diversity.In early summer and in response to the loss of surface ow, the observed decline in taxonomic and functional richness in IRs can be explained by the disappearance of lotic species inhabiting ri e habitats (Bogan et al. 2017;Tonkin et al. 2017).Yet, during the disconnected pool phase in summer, the increase in taxonomic and functional LCBD in IRs can be related to the colonization of species with speci c traits adapted to cope with such conditions (Bonada et al. 2006(Bonada et al. , 2020)).A similar pattern, where sites with higher LCBD are those with lower richness, has been found in other studies of macroinvertebrates (e.g.On the other hand, our results suggest that ow intermittence (i.e.indicated by ZF T ), the number of impacts and the environmental uniqueness were signi cant predictors of taxonomic and functional richness, but only ow intermittence was an important predictor of taxonomic LCBD.Indeed, the permanence of disconnected pools (i.e.indicated by DP i ) was the main predictor explaining the spatiotemporal patterns of IRs.However, there was no signi cant correlation with the number of impacts, despite some impacted IRs also showed higher LCBD values.This could be explained by the fact that resistance and resilience traits developed in response to ow intermittence in IRs (e.g.multivoltism, aerial respiration or mechanisms to tolerate low dissolved oxygen concentrations) may also be useful in coping with anthropogenic impacts (Bonada and Resh 2013; ).Indeed, as the distribution of disconnected pools in IR networks can vary spatially and temporally from one year to another, tools have been developed to even account for their temporal and spatial occurrence (e.g.TRESH software; Gallart et al. 2017), as well as to assess their priority as biodiversity refuges and incorporate them into conservation planning (Hermoso et al. 2013;Yu et al. 2022).Citizen science can also result in a powerful tool to fully understand the hydrological characteristics of these ecosystems (e.g.RiuNet App, DRYRivERS App).The integration of these tools in conservation management of IRs, together with the use of community metrics able to capture their spatiotemporal biodiversity patterns, are key to improve freshwater conservation in the Mediterranean climate region.
al. 2003; Anderson et al. 2011; Socolar et al. 2016; Bush et al. 2016).Taxonomic β diversity can be measured as the compositional dissimilarity in species assemblages, either across space (Baselga 2010; Anderson et al. 2011) or time (Legendre and Gauthier 2014; Legendre 2019; Shimadzu et al. 2015).While links between α, β and γ spatial taxonomic diversity have been widely explored, their temporal dynamics have been overlooked (Anderson et al. 2011; Legendre and Gauthier 2014; McGill et al. 2015).Temporal taxonomic β diversity informs how species composition changes over time (Legendre and Gauthier 2014; Legendre and Condit 2019).Despite the combination of spatial and temporal β diversity having the potential to better represent changes in community composition in response to environmental change, few studies have considered the contribution of both components to regional diversity (e.g.Vander Vorste et al. 2021; Faustino de Queiroz et al. 2022).
al. 2014; Arroita et al. 2018; Arias-Real et al. 2020).In the case of aquatic macroinvertebrates, surface ow cessation and the subsequent formation of disconnected pools imply the disappearance of ri e-dwelling species and the appearance of pool-dwelling species (Bonada et al. 2006; Bogan et al. 2017; Tonkin et al. 2017).With the complete drying of the riverbed, some taxa may emerge, move to other wet habitats, to the hyporheic zone, or enter in a desiccation-resistant life stage (Lytle and Poff 2004; Bogan et al. 2017; Stubbington et al. 2019).Shifts from dry to owing phases following rewetting favour recolonization, contributing to the recovery of local communities in IRs (Leigh et al. 2016; Bogan et al. 2017).IRs communities are highly variable in time, with species and trait composition changing from one period to another (Datry et al. 2014; Bogan et al. 2017).Nevertheless, climate change and increasing human water demand are altering biodiversity patterns and functional processes of IRs (Datry et al. 2014), and thus a better understanding of their contribution to biodiversity is timely (Ruhí et al. 2017; Sánchez-Montoya et al. 2020).In comparison to IRs from other climatic regions, those in Mediterranean climate areas are characterised by being highly predictable in terms of seasonality, resulting in well-known community shifts between owing and non-owing phases (Hershkovitz and Gasith 2013; study area is characterized by Mediterranean climate, with high seasonal and inter-annual variability in precipitation and ow regime (Bonada and Resh 2013; Cid et al. 2017).Sites ranged from 6 to 1100 m.a.s.l., with discharges ranging from 0 to 0.417 m 3 /s.See Soria et al. (2020) for further details of the study area.
(2020)  and estimated as suggested byVilléger et al. (2008).T-LCBD and F-LCBD were estimated for each sample, i.e. for each site and sampling time (spatiotemporal approach;Legendre  and De Cáceres 2013; Legendre and Gauthier 2014).T-LCBD and F-LCBD values indicate the degree of ecological uniqueness of each site at each sampling time in terms of taxa abundance or traits composition, respectively (Legendre and De Cáceres 2013).The T-LCBD was calculated using the taxa abundance matrix (Legendre and De Cáceres 2013).The F-LCBD was calculated using a matrix with the relative abundance of each trait category (columns) across the samples (rows) (Rodrigues-Capítulo et al. 2009).Following the procedures described by Legendre and De Cáceres (2013), T-LCBD and F-LCBD were estimated using the Euclidean distance on Hellinger-transformed data (Legendre and Gallagher 2001) and the signi cance of each LCBD value (i.e. each site in each sampling time) was tested with 999 permutations.Linear mixed-effect models (LME) were used to test the relationship between TRic, FRic, T-LCBD and F-LCBD, and our set of environmental predictors.Speci cally, the variables of ZF T , environmental uniqueness and the number of impacts were used when considering the whole dataset, while DPi, environmental uniqueness and the number of impacts were used when considering data from IRs only.All models included sampling time (t1-t5) as a random factor.All analyses were conducted in R software version 3.6.2(R Core Team 2015), using the packages "ade4"(Dray et al. 2007), "adespatial"(Dray et al. 2018), "car" (Fox and Weisberg 2011), "maptools" (Bivand and Lewin-Koh 2020), "nlme"(Pinheiro et al. 2016), "raster" (Hijmans 2020), "rgdal" (Bivand et al 2020), "rgeos" (Bivand and Rundel 2020), and "vegan"(Oksanen et al. 2013).The code and functions used to run these analyses are available at Data Accessibility Statement.
, NC, LMB and NB conceived the ideas and designed the methodology.NC, PR-L, PF and DV collected the data; RA identi ed the macroinvertebrate data; MS, JO and L-MB analysed the data.MS, JO, LMB, NB and NC led the writing of the manuscript, and CG-C, RA, PR-L, PF, DV, FG and NP contributed to the drafts.All authors gave nal approval for submission.FUNDING AND/OR COMPETING INTERESTSThe study was supported by the LIFE+ TRivers (LIFE13 ENV/ES/000341) project.MS was supported by the Scholarship for studies or projects outside Catalonia funded by "Fundació Universitària Agustí Pedro i Pons" and the Iberoamerican Mobility Scholarship funded by "Banco Santander".CG-C was supported by a Junior Leader Fellowship contract (LCF/BQ/PR22/11920005) funded by "la Caixa" Foundation (ID 100010434).PR-L was supported by a Margalida Comas postdoctoral contract (PD/031/2018) funded by the Government of the Balearic Islands and the European Social Fund and by a Juan de la Cierva-Incorporación fellowship (IJC2019-041601-I).

Figures
Figures
Soria et al. 2020)al regime (i.e.perennial and IRs) and to differentiate river sites affected by natural ow intermittence from those with human-driven ow intermittence.To infer the IRs' phases (i.e.ow, disconnected pools, dry riverbed), two temperature data loggers (UA-002 HOBO) were installed at each river site, which recorded data during the 30-week study period (see details inSoria et al. 2020).Thermal amplitude was then used to calculate the total zero-ow days (i.e.disconnected pools or dry riverbed) Soria et al. 202019)astro et al. (2019), the method described by Legendre and De Cáceres (2013) was used to calculate a measure describing the environmental uniqueness of each site.To do so, all environmental variables (except pH) were log-transformed and then a standardized Euclidean distance matrix was calculated.Samples with high values (closer to one) of the resultant vector (containing the local contribution to environmental variability) are more unique in terms of environmental conditions(Castro et al. 2019).(i.e. from 0 = non-impacted, to 20 = highly impacted, hereafter number of impacts; seeSoria et al. 2020).
(Sánchez-Montoya et al. 2009)riteria (MRC index)(Sánchez-Montoya et al. 2009).The MRC index includes information on invasive species, diffuse pollution sources, land-use intensity, riparian vegetation, river geomorphology, instream habitat conditions and hydrological alterations, and ranges from 0 (highly impacted) to 20 (non-impacted).To facilitate interpretation, the inverse of the MRC index values was used

Table 2
Partial regression coe cients of the explanatory predictors used in linear mixed-effects models for intermittent rivers only.The response variables were the taxonomic and functional local contribution to β diversity (T-LCBD and F-LCBD, respectively), and the local taxonomic and functional richness (TRic and FRic, respectively).DP i : the number of days in the disconnected pool phase since the last sample was taken.Signi cant coe cients (P < 0.05) are indicated in bold.
(Gallart et al. 2017;elationship indicates that less diverse sites tend to hold unique species and trait compositions, highlighting the complementarity of alpha and β diversity indices in describing biodiversity(Heino et al. 2017;da Silva et al. 2018).Incorporating temporal patterns of LCBD is therefore key to capturing the full variation in community composition that exists in these highly dynamic ecosystems(Ruhí et al. 2017;Sánchez-Montoya et al. 2020).Our study also showed that higher values of taxonomic and functional LCBD occur mostly during the disconnected pool phase of IRs.During this phase, taxa inhabiting owing conditions tend to be lost, while lentic taxa progressively colonise the remaining disconnected pools from nearby sites that are drying up, such as Odonata, Coleoptera and Hemiptera(Bogan et al. 2017;Bonada et al. 2020).Indeed, there are studies showing that lentic-dwelling species have higher dispersal abilities than lotic species (e.g.Ribera and Vogler 2000; Hjalmarsson et al. 2015).Disconnected pools can act as refuges for some aquatic taxa such as sh or amphibians during IRs' dry season, which are fundamental to recolonize the river network upon ow resumption(Hermoso et al. 2013;Gallart et al. 2017).In addition, for some species of macroinvertebrates (and also amphibians and sh), disconnected pools are also used as stepping-stones for their dispersal or as key sites for laying eggs and, thus, complete their life cycle(Bonada et al. 2006(Bonada et al. , 2020;; Stubbington et al. 2017;Moidu et al. 2023).Considering that up to 60 per cent of the world's rivers by length are IRs and that they are expected to increase worldwide (Döll and Schmied 2012; Messager et al. 2021), it is expected that disconnected pools will become more abundant and, consequently, management actions will be needed to conserve them(Gallart et al. 2017; Bonada et   al. 2020).
Heino et al. 2017; da Silva et al. 2018; Valente-Neto et al 2020) and other biological (Gallart et al. 2017;Bonada et al. 2020)Villéger et al. 2017;Crabot et al. 2020)oulton et al. 2000).In fact, Mediterranean IRs hold unique species composition adapted to natural ow intermittence, such as a dominance of pooldwelling species during the disconnected pool phase(Bonada et al. 2006(Bonada et al. , 2020;;Cid et al. 2017), which could give them the ability to resist and to recover from drying periods and, at the same time, from anthropogenic impacts.However, persistent and intensifying anthropogenic impacts over time could reduce the ability of these species to cope with ow intermittence(Datry et al. 2017).Considering that IRs have been commonly ignored in conservation planning (Bogan and Lytle 2007; Leigh et al. 2019), it is timely to provide complementary measures to adequately assess their biodiversity.In this sense, the LCBD approach shows a high potential to be used for conservation purposes, as the ecological uniqueness of a site can be compared with other sites sampled in a region(Legendre and De Cáceres 2013; da Silva et al. 2018; Valente-Neto et al 2020).In highly dynamic systems such as IRs, not only spatial but temporal patterns should also be considered to better identify key sites (Ruhí et al. 2017; Rogosch and Olden 2019; Stubbington et al. 2019; Sánchez-Montoya et al. 2020).However, analysing only spatiotemporal species richness and community composition might not be su cient to design conservation plans aimed at protecting the processes that maintain their ecosystem functioning(Leigh et al. 2016(Leigh et al. , 2019;;Villéger et al. 2017;Crabot et al. 2020).Therefore, given the hydrological variability of IRs and the increasing anthropogenic impacts they receive, freshwater conservation planning should consider monitoring the temporal variability of both taxonomic and functional biodiversity in these ecosystems.This might be even more relevant in Mediterranean-climate regions worldwide where IRs constitute one of their predominant freshwater ecosystems(Bonada and Resh 2013;Cid et al. 2017).Special attention should also be given to the disconnected pool phase, as this is key to maintain local and regional aquatic biodiversity(Gallart et al. 2017;Bonada et al. 2020).In this regard, several management-related tools have recently been developed to better predict the ow patterns of IRs, such as wet-dry mapping, remote sensing (e.g.satellite images, xed cameras), eld sensors (e.g.