Assessing the Vulnerability of Aquatic Macroinvertebrates to Climate Warming in a Mountainous Watershed: Supplementing Presence-Only Data with Species Traits

Abstract

Mountainous running water ecosystems are vulnerable to climate change with major changes coming from warming temperatures. Species distribution will be affected and some species are anticipated to be winners (increasing their range) or losers (at risk of extinction). Climate change vulnerability is seldom integrated when assessing threat status for lists of species at risk (Red Lists), even though this might appear an important addition in the current context. The main objective of our study was to assess the potential vulnerability of Ephemeroptera (E), Plecoptera (P) and Trichoptera (T) species to global warming in a Swiss mountainous region by supplementing Species Distribution Models (SDMs) with a trait-based approach, using available historical occurrence and environmental data and to compare our outcomes with the Swiss National Red List. First, we used nine different modelling techniques and topographic, land use, climatic and hydrological variables as predictors of EPT species distribution. The shape of the response curves of the species for the environmental variables in the nine modelling techniques, together with three biological and ecological traits were used to assess the potential vulnerability of each species to climate change. The joint use of SDMs and trait approach appeared complementary and even though discrepancies were highlighted between SDMs and trait analyses, groups of potential “winners” and “losers” were raised out. Plecoptera appeared as the most vulnerable group to global warming. Divergences between current threat status of species and our results pointed out the need to integrate climate change vulnerability in Red List assessments.

1. Introduction

Mountainous running water ecosystems harbour different types of habitats in terms of hydrology regimes and physico-chemical characteristics that are largely determined by the spatial and temporal balance between contrasted water sources (glacier melt, snow melt, groundwater and surface runoff) [1]. This diversity of habitats and their conditions shelters an important biodiversity. These ecosystems are vulnerable to anthropogenic threats (e.g., water abstraction or diversion) but also to climate change, with major changes coming from warming temperatures [2]. Due to warming temperatures and hydrological changes, species assemblages are likely to be modified [3]. Cold-adapted species should suffer from habitat contraction and should show upstream shifts, while warm-adapted species can present an upstream expansion of their altitudinal distribution and generalist species may overall benefit from changes [4]. Species distribution will be affected and some species are anticipated to be winners (increasing their range) or losers (at risk of extinction) [5,6,7,8]. IUCN Red List status [9] based upon population size and geographic distribution criteria reflect notably anthropogenic threats that may have important impacts (through pollution and habitat modification) on species and they are widely used as standards for assessments of species extinction risk. However, climate change vulnerability is seldom integrated when assessing threat status for Red Lists, even though this might appear an important addition in the current context [10,11,12]. In addition, several methods were developed to assess species vulnerability and further work is needed to reach consensus [13].

In this study we propose to assess the potential vulnerability of species to climate changes by supplementing Species Distribution Models (SDMs) with a trait-based approach, using available historical occurrence and environmental data. As far as we know, this kind of approach crossing SDMs and a trait-based analysis has seldom been developed and can provide a valuable contribution to species vulnerability evaluation [14,15].

SDMs relating the occurrence of species to variations in key environmental features are a means of analysing species vulnerability to climate change. SDMs are widespread in ecological modelling [16,17,18,19] and are increasingly applied to freshwater ecosystems [5,20,21,22,23]. Another means of assessing species vulnerability is the use of ecological and biological traits. Different studies demonstrated that species classified as not threatened according to IUCN criteria were, in fact, at high risk of extinction when considering their traits [6,24,25]. Presence-only data typically stored by museums or record centres are typically used in predictive modelling but they often contain scarce information [26], frequently limiting the significance of modelling attempts. Assessment of traits may help to identify the species at greatest risk [27] and could therefore be considered as complementary to SDMs.

We focused on three insect orders: Ephemeroptera (E), Plecoptera (P) and Trichoptera (T). These three orders are worldwide represented in freshwater ecosystems and are an important component of river and stream biota, quantitatively, qualitatively and functionally [28,29,30]. In mountainous regions, they are the major component of aquatic biodiversity as one of the first biotic component that can reach the highest elevations and cope with harsh conditions contrary to fishes that are absent from glacier-fed headwaters [1,31].

The goals of the study were to address the following questions: (i) What are the major environmental drivers of EPT species diversity and distribution in a large catchment (ii) What could be the impact of climate change, specifically an increase in air temperature? (iii) Can a combination of SDMs and trait-based approach lead to more insightful and reliable conclusions on species vulnerability, assuming that species sharing similar suites of traits would tend to show similar trends along key environmental gradients? (iv) Which species or groups of species can be anticipated to be winners or losers under warming temperatures? (v) Are the available threat status of the species provided by the Swiss National Red List [32] in accordance with the outcomes of the SDM and trait analyses?

2. Methods

2.1. Study Area

The selected study area was the Rhone catchment, upstream of Lake Geneva in Switzerland, one of the key case-study areas of the ACQWA project in the Alps [31]. This catchment is at present covered by 11.5% of glaciers [33]. By 2100, it is likely to undergo a winter temperature increase of 4 °C and of more than 6 °C in summer, when compared to 1961–1990 [34]. Winter precipitations are likely to increase and summer precipitations to decrease. Snowfall is anticipated to be more abundant in the higher reaches and reduced at lower levels. Moreover, more frequent geomorphic hazards are expected to occur [35]. Therefore, this catchment is considered as potentially facing major environmental changes including possible hydrological changes in the coming decades [33,36].

The Rhone catchment covers 5338 km² and is located in the South-western part of the Swiss Alps (6°78’–8°48’ E, 45°86’–46°66’ N) (Figure 1a). Its altitude ranges from 370 to 4634 m. It is covered by 71% of forests and pastures, 10% of agricultural lands and 11.5% of glaciers [32]. The population density is about 52 individuals/km² and the highest densities are located in the bottom of the Rhone valley [37]. The Rhone River originates from the Rhone Glacier at 2250 m a.s.l. The natural annual flow regime at the river mouth into Lake Geneva shows a low discharge in winter (about 50–60 m³ s⁻¹ and a higher discharge in summer (about 400 m³ s⁻¹) due to snow and upstream glacier melt [32,38]. However, the Rhone River and its tributaries were highly altered in the 19th and 20th centuries. The Rhone river was subjected to two river channelization (respectively from 1863 to 1894 and from 1930 to 1960) to protect from floods and create space for agriculture and industry [39,40]. Moreover, the hydrology was heavily modified since the fifties due to the implementation of large high-head storage schemes for the production of hydroelectricity [32,38,39]. Consequently an important part of the Rhone network in Switzerland does not show a natural flow regime [32,38,41].

2.2. Species Occurrence Data

The Ephemeroptera, Plecoptera and Trichoptera (EPT) records were provided by the Swiss Biological Records Centre (www.cscf.ch). This centre stores and makes available species occurrence data through a standardized method. In a first step, field observations are sent to the centre or uploaded on its website. Then, a consistency control is performed to detect mistakes in the observation coordinates or encoding transcription. A scientific validation is implemented to confirm the taxonomic identification. The records used in the current study were presence-only data that were recorded between 1991 and 2008. All records were at the species level and concerned larvae as well as adult stages. There were 291 observation points in the catchment that comprised 206 EPT species (32 Ephemeroptera, 57 Plecoptera and 117 Trichoptera) (Figure 1b). Consequently, one observation point could contain several observation dates and species. The number of records per species varied from 1 to 173.

2.3. Species Distribution Models (SDM)

2.3.1. Temporal and Spatial Scales of Environmental Variables

The stream network used in this study was at the 1:25,000 scale (Swiss Federal Office of Topography, Swisstopo). Its total length was 6982 km. In order to associate information at common spatial units, the stream network was divided into 500 m segments. Very small segments (<30 m) were combined to larger ones. Segment lengths varied between 30 and 530 m, with a vast majority of them at 500 m. The total number of river segments was 23,075. EPT occurrences and environmental variables were all associated to a given stream segment. The adult EPT records were associated to the nearest stream segment.

For the stream segments that were associated to a species record, the associated values of time-variable environmental predictors were averaged for the period 1991–2008, except for air temperature, snow and glacier melt, which were associated to the year of observation. For stream segments with pseudo-absences, the environmental variables were calculated for each year and values were randomly associated to stream segments.

2.3.2. Environmental Predictors

An initial set of 19 environmental variables was selected as potential predictors of EPT species distribution. These hydromorphological variables represented the land use, topography and topology, climate, hydrology (including snow and glacier contributions) of the catchment (

Table 1. Summary of the 19 environmental variables at the stream segment scale in the Rhone River catchment. Variables selected for the models *.

	Abbreviation	Variable	Description	Units	Period 1991–2008
Climate					Mean	Range
	AirTemp *	Air temperature	Mean daily air temperature	°C	3.76	−6.77 to 37
Glacier and Snow components
	GlacArea	Glaciated area	Fraction of glaciated area	no	0.13	0–0.69
	GlaMelt *	Glaciermelt	Mean monthly glaciermelt	mm/hr	0.015	0–0.24
	SnoMelt *	Snowmelt	Mean monthly snowmelt	mm/hr	0.09	0–0.57
Hydrology
Magnitude of flow events
	FlowVar	Flow variability	Mean of daily flow variability	m3/s	6.26	0–81.12
	SumFlowVar *	Summer flow variability	Mean of daily flow variability for summer months (June-August)	m3/s	7.93	0–84.2
	Mean	Annual flow	Mean annual flow	m3/s	15.5	0–190.9
	High	High flow conditions	Mean of Q90 of all years/median on whole period	no	14.6	0–419.5
Frequency of flow events	Hispell	High flow spells	Total number of high flow spell on whole period (threshold: 0.9 * mean flow value)	no	127.8	1–410
	Base	Base flow index	Mean annual Ratio of base flow to total flow	no	0.38	0–0.65
Duration of flow events	Zero	Zero flow days	Total number of zero flow days on whole period	days	163.2	0–6575
Timing of flow events	Dure	Duration	Mean annual duration of high spell	days	23.3	1.37–365.2
Rate of change in flow events	Cons *	Constancy	Constancy based on seasonal mean daily flow	no	0.27	0.021–1
	Rise	Rise	Mean annual number of rises	no	49.25	0–80.8
Land use
	Forest *	Forest	Presence of forest in a buffer of 200 m around the river segment	percentage	21.5	0 to 98
	Agric *	Agriculture	Presence of agriculture in a buffer of 200 m around the river segment	percentage	9.5	0 to 95
Topography and topology						372–3332.3
	Sorder *	Stream order	Stream order according to Strahler	class	1–7
	Alti	Altitude	Mean altitude	m	1715.1	372–3332.3
	Slope *	Slope	Mean slope	percentage	33.3	0–240.8

). A major part of environmental variable calculations were made with standard ArcMap 9.3.1 tools (ESRI, Redlands, CA, USA) [42].

• Topography, Land use and Climatic Variables

The stream order [43] was calculated from the 1:25,000 topographic river network using an ArcView script (ArcView 3.3) (ESRI, Redlands, CA, USA) [44] able to cope with braided segments. The stream order ranged from 1 (headwaters) to 7 (downstream part of the Rhone river). Slope and altitude were calculated from the Digital Elevation Model of the Swiss Rhone catchment (Swiss Federal Office of Topography, Swisstopo) and related to the stream network (mean of the variable along each stream segment) with Hawth’s Analysis Tools (Line Raster Intersection Statistics) for ArcGIS [45]. The forest and agriculture variables were extracted from a land use cover layer (Swiss Federal Statistical Office, Land Cover Statistics 92/97). The forest cover combined all types of forested areas (e.g., scattered forests, bushes, undergrowth) and the agriculture cover combined farming, orchards, viticulture and horticulture. Both land use variables were expressed as a percentage coverage in a 200 m-buffer around the stream segments using the Zonal Statistics tool from Hawth’s Analysis Tools [45].

Raw air temperature data came from the ENSEMBLES project [46] using the period of 1991 through 2010 as a control scenario [32]. This data was downscaled (1 km × 1 km) and corrected [47,48] for the Rhone river catchment. This modelled air temperature data was chosen rather than observed data in the purpose of making predictions with climate scenarios in further analyses not covered here. The annual daily air temperature was calculated per raster cell. Raster cell values were then associated to stream segment centroids.

• Hydrological Variables including Snow and Glacier Components

The raw stream discharge, the glaciated area and snow and glacier components were calculated with the distributed rainfall-runoff model TOPKAPI-ETH (Topographic Kinematic Approximation and Integration model-ETH). The version used in this study was a substantially modified version of the original model [49] and included a snow and glacier melt module. The model represented the watershed with regular square grid and required as input data a digital elevation model, soil and land use information as well as air temperature, cloud transmissivity and precipitation time series [32]. The model simulated subsurface flow, overland flow and channel flow using the kinematic wave approximation. As it was not possible to save the discharge for each stream segment, discharge was saved for a number of sub catchments. A total of 298 sub catchment outlets were selected in order to cover the entire Rhone catchment [35], taking into account water withdrawals and river diversions. The discharge was simulated at an hourly scale for the period 1991 to 2008 [32].

Values of simulated snow and glacier melts were expressed as yearly time series for the 298 sub catchments and associated to those of the year of observation of the species records. All the other hydrological variables were expressed for the entire period 1991–2008. The glaciated area was the fraction of glacier calculated in each of the 298 sub catchments.

Two hydrological variables based on the hourly discharge were calculated to express the variation of daily flow (maximum—minimum hourly flow): the range of daily flow (mean of all years for the period 1991–2008) and the range of daily flow for summer months (June, July and August) (mean of all years for the period 1991–2008). Values per sub catchment were also assigned to stream segments as mentioned previously.

Moreover, based on the daily discharges, we selected a number of hydrological variables that expressed the five main aspects of the flow regime: magnitude, frequency, duration, timing and rate of change of flow events [50,51] and that were considered as relevant for various aspects of macroinvertebrate distribution and life-cycle [50,52,53]. Sixteen hydrological variables were initially calculated: the skewness of daily flows, the mean and median annual flow, the mean flow for each season, the baseflow index, the low flow conditions, the frequency of low and high flow spells, the number of zero flow days, the mean duration of high spell, the constancy and the number of rise and falls. TOPKAPI-ETH hourly discharge in each sub catchment was aggregated into daily discharge and then implemented in the River Analysis Package (RAP) [54] to calculate the 16 hydrological variables for each sub catchment for the period 1991–2008. In order to obtain a reduced set of less redundant hydrological predictors, a Principal Component Analysis was performed to eliminate the most correlated variables. This calculation was done using the ade4 package [55] of the R software [56]. Ten hydrological indices were kept (Table 1). For each of the ten variables, the values per sub catchment were assigned to all the stream segments it contained.

• Final Selection of the Environmental Variables

A first selection of environmental variables was made after testing for correlations with a Pearson correlation coefficient (Appendix A) removing variables above with correlation coefficients above 0.8 and selecting the preferred variables based on their ecological relevance (Lehmann et al. 2002). When examining the correlation coefficients and the identity of the environmental variables, a final set of only 9 variables was kept as predictors in the SDMs: air temperature, glacier melt, snow melt, summer flow variation and constancy, forest and agriculture covers, stream order and slope (asterisk in Table 1). Slope was transformed with log(x + 1) to stabilize variances.

2.3.3. Methods included in the Species Distribution Models (SDMs)

EPT species distribution models were calculated using presence-only data supplemented with generated pseudo-absences. Nine different modelling techniques available in the biomod2 Package (version 1.3.5) (UJF, Grenoble, France) [57] were applied. They belong to different families of methods:

-

regression methods: Generalized Linear Models ‘GLM’ [58], Generalized Additive Models ‘GAM’ [59] and Multivariate Adaptive Regression Spline ‘MARS’ [60],

-

machine learning methods: Gradient Boosting Machine ‘GBM’ [61], Artificial Neural Network ‘ANN’ [62], Maximum Entropy ‘MAXENT’ [63] and Random Forest ‘RF’ [64],

-

classification methods (Classification Tree Analysis ‘CTA’ [65] and Flexible Discriminant Analysis ‘FDA’ [66].

Biomod2 was run in the R Software version 2.15.2 [56]. Default settings were used for all methods. SDMs were calculated only for species with at least 10 observations, that is for 63 species out of 206. For each species 10 random partitions of the occurrence localities were made. Each partition was made by randomly selecting 80% of the occurrence localities as training data with the remaining 20% reserved for testing the resulting models. As EPT species records were presence-only data, pseudo-absences (PA) needed to be generated, given that most of the modelling techniques used here only work with presence-absence data. Although PA characteristics (number, location, methodology selection) are suggested to be different according to the method family [67], we used a consensus methodology for all model techniques based on findings of [67]. PA were selected at random and for each species the number of PA was equal to three times the number of observations, reaching therefore a minimum of 40 observations (10 presences, 30 pseudo-absences) [68,69]. PAs could introduce some bias in our analyses and forced us to remain careful regarding our conclusions [67]. Ten different sets of PA were selected and an equal weight was given between presences and PA. To assess the predictive power of each model we used the area under the Receiver-Operating Characteristic (ROC) curve. A ROC value between 0.7 and 0.8 was considered as showing a fair model, a value between 0.8 and 0.9 a good model and higher than 0.9 a very good model. A ROC value under 0.7 indicated a poor model [70].

Explanatory variables contribution was considered. Indeed, biomod2 allowed assessing the environmental variables importance in explaining species distribution [57]. As all models did not rely on the same techniques, biomod2 used a permutation procedure that measured the relative importance of each variable independently of the models. Furthermore, the response curves of the environmental variables for the nine modelling techniques were calculated as in Elith et al. 2005 [71], were used to assess the potential vulnerability of each species to climate change. This method allows to explore the effect of one variable at the time by building a hundred regular steps between its minimum and maximum values, while keeping the other variables constant at their mean values. These resulting predicted response curves allowed an evaluation of the ecological sensibility of our models and consequently of our species.

2.4. Trait-Based Approach

2.4.1. Species Traits

As all the species could not be modelled via SDMs due to low occurrences (70% of the EPT species had less than 10 occurrences), species traits were used as an alternative way to assess vulnerability to climate change for all species. Species traits were extracted from the “freshwaterecology.info” database that gathers information on the autecology of European freshwater organisms [72]. One Plecoptera, Rhabdiopteryx harperi, was not covered by the database. The trait approach was thus developed for 205 species. Three biological and ecological traits were selected: stream zonation, water temperature range and emergence/flight period. The longitudinal distribution of species encapsulates various aspects of their ecology. Stream zonation was therefore used as a surrogate for stream geomorphological characteristics and for temperature/oxygen regime [73]. The water temperature range trait and the emergence/flight period were chosen because of the potential impact of climate change on these features. Stream zonation preferences were coded under 10 categories. Water temperature range preference had three categories and the trait Emergence/Flight period had four categories (

Table 2. Definition and percentage of species documented for each trait and EPT order.

Traits	Categories	Coding System	E (%)	P (%)	T (%)
Stream zonation	eucrenal zone (euc), epipotamal zone (epo) hypocrenal zone (hyc), metapotamal zone (mpo) epirhitrhal zone (erh), hypopotamal zone (hpo) metarhitrhal zone (mrh), littoral zone (lit) hyporhithral zone (hrh), profundal zone (pro)	10 points system	100	100	99.1
Water temperature range	cold stenotherm (cos), water temperature <10 °C warm stenotherm (was), water temperature >18 °C eurytherm (eut), no specific preference	Single category assignment	90.6	96.4	71
Emergence/Flight period	winter (win), spring (spr), summer (sum), autumn (aut)	10 points system	100	100	82.1

).

2.4.2. Trait Analysis

A fuzzy coding analysis [74] was applied to the species x traits matrix in order to calculate trait-based distances between species. It was done with the ade4 package [75] of the R Software [56]. Then a hierarchical cluster analysis using Ward’s method [76] was applied to the between-species distances to group species presenting similar combinations of traits. The simple ’Phenon line’ method [77] was used to retain the number of groups corresponding to the maximal dissimilarity between two successive clustering levels. The incorporation of multiple traits to define functional groups of species allows to account for contrasted life strategies among the species [78]. A vulnerability status was then assigned to each species according to the characteristics of its respective group.

3. Results

3.1. Species Distribution Models

3.1.1. Model Performance

For all modelled species, model performances ranged between 0.66 and 1 according to the ROC evaluation Criteria. Except for two runs, models were from “fair” (0.7 < ROC < 0.8) to very good (ROC > 0.9). Random Forest (RF) ROC distribution was much more homogenous than for the other models.

The average of single species model performances varied between 0.85 and 0.99 (Appendix B). The GBM model did not run for the 30 species with occurrences equal lower than 16. Species with few occurrences did not show worse models than species with higher occurrences.

3.1.2. Contribution of the Variables

Considering all models together (Figure 2), air temperature, stream order and slope were the main variables explaining the distribution of Ephemeroptera, Plecoptera and Trichoptera (EPT) species. Variables importance for each species logically highlighted air temperature, stream order and slope as the main variables explaining species distribution (Appendix C). However, this was not the case for all species. Summer flow variation was the most contributing variable for Ecdyonurus parahelveticus (SumFlowVar: 0.28), Protonemura risi (SumFlowVar: 0.45) and Rhithrogena loyolaea (SumFlowVar: 0.27). Snowmelt appeared as the most explanatory variable for two species: Nemoura obtusa (SnoMelt: 0.46) and Rhabdiopteryx harperi (SnoMelt: 0.53). For Allogamus mendax it was forest cover: Forest (0.32) and for Rhyacophila intermedia it was the agriculture cover: Agric (0.29) (Appendix C).

3.1.3. Assessing Species Vulnerability from SDMs

As air temperature is directly affected by climate change and was one of the most important variables explaining the EPT distribution, we analysed the overall trend of the response curves of the species for this variable across the nine modelling techniques. The vulnerability of the species to increased temperature was then assessed given the shape of their individual response curve.

Three groups of species were proposed on this basis (

Table 3. Grouping of the modelled species according to their response curves for all SDMs and the two overall most contributive environmental variables, Air temperature and Stream order. * indicates species for which air temperature or stream order were not among the most contributive variables. 12 species that had no SDM response to air temperature were left out.

	« Winners ». Occurrence Probability Increases with Air Temperature	« Losers ». Occurrence Probability with a Bell-shaped Response to Air Temperature or Decreases with Air Temperature Increase (Rh. harperi)
Upper course (stream orders 1–4)	Alainites muticus	Melampophylax melampus
	Plectrocnemia geniculata	Philopotamus ludificatus
Median course (stream orders 3–5)		Chloroperla susemicheli
		Cryptothrix nebulicola
		Dictyogenus alpinus
		Drusus biguttatus
		Drusus discolor
		Ecdyonurus picteti
		Epeorus alpicola
		Leuctra teriolensis
		Halesus rubricollis
		Protonemura lateralis
		Protonemura nimborum
		Protonemura brevistyla
Lower course (stream orders 5–7)	Allogamus auricollis	Baetis melanonyx
	Baetis alpinus	Isoperla rivulorum
	Baetis rhodani	Leuctra inermis
	Ecdyonurus venosus	Protonemura nitida
	Leuctra major	Siphonoperla montana
	Perla grandis
	Brachyptera risi	Rhabdiopteryx harperi *
	Protonemura nimborella
	Potamophylax cingulatus
	Rhithrogena hybrida
	Rhithrogena alpestris
	Rhyacophila intermedia *
	Rhyacophila torrentium
No zonation preference	Allogamus hilaris	Ecdyonurus helveticus
	Allogamus mendax *	Leuctra braueri
	Cloeon dipterum	Leuctra rauscheri
	Hydropsyche instabilis	Leuctra rosinae
	Protonemura intricata	Nemoura sinuata
	Rhyacophila dorsalis	Nemurella pictetii
	Rhyacophila pubescens	Perlodes intricatus
	Rhyacophila tristis
	Rhyacophila vulgaris

):

-

“Winners” (24 species): all the species of this group are potentially favoured by an increase in air temperature. Only 2 species in this group were mostly associated with stream upper reaches (stream order ≤ 5). Thirteen species occurred preferably at higher stream order (≥5) and they can be predicted to expand their distribution upstream under warming temperature.

-

“Losers” (27 species): this group contains species presenting either a bell-shaped response or a decreasing response to increasing air temperature. The latter occurred only in the case of the Plecoptera Rhabdiopteryx harperi. Fourteen species in this group showed a preference for low stream orders (≤ 5). These species should either show a contraction (those associated with the lowest stream orders) or an upstream shift of their range, making them vulnerable. In this group, five species associated with the higher stream order (≥5) and seven species with no definite longitudinal preference might be regarded as vulnerable depending on their ability to move upstream in the catchment.

-

Finally, a third group (not presented in Table 3) assembles twelve species that were not significantly related to air temperature in the SDMs.

When looking at the EPT orders separately, 46.7% of Ephemeroptera and 63.2% of Trichoptera species appeared potentially favoured under the anticipated climate change. On the contrary, a majority of modelled Plecoptera species (58.7%) appeared as potential loser under warming temperatures (

Table 4. Number and percentage of species per vulnerability status as identified according to SDMs.

	Ephemeroptera	Plecoptera	Trichoptera	Total
Potential “winners”	7 (46.6%)	5 (17.2%)	12 (63.1%)	24 (38.1%)
Potential “losers”	4 (26.7%)	17 (58.7%)	6 (31.6%)	27 (42.9%)
Not classified	4 (26.7%)	7 (24.1%)	1 (5.3%)	12 (19.0%)

).

3.2. Trait-Based Approach

The cluster analysis allowed gathering species into five groups with different combinations of the selected biological and ecological traits: stream zonation, water temperature and flight period (Figure 3). Group T1 was composed of species that preferred cold waters and were found in the headwaters. Group T2 differed from T1 in the flight period concentrated in autumn. T3 assembled species with a longitudinal preference extended downstream. No data was available concerning the water temperature for the species of this group. T4 was comprised of eurythermal species, also ubiquitous for stream zonation. T5 gathered species associated with warmer waters and the most downstream parts of catchments. Flight periods in T1, T3, T4 and T5 were highly similar, with a peak in summer.

T1 (85 species) and T2 (15 species) appeared the most vulnerable to climate change as they were cold-adapted and found in the headwaters. They can be considered as potential losers under climate change. T3 gathered 31 species found mostly in the median part of the stream zonation. In the absence of temperature preference data, it is impossible to draw conclusions upon this species group. Their distribution and flight period closely resembled those of group T4. The 46 species of T4 appeared as generalists with no restricted traits. The 28 species of T5 could be considered as warm adapted species.

Species of T4 and T5 could be considered as potential winners under climate change. T1, grouping cold adapted species, contained the major part of the Plecoptera species (40 species out of 55) whereas T5, the group of warm adapted species, did not contain any (Appendix D).

In the trait analysis, a majority of Ephemeroptera species (65.6%) was anticipated as potential winners of warming temperature. For the Plecoptera, 78.6% of the species appeared as potential losers and for the Trichoptera, results were less contrasted with 35.9% of winners and 40.2% of losers (

Table 5. Number and percentage of species per vulnerability status as identified in the trait analysis.

Group	Ephemeroptera	Plecoptera	Trichoptera	Total
Potential “winners”	21 (65.6%)	11 (19.6%)	42 (35.9%)	74 (36.1%)
Potential “losers”	9 (28.1%)	44 (78.6%)	47 (40.2%)	100 (48.8%)
Not classified	2 (6.3%)	1 (1.8%)	28 (23.9%)	31 (15.1%)

).

3.3. Confronting SDM and Trait Groupings

The groups of species produced by the SDM and trait analyses were crossed and confronted with the Swiss red list status [36] (

Table 6. Distribution of species according to their taxonomy (E: Ephemeroptera, P: Plecoptera, T: Trichoptera), threat status (RL: Swiss Red List) and grouping obtained from the SDM and trait analyses. Appendix D contains the raw data.

Group	Red List Status				No Red List Status				TOTAL
	Total	E	P	T	Total	E	P	T
“losers” SDM & traits agreement	4	-	3	1	15	2	10	3	19
“losers” by traits (no SDM data)	37	2	14	21	34	4	15	15	71
SDM / traits contradiction “losers” by SDM	1	1	-	-	6	1	3	2	7
“winners” SDM & traits agreement	-	-	-	-	13	5	3	5	13
“winners” or indifferent by traits (no SDM data)	22	4	1	17	63	12	5	46	85
SDM / traits contradiction “winners” by SDM	3	-	1	2	7	1	1	5	10
TOTAL	67	7	19	41	138	25	37	76	205

, derived from the detailed data in Appendix D). For simplification, the four red list status “Critically Endangered” (3 species), “Endangered” (6 species), “Near Threatened” (37 species) and “Vulnerable” (21 species) were amalgamated.

Overall, 90 out of the 205 species considered (i.e., 44%) were considered as potential “losers” under climate change, either by the trait analysis alone or conjunctly by the traits and SDM (Table 6). These 90 species included 41 species already having a threat status on the Swiss red list.

4. Discussion

4.1. Environmental Variables in the SDM

The types of water sources in mountainous regions and their mix driven by hydrological processes and the resulting water temperature determine aquatic macroinvertebrate assemblages [1,79]. In our models, the three main variables explaining species distribution were air temperature, stream order and slope. Stream order and slope were demonstrated to play a major role in the macroinvertebrate assemblage structure [80,81] but they are not affected by climate change, unlike the remaining most important variable, namely air temperature. For aquatic insects, growth and biological processes are particularly dependent on water temperature [82,83] and most of large scale studies on running waters used air temperature as a surrogate for water temperature as both of them tend to be strongly correlated [84]. Indeed, the implementation of large-scale water temperature monitoring in large catchment is generally difficult to carry out. In the case of headwaters, temperature is dependent on the contribution of water type (snow melt, glacier melt, groundwater and surface runoff) although local stream reach conditions can also play an important role. Mean groundwater temperature generally reflects the annual mean air temperature and is consequently less correlated to current air temperature [85,86]. As noticed by Domisch et al. [5], this could be misleading when modelling the distribution of EPT species associated with the crenal zone such as the Trichoptera Plectrocnemia geniculata McLachlan, 1871 [87].

Contrary to our expectations, hydrological variables appeared only secondary. Even though similar findings were made [5] and showed that climatic and topographic variables were the most contributive in explaining aquatic macroinvertebrate distribution compared to hydrological variables, our results may be biased because of the resolution of the data used for hydrological models. Aquatic macroinvertebrates achieve either active or passive dispersal but travel distances have rarely been estimated [88]. Elliott [89] provided ranges of distance travelled between 3.5 to 13.5 m/day. With such distances, we can assume that EPT are influenced by very local environmental conditions. However, data such as the glacier melt component at fine resolution for large catchments is scarce and challenging to obtain from hydrological models which typically use coarser spatial resolutions [90,91].

4.2. Can a Trait-Based Approach Supplement SDMs?

The choice was made here to supplement SDMs with a trait-based approach in order to consider all the available species. Two groups of “winners” and “losers” were identified in each type of analysis. A majority of winners identified in the SDMs were also found in the winner group according to the trait-based approach. Similarly, a majority of losers identified in the SDMs were found in the loser group of the trait analysis. Therefore, both techniques broadly speaking converged. However, some discrepancies (notably for Trichoptera species) were observed. For instance, the Trichoptera Plectrocnemia geniculata appeared as winner according to SDMs but as loser according to the traits. For this crenal species [87], the use of air temperature as surrogate for water temperature in SDMs may have led to erroneous modelling. As the stream zonation preference of this species is known without confusion, confidence should be put on the trait results and this species should be considered as a loser of climate change. Such an outcome can only be decided on a case-by-case basis. Mismatches can also be explained by the fact that precise ecological traits of aquatic macroinvertebrates and especially thermal preferences are missing [4]. In the European freshwater organisms database [92], the temperature range for EPT species is coded with three categories: cold stenothermic, warm stenothermic and eurythermal. No category exists for cool-adapted species. Consequently, mismatches could be also due to the fact that some species have thermal preference in cool waters (between 10 and 18 °C) or for species that have either cold or warm preferences but can also be found in cool waters. This can highlight the difficulty to categorize the ecological preferences of tolerant species showing large distribution patterns. For example, the Ephemeroptera Baetis alpinus is considered as eurythermal (cold stenothermic) [93]. In some Mediterranean islands, colonization of this species seems restricted by a maximal summer temperature just below 20 °C [94]. In the European freshwater database, Baetis alpinus is considered as a cold stenothermic. Similar observations can be done for other species. Trichoptera Rhyacophila tristis and Rhyacophila vulgaris appear as mesostenothermic (preference between 0 and 18 °C) [95]. In the European freshwater database, they are also considered as cold stenothermic. These three species are highlighted as winners in the SDM analysis but as losers in the trait analysis. Implementing a cool-adapted category in the thermal preferences of species could discriminate those species as winners of climate change in a trait-based approach. In essence, more investigations have to be carried out to discriminate more precisely species thermal preferences in order to obtain a finer analysis on species vulnerability to climate change.

According to SDMs results, twelve species were not influenced by air temperature. In the trait-based analysis, these species were distributed heterogeneously in the different groups. For four of them (Ecdyonurus parahelveticus, Nemoura obtusa, Protonemura risi and Rhyacophila loyolea), the probability of occurrence increased with the increase of summer flow variation or snowmelt. Climate simulations for the Alps in 2100 demonstrated that glacier mass could be reduced from 50 to 90% and snow volume from 35 to 90% (depending on the altitude) leading to a reduction in glacier and snow melt [34,96]. These species should consequently be disadvantaged by the end of the 21th century. For the other species of this set, complementary investigations need to be performed to identify the main drivers of their distribution.

Stream zonation encapsulates stream characteristics from the catchment headwaters to the lowlands. However, in aquatic macroinvertebrate trait databases, glacier-fed rivers—the kryal component of the catchment—is not represented. As already mentioned, glacier-fed rivers are very harsh environments where aquatic life is limited and biodiversity increases with an increase of the distance to the glacier [3,31]. However, some species, such as the Trichoptera Rhyacophila angelieri endemic of the French Pyrenees [3] adapted to glacier-fed rivers, are at risk of extinction due to shrinking glaciers. For a better assessment of mountainous aquatic insects vulnerability to warming temperatures, the kryal component of the discharge should be taken into account in stream zonation as some aquatic insects are predicted to disappear with the complete melt of glaciers [31].

Even though there are insights to gain when using SDMs for the modelling of aquatic species (more specific environmental variables and a finer resolution) and when using a trait-based approach (more accurate data), the combined use of both techniques seems beneficial since it can enable confirmation of the level of vulnerability of a species or a group of species to climate change.

4.3. Winners and Losers

When considering species submitted to both analyses and species only used in the trait based approach, it appeared that the majority of loser species belong to the Plecoptera. Mountainous regions hold an important diversity of Plecoptera with a large proportion of endemic species restricted to spring-fed headwaters [85]. But Plecoptera, with narrow ecological requirements and limited dispersal abilities [29] are also one of the most endangered groups of insects due to pollution and habitat disturbance. Within the three studied orders, Plecoptera seem to be at the front line of groups threatened by warming temperatures. EPT species dependent on the glacier-melt contribution should also be included among the most threatened as climate change could lead to an almost total disappearance of glaciers [31,32]. None of the species of the T5 group in the trait analysis were modelled by SDM. Indeed, species in this group—most of them belong to the Trichoptera order—are warm-adapted and are found rather in the potamal stream zonation or the littoral zone. Consequently, their occurrences in the Swiss Rhone stream network are very scarce as those species are in fact lowland and/or almost limnephilous species. For instance, some of the species belonging to this group were found in the French part of the Rhone River [97,98,99]. Such warm-adapted species as well as generalist species (group T4) will benefit from climate change and might progressively shift upstream. They can be considered as part of the future possible “colonizers” in the Swiss Rhone catchment.

On a long term, we can anticipate that due to the homogenization of habitats (progressive disappearance of glaciers) and as common lower altitude species will replace more specialized headwater species, the regional (γ) diversity should be reduced [100,101].

Species ability to adapt to climate change is not only dependent on their thermal preferences. Species intrinsic characteristics such as dispersal abilities, degree of habitat specialization and population trend will interfere [5,6]. Moreover, hydrological changes (directly induced by climate change or due to water abstraction and hydropower management) will likely change the structure of physical and biological gradients in the Swiss Rhone catchment and add additional uncertainties.

4.4. Conservation Perspectives

At present, species extinction risk is commonly assessed through IUCN Red List Status [9]. IUCN criteria focus on population size and geographic distribution reflecting threats on habitat area and/or quality. In running water ecosystems, 65% of global river discharge and the associated aquatic habitats are under moderate to high threat [102]. In the Swiss Rhone Catchment, major threats to habitat area and/or quality are changes to natural hydrological regimes due to channelization, water abstraction and hydroelectricity production [36,41]. These modifications will remain in future years as well as the vulnerability of species to habitat degradation. However, climate change and notably warming temperatures will have an impact on species distribution [103,104,105] and also on their vulnerability due the combined effect of anthropogenic pressures and climate change [106]. The mixed use of SDMs and trait-based approach allowed us to gain additional information on the spectrum of vulnerability of species to global warming. Some groups of species can be highlighted as a result of this analysis (Table 6):

-

Four species were at the same time considered “Near Threatened” in the Swiss red list and identified as potential “losers” in both the SDM and trait analyses: Leuctra rauscheri (P) Nemoura sinuata (P), Siphonoperla montana (P), Cryptothrix nebulicola (T)

These species can be identified as the most at risk given their pre-existing status and the agreement of both analytical methods.

-

Fifteen species, among which are ten Plecoptera, do not currently have a threat status in Switzerland but were identified according to the trait and SDM analyses as being potential “losers” under climate change: Ecdyonurus picteti (E), Epeorus alpicola (E), Chloroperla susemicheli (P), Dictyogenus alpinum (P), Isoperla rivulorum (P), Leuctra braueri (P), Leuctra rosinae (P), Leuctra teriolensis (P), Perlodes intricatus (P), Protonemura brevistyla (P), Protonemura lateralis (P), Protonemura nimborum (P), Drusus discolor (T), Halesus rubricollis (T), Melampophylax melampus (T)

These species certainly deserve consideration as potentially at threat under climate change and care must be taken that no additional pressure is put upon their natural habitat.

-

Thirteen species (among which only three Plecoptera) have no current threat status in Switzerland and were identified as potential “winners” under climate change by both methods: Baetis muticus (E), Baetis rhodani (E), Cloeon dipterum (E), Rhithrogena alpestris (E)m Rhithrogena hybrida (E), Brachyptera risi (P), Leuctra major (P), Protonemura intricata (P), Allogamus auricollis (T), Hydropsyche instabilis (T), Potamophylax cingulatus (T), Rhyacophila torrentium (T), Rhyacophila dorsalis (T)

In addition to the groups listed above, for which both the SDM and trait analyses provided concordant results, other sets of species can be put forward:

-

Seventy-one species were considered as potential “losers” under climate change by the trait analysis (see Appendix D). Thirty-seven of them also had a threat status in the Swiss red list (including two of the three “Critically Endangered” species and three of the six “Endangered”). As no SDM were available in this set, particular attention should be paid to these species, including efforts to gather more precise information about their ecological requirements.

-

Twenty-two species with a threat status in Switzerland were not identified as potential losers, neither by the trait, nor by the SDM analyses. Some of them could even be considered as potential “winners” in the trait analysis (see Appendix D). This seems to indicate that these species are vulnerable because their habitat is generally under threat but potentially not because of climate change.

Among the 56 Plecoptera species recorded for the Swiss Rhône catchment, 42 (i.e., 75%) could be considered as potential “losers” with the same criteria, making this group the most vulnerable of the three under consideration.

5. Conclusions

Besides stream order and slope, air temperature, which is directly related to climate change, played a major role in explaining the distribution of EPT species. Therefore, we confirm the very likely important impact of global warming upon this key component of freshwater biodiversity.

The joint use of SDMs and trait approach appeared complementary. First, it enabled to confirm the level of vulnerability of EPT species to climate change. Secondly, for species with too few occurrences to be modelled in the SDMs, the trait analyses allowed identification of potential colonizers of the Swiss Rhone catchment under warming temperatures.

Even though discrepancies were highlighted between SDMs and trait analyses, groups of potential “winners” and “losers” were raised out and Plecoptera appeared as the most vulnerable group to global warming. We can expect a loss of EPT species at the scale of the Swiss Rhone catchment. An important complementary analysis would be to assess impacts of climate change upon the functional diversity of the biota.

An important point for re-assessing species vulnerability would be to identify, between habitat alteration and warming temperatures, which one is the more threatening [6]. Thus, an improvement to Red List Status would be to account for species vulnerability to climate change [10,12]. Species disadvantaged by climate change would then appear as threatened and species both threatened by habitat destruction and climate change would appear even more vulnerable.