Large-scale factors controlling biological communities in the Iberian Peninsula: an insight into global change effects on river ecosystems

The ongoing global environmental change poses a serious threat to rivers. Comprehensive knowledge of how stressors affect biota is critical for supporting effective management and conservation strategies. We evaluated the major gradients influencing spatial variability of freshwater biodiversity in continental Spain using landscape-scale variables representing climate, land use and land cover (LULC), flow regime, geology, topography, and diatom (n = 117), macroinvertebrate (n = 441), and fish (n = 264) communities surveyed in minimally impacted streams. Redundancy analysis identified the environmental factors significantly contributing to community variability, and specific multivariate analyses (RLQ method) were used to assess trait–environment associations. Environmental variables defined the major community change gradients (e.g., mountain–lowland). Siliceous, steep streams with increased precipitation levels favored stalked diatoms, macroinvertebrates with aquatic passive dissemination, and migrating fish. These traits were replaced by adnate diatoms, small macroinvertebrates, and non-migratory fish in lowland streams with warmer climates, calcareous geology, agriculture, and stable flow regimes. Overall, landscape-scale environmental variables better explained fish than diatom and macroinvertebrate community variability, suggesting that these latter communities might be more related to local-scale characteristics (e.g., microhabitat structure, substrate, and water physicochemistry). The upslope environmental gradient of river networks (e.g., slope, temperature, and LULC changes) was paralleled to the observed taxonomy-based and trait-based spatial variability. This result indicates that global change effects on riverine biodiversity could emerge as longitudinal distribution changes within river networks. Implementing management actions focusing simultaneously on water temperature, hydrological regime conservation (e.g., addressing LULC changes), and river continuity might be the best strategy for mitigating global change effects on river biodiversity.


Introduction
Rivers provide essential goods and services for human and societal well-being (Baron et al. 2002).However, over the past century, accelerated human-driven environmental changes, referred to as "global change", have been impairing both biotic and abiotic components and processes on Earth (Steffen et al. 2005).This has been especially relevant in freshwater ecosystems (Vörösmarty et al. 2010;Dudgeon 2019).Multiple anthropogenic pressures on rivers [e.g., water pollution, water abstraction, and hydromorphological changes (Vörösmarty et al. 2010;Dudgeon, 2019)] have led to alarming rates of biodiversity loss worldwide (Geist 2011;Baker et al. 2020).Given this context, it is fundamental to gain knowledge on ecological responses to global changerelated factors so that effective management actions that are 95 Page 2 of 17 in line with the post-2020 global biodiversity framework (Convention on Biological Diversity 2021; Nicholson et al. 2021) can be implemented to protect freshwater resources and biodiversity (Grizzetti et al. 2017;van Rees et al. 2021).
Changes in climate and land use and land cover (LULC) are the global changes that most compromise riverine ecosystems worldwide, in addition to dam construction (Liermann et al. 2012;Dudgeon 2019;Reid et al. 2019;Tickner et al. 2020;Caro et al. 2022).Climate change produces alterations in historical precipitation and temperature regimes, directly affecting water and sediment flows in rivers (Arnell and Gosling 2013;Schneider et al. 2013;Peng et al. 2020;Syvitski et al. 2022).Moreover, thermal regime changes are critical to riverine biota and reduce species richness, simplify assemblages, and change primary productivity (Durance and Ormerod 2007;McDonald et al. 2011;Song et al. 2018;Mouton et al. 2022).LULC changes (e.g., deforestation, intensive agriculture, and urbanization) across catchment and riparian areas frequently result in habitat degradation (Allan 2004;rösmarty et al. 2010) and changes in species distributions (Guo et al. 2018;Zeni et al. 2019;Liu et al. 2021).At the same time, hydropower expansion and water supply reservoirs constructed to meet increasing demands have led to flow, thermal, sediment, and nutrient regime alterations and the loss of habitat connectivity, among other remarkable impacts (Sauer et al. 2010;rösmarty et al. 2010;Zarfl et al. 2019).All these pressures profoundly alter habitat characteristics and energy fluxes in rivers, producing effects from individual to ecosystem levels, e.g., affecting organismal growth and life cycles (Larsen et al. 2016;Lourenço et al. 2023) or river primary production and ecosystem respiration (Bernhardt et al. 2022).
In recent years, a growing number of large-scale studies have investigated global change-related impacts (i.e., climate change, LULC, and flow alteration) on riverine biodiversity (Soininen et al. 2016;Comte et al. 2021;Mouton et al. 2022;Manfrin et al. 2023).However, most have focused on specific communities or populations, and the use of multiple and diverse assemblages to assess global change effects on biodiversity is still rare (but see Pound et al. 2021;Alric et al. 2021;Tison-Rosebery et al. 2022).Moreover, timeseries data from multiple biological communities and river reaches covering large spatial scales are very scarce (Blois et al. 2013;Lento et al. 2022).Thus, space-for-time substitution might be an appropriate technique to gain insights into the future impacts of global change on river ecosystems (Wogan and Wang 2018).
Freshwater biological communities (e.g., diatoms, macrophytes, benthic macroinvertebrates, and fish) have shown asymmetrical responses to stressors (Hering et al. 2006;Guse et al. 2015;Tison-Rosebery et al. 2022), suggesting that stress impact direction and intensity are highly dependent on food web trophic level or biological and ecological characteristics of the communities under study.Thus, research encompassing large spatial extents, highly variable environmental conditions, and multiple biological communities (e.g., primary producers to apex predators) is paramount for improving the understanding of the effects that global change might have on river ecosystems (e.g., De Castro-Català et al. 2020;Lento et al. 2022).
Significant biodiversity variability is expected when carrying out studies covering large environmental gradients (Liu and Wang 2018;Troia and McManamay 2019), as organisms are filtered according to their ecological requirements (Poff 1997;McGill et al. 2006).However, even across wide environmental gradients, general patterns underlying community assembly are not always evident due to the historic, dispersal-related, and evolutionary processes controlling regional species pools (Lessard et al. 2012;Soininen et al. 2016).In consequence, integrative approaches (e.g., the use of organisms' biological and environmental characteristics) that unify the species pool, reduce the influence of biogeographical processes, and that focus on biodiversity-environment relationships are essential (Lessard et al. 2012;Alric et al. 2021).In this context, trait-based approaches are a valuable alternative to taxonomy-based approaches, as traits consistently respond to habitat filters (Mouillot et al. 2013) and reduce biogeographic constraints imposed by single species at large spatial scales (Bonada et al. 2006;Menezes et al. 2010;Chen and Olden 2018;De Castro-Català et al. 2020).Furthermore, traits provide a mechanistic perspective for examining the links between biodiversity and environmental gradients or stressors (McGill et al. 2006;Olden et al. 2010;Culp et al. 2011).
The main goal of this exploratory study was to identify the key dynamic (i.e., climate, LULC, and hydrology) and static (i.e., geology and topography) environmental variables influencing the spatial variability in riverine biodiversity.This information will serve as a baseline for monitoring riverine biological community changes and helping to track anthropogenic effects related to global change.Biological datasets covering diatom, macroinvertebrate, and fish communities were originally collected to address the European Water Framework Directive (WFD) (European Commission 2000) water quality objectives in Spain, a country with highly variable environmental conditions (e.g., semiarid, Mediterranean, and temperate climates) in a relatively small area (over 0.5 million km 2 ).Biological sampling was intentionally limited to streams with good ecological status and without local anthropogenic pressures, so community shifts were mainly restricted to large-scale environmental gradients.
In this study, we expected that (1) the main ecological responses associated with climate, LULC, and hydrology would be identified since these variables can act as filters selecting taxonomic groups and traits.We also expected that (2) large-scale climate, hydrology, LULC, geology, and topography gradients would contribute differently to structuring diatom, macroinvertebrate, and fish communities, as these organisms have contrasting biological and ecological characteristics.Finally, we expected that (3) biological and ecological traits would be useful in determining these ecological relationships, as trait-based variability should be highly related to environmental gradients compared to taxonomy-based variability.

Study area
The study area corresponds to continental Spain (40°23′ N, 3°33′ E; Fig. 1).The relief is characterized by an extensive high inland plateau [average of 650 m above sea level (a.s.l.)] surrounded by mountain ranges, with peaks reaching 3400 m a.s.l.(Rivas-Martínez et al. 2004).Different climatic zones are distributed throughout the country: a mountain climate at high altitudes, with cold winters and abundant precipitation (annual averages reaching 3000 mm); an oceanic climate in the northwestern region; a hot steppe climate in the southeast, with minimum rainfall levels for Spain (annual averages < 150 mm); and Mediterranean climates with dry and hot summers (AEMET and IMP 2011;Serrano-Notivoli et al. 2018).Relief-climate interactions generate significant variability in environmental conditions and hydrological regimes, influencing riverine biodiversity across continental Spain (Peñas et al. 2014).

Hydrological and environmental characteristics
To integrate biological, hydrological, and environmental information into a spatial framework, virtual watersheds (Barquín et al. 2015;Benda et al. 2016) with flow directions inferred from a 10-m digital elevation model (Peñas and Barquín 2019) were used.At the selected sites, 26 catchment-scale environmental variables (Table 1), obtained from national and regional databases, were divided into two groups (Table 1 and Fig. 2): (1) dynamic environmental variables: climate (average annual temperature, precipitation, and evapotranspiration), LULC (mean fractions in the catchment occupied by agriculture, urbanization, forests, etc. ), and hydrology (non-correlated synthetic hydrological indices); and (2) static environmental variables: geology (mean catchment fraction occupied by different geological structures) and topography (slope, altitude, distance to river outlet, and catchment area, i.e., variables describing the position throughout the river network).For practical purposes, changes in static variables are considered imperceptible over time (e.g., soil type changes occurring over geologic time), whereas dynamic variables (i.e., LULC, climate, and hydrology) are expected to undergo relevant alterations over shorter timescales, especially when linked to human-driven activities such as climate change or land use change (Stanton et al. 2012).
The hydrological variables were composed of 85 hydrological indices (HIs), which were calculated from the normalized daily flow series recorded by 282 natural flow gauges throughout the country (Peñas and Barquín 2019).The HIs were reduced to a set of non-correlated synthetic indices (SI NAT ) and predicted for all river reaches using  S1 in Supplementary Information; further details in Peñas and Barquín, 2019): high SI NAT 1 scores corresponded to low intraannual variability in daily flows, high SI NAT 2 scores were associated with longer but less frequent high flow events, high SI NAT 3 scores corresponded to low interannual variability in the minimum flow timing, and high SI NAT 4 scores were associated with an increased magnitude and variability in April flows, i.e., higher and variable spring flows.

Characterization of biological communities and traits
In this study, national biological monitoring databases (NABIA; Spanish Royal Decree 817/2015 on Water Policy; BOE, 2015) were used.These databases have been compiled since 2008 by the Spanish Ministry for the Ecological Transition (MITECO), which gathers biological surveys carried out at the national level to assess water body ecological status in compliance with the European WFD (European Commission, 2000).
From these databases, only biological community samples collected in wadeable streams between June and October were selected to reduce the intraannual variability (Leathwick et al. 2005).Then, to reduce the effects of local anthropogenic pressures commonly found in rivers (see red boxes in Fig. 2), we retained only samples surveyed in the river reaches cataloged as "good" or higher in terms of ecological and chemical status (River Basin Management Plans 2015-2021;BOE, 2015) and those unaffected by local anthropogenic pressures, such as dams, embankments, or significant water abstractions upstream.This site selection process, covering only minimally disturbed streams, allowed us to examine the effects of dynamic variables over the selected riverine biological communities without the confounding effect of other stressors, such as water pollution, hydromorphological pressures, or the presence of dams and reservoirs (Fig. 2).
A total of 177 diatom, 441 macroinvertebrate, and 264 fish sampling sites were retained (Fig. 1).Physicochemical measurements were not available in all sampling sites and therefore were not included in the study.
Diatoms and fish were reported at the species level.Macroinvertebrates were mostly identified at the family level, as required by the Iberian Biomonitoring Working Party (IBMWP) scoring system (Alba-Tercedor et al. 2002), except for Acariformes and Oligochaeta.Ostracods were removed from the taxa list, as they are considered microcrustaceans and were not accounted for in all surveys.At the sites sampled for multiple years, the abundance of taxa was averaged to obtain a site-specific assemblage structure, following previous studies (e.g.Paavola et al. 2003;Filker et al. 2016).Given the large number of taxa retained, the rare taxa in the diatom and macroinvertebrate communities were eliminated (those representing less than 2% of total sampled individuals when considering the total number of occurrences; Lavoie, Dillon & Campeau, 2009).
Diatom species were assigned to ecological guilds, i.e., growth morphologies: low profile (including species of short stature, adhered to the substrate, and slow-moving species, adapted to high current velocities and low nutrient concentrations), high profile (tall-statured, erect, filamentous, branched, stalked, and colonial diatoms, adapted to high nutrient concentrations and low physical stress), motile (fast moving species, free of both resource limitation and disturbance stress), and planktonic (species adapted to lentic environments with morphological adaptations to resist sedimentation; Passy 2007;Rimet and Bouchez 2012).Diatom life forms (colonial, stalked, adnate, and pioneer) and cell sizes (biovolume) were also selected as traits, as seen in Table S2 in Supplementary Information.Macroinvertebrate families were assigned to biological traits using a fuzzy-coding approach (Chevenet et al. 1994), which uses positive scores (0 to 5 scale) to describe the affinity of taxa for different biological characteristics, e.g., feeding habits, food types, maximum body size, lifespan, the annual number of reproductive cycles,  S3 in Supplementary Information).Since macroinvertebrate traits were mainly reported at the genus level (Tachet et al. 2010), each trait affinity was averaged and rescaled to the family level using affinities of all genera calculated within a family.Although the use of family-level identification can result in loss of ecological information, previous studies (Dolédec et al. 2000;Gayraud et al. 2003;García-Roger et al. 2014) found that this approach had minor impacts on the trait-based description of macroinvertebrate communities when compared with the genuslevel approach.Finally, fish species were assigned to traits describing maximum body length, overall tolerance to stressors, habitat use, feeding habits and habitat, and migration patterns based on available information for Iberian species (Table S4 in Supplementary Information).

Data analyses
Relationships among environmental variables and biological community spatial patterns were explored with two different methods.First, using redundancy analysis (RDA) (Legendre and Legendre, 1998), we identified the environmental factors that best explained the communities' spatial variability.RDA establishes a model measuring the variance in biological datasets in relation to environmental characteristics (Legendre and Legendre 1998).For each biotic group, using both taxonomic and trait-based datasets [i.e., community-weighted mean trait values (CWM)], a subset of explanatory variables was obtained using a stepwise forward selection procedure with 9999 permutations and a double-stopping criterion (α < 0.05 and the global model adjusted coefficient of determination) to avoid type I error inflation (Blanchet et al. 2008).The CWM trait approach calculates an aggregated trait value based on the  3 ± 327.4 329.3-1796.3 1162.4 ± 354.3 329.3-1991.7 1132.1 ± 342.6 455.9-1947 relative abundance of species in a given community (Ackerly and Cornwell 2007).Then, the adjusted coefficient of determination (R 2 adj ) of each RDA model was calculated accounting for the number of predictors and samples, which allowed comparisons between groups (Peres-Neto et al. 2006).For the taxonomy matrices, Hellinger transformation was used to reduce the influence of abundant taxa in the community structure (Legendre and Gallagher 2001); for the trait matrices, the average trait expression of all Hellinger-transformed taxa was weighted by taxa abundance at each site (Kleyer et al. 2012).
The second analysis investigated the overall and individual trait-environmental relationships using RLQ and fourth-corner analyses (Dray et al. 2014).These two methods constitute a powerful approach for investigating trait-environment relationships in community ecology (Peres-Neto et al. 2017), with numerous examples in the literature (e.g., Díaz et al. 2008;Brown et al. 2019;Zeni et al. 2019).The RLQ analysis (Dolédec et al. 1996), an exploratory multivariate method, summarizes the joint structure among three tables: Table R (n × m), containing the measurement of m environmental variables at n sites; Table L (n × p), with the abundance of p taxa at n sites; and Table Q (p × s), describing s traits related to the same p taxa.The method generates a fourth matrix representing the cross-variance structure explained by environmental variables and traits, providing a broad overview of how they are associated (Brown et al. 2014;Thioulouse et al. 2018).The overall significance of the trait-environment relationship can be assessed with 9999 permutations using a two-step procedure (model 6) (Dray et al., 2014), which simultaneously tests the permutation of sites (rows) and species (columns) while controlling the type I error (ter Braak et al. 2012).The fourth-corner method, in contrast, calculates a matrix (s × m) with the one-to-one trait-environment correlations (Dray et al. 2014).The combination of these two methods allowed us to identify the main patterns and individual trait-environment relationships.In addition, we evaluated the responses of individual taxa to environmental variables using Spearman rank correlations for the three biotic groups.
Prior to the analyses, the environmental variables were standardized (zero mean and unit standard deviation) to reduce the effects of different measurement units and improve normality (Quinn and Keough 2002).All analyses were performed in the R environment (R Core Team 2020), with the packages ade4 (Dray and Dufour 2007), adespatial (Dray et al. 2020), and vegan (Oksanen et al. 2019).

Contribution of environmental variables to biological variability
The selected environmental variables obtained with the RDA analyses (Table 2 and Table S5 in Supplementary Information) accounted for 14.9%, 19.0%, and 35.9% of the explained taxonomic variation in the diatom, macroinvertebrate, and fish communities, respectively.The dynamic variables (i.e., linked to climate, LULC, and hydrology) showed significant contributions in the three biotic groups.The climate-related variables were clearly the most important variables explaining the taxonomic variation in diatoms (TEMP with R 2 adj of 0.047, which corresponded to 31.5% of the total explained variability of the final RDA model, as seen in Table 2 and Table S5 in Supplementary Information).In the macroinvertebrate communities, climate also played a major role in explaining spatial variability (PREC and TEMP accounting together for an R 2 adj = 0.079, which corresponded to 42% of the total explained variability).In the fish communities, hydrology accounted for a major contribution (SI NAT 2, representing the magnitude and frequency of high flow events, accounted for an R 2 adj of 0.075, which corresponded to 21% of the total explained variability).Among the static variables, topography (slope and altitude; SLP and ELEV, respectively) and geology (siliceous and calcareous bedrocks; SLIC and CALC, respectively) also contributed significantly to the variation in the three biotic groups, especially in macroinvertebrates.
When using community-weighted mean traits, the explanatory variables selected by the CWM-RDA (Table 2 and Table S6 in Supplementary Information) accounted for 12.9%, 19.2%, and 54.6% of the diatom, macroinvertebrate, and fish community variability, respectively.LULCrelated variables (agriculture, broadleaf forests, and urban areas; AGR, BLF, and URB, respectively) accounted for most of the explained variation in diatom trait variability (Table S6 in Supplementary Information).In contrast, calcareous bedrock (CALC) had an R 2 adj of 0.043, which corresponded to 22% of the total explained variability in the macroinvertebrate trait variability.Fish trait variability was mainly related to precipitation (PREC with R 2 adj = 0.241, which corresponded to 44% of the explained variation).
In both the taxonomic and trait-based analyses, fish communities had the highest amount of variation explained by the catchment-scale environmental variables (on average, 45.3%, in contrast to 13.9% for diatoms and 19.1% for macroinvertebrates).Moreover, when using the trait-based approach, the environmental variables accounted for larger amounts of explained variation in fish (an increase of 19% in relation to taxonomy), whereas in diatom and macroinvertebrate communities, the contribution of environmental variables scarcely changed (decrease of 2% and an increase of < 1%, respectively, in relation to taxonomy).

Diatom assemblages
In the 177 diatom samples, 95 diatom species were retained.The most abundant genera were Achnanthidium (corresponding to 42.2% of all sampled diatoms), Cocconeis (14.8%),Gomphonema (9.4%), and Nitzschia (5.4%).Achnanthidium minutissimum was the most ubiquitous species, with an average relative abundance of 20.1% at all sampling sites.The high-profile guild, stalked life forms, and small diatoms dominated most locations.
The first two axes of the RLQ ordination representing diatom assemblages (Fig. 3) accounted for 78.2% of the crossvariance between environmental variables and diatom traits.A global RLQ permutation test indicated a nonsignificant overall association between the environment and traits (permutation across sites p < 0.001; permutation across species p = 0.106).
The first axis (64.0% of the cross-variance) showed a marked gradient of stream size (i.e., longitudinal variation) associated with topography, climate, LULC, and hydrology.Mountain and headwater streams, characterized by high slopes (+ SLP) and precipitation intensities (+ PREC), were positively correlated with species such as Achnanthidium pyrenaicum (SLP r = 0.26, p < 0.001; PREC r = 0.26, p < 0.001) and Gomphonema pumilum (SLP r = 0.15, p = 0.041; PREC r = 0.21, p = 0.006).These sites were segregated from lowland streams, which had larger catchment areas (+ AREA), agricultural (+ AGR) and urban activities (+ URB), and streams with more stable flow regimes (+ SI NAT 2) and were positively linked to species such as Amphora pediculus (AREA r = 0.32; AGR r = 0.46; URB r = 0.44; SI NAT 2 r = 0.19; p < 0.01 in all cases) and Nitzschia Fig. 3 Ordinations showing the first two RLQ axes for diatom communities.a Site scores according to diatom genus dominance; b coefficients for key environmental variables (see Table 1 for codes) with significant links to traits (see Fig. S2 in Supplementary Information); c coefficients for traits (see Table S2 for codes).The "d" values represent the grid size for scale comparison across the graphs 95 Page 8 of 17 dissipata (AREA r = 0.21; AGR r = 0.34; SI NAT 2 r = 0.22; p < 0.005 in all cases).This gradient was also reflected in the diatom traits, with a higher dominance of stalked diatoms (Stalk, e.g., Achnanthidium lineare) in the headwater streams being replaced by adnate life forms (Adnate, e.g., Cocconeis euglypta), pioneer (e.g., A. pediculus), and motile diatoms (e.g., Nitzschia fonticola) downstream (see the fourth-corner analysis in Fig. S2 in Supplementary Information).
The first two axes of the RLQ ordination (Fig. 4) accounted for 81.6% of the cross-variance between the environmental variables and macroinvertebrate traits.A global RLQ permutation test did not show a significant overall trait-environment association (permutation across sites p < 0.001; permutation across species p = 0.763).

Fish assemblages
Thirty-eight species of fish were retained for the community analyses.The most abundant taxa were Salmo trutta (corresponding to 29.8% of all sampled fish and present in most samples), Phoxinus bigerri (28.7%), and Parachondrostoma miegii (6.7%).
The first two axes of the RLQ ordination (Fig. 5) accounted for 86.2% of the cross-variance between the environment and fish traits.A global RLQ permutation test indicated a significant trait-environment relationship (permutation across sites p < 0.001; permutation across species p = 0.041).

Discussion
In this study, we investigated the main environmental gradients related to the spatial variability in diatom, macroinvertebrate, and fish communities in Iberian rivers.As expected, the power of the landscape-scale environmental variables to explain the spatial variability of the studied biological communities differed significantly among biotic groups and taxonomic-based or trait-based approaches.Both dynamic (i.e., temperature, precipitation, extreme flow events, agricultural, and urban areas) and static (i.e., slope, altitude, calcareous, and siliceous geology) variables contributed significantly to the observed biological variability.The most consistent spatial pattern observed on biological community changes was associated to the longitudinal variation on environmental conditions from mountain streams to lowland rivers for the three studied riverine communities.This result is particularly relevant for establishing past or current baselines for monitoring current or future changes, respectively, on diatom, macroinvertebrate, and fish community composition along river networks.The approach developed in this study is not valid for demonstrating specific cause-effect relationships, but the obtained results are highly useful for proposing specific hypotheses and future research directions regarding the potential effects of global change on riverine ecosystems.

Role of dynamic environmental factors
Climate, LULC, and hydrology contributed significantly to the spatial variability in the three biotic groups, suggesting that these communities might be highly sensitive to global change effects (see Table 3).
Our findings showed a significant influence of hydrology and climate on diatom community spatial variability, in agreement with the results of previous studies (Wu et al. 2019;Guo et al. 2020).Pioneer diatoms, Nitzschia, and other genera dominances (e.g., Amphora and Navicula) were positively linked to longer but less frequent high flow events but negatively linked to precipitation.This pattern was expected because loosely attached diatoms (e.g., Nitzschia and Navicula) are more susceptible to drift and are easily displaced from periphyton after floods (Biggs and Thomsen 1995;Passy 2007;Schneck and Melo 2012).Moreover, mountain streams (i.e., steep headwaters, where agriculture or urbanization are usually less intense than in lower areas) were associated with high profile and stalked diatoms, whereas lowland rivers with more intensive human land uses were linked to motile diatoms (Table 3).The lower grazing pressures usually encountered in headwaters may explain the dominance of diatoms growing in the biofilm upper layers (e.g., stalked and high profile; Steinman, 1996;Marcel et al., 2017).The motile guild increased along nutrient gradients: Nitzschia (a genus typical of nutrient-rich and organically polluted waters; van Dam, Mertens & Sinkeldam, 1994; Table 3 Main significant community shifts (increases and decreases in abundance) in response to dynamic (i.e., global change-related) variables identified with the fourth-corner analysis in the diatom, macroinvertebrate, and fish assemblages in continental Spain Decrease in precipitation levels and increase in temperature Decrease of natural forest by agriculture, urbanization, or others Increase in the number of high flow events and reduction in their length • Increase in pioneer diatoms (e.g., Amphora pediculus) • Increase in nonmigratory species (e.g., Phoxinus bigerri) • Increase in motile diatoms (e.g., Nitzschia fonticola) • Increase in macroinvertebrates with small body sizes and no resistance strategies • Increase in nonmigratory fish species (e.g., Phoxinus bigerri) • Increase in rheophilic fish (e.g., Parachondrostoma turiense) • Decrease in scrapers (e.g., mayflies and caddisflies) • Decrease in stalked diatoms (e.g., Cymbella excise and Rhoicosphenia abbreviata) • Decrease in macroinvertebrates with active aquatic dissemination (e.g., Nemouridae and Ephemerellidae) • Decrease in intolerant fish species (e.g., Alosa alosa) • Decrease in pioneer diatoms (e.g., Amphora pediculus) • Decrease in piercers (e.g., Athericidae and Gerridae) Page 11 of 17 95 Biggs & Smith, 2002;Passy, 2007) and other dominant genera (e.g., Amphora and Navicula) were replaced by the Achnanthidium genus, which is adapted to low nutrient conditions (Berthon et al. 2011;Rimet and Bouchez 2012).Hydrology played a significant role in the spatial variability in the macroinvertebrate community.The presence of macroinvertebrates with active dissemination modes (aquatic or aerial) was associated with the increased magnitude and variability in spring flows.Spring flow disturbances are known to be crucial for univoltine and seasonal migrating insects (e.g., Coleoptera), as these disturbances may occur during a sensitive part of their life cycle (Williams 1996;Lytle 2008).Similarly, Belmar et al. (2019) reported an increase in small-sized and multivoltine macroinvertebrates in regulated Mediterranean rivers in Spain.
Moreover, climate and LULC influenced macroinvertebrate spatial variability, as increases in the mean temperature and urban and agricultural land uses favored small and multivoltine macroinvertebrates, while the relative abundance of scrapers (e.g., some mayflies and caddisflies) decreased.These results agree with Lourenço et al. (2023), who reported positive correlations among multiple stressors (i.e., nutrient enrichment, oxygen depletion, and thermal stress) and fast-living macroinvertebrate strategies.Temperature, daylight duration, and other climatic factors are known to control the development, life history, spawning time, and distribution of macroinvertebrates and other organisms at the catchment and larger spatial scales (Lancaster and Downes 2013).Therefore, climate change, as well as human land-use intensification and habitat loss, which are also major drivers of biodiversity losses in rivers (Pereira et al. 2012;Caro et al. 2022), can cause profound changes in macroinvertebrate communities (Durance and Ormerod 2007;Flenner et al. 2010).In line with these results, Mouton et al. (2022) reported increases in the populations of multivoltine species and longer adult lifespans, which were linked to climate and land use changes during a 26-year period in New Zealand, and Manfrin et al. (2023) found a significant increase in taxa richness adapted to warm temperatures in small low mountain streams in Central Europe over 25 years.Other studies, in agreement with our findings, also showed that global warming had a significant impact on macroinvertebrates, with increasing water temperatures leading to the homogenization of EPT assemblages (Timoner et al. 2020) and the predominance of resistance and resilience traits (small organisms and fast and seasonal development) in streams affected by urban and agricultural activities (Liu et al. 2021).
Fish assemblages also demonstrated important responses to dynamic variables.Decreases in precipitation and frequency of high flow events seemed to favor generalists (e.g., tolerant and nonmigratory) over specialists (e.g., migratory species), as also found in the results of other studies (Zampatti et al. 2010;Freitas et al. 2013;Whiterod et al. 2015).
Drought events and changes in hydrological patterns can be especially harmful to migrant species, as their movements may be constrained by the loss of connectivity among habitats (Milton 2009;Freitas et al. 2013).

Role of static environmental drivers and interactions among environmental variables
Static variables, in combination with dynamic variables, played a key role in the main community shifts identified in this study.Topography (slope and catchment area as surrogates for position in the river network), climate (precipitation and temperature), LULC (agricultural and urban land uses), and hydrology (magnitude and frequency of high flow events) were related to community changes across the upstream-downstream gradient.Stalked diatoms, EPT, and long-distance migratory fish in headwaters were replaced by pioneer and motile diatoms, noninsects, and nonmigratory fish in lowland streams.Previous research reported similar findings, with temperature and position in a river network as the major factors controlling the distribution of diatom (Potapova and Charles 2002;Tornés et al. 2022), macroinvertebrate (Griffith et al. 2001;Verdonschot 2006;Díaz et al. 2008), and fish communities (Buisson and Grenouillet 2009;Troia and McManamay 2019;Herrera-R et al. 2020).In the context of ongoing global change, this altitudinal gradient emphasizes the vulnerability of riverine biodiversity.Temperature increases and land use changes (e.g., forest loss) accelerate the upslope movement of species seeking suitable habitats (Guo et al. 2018); at the same time, the rate of warming is enhanced with altitude, increasing pressure on mountain ecosystems (Pepin et al. 2015).
Our findings also highlight that dynamic variables such as temperature, precipitation, and agricultural land use are often correlated in strong hierarchical relationships with natural static (e.g., altitude, latitude, position along the river network, and geomorphology; Vannote et al. 1980;McCain and Grytnes, 2010) and anthropogenic (e.g., land use; Allan, 2004) gradients, which could lead to an overestimation of variables' contribution to impacts (Snelder and Lamouroux 2010).Thus, the results emphasize the need to consider static and dynamic variables, as well as their interactions, when assessing global change impacts.

Study and method limitations
As expected, landscape-scale variables had different contributions to diatom, macroinvertebrate, and fish community variability.Relatively low amounts of explained variation and nonsignificant overall environmental-trait relationships were observed in diatoms and macroinvertebrates.These results could be attributed to the lack of relevant local-scale environmental descriptors (e.g., water physicochemistry, 95 Page 12 of 17 sediments, and other physical habitat characteristics at the reach level), which were not available at all sampling sites and hence not included in the analyses.Although the study focused on large-scale environmental gradients, the localscale descriptors are known for controlling diatom and macroinvertebrate community structure and composition (Dong et al. 2012;Soininen 2016;Álvarez-Cabria et al. 2017).In contrast, when considering fish assemblages, our environmental predictors explained relatively higher amounts of biological variation, and the global relationship between traits and environmental variables was significant.These findings suggest that landscape-scale environmental variables have an important influence on fish taxonomy-based and trait-based spatial variability, as also seen in the results of other studies (Esselman and Allan 2010;Pompeu et al. 2022).
The noisiness of the macroinvertebrate datasets, which were mainly reported at the family level, could also reflect their low explained variation.Although macroinvertebrate family-level taxonomy is considered suitable for cost-effective water quality assessment (Alba-Tercedor et al. 2002;Chessman et al. 2007), it can mask some of the relationships with environmental and hydrological factors (Monk et al. 2012;Jiang et al. 2013;Powell et al. 2023) and affect biodiversity patterns in regions with high species diversity (Heino and Soininen 2007).Some species exhibit greater heterogeneity in their habitat preferences and responses to disturbances in comparison to a family's average characteristics (e.g., water quality and flow preferences of species belonging to the Hydropsychidae family; Lenat and Resh, 2001).Thus, finer-level taxonomy (e.g., identification at the genus level only in families with a wide intrafamilial variation in environmental tolerance (Chessman et al. 2007) may yield more accurate results in macroinvertebrate analyses.
The CWM-RDA reported a slightly lower total explained variation in diatom communities when compared with the taxonomy-based RDA, contradicting hypothesis 3, which considered traits more advantageous in comparison to taxonomy when exploring the ecological responses to environmental gradients.This could have been related to the number of diatom traits (9 traits divided into 3 categories), which was much lower than the number of taxa (95 diatom species), resulting in considerably lower total inertia in the trait-based analysis relative to that in the taxonomy-based analysis (Heino et al. 2007).This difference was less pronounced in the macroinvertebrate (11 trait groups, 60 taxa) and fish (6 trait groups, 38 taxa) assemblages, in which traitbased datasets revealed higher total explained variation than the taxonomy-based datasets.
It must also be remarked that the site selection along the study area was not completely balanced.A large percentage of the study sites were concentrated in the northern catchments, whereas the proportion of sites in the southern ones was significantly lower.The heterogeneity in the site distribution along the country is due to the data availability and the site selection criteria.In this regard, the northern catchments have a larger proportion of rivers in good and very good ecological status and free of human pressures than the southern catchments.This situation can lead to the underrepresentation of rivers subjected to specific environmental conditions (e.g., arid climates).Nonetheless, the site selection procedure ensured that anthropogenic disturbances are minimized in our dataset and that the observed ecological patterns were arising from natural environmental gradients.

Implications for investigating global change effects on rivers and mitigation actions
This study has shown that shifts in biological communities were significantly linked to temperature, precipitation, hydrological indices, agriculture, geology, and slope (i.e., a proxy for position in the river network).These findings are promising for forecasting and anticipating the potential effects of global change on river biota since these dynamic variables are greatly related to the main sources of change, i.e., climate, LULC, and hydrological changes, and will help focus future research priorities.Specific study designs and ad hoc experiments accounting for the effect of these dynamic variables are paramount for determining the causeeffect relationships and cascade effects on ecosystem functioning and, ultimately, on the ecosystem services that rivers provide.
The three biotic groups, which are major components of riverine food webs, presented asymmetric responses to environmental variables, as also reported by Tison-Rosebery et al. (2022) and Lento et al. (2022).Diatoms and macroinvertebrates were more sensitive to climateand LULC-related factors than to other factors, indicating that these communities may be good indicators of the ecological effects of land use (e.g., water quality deterioration, organic pollution, and nutrient enrichment; Guse et al. 2015) and climatic changes (e.g., changes in temperature and precipitation trends; Timoner et al. 2020).Hydrology had an important effect on fish community variability, suggesting that fish communities may be good indicators of flow alteration and habitat connectivity (Geist and Hawkins 2016;Lento et al. 2022).Finally, environmental variables related to the position in the river network (i.e., upstream-downstream gradient) contributed substantially to community shifts in the three biotic groups.The upslope gradient can also be useful for delimiting climate-sensitive zones in which species or traits rapidly respond to climatic conditions and, therefore, for efficiently monitoring climate change impacts (Bässler et al. 2010;Shah et al. 2015).Furthermore, the changes observed across the longitudinal patterns (i.e., Page 13 of 17 95 from mountain streams to lowland rivers) in the three riverine communities reinforce the need to restore river connectivity to mitigate the effects of global change on these ecosystems.Finally, the results obtained in this study are important for not only increasing the understanding of the effects of global change on fluvial ecosystems but also informing integrated catchment management and guiding policy development aimed at mitigating these effects.For instance, many authors (e.g., Lammert and Allan 1999;Álvarez-Cabria et al. 2017;Segurado et al. 2018;Silva et al. 2018;Ladrera et al. 2019) have already reported significant changes in riverine biological communities due to agriculture and its impacts (e.g., nutrient enrichment, sediment deposition, water abstraction, and riparian habitat degradation).Thus, strategies such as conserving wetlands and riparian zones, which are capable of regulating the energy and material transfer between terrestrial and aquatic ecosystems, can restore the natural capacity of rivers to buffer impacts caused by LULC, hydrological alteration, and climate change (Pusey and Arthington 2003;Palmer et al. 2008;Lind et al. 2019).Relatedly, Spain, a climate change hotspot (IPCC 2021), dedicates approximately 17 million hectares to croplands, of which 23% are irrigated (MAPAMA 2021), making this activity a key economic engine of the country.Hence, to prevent the long-term decline in riverine biodiversity, adaptation strategies and water policy should account for the actual stressors affecting river ecosystems and avoid partial or highly local solutions.Nature-based solutions integrated into a landscape network should be strongly encouraged to help buffer changes in water temperature, hydrology, and diffuse pollution from agroforestry activities.

Fig. 1
Fig.1Study area (505,000 km 2 ) and selected sampling sites: a total of 882 locations (177 for diatoms, 441 for macroinvertebrates, and 264 for fish).See Fig.S1in the Supplementary Material for sampling site distribution in each biotic group

Fig. 2
Fig. 2 Interactions among environmental factors at large-and localspatial scales (gray boxes), anthropogenic pressures (red boxes), and riverine biological communities (Color figure online)

Fig. 4
Fig. 4 Ordinations showing the first two RLQ axes on macroinvertebrate communities.a Site scores according to macroinvertebrate taxa dominance: Ephemeroptera, Plecoptera and Trichoptera (EPT), Diptera, and Odonata, Coleoptera and Heteroptera (OCH), Diptera, and non-insects; b coefficients for key environmental variables (see

Fig. 5
Fig. 5 Ordinations showing the first two RLQ axes on fish communities.a Site scores according to fish subfamily dominance; c coefficients for key environmental variables (seeTable 1 for codes) with

Table 1
Mean, standard deviation, and range of environmental variables describing dynamic (climate, LULC, and hydrology) and static (topography and geology) variables at the diatom, macroinvertebrate, and fish sampling sites.Only uncorrelated variables in each group (Pearson's |r| ≤ 0.7) are shown.Variables refer to the mean catchment values upstream of the sampling site

Table 2
Subset of environmental variables in the RDA models (listed in their decreasing order of contribution) explaining taxonomic and trait-based diatom, macroinvertebrate, and fish community variability.The adjusted coefficient of determination (Radj2) values are shown.See Table1for environmental variable codes.The contribution of each independent variable is shown in TablesS5 and S6in the Supplementary Information NAT 1, CNF, CALC, AGR, CONG, TEMP, OUT, PLT, SI NAT 3, AREA 0.546 Page 7 of 17 95