Romani Am

Multi-model assessment of hydrological and Environmental impacts on streambed microbes in Mediterranean catchments

Microbes inhabiting intermittent streambeds are responsible for controlling and developing many biogeochemical processes essential for the ecosystem functions. Although streambed microbiota is adapted to intermittency the intensification of water scarcity and prolonged dry periods may jeopardise their capacity to cope with hydrological changes. This study aims to evaluate whether, and to what extent, the duration of dry periods affects streambed microbial density, diversity, composition (16S rRNA gene diversity) and functions (extracellular enzyme activities and respiration). Our results highlight the fact that hydrology modulates the community composition and, to some extent, the functions carried out under different environmental conditions. The relative abundance of certain taxa inhabiting the driest intermittent communities differs significantly from those found at sites with continuous flow.

Microbial functional metrics revealed a progressive increase in recalcitrant carbon degradation activity at sites with an extended dry phase. In contrast, bacterial density and diversity were mainly influenced by the catchment land use, agriculture enhanced density but reduced diversity, and the presence of riparian vegetation supported greater streambed bacterial diversity. From this perspective, a combination of prolonged dryness with reduced riparian vegetation and increased agricultural land cover could compromise the ecosystem functioning by threaten microbially mediated processes linked to the carbon cycle.

INTRODUCTION

Dry streambed sediments are dynamic ecotones where terrestrial and aquatic ecosystems interact, their functioning is largely driven by microorganisms hosted in the sediment and shaped by temporal environmental fluctuations [START_REF] Datry | Intermittent Rivers and Ephemeral Streams -Ecology and Manage-ment[END_REF][START_REF] Arce | A conceptual frame-work for understanding the biogeochemistry of dry river-beds through the lens of soil science[END_REF].

Currently, dry streambeds are expanding worldwide due to increased flow intermittency.

Reduced annual precipitation, increased water extraction and the global temperature rise are the major factors intensifying flow intermittency, especially in Mediterranean areas [START_REF] Cipriani | Impact of climate change on aquatic eco-systems along the Asse River network[END_REF][START_REF] Schewe | Multimodel assessment of water scarcity under climate change[END_REF]. Most Mediterranean intermittent streams, such as those in the northeastern Iberian Peninsula, usually cease flowing for between 1 and 2 months each year in the summer period. It has been reported that, currently, intermittent streams may exhibit extended dry-phase episodes also coupled with intense and sporadic precipitation [START_REF] Duran | Climate change in Catalonia: executive summary of the third report on cli-mate change in Catalonia Institute of Catalan Studies and the Government of Catalonia[END_REF]. The increase of streambed aridity in intermittent streams and the consecutive reduction of sediment-water content due to extended exposure to the air may affect the hosted microbiota and compromise microbial metabolism [START_REF] Rees | Variability in sediment microbial communities in a semipermanent stream: impact of drought[END_REF][START_REF] Marxsen | Microbial communities in streambed sediments recovering from desiccation[END_REF][START_REF] Febria | Microbial responses to changes in flow status in temporary headwater streams: a cross-system comparison[END_REF][START_REF] Gionchetta | Key role of streambed moisture and flash storms for microbial resistance and resilience to long-term drought[END_REF]. Microbial streambed communities are responsible for key ecosystem processes such as nutrient cycling and biogeochemical processes, and their metabolic capabilities are fundamental for maintaining fluvial ecosystem functioning [START_REF] Datry | Intermittent Rivers and Ephemeral Streams -Ecology and Manage-ment[END_REF]. Studies report that a long period of water stress reduces specific microbial activities linked to carbon degradation, probably due to: (i) the transition of microbes to a dormant stage [START_REF] Schimel | Microbial stress-response physiology and its implications for eco-system function[END_REF][START_REF] Jones | Dormancy contributes to the maintenance of microbial diversity[END_REF];

(ii) the allocation of a large fraction of cell metabolic costs for developing survival strategies to cope with desiccation (such as the production of extracellular polymeric substances) [START_REF] Lake | Ecological effects of perturbation by drought in flowing waters[END_REF][START_REF] Frossard | Disconnect of microbial structure and function: enzyme activities and bacterial communities in nascent stream cor-ridors[END_REF][START_REF] Romaní | The biota of intermittent rivers and ephemeral streams: prokaryotes, fungi and protozoans[END_REF]; and/or (iii) reduced labile organic matter available for microbes during droughts and increased accumulation of recalcitrant compounds on the streambed surface [START_REF] Ylla | Organic matter availability during pre-and post-drought periods in a Mediterranean stream[END_REF][START_REF] Freixa | Sediment microbial communities rely on different dissolved organic matter sources along a Mediterranean river continuum[END_REF]. On the other hand, and in a similar way to arid soil systems, frequent or sporadic precipitations interrupting dry periods can cause sharp exceptional increases in streambed sediment respiration activity, as a result of the awakening of dormant microbes [START_REF] Bérard | Resilience of soil microbial communities impacted by severe drought and high temperature in the context of Mediterranean heat waves[END_REF]Evans et al., 2014;[START_REF] Meisner | Partial drying accelerates bacterial growth recovery to rewetting[END_REF][START_REF] Marcé | Emissions from dry inland waters are a blind spot in the global carbon cycle[END_REF]. Streambed microbial assemblages have been shown to be affected by drying-rewetting cycles, especially when these include unusually extended dry phases [START_REF] Pohlon | Desiccation affects bacterial community structure and function in temperate stream sediments[END_REF][START_REF] Meisner | Drought legacy effects on the composition of soil fungal and prokaryote communi-ties[END_REF]. Their taxonomic composition can shift towards a community able to cope with more extreme environmental conditions, resulting in a more specialized community able to develop specific life strategies for coping with drought [START_REF] Jones | Dormancy contributes to the maintenance of microbial diversity[END_REF][START_REF] Shade | Fundamentals of microbial community resistance and resilience[END_REF]Evans and Wallenstein, 2014). In addition to hydrology, streambed microbes in intermittent systems may be further affected by cooccurring environmental factors related to the characteristics of the catchment acting at different spatial scales. Physical characteristics, such as the type of streambed substrata (e.g. cobbles/bedrock, fine/coarse sediment, accumulation of leaf litter), can influence the microbial response to hydrological intermittency [START_REF] Duarte | Responses of microbial decomposers to drought in streams may depend on the environmental context[END_REF][START_REF] Arce | A conceptual frame-work for understanding the biogeochemistry of dry river-beds through the lens of soil science[END_REF].

Headwaters are usually more heavily affected by natural intermittency than downstream areas, due to their hydro-morphological features (e.g. pronounced slope, distinct geomorphology reducing the hyporheic zone); also, these areas have been described as microbial biodiversity reservoirs, reflecting the high variability of habitats and resources (varying quantity and quality of organic matter), and/or higher levels of microbial colonization from the adjacent soils than that observed downstream [START_REF] Besemer | Headwaters are critical reservoirs of microbial diversity for fluvial networks[END_REF]. In headwater sections, a greater accumulation of leaf litter combined with a reduced percentage of sand and gravel sediment particles could lessen the loss of organic matter and water content during dry periods, limiting the reduction of microbial activity [START_REF] Freixa | Sediment microbial communities rely on different dissolved organic matter sources along a Mediterranean river continuum[END_REF]. Other studies have reported that leaf packs buried in the sediment maintain high humidity under extreme dry conditions, becoming temporal shelters for microbial survival [START_REF] Gionchetta | Key role of streambed moisture and flash storms for microbial resistance and resilience to long-term drought[END_REF]. Similarly, the presence of riparian vegetation, able to curtail water loss from sediments under dry conditions, has been considered a hydrological refuge for microbiota inhabiting the streambed [START_REF] Mclaughlin | Hydrologic refugia, plants, and climate change[END_REF]. At a larger scale, catchment land uses can also influence the stream and streambed biota and microbiota responses, especially through factors relating to the use of the surrounding soils (e.g. agriculture or industrial activities), or a reduction in the natural riparian vegetation [START_REF] Allan | Landscapes and Riverscapes: the influ-ence of land use on stream ecosystems[END_REF][START_REF] Steward | When the river runs dry: human and ecological values of dry riverbeds[END_REF][START_REF] Steward | Assessing the ecological health of rivers when they are dry[END_REF][START_REF] Bruno | Impacts of environ-mental filters on functional redundancy in riparian vegetation[END_REF]. To date, most studies related to the consequences of intermittency on streambed microbes have been mainly focused on the effects provoked by wet (and rewetting) periods [START_REF] Timoner | Biofilm functional responses to the rehydration of a dry intermittent stream[END_REF]. Less attention has been paid to the drought period, with only a few field studies focused on a single or limited number of stream ecosystems [START_REF] Timoner | Functional responses of stream biofilms to flow cessation, desiccation and rewetting[END_REF][START_REF] Zeglin | Stream microbial diversity in response to environmental changes: review and synthesis of existing research[END_REF][START_REF] Freixa | Sediment microbial communities rely on different dissolved organic matter sources along a Mediterranean river continuum[END_REF].

However, to better understand streambed microbial responses to intermittency, field studies embracing the multiple pressures (including the hydrological history and other environmental factors) are needed, as shown by recent studies demonstrating the interactive and joint effect of multiple factors on biological communities in natural ecosystems [START_REF] Besemer | Headwaters are critical reservoirs of microbial diversity for fluvial networks[END_REF][START_REF] Dickie | Soil microbial community structure explains the resistance of respiration to a dry-rewet cycle, but not soil functioning under static conditions[END_REF][START_REF] Bruno | Impacts of environ-mental filters on functional redundancy in riparian vegetation[END_REF][START_REF] Gieswein | Additive effects prevail: the response of biota to multiple stressors in an intensively monitored watershed[END_REF]. Accordingly, this work has attempted to fill this gap by evaluating the functional and structural responses of streambed microbial communities to a natural dryness gradient, by studying intermittent streams located in the Mediterranean basin. As the duration of dry and wet periods is highly unpredictable, the hydrology of each of the 37 selected sites was monitored in situ and later classified based on the natural gradient of their dryness history. From the 8 months of hydrology recorded, several parameters were calculated, including the duration of the total dry phase, as well as the length of the last wet event before our sampling. In parallel, a range of environmental factors

were measured and included as potential predictors for the responses of the streambed microbial community (e.g. riparian vegetation, organic matter content and land use).

Specifically, this investigation aimed to (i) elucidate whether, and to what extent, the hydrological history has influenced the microbial diversity, composition and activity of sediment microbial communities in the selected Mediterranean intermittent streams; (ii) determine the relative contribution of catchment and sediment characteristics (riparian vegetation, organic matter and water content, and land use) to modulating the microbial diversity, composition and activity according to the various hydrological factors; and (iii) define patterns of, and linkages between, the composition of the bacterial community and its metabolic activities along the dryness gradient. At sites with long dry-phases, we expected decreased microbial labile carbon degradation capacities but increased recalcitrant material degradation capacities, together with a transition of the microbial community composition towards a soil-like assemblage. We hypothesized to observe clearer microbial functional responses to the changes of sediment-water content compared with those related to the microbial composition and density.

RESULTS

Site classification and catchment characteristics

The principal component analysis (PCA) clearly distinguished the hydrological and catchment characteristics along the first and second axes respectively (Fig. 1). Factors with correlation coefficients >0.80 (Pearson correlation, rp) with the first two PCA components (PC1 and PC2) were selected to define gradients between sites. The first axis accounted for 40% of the total variability, describing the increase in both the total dry duration (T_DRY, rp: 0.941, Fig. 1) and the percentage of dry phase duration (P_DRY, rp: 0.941, Fig. 1), as well as the concurrent decrease (in the opposite PC1direction) in the number of consecutive wet days prior to the sampling (REW, rp: -0.946, Fig. 1). The variability along the second axis, explaining 17.3% of the total variation, mainly corresponded to the increase in forested land cover (FORE, rp: 0.861, Fig. 1) and the decrease (in the opposite PC2direction of the) in agricultural land-use (AGRI, rp: -0.888, Fig. 1). The site distribution was clearly ordinated along PC1, which identified four hydrological history categories (FL, flow sites; SD, short-dry phase sites; MD, medium-dry phase sites; LD, long-dry phase sites see Experimental procedures section) based on the duration of their dry phase. The PC2 axis ordered the sites according to their forested or agricultural land cover (Table 1). Hydrology and catchment characteristics were the main drivers shaping the site distribution, where LD-phase sites had significantly lower riparian vegetation (RIP_VEG) and higher organic matter content (ORG_MAT) compared with FL sites (p < 0.05, Table 1). All the sites showed similar percentages of sediment water content since the sampling was conducted over the first days of the wet season (WAT_CON, Table 1). As expected, LD sites showed the longest dry period (L_DRY) prior to the sampling and the shortest rewetting periods (REW) in comparison to the other sites (p < 0.05, Table 1). A similar frequency of dry/rewetting cycles was observed between the sites, and only MD sites presented significantly higher frequencies than those of FL sites (p < 0.05, Table 1).

Influence of hydrological history on streambed microbial functions

The streambed extracellular enzymes and community respiration were not significantly different between the four hydrological categories (Fig. S2), but their variability was particularly high for SD sites, in the case of hydrolytic activities (β-glucosidase, GLU and βxylosidase, XYL), and for MD and LD sites, in the case of sediment respiration (Fig. S2).

Both variance partitioning (Fig. S3) and averaged models (Table 2) indicated that hydrological variables were important for phenoloxidase (PHE) and community respiration (RESP), while catchment characteristics explained more of GLU and XYL variability (Fig. S3, Table 2). Specifically, the total duration of the dry phase explained the largest variance in the case of PHE activity (PHE, VAR. PAR. T_DRY = 13%, Fig. S3), while the interaction between the total dry duration (T_DRY) and the length of the last wet period (REW) influenced the variation in community respiration (RESP, VAR. PAR. T_DRY*REW = 13%, Fig. S3). In contrast, the organic matter content and riparian vegetation cover better explained the variance of XYL (VAR. PAR. ORG_MAT = 22%, RIP_VEG = 25%) while both organic matter and water content of sediments were relevant for GLU variance explanation (VAR.

PAR. ORG_MAT = 17%, WAT_CON = 15%, Fig. S3). Microbial function plots showed the relationship (significant or not) between each functional variable and the total dry duration as well as the effect (e.g. additive, synergistic or antagonist) of the duration of the last wet period (Fig. 2). According to the model plots, despite the high variability seen between the intermittent sites, the hydrolytic GLU and XYL activities decreased slightly when T_DRY increased, and this trend was significant in the case of XYL activity (Fig. 2A andB, Table 2).

In contrast, PHE activity increased significantly with longer dry periods, with mean PHE activities doubling at LD sites compared with those at SD sites (126.63 13.64 and 251.13 55.19 nmol DIQC gDW-1 h, at SD and LD sites respectively), while no interaction with REW was found (Fig. 2C, Fig. S2 and Table 2). With regard to catchment characteristics, sediment organic matter content was positively related to GLU and XYL hydrolytic activities (Table 2), while the riparian vegetation cover was negatively related to XYL activity (Table 2). Community respiration presented a more complex response to the hydrological conditions, with peaks occurring with longer rewetting periods in largely intermittent and extremely dry sediments (MD and LD sites in Fig. 2D, Table 2).

Streambed microbial diversity responses to hydrological history and relationship to

function Remarkably, the sediment microbial community diversity metrics [e.g. Shannon-Wiener diversity (DIV), richness (RICH)], and bacterial density (BAC), showed stronger relationships with catchment characteristics than hydrology (Table 2, Fig. S3). Riparian vegetation cover explained the largest part of the variance in the case of taxonomic diversity (DIV, VAR. PAR. RIP_VEG = 28%, Fig. S3) and taxonomic richness (RICH, VAR. PAR. RIP_VEG = 33%, Fig. S3), but only 13% for bacterial density within the sediment communities (Table 2 and Fig. S3). A small fraction of the variability was explained by hydrological parameters, such as T_DRY, REW and T_DRYxREW, accounting for 5% of the diversity, 9% of the richness and 5% of the bacterial density (Fig. S3). Significant positive relationships were found between both the diversity and richness indices (DIV, RICH) and the riparian vegetation coverage (RIP_VEG, Fig. 3, Table 2). Specifically, linear models showed that abundant riparian vegetation determines a positive and significant increase in both the microbial diversity and richness indices, particularly at MD sites (Fig. 3A andB; Table 2).

With regard to land-use, agricultural land cover explained the largest part of the bacterial density variance (BAC, VAR. PAR. AGRI = 33%) and it was also positively and significantly related to the bacterial density, but negatively related to community diversity (Table 2). The microbial composition of the sediments changed in response to the duration of the dry phase; sites with continuous flow had a significantly different community composition to that at intermittent sites [permutation analysis of variance (PERMANOVA), Fig. 4A]. This result was confirmed by the relation existing between NMDS1 and PC1 (Fig. 4B). In terms of taxonomy, the phylum Proteobacteria dominated the streambed communities at all sites, although its relative abundance decreased at intermittent sites (Fig. S4). Specifically, there was a significant reduction in the classes Delta-proteobacteria and Epsilon-proteobacteria between FL and LD sites [analysis of variance (ANOVA) p < 0.05; Fig. S4, Fig. 4A].

Similarly, the relative abundance of the phylum Planctomycetes decreased along the dry phase gradient, with a significant decrease in the relative abundance of the class Planctomycetia at MD and LD sites compared with FL sites (ANOVA p < 0.05; Fig. S4, Fig. 4A). In contrast, the relative abundance of the phyla Actinobacteria and Bacteroidetes doubled at sites with the longest dry phase (Fig. S4). More specifically, of the classes Actinobacteria (rs = 0.76), Thermoleophilia (rs = 0.56) and Sphingobacteria (rs = 0.65) correlated strongly with the MD and LD communities (Fig. 4A); the relative abundance of Actinobacteria increased significantly at MD and LD sites (ANOVA, p < 0.05; Fig. S4, Fig. 4A), and Sphingobacteria at LD sites, compared with FL sites (ANOVA, p < 0.05; Fig. S4, Fig. 4A). The global correlation between the community composition and the functions was not significant (p = 0.182, rs = 0.10), suggesting that these are uncoupled (Supplementary Information S4).

DISCUSSION

Overall, our findings confirm the hypothesis that microbial functions respond more strongly to hydrological factors than bacterial density and diversity do. The microbial community composition differs between flowing and intermittent sites whereas bacterial diversity and density respond to the environmental features of the catchment area, such as riparian vegetation cover and agricultural land use. Generally, reduced microbial activity, diversity changes and shifts in community composition are expected as the dry period lengthens [START_REF] Rees | Variability in sediment microbial communities in a semipermanent stream: impact of drought[END_REF][START_REF] Amalfitano | Recovery and quantification of bacterial cells associated with streambed sedi-ments[END_REF][START_REF] Marxsen | Microbial communities in streambed sediments recovering from desiccation[END_REF], although co-occurring environmental factors are also often cited as important modulators for microbial responses to drought [START_REF] Drenovsky | Land use and climatic factors structure regional patterns in soil microbial communities[END_REF]RuizGonzález et al., 2018). As previously reported, streambed sediments are completely exposed to the atmosphere under persistent dry conditions and may acquire characteristics similar to those of the adjacent soils [START_REF] Steward | When the river runs dry: human and ecological values of dry riverbeds[END_REF][START_REF] Morandi | How is success or failure in river restoration pro-jects evaluated? Feedback from French restoration projects[END_REF][START_REF] Arce | A conceptual frame-work for understanding the biogeochemistry of dry river-beds through the lens of soil science[END_REF]. A longer dry phase and modifications in the frequency of the dry-wet cycle in river ecosystems often determine the degree of similarity between streambed and terrestrial habitats [START_REF] Harms | Responses of trace gases to hydrologic pulses in desert floodplains[END_REF][START_REF] Mori | Spatiotemporal heterogeneity of actual and potential respiration in two contrasting floodplains[END_REF].

Drying gradient effects on streambed microbial function

Similar to previous studies, we found that the presence of specific extracellular enzyme activities (EEA) is influenced by the extended duration of the dry phase [START_REF] Zoppini | Organic priority substances and microbial processes in river sediments subject to contrasting hydrological conditions[END_REF][START_REF] Freixa | Sediment microbial communities rely on different dissolved organic matter sources along a Mediterranean river continuum[END_REF][START_REF] Gionchetta | Key role of streambed moisture and flash storms for microbial resistance and resilience to long-term drought[END_REF]. With increasing dryness, PHE activity is enhanced indicating a greater capability to degrade lignin compounds, while the tendency of β-glucosidase and β-xylosidase to decrease suggested a reduced capability to degrade simple polysaccharides. This shift in enzyme capabilities suggests a gradual modification of the quality of organic matter available in streambed sediments when desiccating (Figs. 2 and5, PHE, GLU and XYL). According to previous studies, reduced flow connectivity, accumulation of old plant litter and the development of terrestrial vegetation on the streambed [START_REF] Battin | Contributions of microbial biofilms to ecosystem processes in stream mesocosms[END_REF][START_REF] Fazi | Stream hydrological fragmentation drives bacterioplankton community composition[END_REF][START_REF] Romaní | Microbial biofilm structure and organic matter use in mediterranean streams[END_REF] may enrich the sediment with lignin-like compounds. At the same time, desiccation may reduce the availability of simple polysaccharides due to the decreased input of freshly available organic compounds, such as those produced by algae in wet conditions [START_REF] Acuña | Meteorological and riparian influences on organic matter dynamics in a forested Mediterranean stream[END_REF][START_REF] Romaní | Microbial biofilm structure and organic matter use in mediterranean streams[END_REF][START_REF] Freixa | Sediment microbial communities rely on different dissolved organic matter sources along a Mediterranean river continuum[END_REF]. Surprisingly, β-xylosidase activity was shown to be negatively related to riparian vegetation, even though this enzyme is usually involved in the degradation of hemicellulose from plants [START_REF] Romaní | Microbial biofilm structure and organic matter use in mediterranean streams[END_REF][START_REF] Romaní | Interactions of Bacteria and Fungi on Decomposing Litter: Differential Extracellular Enzyme Activities[END_REF]. The large variability in enzyme activity observed between the intermittent sites may indicate a high level of plasticity in the microbial community, enabling them to select the source of organic matter. The overall expression of microbial enzyme activity is modulated by both the availability and the type of resources. For instance, as simple polysaccharides from algal activity and allochthonous plant materials from riparian vegetation are available, the production of the β-glucosidase enzyme would be favoured rather than β-xylosidase, since the end-product, glucose, is more energetically important than xylose. Nevertheless, the activity of β-xylosidase could be favoured when low riparian cover and scarce readily degradable algal compounds are present (Fig. 5). These results suggest that the coupled quantity and quality of organic matter influences microbial capacity to use organic materials, and they reveal adaptation to the sources available resulting from changes in intermittency and co-occurring factors, like riparian vegetation cover (Fig. 5). Long dry-phase sites are clearly characterized by the use of more terrestrial-origin materials (shown, in particular, by the increased PHE activity) despite presenting much less riparian vegetation (Table 2, Fig. 5). Besides the effect of the duration of the dry period, the occurring of wet episodes may enhance the microbial activity. This was observed in the case of respiration, which was particularly enhanced when prolonged rewetting occurred at hardly dry sites. In addition, the activity of the two hydrolytic enzymes varied greatly between the sites with short and intermediate length dry periods (SD and MD, Fig. S2). The large variability in these activities suggests that greater intermittency (i.e. alternation of dry and wet episodes) may determine larger fluctuations in microbial activity, which in turn could mask the expected decrease in β-glucosidase activity at the driest sites. Indeed, wet events make simpler and more labile compounds available, mobilizing these into the sediment [START_REF] Birch | Microbiological changes in freshly moistened soil[END_REF][START_REF] Fierer | A proposed mechanism for the pulse in carbon dioxide production commonly observed following the rapid rewetting of a dry soil[END_REF][START_REF] Rees | Variability in sediment microbial communities in a semipermanent stream: impact of drought[END_REF]. In this study, respiration was not measured directly after each wet event, but the highest values of βglucosidase activity were observed in the frequently rewetted streambeds (MD sites), suggesting that the short wet events may reactivate a large number of dormant microbes [START_REF] Timoner | Functional responses of stream biofilms to flow cessation, desiccation and rewetting[END_REF] or still-active oxidative enzymes from damaged or dead cells [START_REF] Burns | Soil enzymes in a changing environment: current knowl-edge and future directions[END_REF][START_REF] Meisner | Prolonged drought changes the bacterial growth response to rewetting[END_REF].

Drying gradient effects on streambed community composition

The microbial community composition reflects clear differences in the relative abundance of certain taxa inhabiting the flowing and intermittent sites, suggesting that the presence of no flow episodes exerts an important effect, irrespective of the degree of intermittency. Microbes are capable of adapting rapidly to changing environmental conditions through structural and functional changes [START_REF] Allison | Resistance, resil-ience, and redundancy in microbial communities[END_REF][START_REF] Wallenstein | A trait-based frame-work for predicting when and where microbial adaptation to climate change will affect ecosystem functioning[END_REF]. Consequently, the stress provoked by a single dry period may determine the changes in intermittent microbial communities through the rapid osmotic adaptation of cells to cope with alternating dry and wet periods [START_REF] Shade | Fundamentals of microbial community resistance and resilience[END_REF][START_REF] Attermeyer | Ecosystem-level studies of terrestrial carbon reveal contrasting bacterial metabolism in different aquatic habitats[END_REF][START_REF] Meisner | Drought legacy effects on the composition of soil fungal and prokaryote communi-ties[END_REF]. As expected, in our study, increased dryness tended to modify the microbial communities' composition with respect to those of the flowing sites, involving higher relative abundances of drought-tolerant taxa and fewer drought sensitive taxa. This microbial community transition is characterized by a decrease in Gram-negative bacteria compared with an increase in the proportion of Grampositive bacteria, which are more resistant and adapted to live in soil systems submitted to extreme drought (de Vries and [START_REF] De Vries | Controls on soil micro-bial community stability under climate change[END_REF][START_REF] Acosta-Martínez | Pre-dominant bacterial and fungal assemblages in agricultural soils during a record drought/heat wave and linkages to enzyme activities of biogeochemical cycling[END_REF][START_REF] Naylor | The role of drought on root-associated bacterial communities across diverse cereal grass species and over a developmental gradient[END_REF]. Furthermore, the microbial composition changes were driven by either just one or a few members of a microbial phylum, as also observed by Barnard and colleagues (2013) in soil systems. Specifically, the assemblage shifts in these streambeds included a significant decrease in the Gram-negative classes of Epsilon-proteobacteria and Deltaproteobacteria (phylum Proteobacteria) and Planctomycetia (phylum Planctomycetes).

Conversely, we found significant increases in the Gram-positive classes Actinobacteria and Thermoleophilia (phylum Actinobacteria) as observed in extremely desiccated soils [START_REF] Bouskill | Pre-exposure to drought increases the resistance of tropical forest soil bacterial com-munities to extended drought[END_REF]. The spread of the class Sphingobacteria (phylum Bacteroidetes) in the driest intermittent communities could be due to rewetting effects (last rewetting event prior to sampling). Sharp increases in the relative abundance of the class of Sphingobacteria have often been associated with abrupt wet events in arid and desert soil communities [START_REF] Aslam | Soil compartment is a major determinant of the impact of simulated rainfall on desert microbiota[END_REF][START_REF] Meisner | Drought legacy effects on the composition of soil fungal and prokaryote communi-ties[END_REF]. In contrast to community composition, richness and diversity were not significantly modified by hydrological history, suggesting that there was taxadominance alternation or substitution without increased dryness significantly affecting the overall diversity. This suggests that desiccation may also play a role in decreasing diversity and richness with a concomitant decrease in riparian vegetation cover.

Relevance of environmental catchment characteristics

Riparian vegetation and land use (agricultural cover) were significant drivers of streambed microbial community richness and diversity. While riparian vegetation positively influenced the diversity of the microbial community, agricultural cover negatively affected bacterial diversity but enhanced bacterial density. Riparian vegetation may help conserve stream biodiversity and ecosystem functioning [START_REF] Elosegi | Effects of hydro-morphological integrity on biodiversity and functioning of river ecosystems[END_REF]. The greater bacterial diversity and richness associated with the greater riparian vegetation cover could be related to more microhabitats available for microbes (e.g. through leaf fall, pool persistence, colonization by terrestrial plants) [START_REF] Romaní | The biota of intermittent rivers and ephemeral streams: prokaryotes, fungi and protozoans[END_REF][START_REF] Arce | A conceptual frame-work for understanding the biogeochemistry of dry river-beds through the lens of soil science[END_REF]. On the other hand, although we found no direct effect of hydrological history on microbial diversity in our study, reduced riparian vegetation co-occurred with longer dry periods. With regard to land use, agricultural cover was the major force driving reduced streambed bacterial diversity while at the same time increasing bacterial density. Even though most of the sites studied presented similar catchment characteristics, some had important agricultural field cover in their surroundings.

These microbial community responses might be a consequence of enhanced nutrient input from agricultural practices in the fields around the sites; this may induce bacterial growth and favour the spread of certain dominant species, thus reducing diversity, as expected in eutrophic conditions or nutrient-contaminated systems [START_REF] Zeglin | Stream microbial diversity in response to environmental changes: review and synthesis of existing research[END_REF][START_REF] Li | Factors affecting soil microbial biomass and func-tional diversity with the application of organic amendments in three contrasting cropland soils during a field experi-ment[END_REF]. The negative effect of agricultural cover on bacterial diversity may be further enhanced by longer droughts (negative interaction T-DRYxAGRI, Table 2), as also suggested in previous studies on leaf litter decomposition and soils [START_REF] Duarte | Responses of microbial decomposers to drought in streams may depend on the environmental context[END_REF][START_REF] Meisner | Drought legacy effects on the composition of soil fungal and prokaryote communi-ties[END_REF].

CONCLUDING REMARKS

Considering both microbial function and the response of streambed community composition, these findings highlight the fact that hydrological history has a significant effect although functional, composition and diversity responses did not occur in parallel (absence of correlation in Mantel test). Interestingly, the response of function, in terms of organic matter use capability (i.e. extracellular enzymes) presented a gradient with increased use of PHE (related to lignin compounds) and decreased use of β-glucosidase and β-xylosidase (linked to simple polysaccharides degradation) in streambeds suffering longer dry periods, as well as a high degree of variability at intermittent sites. In contrast, streambed community composition differed between sites with flowing waters and intermittent conditions, irrespective of the drying intensity, while richness and diversity were strongly modulated by the characteristics of the catchment area. Contradictory results are also reported in the literature, where some studies present coupled compositional, diversity and functional microbial responses to water stress, while others show activity responses to desiccation decoupled from changes in microbial community diversity and/or composition [START_REF] Frossard | Disconnect of microbial structure and function: enzyme activities and bacterial communities in nascent stream cor-ridors[END_REF][START_REF] Pohlon | Desiccation affects bacterial community structure and function in temperate stream sediments[END_REF][START_REF] Febria | Microbial responses to changes in flow status in temporary headwater streams: a cross-system comparison[END_REF][START_REF] Zoppini | Microbial responses to polycyclic aromatic hydrocarbon contamination in temporary river sediments: experimental insights[END_REF]. The intrinsically high variability of a field assessment could explain the unexpectedly weak correlation observed between the microbial community composition and functions, bearing in mind the fact that the activities and diversity measured in this work do not exhaustively explain the entire system. Overall, the maintenance of microbial diversity under extreme dry streambed conditions potentially ensures reduced losses in system functioning. Nevertheless, further studies are needed to confirm whether longer dry periods combined with reduced riparian vegetation and increased agricultural land cover could threaten microbially mediated processes linked to the carbon cycle.

EXPERIMENTAL PROCEDURES

Study area, site selection and sampling strategy

A total of 37 upstream sites were selected and surveyed in streams situated in the northeastern Iberian Peninsula (Catalonia). These sites belong to nine watersheds corresponding to the main rivers in that area, flowing into the Mediterranean Sea, i.e. the Muga, Fluvià, Ter, Tordera, Llobregat, Francolí, Foix, Besós and Ebro rivers (Fig. S1). The sampling sites include lowland and midmountain watercourses (altitude 80% totally shaded, considered length 100 m) and riparian vegetation cover (e.g. the presence of trees, shrubs or small bushes was considered as a percentage, from 75%, for the first 10 m width of riparian vegetation) were measured in situ following a fieldassessment procedure. Once at the laboratory, sediment subsamples of about 1 ml volume (1 cm2 projected surface area) were created and collected from each sediment corer, using an uncapped syringe, for EEA, respiration, bacterial community composition and density analyses. Community respiration and EEA were measured within 24 h of sampling. Bacterial density samples were fixed with 10 ml of filtersterilized (0.2 μm pore size nylon filters Whatman, Kent, UK) synthetic water (13 mg L-1 Na2SO4, 16.1 mg L-1 Na2SiO3, 29.4 mg L-1 CaCl2*2H2O, 0.6 mg L-1 KCl, 3 mg L-1 MgSO4*7H2O, 26.5 mg L-1 Na2CO3, 0.6 mg L-1 NH4H2PO4 and 7.3 mg L-1 (NH4) (NO3) in MQwater) and formaldehyde (2% final concentration), and then stored at room temperature. For the microbial diversity analysis of the 16S rRNA gene, the sediments were stored at -20C. The remaining sediment was used to determine organic matter and water content.

Hydrological history and catchment characteristics of the study sites

One year prior to the sampling, temperature and pressure transducers (Solinst Levelogger Gold Model 3001 and Solinst Barologger Gold Model 3001,Solinst,Georgetown,ON,Canada), and temperature data loggers (ACR SmartButton Logger, MicroDAQ) were installed at the study sites to characterize their hydrology. All the temperature loggers recorded data every 60 min. Once the entire dataset for each sampling site had been obtained, the daily variations in the streambed and air temperature were estimated. To characterize the stream hydrology in the intermittent stream, we used the daily variation of the streambed temperature corrected for the barometric pressure and air temperature. This daily variation was determined as the difference between the maximum and minimum temperature per day and the daily higher rate of change per hour (Colls et al., 2019) (see Supplementary Information S1, Arias-Real et al., 2018). In addition, the physical and chemical characteristics of the sites were recorded (Supplementary Information S2). In total, the hydrological data collected covered an 8 month period (245 days) of hydrology history for each stream site, and the related parameters calculated from this database were: the total number of dry days in the 8-month period (T_DRY); the percentage of dry days in the 8-month period (P_DRY); the duration (days) of the last dry event prior to the sampling (L_DRY); the duration of the last rewetting period (i.e. the number of consecutive wet days prior to the sampling) (REW) and the number of dry/rewetting cycles over the 8-month period (FREQ). As well as hydrological descriptors, in situ data corresponding to sediment organic matter, water content, shade cover and riparian vegetation were collected, and land use in the catchment area of the sites was determined. The organic matter content was measured as AFDW (Ash-Free Dry Weight) and 1 ml of sediment sample was dried at 70C for 72 h and then burnt for 4 h at 450C using a muffle furnace (AAF 1100, Carbolite, UK). The results were expressed as grams of AFDW per gram of sediment dry weight (DW). In order to estimate the sediment saturation, sediment water content was calculated as the percentage of water loss (% difference) between the fresh and DW (1 ml for each sample). Major land covers (urban zones, agricultural and forested areas, Corine Land Cover 2012, https://www.eea.europa. eu/data-and-maps/data/clc-2012-raster) were calculated as percentage considering the adjacent 1 km upstream of the sampling point, within the catchment area.

Microbial response variables

Bacterial density and community respiration. Bacterial density in sediments was determined through flow cytometry (FACSCalibur, Becton Dickinson), following a protocol adapted from [START_REF] Amalfitano | Recovery and quantification of bacterial cells associated with streambed sedi-ments[END_REF]. Sediment samples were sonicated for 1 min, shaken for 30 s, and sonicated again for 1 min to extract the biofilm from the sediment grains (Ultrasons, Selecta).

A subsample of the obtained extract (1 ml) was pipetted into a glass vial and 9 ml of detaching solution (distilled water, NaCl at 0.85% final concentration, Tween 20 at 0.5% final concentration and sodium pyrophosphate at 0.1 M, final concentration) was added to promote the separation of the cells. The samples were shaken for 30 min (150 rpm, in the dark and at room temperature), left for 10 min at 4C, and then sonicated with ice for two 1 min cycles.

Next, the samples were shaken again for 1 min and left for 5 min, to allow the sedimentation of the larger particles, and, finally, 1 ml of supernatant was transferred into a clean sterile Eppendorf tube. Once the sediment extracts had been obtained, each sample was purified through Nycodenz (1 ml), added to the bottom of each Eppendorf tube, and centrifuged (14 000 rpm) for 90 min at 4C. The purified sediment extracts (400 μl) were stained with Syto13 (4 μl Fisher, 5 μM solution) and incubated in the dark for 30 min. The stained samples were counted using flow cytometry (FACSCalibur, Becton Dickinson). To normalize fluorescence data, a known concentration of bead solution (10 μl of 106 beads ml-1 , Fisher 1.0 μm) was added to the samples. The bacterial density (BAC) results were expressed as cells × 108 per gram of sediment DW. Resazurin assay was applied to determine community respiration (RESP), using a modified version of the technique in Haggerty and colleagues (2010). The resazurin (7-Hydroxy-3H-phenoxazin3-one-10-oxide sodium salt, RAZ) method consists of using a redox dye (blue in the oxidized state) that indicates the respiratory activity of microorganisms (it turns pink when reduced). Before applying the resazurin dye, we prepared the sediment extract for each sample replicate from 1 ml of sediment, placed in a sterile vial with 5 ml of pyrophosphate solution 50 mM [START_REF] Manini | Synoptic determination of living/dead and active/dormant bacterial fractions in marine sediments[END_REF]. Each sample was sonicated for 1 min in ice in order to preserve the cell structure [START_REF] Amalfitano | Recovery and quantification of bacterial cells associated with streambed sedi-ments[END_REF], shaken for 30 s, and sonicated again for 1 min to extract the biofilm from the sediment grains (Ultrasons, Selecta). One milliliter subsamples of the extract was diluted in 9 ml of filtered water (0.2 μm pore size nylon filters Whatman, Kent, UK). Aliquots of each sediment extract (166 μl) were placed in a transparent microplate (MicroWell™ Nunc™, ThermoFisher Scientific) and mixed with ringer solution (28 μl) and resazurin stock substrate (6 μl of 0.35 mM Sigma-Aldrich, USA). The plate was closed with a lid, then incubated (target temperature 20C), and each cycle started with shaking for 90 s (shaking amplitude 3 mm), with a pause of 10 s before each absorbance reading. The incubation period lasted 24 h and the absorbance was read every 4 min (603 and 570 nm resazurin and resorufin wavelengths). Subsequently, we used the standard resorufin (pink compound that indicates the conversion of resazurin to resorufin) curve and the absorbance values obtained after 18 h of incubation to define the community respiration as micromoles of resazurin respired per grams of DW h-1 (μmol RAZ gDW-1 h-1).

Extracellular enzyme activities

Fresh sediment subsamples (1 ml) were analysed for the activity of two hydrolytic enzymes, β-D-1, 4-glucosidase (EC 3.2.1.21, GLU), β-xylosidase (EC 3.2.1.37, XYL), as well as for the oxidative enzyme PHE (EC 1.14.18.1, PHE), as indicators of the cellulose, hemicellulose and lignin compound decomposition capacity respectively. For the EEA, we followed the procedure from Gionchetta and colleagues (2019). The hydrolytic enzyme (GLU and XYL) values were expressed as nmol MUF gDW-1 h-1, and the PHE enzyme values were expressed as nmol 2,3-dihydroindole-5,6-quinone-2-carboxilate (DIQC) gDW-1 h-1.

Community composition and diversity

DNA extraction was performed from 0.5 g of sediment for each composite sample (n = 5 for each stream site) using the FastDNA™Spin Kit for Soils (MP Biomedicals, Irvine, CA) according to the manufacturer's instructions. The DNA concentration and quality in each sample were measured spectrophotometrically (Nanodrop2000, Thermo Scientific™, Waltham, USA) and further assessed with PCR (95C × 5 min; 35 cycles: 95C × 1 min; 55C × 1 min; 72C × 1 min; 72C × 10 min) and 1% agarose gel electrophoresis. The DNA extracts were stored at -20 °C until sequencing analysis at the Roy J. Carver Biotechnology Center (University of Illinois, IL, USA). The diversity of prokaryotes (including bacteria and archaea) was analysed utilizing the primer pair V4_515F/V4_806R61 [START_REF] Caporaso | Global pat-terns of 16S rRNA diversity at a depth of millions of sequences per sample[END_REF], targeting the V4 region of the 16S rRNA gene, using Illumina MiSeq technology. We analysed both bacteria and archaea species in sediments, but the latter was excluded from the analysis due to extremely low read values. IM-TORNADO pipeline version 2.0.3.2 [START_REF] Jeraldo | IM-TORNADO: a tool for compari-son of 16S reads from paired-end libraries[END_REF] was used to generate the operational taxonomic unit (OTU) table and assign taxonomy to the corresponding OTUs (≥97% sequence similarity). Bacterial taxonomy was assigned through the Ribosomal Database Project (RDP version 10, [START_REF] Cole | Ribosomal database project: data and tools for high throughput rRNA analysis[END_REF].

Diversity analyses were performed using the IMTORNADO generated output BIOM table.

The range of sequencing depth was between 9668 and 127 452 of DNA for total prokaryotes. OTU tables were generated from the BIOM table using the phyloseq package for analysing microbiome data, by R software (R version 3.4.1). Rarefaction curves were generated (rarefy_even_depth function, phyloseq package) and the rarefied OTU table was set to the smallest count of bacterial reads (9211 reads/sample) in order to minimize bias due to different numbers of total sequences present in the samples. Two alpha-diversity metrics were calculated from the OTU table, i.e. the Shannon-Wiener index (DIV) and the Chao1 richness index (RICH), using the phyloseq package by R software (R version 3.4.1). The raw sequence dataset was deposited in the NCBI Sequence Read Archive database under Accession No.

PRJNA557375.

Data analyses

To summarize the hydrological history and catchment characteristics, and assess the existence of gradients between the study sites and key driving factors, we performed a PCA including all the parameters analysed (total number of dry days in the 8-month period, T_DRY; the duration (days) of the last dry event prior to the sampling, L_DRY; the number of consecutive wet days prior to the sampling, REW; frequency of dry/rewetting cycles in the 8-month period, FREQ; sediment water content, WAT_CON; organic matter content, ORG_MAT; riparian vegetation cover, RIP_VEG; shade cover, SHAD; forested land use, FORE; agricultural land use, AGRI; urban land use, URBA). Although the hydrological history was mainly studied as a continuous factor, in order to allow some general analyses, the streambed sites were classified according to the duration of their dry phase in a similar way to the proposal of [START_REF] Costigan | Flow regimes in inter-mittent rivers and ephemeral streams[END_REF], and four classes were defined for our study sites: flowing sites (FL, continuous flow); short-dry phase sites (SD, from 1% to 15% dry days); mediumdry phase sites (MD, 15%-85% dry days over the 8-month period); and longdry phase sites (LD, average 85%-100% dry days) (Fig. S1, Table 1). The differences between the hydrological categories (FL, SD, MD and LD sites) were analysed for EEA, community respiration, bacterial density and microbial diversity and richness indices, using one-way ANOVA. To analyse the response of microbial functions (extracellular enzymes and respiration), and diversity and bacterial density, to the hydrological and catchment characteristics, we followed the multiple-predictor analytical procedure using survey data (see [START_REF] Feld | Analysing the impact of multiple stressors in aquatic bio-monitoring data: a 'cookbook' with applications in[END_REF]. First, we ran exploratory analyses to reduce skewness when necessary: GLU (β-glucosidase activity) and XYL (β-xylosidase activity) were square-root transformed; PHE (PHE activity), RESP (community respiration), BAC (bacterial density), DIV (Shannon diversity index), RICH (Chao1 richness index) and RIP (riparian vegetation) were logtransformed). Consequently, we standardized and checked the collinearity between the variables [START_REF] Feld | Analysing the impact of multiple stressors in aquatic bio-monitoring data: a 'cookbook' with applications in[END_REF]. Second, using Random Forest (RF) analysis (package ForestSRC, R software) [START_REF] Ishwaran | Random survival for-ests for competing risks[END_REF] we ranked, identified and visualized all the predictor interactions for each response variable: hydrological history parameters (T_DRY, L_DRY, REW, FREQ); sediment water and organic matter (WAT_CON, ORG_MAT); shade cover (SHAD); riparian vegetation (RIP_VEG); and agricultural (AGRI), forested (FORE), and urban (URBA) land use. The number of trees in the RF was 3000 and the variables used in each split was set to three [START_REF] Feld | Analysing the impact of multiple stressors in aquatic bio-monitoring data: a 'cookbook' with applications in[END_REF]. Afterwards, we applied generalized linear models to model the response of each variable to single and combined predictors [START_REF] Feld | Analysing the impact of multiple stressors in aquatic bio-monitoring data: a 'cookbook' with applications in[END_REF][START_REF] Grueber | Multimodel inference in ecology and evolution: challenges and solutions[END_REF][START_REF] Barton | MuMIn: Multi-Model Inference[END_REF]. For each model, we used a maximum of four predictors, selected considering the RFs [START_REF] Feld | Analysing the impact of multiple stressors in aquatic bio-monitoring data: a 'cookbook' with applications in[END_REF], including two hydrological predictor terms (T_DRY and REW), and at least two environmental predictors for each response variable. Subsequently, the multi-model averaging approach (MuMIn R package) was applied. This procedure generated all possible models and ranked them based on model weight < 0.05 for all tests. Details on the analytical procedure are presented in Supplementary Information S3. The differences in microbial community composition between the study sites were analysed using non-metric multidimensional scaling (NMDS) ordination, based on Bray-Curtis dissimilarity performed at the OTU taxonomic level. The OTU relative abundance matrix used for the ordination analyses was previously fourth root transformed. A PERMANOVA (9999 permutations) was applied to test the significant differences in the community composition between the four hydrological history categories (PRIMER 6 version 6.1.11 and PERMANOVA+ version 1.0.1). In order to identify the taxonomic classes significantly related to the OTU distribution, specific vector fitting analysis was applied to the NMDS [START_REF] Legendre | Numerical Ecology[END_REF]. A further one-way ANOVA was applied to determine significant differences between the relative abundances of the taxonomic classes inhabiting the different hydrological history categories (FL, SD, MD and LD sites). Finally, a Mantel test was performed between community composition and the functional matrices to assess whether they responded similarly to flow intermittency (Bray-Curtis dissimilarity for the compositional matrix and Euclidian distance for the functional data, and then RELATE analysis with Spearman rank correlation in Supplementary Information S4, PRIMER 6 version 6.1.11 and PERMANOVA+ version 1.0.1). xylosidase; and PHE, phenoloxidase) and community respiration (RESP) to the hydrology (xaxis, total dry days, T_DRY) and the duration of the last wet period (REW) using linear models. The wet duration was fixed at three levels: the 5% (low), 50% (mean) and 95% (high) percentiles. The predictor selection and interaction for each model were based on the Random Forest analyses (Table 2). The colours of each site define the corresponding hydrological group (legend). [Color figure can be viewed at wileyonlinelibrary.com] Fig. 3. Response of the microbial taxonomic diversity (DIV, Shannon-Wiener index) and richness (RICH, Chao1 index) to the hydrology (x-axis, total dry days, T_DRY) and the abundance of riparian vegetation (RIP_VEG) using linear models. The abundance of riparian vegetation was fixed at three levels: the 5% (low), 50% (mean) and 95% (high) percentiles.

FIGURE LEGENDS

The predictor selection and interaction for each model were based on the Random Forest analyses (Table 2). The colours of each site define the corresponding hydrological group (legend).

[Color figure can be viewed at wileyonlinelibrary.com]

Fig. 1 .

 1 Fig. 1. Principal component analysis (PCA) showing the importance of the environmental and hydrological factors that best distinguished our 37 study sites, represented by arrows. The colours and shape of each site define the corresponding hydrological group (legend). [Color figure can be viewed at wileyonlinelibrary.com]

Fig. 4 .

 4 Fig. 4. A. NMDS ordination plot representing the bacterial community composition (OTUs

Table 2 .

 2 Results of the linear models created for each microbial functional and structural response variable. These models examined the effect of total dry days (T_DRY), the rewetting duration prior to the sampling (REW), their interaction (T_DRY*REW) and at least two environmental predictors. The predictor set for each model was selected considering the results of the Random Forest analyses. Predictors and interaction standardized effect size (SES), standard error (SE) and significance (p value) for each model are shown (significant values are indicated in bold and by: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1). The na interaction values indicate unavailable results that were automatically excluded from the model as their importance was too low. The relative predictor importance of the model (mod.rel.imp.) and the goodness of-fit of the model (r2) are also specified for each predictor. The microbial functional and structural variables considered for modelling were: β-glucosidase, GLU; βxylosidase, XYL; phenol-oxidase, PHE; community respiration, RESP; Shannon diversity index, DIV; richness index, RICH; bacterial density, BAC.

  Fig 1

ACKNOWLEDGEMENTS

We gratefully acknowledge Veronica Granados and Miriam Colls for assistance in the field, and Jordi Bou Manobens for his help with the creation of the GIS map. Special thanks go to the team at Laboratoire Microorganismes: Génome et Environnement for their assistance in the laboratory, and to Cayetano Gutierrez-Canovas for his help with the statistical analyses. Giulia Gionchetta acknowledges the UdG fellow-ship (IFUdG 2016) from the University of Girona. This study was supported by the Spanish Ministry of Economy and Competitiveness through the project FUNSTREAM (CGL2014-58760-C3-R), by the Economy and Knowledge Department of the Catalan Government, through Grant/ Award Number: 2014 SGR 484, by the project DryHarshSal (RTI2018-097950-B-C21) of the Spanish Ministry of Science, Innovation and Universities, and through the Short-Term Scientific Mission (STSM) Grant awarded by the Science and Management of Intermittent Rivers and Ephemeral Streams (SMIRES) Cost-Action (Action number: 40271). We are grateful to the two anonymous reviewers that helped improve the initial version of the manuscript.

Supplementary information: Table of Contents:

 Supplementary Information S1: Detailed hydrological measurements, Page 2  Supplementary Information S2: Physical-chemical streams characteristics (when flow was present in most of the sites), Page 3  Supplementary Information S3: Details of the modelling procedure, Page 4  Supplementary Information S4: RELATE Mantel-like test, Page 5  Supplementary Figure S1: Sampling sites map, Page 6  Supplementary Figure S2: EEAs, RESP, DIV and RICH hydrological groups boxplots, Page 7  Supplementary Figure S2: Variance partitioning plot, Page 8  Supplementary Figure S3: Bacterial composition histograms, Page 9

Supplementary information S1: Detailed hydrological measurements

Daily variation in streambed and air temperatures were used to identify, for each sampling site, the presence of water in the streambed, which is strongly related to its temperature daily variation [START_REF] Constantz | Analysis of streambed temperatures in ephemeral channels to determine streamflow frequency and duration[END_REF][START_REF] Stromberg | Effects of stream flow intermittency on riparian vegetation of a semiarid region river (San Pedro River, Arizona)[END_REF]. The temperature transducers (Solinst Levelogger Gold Model 3001,Solinst Ltd,Georgetown,ON,Canada) and temperature data loggers (ACR SmartButton Logger, MicroDAQ) used for the hydrological characterization were allocated in all streams, on the bottom of riffles areas. Furthermore, the air temperature was obtained from the data loggers (Solinst Barologger Gold Model 3001,Solinst Ltd,Georgetown,ON,Canada) previously installed in the riparian zone of each stream site, or obtained from nearby meteorological stations (Servei Meteorologic de Catalunya; http://www.meteo.cat/). All the applied temperature loggers recorded data every 60 minutes. Once obtained the entire dataset for each sampling site the daily variations in streambeds and air temperatures were estimated by two ratios as presented in Colls et al., (2019).

Once data had been obtained, the daily variations in streambed and air temperatures were characterized by two ratios: the daily streambed-to air temperature amplitude ratio (DA) and the streambed-to-air temperature change rate ratio (that is, heating or cooling; RTC). Both ratios were based on the relationship between the daily variations in streambed and air temperature. To reinforce results, each one assessed a different aspect of the temperature oscillations. DA was determined by the difference between the maximum (max(T)) and the minimum (min(T)) temperatures in the streambed and the air during a whole day. RTC was determined as the ratio between the highest hourly temperature change rate in the streambed and the air, in each sampling site.

A moving average of the DA and RTC ratios was calculated, comprising a total of 5 days to overcome potential erroneous interpretations of the stream's flow or no-flow state. This allowed dampening extreme values and achieving the best fit to water level data. The output value was named hydrological status (HS). Starting from the HS serial datasets we calculated all the hydrological parameters for each streambed.

For more specific information about the hydrological characterization and to visualize the DA, RTC and HS equations we referred to Colls et al. (2019) and Arias-Real et al., (2018), being part of the same multi-sites project.

Supplementary information S2: physical-chemical streams characteristics

Supplementary information S3: Details of the modelling procedure

Specifically, the procedure followed consisted in 4 main steps:

The first step is to ensure the requirements for data quality and consistency for the subsequent analyses. Once the data transformation and standardisation was assessed we proceed addressing the collinearity through the Variance Inflation Factor (VIF), which is a way to quantify collinearity. This method also accounts for non-linear relationships, which may remain undetected using correlation analysis. A VIFs>8 indicates variance-inflated variables [START_REF] Zuur | A protocol for data exploration to avoid common statistical problems[END_REF]. If more than one variable exhibit high VIFs, we exclude them stepwise, starting with the variable that has the highest VIF.

The second step consisted in the application of the exploratory analysis (Random Forest) which helped to decide which variables to include in the model of each response variable, bearing in mind, however that those known to be strongly related to biology should be preferred [START_REF] Feld | Analysing the impact of multiple stressors in aquatic bio-monitoring data: a 'cookbook' with applications in[END_REF]. Random Forest (RF) analysis is a flexible, non-parametric regression tool belonging to the family of Classification and Regression Tree analyses (CART). RF can handle a large number of predictors even in small datasets with a low number of observations. RF is suited to analyse non-linear relationships and complex interactions (Breiman, 2001) and can handle missing values. For this study, we chose to use rfsrc() function in the ForestSRC package (r software), because of its capability to rank, identify and visualise predictors interactions. This step aims to quantify and test the significance of single and combined effects of the explanatory variable candidates identified by RF analysis. To determine if these candidates have significant additive or interactive effects and to compare the effects, we need to generate standardised effect sizes (SES), see [START_REF] Grueber | Multimodel inference in ecology and evolution: challenges and solutions[END_REF]. The third step consists in modelling. Here, we used multi-model inference to select the final model for testing multiple-predictors effects for each response variable. In multiple-predictors modelling, the primary goal is to determine the importance and significance of various additive and combined (interactive) explanatory variables. Multi-model inference is a more objective method of model selection providing consensus models and variable importance estimates across a set of "top models", instead of focussing on a single model [START_REF] Grueber | Multimodel inference in ecology and evolution: challenges and solutions[END_REF]. The R package MuMIm [START_REF] Barton | MuMIn: Multi-Model Inference[END_REF] offers a useful set of functions to perform the steps required from model ranking to averaging. First, a global GLM (full) model including all the variables of interest is run. Second, the function dredge() automatically runs all models possible with different combinations of predictors variables based on a minimum and maximum number of predictors (1 for each 10 data; predictors selection based on RF analysis and coherent with the objectives of the study). Third, we select a set of best models that potentially best describe multiple-predictors patterns using the function get.models(). The argument "subset" allowed restricting the output based on thresholds of cumulative AIC ≤ 0.90, in order to obtain the best model in the set. Thus, we derived an averaged model for each response variable using the function model.avg() based on AIC values [START_REF] Grueber | Multimodel inference in ecology and evolution: challenges and solutions[END_REF]. The last step consists in validate the final model by checking model assumptions in order to represent a statistically valid linear model, its residuals need to exhibit normality, homoscedasticity (i.e.

homogeneous variability of residuals along the response gradient), without strong influential observations and any spatial or temporal autocorrelation in the residuals [START_REF] Zuur | A protocol for data exploration to avoid common statistical problems[END_REF]. References:


Silke Feuerriegel

Thomas A Kohl

Christian Utpatel

Sönke Andres

Florian P Maurer

Jan Heyckendorf

Agathe Jouet

Nelly Badalato

Lynda Foray

Rashidatu Fouad Kamara

Osman S Conteh

Philip Supply

Stefan Niemann
email: sniemann@fz-borstel.de

Rapid genomic first-and second-line drug resistance prediction from clinical Mycobacterium tuberculosis specimens using Deeplex®-MycTB

Targeted next generation sequencing using the Deeplex-MycTB assay can rapidly generate comprehensive drug resistance profiles from Mycobacterium tuberculosis complex cultures and directly from tuberculosis patient samples to guide personalized treatment.

To the Editor:

The emergence of multidrug-resistant (MDR) and extensively-drug-resistant (XDR-TB) Mycobacterium tuberculosis complex (Mtbc) isolates challenges tuberculosis (TB) control worldwide [1] and requires the rapid determination of extensive resistance profiles [START_REF] Dheda | The epidemiology, pathogenesis, transmission, diagnosis, and management of multidrug-resistant, extensively drug-resistant, and incurable tuberculosis[END_REF], enabling prompt initiation of effective treatment regimens. Phenotypic drug susceptibility testing (pDST) takes up to six weeks [START_REF] Köser | Whole-genome sequencing for rapid susceptibility testing of M. tuberculosis[END_REF], and is unreliable and/or not standardized for several drugs according to WHO guidelines [START_REF] Andres | Bedaquiline-resistant Tuberculosis: Dark Clouds on the Horizon[END_REF]. Molecular DST (mDST) assays, like Xpert MTB/RIF (Cepheid) or MTBDRplus/sl line probe assays (LPAs, Hain Lifesciences) can be performed directly from clinical specimens, but only target a limited number of resistance variants [START_REF] Makhado | Outbreak of multidrug-resistant tuberculosis in South Africa undetected by WHO-endorsed commercial tests: an observational study[END_REF].

This limitation can be overcome by whole-genome sequencing (WGS) allowing for the most comprehensive analysis for resistance variants [START_REF] Meehan | Whole genome sequencing of Mycobacterium tuberculosis: current standards and open issues[END_REF]. However, direct sequencing of sputum samples is challenging due to low amounts of Mtbc DNA [START_REF] Brown | Rapid Whole-Genome Sequencing of Mycobacterium tuberculosis Isolates Directly from Clinical Samples[END_REF].

To overcome this difficulty, an all-in-one targeted deep-sequencing assay of a 24-plexed amplicon mix, named Deeplex ® -MycTB (Genoscreen) has been developed. The assay covers 18 regions associated with drug resistance (rpoB, katG, fabG1, ahpC, inhA, pncA, embB, gidB, rpsL, gyrA, gyrB, ethA, eis, rrs, tlyA, rplC, rrl, rv0678), and enables species identification and genotyping of Mtbc isolates. However, very limited data are available on the performance of this assay for predicting resistance [START_REF] Makhado | Outbreak of multidrug-resistant tuberculosis in South Africa undetected by WHO-endorsed commercial tests: an observational study[END_REF][START_REF] Tagliani | Culture and Next-generation sequencing-based drug susceptibility testing unveil high levels of drug-resistant-TB in Djibouti: results from the first national survey[END_REF]. Therefore, we evaluated the performance of the Deeplex ® -MycTB assay for resistance and susceptibility prediction from DNA extracts from 81 cultured isolates and 50 clinical samples, by comparison with results from WGS, pDST and/or MTBDRplus/sl LPAs.

DNA was obtained from cultured isolates from an Mtbc isolate collection from Sierra Leone [START_REF] Feuerriegel | Sequence analysis for detection of first-line drug resistance in Mycobacterium tuberculosis isolates from a high-incidence setting[END_REF]. Clinical specimens were collected at the German National Reference Laboratory for Mycobacteria from January-November 2018 by convenience sampling, stored at -80°C and processed as follows: after decontamination/homogenization (2-5 ml sample) with 3% NaOH/NALC, heat inactivation (20 min at 95°C), DNA was purified by using the QIAamp DNA-Mini-Kit (Qiagen).

for Deeplex ® -MycTB PCR. Next-generation sequencing (NGS) was performed using 1 ng of prepared amplicons as input for the Nextera XT library preparation kit (Illumina) and sequencing on the MiniSeq and NextSeq 500 platforms. FASTQ files were uploaded to the cloud-based analysis software (app.deeplex.fr) or analysed with the MTBseq pipeline when performing WGS of genomic DNA from cultured isolates [START_REF] Kohl | MTBseq: a comprehensive pipeline for whole genome sequence analysis of Mycobacterium tuberculosis complex isolates[END_REF].

Using the Deeplex ® -MycTB assay, an average read coverage depth between 2,617 and 7,593 was obtained on individual samples with a mean average coverage depth of 4,476 over the reference targets across the 81 samples. On DNA samples from Sierra Leone [START_REF] Feuerriegel | Sequence analysis for detection of first-line drug resistance in Mycobacterium tuberculosis isolates from a high-incidence setting[END_REF], 466 variants were detected by Deeplex ® -MycTB in the 18 resistance genes analysed, while WGS analysis identified 464 of these 466 variants. Of these, 100 and 98 were annotated as resistance mediating mutations by Deeplex ® -MycTB and WGS, respectively. Accordingly, the overall concordance between both assays was 99.6% for all SNPs with 98% for resistance mutations. The discordance was 2 minority resistance variants (frequency below 50%) in rpoB and embB with a frequency of 6.8% and 4.9%, determined by the Deeplex ® -MycTB assay. These mutations were detected by Deeplex ® -MycTB with amplicons at 2,222fold and 11,901-fold coverage depths, respectively, but were undetected by WGS due to comparatively lower coverage depths on the corresponding gene positions.

When testing the performance of Deeplex ® -MycTB on 50 clinical samples, complete resistance predictions with Deeplex ® -MycTB could be made for 39. All smear positive (+, 2+, 3+) samples resulted in a sequencing acceptability grading of +, 2+ or 3+, by the Deeplex ® -MycTB software, and allowed complete resistance predictions. In addition, five out of seven (71.4%) smear +/-samples and two out of nine (22.2.%) smear negative samples also generated fully interpretable sequencing results. Overall, seven out of eleven samples that failed were microscopically negative, two were +/-, and for two samples no microscopy results were available. Overall, this resulted in the detection of 45 resistance mediating variants in rpoB, katG, inhA, pncA, embB, gidB, rpsL, gyrA, eis, rrs, and ethA. Of those, 26 variants were fixed variants (frequency of at least 97%), five were unfixed variants (frequencies below 97%) in genes rpoB, pncA, embB, rrs and ethA and 14 were minority variants (below 50%) in genes rpoB, katG, ahpC, pncA, embB, gidB, rrs, gyrA (n=5) and ethA (n=2).

All of the 16 fixed or majority resistance-mediating mutations (frequency of at least 97%), detected by Deeplex ® -MycTB in the regions covered by MTBDRplus/sl LPAs were concordantly detected by LPA testing, also performed on clinical samples (Tab.1). However, three instances of minority variants in rpoB (S450L; 27.8%), katG (S315G; 5.5%) and gyrA (A90V, 10.5%, 19.3%, 25.9%, 28%) detected by Deeplex-MycTB were undetected by LPAs also in repeated assays. This is somewhat surprising as previous papers suggest a lower mixed population threshold of ≥5% for LPAs [START_REF] Rigouts | Fluoroquinolone heteroresistance in Mycobacterium tuberculosis: detection by genotypic and phenotypic assays in experimentally mixed populations[END_REF]. However, this may be variable and depend e.g. on the nature of the sample and the binding kinetics of particular probes.

Concordances between predictions of resistance and susceptibility by Deeplex ® -MycTB and pDST data for the first-line drugs obtained after culture of these clinical samples ranged between 94.9% (INH) and 97.4% (RIF, PZA, EMB) (Tab.1). For second-line drugs, the concordances ranged from 66.7% (FQ) to 75% (PTH) and 100% (aminoglycosides, LIN, BDQ).

The discordance was observed for samples with minority variants known to confer low-level resistance. One sample contained Mtbc DNA with a katG S315G mutation reportedly associated with low-level resistance to INH [START_REF] Kandler | Validation of Novel Mycobacterium tuberculosis Isoniazid Resistance Mutations Not Detectable by Common Molecular Tests[END_REF] at a frequency of 5.5%. The mutation was not detected by LPA and the culture isolate subsequently tested susceptible at the INH critical concentration (CC). Four samples contained minority gyrA A90V subpopulations reported to confer low-level resistance to MXF [START_REF] Sirgel | gyrA mutations and phenotypic susceptibility levels to ofloxacin and moxifloxacin in clinical isolates of Mycobacterium tuberculosis[END_REF], which are particularly difficult to detect phenotypically when present at low frequency only. Likewise, these subpopulations were not detected by LPAs and showed susceptible phenotypes at the MXF CC. The single discordance among 39 samples phenotypically tested for PZA carried a pncA Q10P mutation at a frequency of 10.6%. Moreover, pDST for PZA is considered less reliable [START_REF] Miotto | A standardised method for interpreting the association between mutations and phenotypic drug resistance in Mycobacterium tuberculosis[END_REF].

For RIF, the single discordance consisted of a sample carrying a low-level resistance rpoB H445D (99.7 %) mutation, which tested RIF susceptible by pDST at the CC. This mutation is often linked with susceptible pDST results but is known to confer clinical resistance to RIF [START_REF] Rigouts | Rifampin resistance missed in automated liquid culture system for Mycobacterium tuberculosis isolates with specific rpoB mutations[END_REF].

The remaining discordance included a sample without detected INH resistance-associated mutation that tested INH resistant by pDST, two samples without ETH/PTH resistanceassociated mutations tested PTH resistant by pDST, and one sample carrying an embB G406S mutation tested EMB susceptible by pDST (Tab.1). The embB G406S mutation is associated with resistance to EMB. However, the pDST is known to be subjected to reproducibility problems [START_REF] Walker | Whole-genome sequencing for prediction of Mycobacterium tuberculosis drug susceptibility and resistance: a retrospective cohort study[END_REF].

The analytical range of Deeplex ® -MycTB is superior to other mDST assays due to Deeplex ® -MycTB covering resistance-mediating genes for all WHO group A and B drugs and covering much larger regions (rpoB, katG, embB, rrs). This advantage has been illustrated in a study performed in South Africa, reporting MDR outbreak isolates carrying a I491F mutation outside the RIF resistance-determining region [START_REF] Makhado | Outbreak of multidrug-resistant tuberculosis in South Africa undetected by WHO-endorsed commercial tests: an observational study[END_REF].

Comprehensive interrogation of genomic regions involved in resistance development is crucial for the prediction of susceptibility and resistance, which is potentially allowing to replace pDST for a large fraction of diagnostic samples. Here, all of the 98 resistanceassociated mutations detected by WGS were also identified by Deeplex®-MycTB in the tested isolates. However, data included mainly mutations conferring resistance to first-line drugs.

Further, another low-level frequency variant (18.4%) detected in ahpC (G-48A) by Deeplex ® -MycTB and missed by the LPA in a sample, was confirmed as INH resistant by pDST, thus supporting the superiority of Deeplex ® -MycTB to detect heteroresistance.

In conclusion, our results suggest that the Deeplex ® -MycTB assay is a valuable tool for rapid resistance and drug susceptibility testing in TB diagnostics, with the potential to replace pDST for a large fraction of patients with drug susceptible, MDR, and XDR-TB. Especially when performed on clinical specimens, this test will reduce the load of pDST in TB diagnostic laboratories and the turnaround time for reporting. The more comprehensive coverage of resistance targets, and its capability to detect bacterial subpopulations, are important advantages compared to other mDST assays [START_REF] Nikam | Redefining MTBDRplus test results: What do intermediate results actually mean?[END_REF]. Additional prospective studies with larger sample numbers are under way to further assess the challenges in implementation and performance characteristics of Deeplex ® -MycTB in settings with different epidemiological backgrounds and infrastructure. RIF, rifampicin (critical concentration 1.0 µg/ml); INH, isoniazid (0.1 µg/ml); PZA, pyrazinamide (100 µg/ml); EMB, ethambutol (5 µg/ml); SM, streptomycin; FQ, fluoroquinolones (moxifloxacin 0.25 µg/ml and 1 µg/ml); ETH, ethionamide; PTH, prothionamide (2.5 µg/ml); KAN, kanamycin (2.5 µg/ml); AMI, amikacin (1.0 µg/ml); CAP, capreomycin (2.5 µg/ml); LIN, linezolid (1.0 µg/ml); BDQ, bedaquiline (1.0 µg/ml); CFZ, clofaximine 3 One sample with a minority mutation pncA Q10P (10.6%) was tested PZA S by pDST. 4 One sample with a resistance-associated embB G406S mutation (99.7%) was tested EMB S by pDST. 5 Four samples with minority, low level moxifloxacin-conferring mutation gyrA A90V (10.5%, 19.3%, 25.9%, 28%) tested moxifloxacin S [START_REF] Sirgel | gyrA mutations and phenotypic susceptibility levels to ofloxacin and moxifloxacin in clinical isolates of Mycobacterium tuberculosis[END_REF]; eight without mutation were concordantly tested FQ S. 6 Two samples without resistance-associated mutations tested PTH R by pDST; three other samples with such mutations and three without such mutations concordantly tested PTH R and PTH S by pDST, respectively. 7 Four samples without mutations tested KAN S; the sample with a KAN R-associated mutation was not tested. 8 One sample with a resistance-associated mutation and seven without such mutations concordantly tested KAN, AMI and/or CAP R and S by pDST, respectively. 9 Five samples without resistance-associated mutation concordantly tested CAP S.

1o rpoB S450L (27.8%; in combination with rpoB H445D at 75.3% in a same sample) was not detected by Hain MTBDRplus.

11 katG S315G (5.5%) was not detected by Hain MTBDRplus.

  Tab. 1: Prediction of resistance and susceptibility by the Deeplex ® -MycTB assay for 39 clinical samples compared to phenotypic DST results and Hain MTBDRplus and MTBDRsl assay results. R, number of isolates tested resistant; S, number of isolates tested susceptible; nd, not done; second line phenotypic DST testing was only performed on MDR isolates. pDST was performed in the BACTEC MGIT960 system. From the 11 specimens that failed to give an interpretable sequencing result with the Deeplex ® -MycTB assay, 7 were microscopy negative, 2 were +/-and for 2 no microscopy result was

	(1.0 µg/ml); available. 1 A sample carrying a "disputed" rpoB H445D (99.7%) was tested RIF S by pDST [15]. Deeplex ® -MycTB assay pDST		Hain MTBDRplus/sl
	drug 2 One sample without mutation tested INH R by pDST; one sample with a minority, low level INH resistance-conferring mutation katG S315G locus R R S R S concordance R S concordance
	(minority mutation) (5.5%) was tested INH S by pDST [12].
	RIF	rpoB	8	0	31	7 1	32 1	38/39 1 (97.4%)	8 10	31	39/39 (100%)
	INH	katG/fabG1/	6	2	31	8 2	31 2	37/39 2 (94.9%)	4 11	34	37/38 (97.4%) (katG)
		ahpC/inhA							3	34	37/37 (100%) (inhA)
	PZA	pncA	4	1	34	4	35 3	38/39 3 (97.4%)	-	-	-
	EMB	embB	4	0	35	3 4	36 4	38/39 4 (97.4%)	1	2	3/3 (100%)
	SM	gidB/rpsL/rrs	5	2	32	nd	nd	nd	-	-	-
	FQ	gyrA/gyrB	0	4	35	0 5	12 5	8/12 5 (66.7%)	0 12	39	35/39 (89.7%) (gyrA)
									0	39	39/39 (100%) (gyrB)
	ETH/	inhA/fabG1/	3	0	36	5 6	3 6	6/8 7 (75.0%)	3	34	37/37 (100%)
	PTH	inhA
	KAN	eis	1	0	38	0 7	4 7	4/4 8 (100%)	-	-	-
	only
	KAN	rrs	1	0	38	1 8	7 8	8/8 8 (100%)	1	38	39/39 (100%)
	AMI
	CAP
	CAP	tlyA	0	0	39	0	5	5/5 9 (100%)	-	-	-
	only
	LIN	rplC/rrl	0	0	39	0	2	2/2 (100%)	-	-	-
	BDQ	Rv0678	0	0	39	0	2	2/2 (100%)	-	-	-
	CFZ

  12 Minority gyrA A90V variants (10.5%, 19.3%, 25.9% in combination with gyrA D94G minority variant at 8.3% in a same sample, 28%) not

	detected	by	Hain	MTBDRsl.

Acknowledgement

We thank Tanja Niemann, Vanessa Mohr and Fenja Boysen, Molecular and Experimental Mycobacteriology, Research Centre Borstel, Borstel, Germany, for excellent technical assistance.