# Climate and land-use driven reorganisation of structure and function in river macroinvertebrate communities

Théophile Mouton, Fabien Leprieur, Mathieu Floury, Fabrice Stephenson, Piet Verburg, Jonathan Tonkin

<figure><img src="image_1.png" /></figure>

## To cite this version:

Théophile Mouton, Fabien Leprieur, Mathieu Floury, Fabrice Stephenson, Piet Verburg, et al.. Climate and land-use driven reorganisation of structure and function in river macroinvertebrate communities. Ecography, 2022, 2022 (3, e06148), 10.1111/ecog.06148. hal-03601039

HAL Id: hal-03601039

https://hal.science/hal-03601039

Submitted on 2 Jun 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
ECOGRAPHY

Research

Climate and land-use driven reorganisation of structure and function in river macroinvertebrate communities

Théophile L. Mouton, Fabien Leprieur, Mathieu Floury, Fabrice Stephenson, Piet Verburg and Jonathan D. Tonkin

T. L. Mouton (https://orcid.org/0000-0003-2571-5777) ✉ (theophilem@outlook.fr) and F. Leprieur, MARBEC, Univ. Montpellier, CNRS, Ifremer, IRD, Montpellier, France. TLM also at: FRB – CESAB, Inst. Bouisson Bertrand, Montpellier, France. FL also at: Inst. Universitaire de France, Paris, France. – M. Floury, Univ. Lyon, Univ. Claude Bernard Lyon 1, CNRS, ENTPE, UMR 5023 LEHNA, Villeurbanne, France. – F. Stephenson and P. Verburg, National Inst. of Water and Atmospheric Research, Hamilton, New Zealand. – J. D. Tonkin (https://orcid.org/0000-0002-6053-291X), School of Biological Sciences, Univ. of Canterbury, Christchurch, New Zealand.

Ecography

2022: e06148

doi: 10.1111/ecog.06148

Subject Editor: Luis Mauricio Bini

Editor-in-Chief: Miguel Araújo

Accepted 20 December 2021

<figure><img src="image_3.png" /></figure>

Understanding temporal changes in the composition of species communities over spa-
tial and temporal scales relevant to conservation management is crucial for prevent-
ing further biodiversity declines. Here, we assessed patterns and potential drivers of
taxonomic and functional temporal β diversity over 26 years (1991–2016) of 64 river
macroinvertebrate communities, and the length of New Zealand (37°00′N, 46°00′S).
We further examined changes in population size and range shifts of species pools,
and related these to taxonomy and functional traits. We found increasing climate
and land-use driven differences in both the taxonomic and functional composition of
communities over time, coupled with poleward species colonisations and increasing
extirpations in northern locations. Increases in population and species range size were
more prevalent than decreases in population and range size. Species shifted their ranges
towards higher latitudes on average by 50 km per decade. Despite little to no relation-
ship with taxonomy, we uncovered distinct relationships between functional traits and
population trends and latitudinal species range shifts. Species with a high number of
reproductive cycles per year and long-life duration of adults tended to increase their
population size, while larger size species with a high number of descendants per repro-
ductive cycle tended to shift their range towards more southern latitudes. Our results
suggest that the intensity of disturbances, the geographic location of individuals and
communities, and species ecological and functional characteristics, are major determi-
nants of riverine biodiversity reorganisation in the Anthropocene.

Keywords: biodiversity, climate change, freshwater macroinvertebrates, functional diversity, land use change, population trends, species range shifts, temporal β diversity

<figure><img src="image_4.png" /></figure>

www.ecography.org

© 2022 The Authors. Ecography published by John Wiley & Sons Ltd on behalf of Nordic Society Oikos

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.
**Introduction**

Current rates of biodiversity change correlate with recent human-induced climate change and land-use activities (IPCC 2014, Brondizio et al. 2019). These biodiversity changes are predicted to increase as human activities intensify (Urban 2015, Johnson et al. 2017), reflecting local extirpations and/or colonisations of species and populations across landscapes (Tilman et al. 1994, Maclean and Wilson 2011, Brondizio et al. 2019). Such demographic shifts could potentially lead to major alterations to the functioning of ecosystems. However, the responses of species likely vary in magnitude among regions of the world and across latitudes (Heino et al. 2009, Lenoir et al. 2019). Biodiversity trends are, however, overly complex, with observations of local biodiversity increases sometimes contradicting reports of a global climate-driven biodiversity crisis (Antão et al. 2020, Li et al. 2020, Outhwaite et al. 2020). There is therefore an imperative need for biodiversity time-series studies over wide and understudied geographical regions to inform conservation planning and policy (Olden et al. 2018, Magurran et al. 2019).

Uncovering the drivers of biodiversity change requires tools that can deconstruct the mechanisms underpinning such change. Research on changes in the composition of species communities (changing β diversity) has typically focused on temporal changes in spatial β diversity (Olden et al. 2018). However, ongoing temporal change in species composition is also a pressing (Kuczynski et al. 2018, Blowes et al. 2019, Antão et al. 2020), yet less investigated aspect of biodiversity change. Temporal β diversity can capture the biogeographical, ecological, functional and ecological processes involved in temporal changes in communities (Magurran et al. 2019), including local extirpations or colonisations (Legendre and Gauthier 2014, Shimadzu et al. 2015, Gotelli et al. 2017, Legendre 2019, Magurran et al. 2019).

The responses of species to environmental change, including local losses and gains associated with temporal β diversity, are underpinned by their ecological and life-history traits (Dawson et al. 2011, MacLean and Beissinger 2017, Pacifici et al. 2017, McLean et al. 2018a). For example, populations of ecological specialists, with long generation times, and low reproductive output are predicted to be highly vulnerable to climate change (Vié et al. 2009, Chin et al. 2010, Conti et al. 2014). Similarly, the degree to which species reduce, expand and/or shift their ranges may depend on their ability to disperse and the geographic location of populations (Comte et al. 2014, Lenoir et al. 2019). Theoretically, the lowest latitudinal margins of continents or islands are more likely to experience loss of taxa under climate change, because colonisations from lower latitudes are not possible (Heino et al. 2009, MacLean and Beissinger 2017).

In running waters, temporal β diversity is increasingly being used as a tool to understand the dynamics of change in macroinvertebrate communities, including identifying sites that are pivotal for maintaining biodiversity at the landscape scale (Ruhí et al. 2017). Changing climate or flow regimes has been identified as a common driver of change in community

structure. For instance, Cañedo-Argüelles et al. (2020) found a strong relationship between structural changes in macroinvertebrate communities and changes in precipitation regimes. Similarly, Crabot et al. (2020) found that the temporal variability of community structure was related to the frequency and duration of drying events. However, few studies have examined temporal β diversity of stream invertebrate communities over large spatial and temporal scales.

Here, we complement recent findings from Mouton et al. (2020), who reported decreasing spatial taxonomic β diversity (taxonomic homogenisation) but increasing spatial functional β diversity (functional differentiation) among river macroinvertebrate assemblages. We capitalise on the same high-resolution time-series datasets, comprising macroinvertebrate communities collected annually from 1991 to 2016, at 64 mainstem river sites across New Zealand's two mainland islands. The two studies are complementary in the sense that Mouton et al. (2020) explored how climate change affects the spatial organization of assemblages while the present study aims at evaluating how both climate and land-use changes influence temporal differences in community composition within sites.

We measured taxonomic and functional temporal β diversity over the period 1991–2016, which we decomposed into indices of colonisations and extirpations. We tested for the effects of climate and land-use change, nested within hierarchically organized environmental spatial scales, in mediating taxonomic and functional temporal β diversity of these communities. Finally, we examined changes in species' population and range size and latitudinal range shifts which we related to their taxonomic affiliation and a set of functional traits.

Given widespread evidence of freshwater macroinvertebrates response to ongoing climate change (Pyne and Poff 2017, Floury et al. 2018, Mouton et al. 2020), we first expected a latitudinal pattern in temporal β diversity ($E_1$), owing to climate-driven species range shifts. Specifically, we expected ($E_{1a}$) increasing species colonisations at the leading edges (southern margins of New Zealand) but ($E_{1b}$) increasing extirpations at the rear edges (northern margin of New Zealand). Second, given globally observed species population declines and poleward range shifts (Urban 2015, Olden et al. 2018, Lenoir et al. 2019), we expected ($E_2$) greater decreases in species' population and range sizes than increases, and poleward range shifts. Finally, given the relationship between the taxonomic identity of species, their functional traits and demographic responses to environmental change (MacLean and Beissinger 2017, Pacifici et al. 2017, Daskalova et al. 2020), we expected a relationship between trends in population sizes and species range shifts with the taxonomic identity of species and with their functional traits ($E_{3a}$ and $E_{3b}$, respectively).

**Methods**

**Data acquisition**

**Biodiversity data**

Macroinvertebrate communities were sampled from mainstem rivers, at 64 wadeable sites (mean Strahler stream
order=6; min=3, max=8), located in 35 catchments of New Zealand, between latitudes 46 and 35°S (Supporting information). Surveys were conducted once a year, during late austral summers (February–April) from 1991 to 2016. These surveys were conducted for New Zealand's National Rivers Water Quality Network (NRWQN, Smith and McBride 1990), which is operated and maintained by the National Institute of Water and Atmospheric Research (NIWA). Samples were collected following a standardized protocol (Smith and McBride 1990) and under baseflow conditions ($Q < Q_{median}$). Seven Surber samples (0.1 m² and 250 µm mesh net) were collected on all sampling occasions during which macroinvertebrates were removed from a 0.1 m² area in the sampler down to a depth of ca 10 cm and from as many substrate types as possible. Individuals were later identified in the laboratory, to the lowest practicable taxonomic level (species=26%; genus=47%, family=21%; Quinn and Hickey 1990). The same taxonomic resolution was maintained throughout the entire period to allow analyses of long-term changes (Scarsbrook et al. 2000).

We described each macroinvertebrate taxon (n=113) using sixteen functional traits related to morphology, life-history, dispersal strategies and resource acquisition methods extracted from the New Zealand freshwater macroinvertebrate trait database (NIWA, <https://niwa.co.nz/sites/niwa.co.nz/files/nz trait_database_v19_2_18.xlsx>). This trait database has been explicitly developed for New Zealand's standardised freshwater macroinvertebrate sampling protocols (Dolédec et al. 2006, Doledec et al. 2011). Functional traits were divided into 59 modalities and fuzzy-coded from 0 to 3 (Chevenet et al. 1994).

## Environmental predictors

We defined a set of environmental predictors at four different spatial scales, commonly identified as prominent scales operating in river networks (Poff 1997, Allan 2004): the global, regional, catchment and reach scales.

Global-scale predictors were defined as temporal changes in air temperature and precipitation (Brown et al. 2013). At each site, we extracted daily values for the period 1991–2016, from 5 km² gridded layers of New Zealand, using NIWA's Virtual Climate Station Network (VCSN; <https://data.niwa.co.nz/>). We used the slopes of linear regressions between climate variables and years as estimates of rates of change (OLS; the *lm* function in R ver. 4.0.2; <www.r-project.org>). The following variables were used: annual and seasonal (for the winter, spring and summer seasons) mean precipitation, precipitation seasonality (i.e. the coefficient of variation × 100; Fick and Hijmans 2017), mean air temperature and air temperature seasonality (i.e. the standard deviation of the mean). We applied a principal component analysis (PCA; the *dudi*.pca function in the Ade4 package ver. 1.7-15 Dray and Siberchicot 2020) individually to our sets of trends in air temperature and precipitation variables. Based on the correlations among each trend in climate variable and the axes of the PCAs, we created synthetic indices of climate change using the first two axes of each climate PCA. In the PCA of

trends in air temperature variables, the first axis (35.9% of variation explained) described a gradient of increasing mean air temperature (TMean) and the second (33.0%) described changes in air temperature seasonality (TSeas). For precipitation variables, the first axis (28.2%) described changes in precipitation seasonality (Prec CV) and the second (22.7%) described changes in mean precipitation (Prec).

Regional-scale predictors were 1) altitude (metres above sea level) of the sampling site, 2) phosphorus, 3) calcium concentrations and 4) mean hardness (induration) of surface rocks of the upstream catchment (respectively the variables USPhosphorus, USCalcium and USHardness extracted from Leathwick et al. 2010). The latter three variables are descriptors of catchment geology, hence likely reflect regional bio-geochemical characteristics.

For catchment-scale predictors, we used descriptors of changes in catchment land-use and catchment hydro-morphology. For land-use, we used changes (1990–2012) in 4 land-cover types (defined as the proportion of catchment occupied by combinations of 1) high producing grassland, 2) shrub/grassland, 3) plantation forest and 4) non-plantation forest (Landcare Research 2015, Julian et al. 2017). These data were available for the year 1990 and 2012, we used the difference in land-cover between these two dates as estimates of temporal change, following Julian et al. (2017). Similarly, we also used changes (1990–2012) in catchment stock unity density (SUD) of dairy, beef, sheep and deer (SUhای-1; Statistics NZ (territorial authority), Julian et al. 2017). Changes in land-cover and changes in stock unity density were synthesised using the first two axes of two PCAs, which we interpreted as synthetic predictors of changes in land-cover (LC PC1, 61.4% and LC PC2, 23.9%; Supporting information) and changes in stock unity density (SUD PC1, 44.3% and SUD PC2, 35.2%; Supporting information), respectively. Catchment hydro-morphology was given by the area of the upstream catchment (m²; USCatchArea) and the average slope of the upstream catchment (degrees; USAvgSlope) from each river segment (Leathwick et al. 2010).

For reach-scale descriptors, we used: 1) the slope (degrees) of the stream segment at each sampling site (SegSlope; Leathwick et al. 2010), 2) the predicted wetted river width (m) at the 7-day mean annual low flow (WidthMALF; Booker 2015), 3) the estimated proportion of riparian shading at each river segment (SegRipShade; measured from satellite imagery by Leathwick et al. 2010), 4) temporal changes in water quality, 5) temporal changes in flow and 6) temporal changes in substrate size. Temporal changes in water-quality, flow and substrate size were all estimated as slopes of variables (described below) and years. For temporal changes in water-quality, we compiled data from 1991 to 2016 of median annual values (mg m⁻³) of nitrate (NO₃-N), ammoniacal-nitrogen (NH₄-N), dissolved reactive phosphorus (DRP) and dissolved oxygen (DO) concentrations and water clarity (metres; CLAR). These data originate from samples collected monthly, at the same sites as those sampled for macroinvertebrate communities (Davies-Colley et al. 2011). We log₁₀-transformed water-quality variables and calculated temporal
changes for each variable at each site. We ordinated temporal
changes in water-quality using PCA and retained the first
axis (33.3%) as a synthetic indicator of temporal changes in
water quality (WQ). For changes in flow, we compiled mean
daily flow data from hydrological gauges located at each site,
which we divided by upstream catchment area to obtain a
measure of run-off per day (Vaughan and Gotelli 2019). We
then calculated mean annual and seasonal (for winter, spring
and summer) flow and flow coefficient of variation, as was
done for precipitation. We ordinated temporal changes using
PCA and kept the first two axes (Flow PC1, 54.1% and Flow
PC2, 20.6%; Supporting information) as indicators of tem-
poral changes in flow. Finally, for changes in substrate size
we used temporal changes in the substrate size index (SSI;
Jowett et al. 1991). Substrate composition was measured by
randomly selecting 100 particles at 1-m intervals along a path
of 45 degrees to the riverbank in a zig-zag manner. Particles
were assigned to each of 8 size classes: bedrock, boulders (>300 mm), large cobles (300-128 mm), small cobles (128-64 mm), large gravel (32-64 mm), small gravel (2-32 mm), sand (62.52 µm-2 mm) and silt (<62.52 µm). We measured the substrate size index (SSI) for each sampling occasion fol-
lowing Jowett et al. (1991) (1=silt only, 8=bedrock only).
Temporal changes in SSI were then measured for each site.
Detailed description of the variables and statistical outputs
in regards to the environmental descriptors used for analy-
ses are given in Mouton et al. (2020) and in the Supporting
information.

Statistical analyses

Temporal β diversity indices

We used the temporal β diversity index from Legendre
(2019) (TBI; the *TBI* function in the *adespatial* package ver.
0.3-8) to measure temporal changes in the taxonomic and
functional composition of each of the 64 communities. This
index has been specifically developed 1) to calculate differ-
ences in assemblage composition (dissimilarity) between two
time periods from pairwise distances, and 2) to decompose
these temporal differences into indices of local colonisa-
tions and extirpations. Macroinvertebrate abundances were
log(x + 1)-transformed and pairwise distances in community
composition were calculated using the percentage difference
index of dissimilarity (%diff; Odum 1950, also known as the
Bray–Curtis index). To estimate temporal changes in func-
tional composition (i.e. the composition of species functional
traits within each community), the species-by-site matrix was
replaced by a trait-by-site matrix using the community-level
abundance weighted means of functional trait values (CWM;
Lavorel et al. 2008; using the *dbFD* function in the FD pack-
age ver. 1.0-12; Laliberté et al. 2014).

We used the first year of sampling (1991) as a baseline for
each time-series, and then successively compared it to each of
the following years (1992, onwards (Magurran and Henderson
2010, Dornelas et al. 2014, Antão et al. 2020)). This compu-
tation therefore yielded a value of dissimilarity, colonisations
and extirpations for each year (except the baseline year) and

site. For each site, we regressed values of each index (taxo-
nomic and functional dissimilarity, extirpations and coloni-
sations, respectively) against years and interpreted the slopes
from the regression models (multiplied by ten) as a measure
of trends per decade. This method allows to examine whether
changes in composition of each assemblage (relative to the
baseline year) consistently increase (positive slope) or decrease
(negative slope) over time (Antão et al. 2020).

To test our first expectations ($E_{1a}$ and $E_{1b}$), we tested for
relationships between trends in temporal $\beta$ diversity, includ-
ing colonisations and extirpations with latitude using linear-
mixed effects models (the *lme* function in the *nlme* package
ver. 3.1-152 (Pinheiro et al. 2017)) setting island (North
Island versus South Island) as random effect.

Drivers of temporal β diversity

We performed hierarchical generalised additive mixed effects
models (GAMMs; the *gamm* function in the *mgcv* pack-
age ver. 1.8-33 (Wood and Wood 2015)) to relate temporal
changes in each of the six indices of trends in temporal $\beta$
diversity to our set of environmental descriptors. Predictor
variables were Box–Cox transformed prior to analyses (Box
and Cox 1964; the *BoxCoxTrans* function in the caret pack-
age ver. 6.0-84 (Kuhn et al. 2020)), and standardised to zero
mean, one unit variance (using the *decastand* function in the
vegan package ver. 2.5-6 (Oksanen et al. 2014)).

GAMMs were fitted with catchments nested within
islands as random effects to account for spatial structure
in data (Dormann et al. 2007). We used regression splines
to account for potential non-linear relationships, restricted
maximum likelihood (REML) to optimize the parameter
estimates and assumed Gaussian-type distribution errors.

Each response variable was first modelled against global-
scale descriptors. We kept only the global-scale descriptor(s)
that maximised the coefficient of determination (adjusted-
R²) as the best model following Van Looy et al. (2017) and
Floury et al. (2018). The same step was then repeated in a
descending way, using successively regional, catchment and
reach-scale descriptors as predictor variables and the residuals
from the previous model as response variable. At each step,
we quantified the percentage of relative importance of each
variable in the model following methods described in Kuhn
(2008) (the *varImp* function in caret). We tested for resid-
ual spatial autocorrelation in the residuals of the reach-scale
model using Moran's index of spatial autocorrelation (the
*Moran.I* function in the ape package ver. 5.5 (Paradis et al.
2019)). Residual spatial autocorrelation was not observed (p
> 0.05 for every model).

Population trends and species range shifts

To test our second expectation ($H_2$), we assessed tempo-
ral trends in the abundance, range size and the latitudinal
distribution of 83 taxa that were recorded for at least 10
years (not necessarily consecutive; following Dornelas et al.
2019). For population trends (hereafter referred to as tem-
poral trends in the abundance of each taxa), we did not
include the time when a taxon was absent in the time series,
because this would tend to flatten the slope towards zero. Macroinvertebrate abundances were log(x + 1) transformed and then standardised, so that each time-series had a mean of 0 and a standard deviation of 1. This transformation put all time-series into common units that are more appropriate for comparisons of taxa with disparate population sizes. We estimated population trends by fitting a linear regression model of abundances for each taxon individually, against years.

Second, we measured the range size of each taxon as the convex hull area encompassing the sites where the taxon was present each year. We regressed range sizes against years using ordinary least squares models. Finally, we assessed latitudinal distribution shifts of each taxon, by regressing the centroid of its range against years. We used the slope of each linear regression model as an estimate of changes in population and range size and of latitudinal range shifts, for each taxon.

## Relationships with taxonomy and functional traits

To test our third expectations ($H_{3a}$ and $H_{3b}$), we started by relating changes in community composition to the taxonomic classification of our taxa ($H_{3a}$). For this purpose, we compiled a taxonomic classification of the studied taxa (using phylum, subphylum, class, subclass, family, subfamily, order, genus and species names), from the National Center for Biotechnology Information (using the taxonomy function in the myTAI package Drost et al. 2018). From this classification, we created a taxonomic tree by calculating taxonomic distances between taxa (using Gower's distance and the taxa-2dist function in vegan). We used hierarchical clustering (the hclust function in R) to produce a tree from the taxonomic dissimilarity matrix. We checked for potential relationships between taxonomic distance and Euclidean dissimilarity in population size and/or species range shifts using a Mantel test (Mantel 1967) with 9999 permutations (the mantel.rtest function in Ade4).

To relate changes in community composition to functional traits of macroinvertebrates ($H_{3b}$), we tested for relationships among population trends, changes in range size and latitudinal range shifts with the axes of a functional trait space (Loiseau et al. 2020). For this purpose, we built a functional trait space of the macroinvertebrate taxa by ordering the species-by-trait matrix in a multidimensional space, using fuzzy correspondence analysis (the dudi.fca function in Ade4). We ensured that equal weights were given to each of the trait categories (e.g. for the trait 'maximum potential size' which contained 5 categories, each category was given a weighting of 1/5). We retained the first three axes of the functional space as synthetic traits of macroinvertebrates (cumulative percentage of variation explained = 31.8%). We used non-parametric Spearman's rank correlation tests to relate changes in population and/or range size and latitudinal range shifts, respectively, to each axis of the trait space and mapped significant relationships on trait space to provide a better visual assessment of potential relationships. Statistical analyses were all performed under the R environment (<www.r-project.org>).

# Results

## Spatial patterns of changes in temporal β diversity

Trends in taxonomic and functional temporal β diversity increased on average across New Zealand over the 26-year time-series (taxonomic dissimilarity: mean rate = 0.03, minimum = -0.08, maximum = 0.17; functional dissimilarity: mean = 0.02, min = -0.07, max = 0.14; Fig. 1). This increasing trend was driven by increasing taxonomic and functional colonisations and functional extirpations (taxonomic colonisations: mean = 0.04, min = -0.01, max = 0.2, functional colonisations: mean = 0.01, min = -0.03, max = 0.06, functional extirpations mean = 0.01, min = -0.04, max = 0.08) but decreasing taxonomic extirpations (mean = -0.01, min = -0.15, max = 0.14, Fig. 1).

Trends in taxonomic temporal dissimilarity and colonisation decreased with latitude, while those of extirpations increased (Table 1, Fig. 1). However, trends in functional temporal β diversity, colonisations and extirpations showed no or weak relationships with latitude (Table 1, Fig. 1).

## Drivers of observed changes

Hierarchical GAMMs had good ($R^2$=0.34–0.44) and fair ($R^2$=0.19–0.24) fits for trends in taxonomic and functional indices, respectively (Table 2). Changes in taxonomic dissimilarity were mostly influenced by reach and catchment-scale descriptors (Table 2). However, trends in the five other indices were mostly influenced by global scale predictors (Table 2), followed by reach-scale predictors for taxonomic colonisations and catchment-scale predictors for all other indices (Table 2; Supporting information).

Changes in taxonomic dissimilarity were mostly influenced by changes in flow regimes (PC1: % of relative influence = 24.10), changes in stock unit density (PC2; 16.4%), upstream catchment area (9.2%) and changes in flow PC2 (10.8%; Fig. 2 $T_{Dis}$). Changes in taxonomic colonisations were mostly influenced by changes in flow (PC2; 22.2%), changes in water quality (15.1%), upstream catchment hardness (13.7%) and changes in precipitation seasonality (12.8%; Fig. 2 $T_{Col}$). Changes in taxonomic extirpations were mostly influenced by changes in land-cover (PC2; 18.5%), changes in stock unit density (PC2; 16.1%), segment slope (12.8%) and changes in precipitation seasonality (10.7%; Fig. 2 $T_{Ext}$).

By contrast, changes in functional dissimilarity (Fig. 2 $F_{Dis}$) were mostly influenced by changes in precipitation (22.4, 27.9 and 21.5%, respectively), followed by changes in flow (PC1; 19.6%), changes in land-cover (PC1; 17.4%) and changes in stock unit density (15.8%). Changes in functional colonisations (Fig. 2 $F_{Col}$) were mostly influenced by changes in changes in land-cover (PC1 and PC2; 17.7% and 12.9%) and changes in flow (PC1; 12.3%). Changes in functional extirpations (Fig. 2 $F_{Ext}$; Supporting information) were mostly influenced by changes in land-cover (PC1; 20.5%), changes in stock unit density (PC1; 15.9%) and changes in flow (PC1; 13.3%).
T<sub>Dis</sub>

F<sub>Dis</sub>

T<sub>Col</sub>

F<sub>Col</sub>

T<sub>Ext</sub>

F<sub>Ext</sub>

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

Figure 1. Maps of New Zealand illustrating rates of changes per decade (Trend decade⁻¹) in taxonomic and functional temporal β diversity indices of macroinvertebrate assemblages: dissimilarity (T<sub>Dis</sub> and F<sub>Dis</sub> respectively), colonisations (T<sub>col</sub> and F<sub>col</sub>) and extirpations (T<sub>Ext</sub> and F<sub>Ext</sub>). Temporal β diversity was measured by comparing the first year of sampling (1991) to each of the following years (1992–2016). Trends per decade are the slopes of linear regression models used to regress each temporal β diversity index against year. Point sizes are proportional to the absolute rate of change per decade. Red dots indicate decreasing trends, whereas blue dots indicate increasing trends. Density curves at the bottom right corner of each map illustrate the density distribution of each index. The x axis of each density curve corresponds to the limits of the legend.

Population trends and species range shifts

Population sizes increased on average over the time series (mean slope = 0.002 ± 0.009), and this trend was driven by 16% of the taxa which exhibited significant population trends. Significantly increasing populations (12%; mean slope = 0.011, maximum = 0.022, minimum = 0.003)

exceeded significantly decreasing populations (4%, mean slope = -0.015, min = -0.012, max = -0.018). Range sizes also increased on average over the 25-year time series (mean slope = 7536 ± 32 952 km² decade⁻¹), but only 14% of taxa exhibited significant trends. Range expansions (11%; mean slope = 27 270, min = 10 380, max = 63 330 km² decade⁻¹) exceeded range contractions (3%; mean slope = 14 640,
<table><caption>Table 1. Coefficients of determination (R²), intercept (standard error), F-value (F) and p-value (p) of latitude in each generalised linear mixed effect model relating temporal β diversity indices to latitude with island as random effect. Values are the estimated mean values of the random intercept and slope.</caption><thead><tr><th>Temporal β diversity index</th><th>R<sup>2</sup></th><th>Intercept (SE)</th><th>Slope (SE)</th><th>F</th><th>p</th></tr></thead><tbody><tr><td><strong>Taxonomic</strong></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Dissimilarity</td><td>0.06</td><td>0.112 (0.044)</td><td>0.111 (0.044)</td><td>3.72</td><td>0.06</td></tr><tr><td>Colonisations</td><td>0.41</td><td>-0.238 (0.064)</td><td>-0.238 (0.064)</td><td>13.8</td><td>&lt; 0.001</td></tr><tr><td>Extirpations</td><td>0.47</td><td>0.249 (0.034)</td><td>0.249 (0.034)</td><td>54.14</td><td>&lt; 0.001</td></tr><tr><td><strong>Functional</strong></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Dissimilarity</td><td>0.02</td><td>0.058 (0.046)</td><td>0.058 (0.046)</td><td>1.59</td><td>0.21</td></tr><tr><td>Colonisations</td><td>0.06</td><td>0.088 (0.045)</td><td>0.088 (0.045)</td><td>3.7</td><td>0.06</td></tr><tr><td>Extirpations</td><td>&lt; 0.01</td><td>0.028 (0.047)</td><td>0.028 (0.047)</td><td>0.36</td><td>0.55</td></tr></tbody></table>

min = 11 570, max = 17 700). Species shifted south on average by 56 km per decade (mean = 56.31, min = 1.05, max = 242.98 km decade⁻¹), with 22% of taxa having significant poleward range shifts (mean slope = 108.14, min = 27.37, max = 242.98 km decade⁻¹).

## Relationship with taxonomy and functional traits

Dissimilarity in population trends was significantly, albeit weakly, related to taxonomic distance (Mantel test: r = 0.17; p = 0.03, Supporting information; Fig. 3d). However, there was no relationship between taxonomic distance and dissimilarity in changes in range size, nor with dissimilarity in latitudinal range shifts (Mantel's r = 0.06; p = 0.27 and Mantel's r = 0.01; p = 0.40, respectively; Fig. 3e–f).

Population trends were positively correlated with the first axis of the functional trait space (r = 0.28, p < 0.01; Supporting information). This axis was positively correlated with species with several reproductive cycles per year (two or more) and long-life duration of adults (> 30 days; Supporting information). Changes in range size were not significantly related to functional trait space axes (p > 0.05). However, latitudinal range shifts were negatively correlated (r = 0.30, p < 0.01) with the third axis of the functional space, albeit with a large amount of scatter (Supporting information; this included large size taxa (maximum potential size > 40 mm), with a high number of descendants per reproductive cycle (> 1000), terrestrial oviposition sites and aerial respiration of aquatic stages; Supporting information).

Changes in functional composition were evident on the ordinations (Fig. 4a and b): for changes in population size (Fig. 4a), most winners and losers occupied exclusive areas of

trait space, signifying a temporal turnover of species populations and their traits over the 26-year time-series. However, taxa with the greatest latitudinal range shifts were mostly positioned in the centre of the functional space (Fig. 4b) indicating that these were rather generalists than specialist taxa (Mouillot et al. 2013).

## Discussion

We observed climate and land-use driven changes in temporal β diversity of taxa and functional traits across New Zealand rivers over a 25-year period relative to the first year of sampling. Trends in taxonomic colonisations increased with latitude while those of extirpations showed opposite patterns. Functional colonisations and extirpations showed weak or non-existent relationships with latitude (accepting our expectation E₁a but rejecting E₁b). Discrepancies among spatial patterns of taxonomic and functional diversity in river macroinvertebrate communities have also been reported elsewhere (Crabot et al. 2020, Sarremejane et al. 2020). This may suggest that communities are characterized by different degrees of functional redundancy across the landscape, allowing maintenance of functional diversity despite species losses, which may have important implications for ecosystem functions and services (McLean et al. 2019, Crabot et al. 2020).

Hierarchical models indicated a role of climate and land-use change in driving taxonomic and functional temporal β diversity. Increasing mean temperature, temperature variability and precipitation seasonality, tended to increase temporal turnover (sensu Anderson et al. 2011) in taxonomic composition. Many of the species that exhibited negative

<table><caption>Table 2. Coefficients of determination of each generalised additive mixed effect model for trends in each taxonomic and functional temporal β diversity index, spatial scale and the total of the four spatial scales.</caption><thead><tr><th>Temporal β diversity index</th><th>Global</th><th>Regional</th><th>Catchment</th><th>Reach</th><th>Total</th></tr></thead><tbody><tr><td><strong>Taxonomic</strong></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Dissimilarity</td><td>0.02</td><td>0.00</td><td>0.18</td><td>0.31</td><td>0.44</td></tr><tr><td>Colonisations</td><td>0.22</td><td>0.05</td><td>0.00</td><td>0.12</td><td>0.34</td></tr><tr><td>Extirpations</td><td>0.21</td><td>0.00</td><td>0.21</td><td>0.08</td><td>0.42</td></tr><tr><td><strong>Functional</strong></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Dissimilarity</td><td>0.10</td><td>0.00</td><td>0.09</td><td>0.04</td><td>0.21</td></tr><tr><td>Colonisations</td><td>0.14</td><td>0.00</td><td>0.10</td><td>0.02</td><td>0.24</td></tr><tr><td>Extirpations</td><td>0.09</td><td>0.00</td><td>0.09</td><td>0.02</td><td>0.19</td></tr></tbody></table>
<figure><img src="image_11.png" /><figcaption>Figure 2. Relative influence of the environmental variables selected in each generalised additive mixed-effects model to explain rates of changes in taxonomic and functional temporal β diversity of macroinvertebrate assemblages. Bars are coloured by their respective spatial scale (as is indicated by the legend). T<sub>Dis</sub>: taxonomic temporal β diversity; F<sub>Dis</sub>: functional temporal β diversity; T<sub>Col</sub>: taxonomic colonisations; F<sub>Col</sub>: functional colonisations; T<sub>Ext</sub>: taxonomic extirpations; F<sub>Ext</sub>: functional extirpations. Acronyms: Prec: changes in precipitation; Prec CV: changes in precipitation coefficient of variation; TMean: changes in mean air temperature; TSeas: changes in temperature seasonality; USPhosphorus: phosphorus concentration of upstream surface rocks; USCalcium: calcium concentration of upstream surface rocks; USHardness: mean hardness of upstream surface rocks; Altitude: altitude of the sampled site; USAvgSlope: upstream average slope; USCatchArea: upstream catchment area; LC PC1: changes in land-cover (PC 1); LC PC2: changes in land-cover (PC2); SUD 1: changes in stock unit densities of sheep, dear, dairy and beef (PC1); SUD 2: changes in stock unit densities of sheep, dear, dairy and beef (PC2); Flow PC1: changes in flow regimes (PC1); Flow PC2: changes in flow regimes (PC2); SegSlope: average river segment slope; SegRipShade: estimated river segment riparian shading; WQ: changes in water-quality; WidthMALF: river width at mean annual low flow; SSI: changes in substrate size.</figcaption></figure>

trends were those that specialise in cold-water pristine con-
ditions, including the stonefly *Stenoperla prasina* and the
mayfly *Amelotopsis perscitus*. One of the key mechanisms of
climate change-related temporal turnover in stream com-
munities globally has been a replacement of cold-dwellers
with warm-dwellers, including in streams of northwestern
Europe (Haase et al. 2019) and New South Wales, Australia
(Chessman 2009, Haase et al. 2019). By contrast, trends
in functional temporal β diversity responded to changes in
mean precipitation only (Supporting information). Changes
in patterns of precipitation has been shown as a key structur-
ing mechanism for freshwater temporal β diversity in many
locations (Tonkin et al. 2017, Cañedo-Argüelles et al. 2020).

Moreover, land-use change (increasing intensively man-
aged land and/or stock unit density of dairy and beef) was
a predominant driver of temporal β diversity in our analyses.
Increasing human land-use tended to increase taxa colonisa-
tions in these river systems (Supporting information). Human
land-use intensification has been found to hasten biodiversity
change in streams worldwide (Allan et al. 1997, Petsch et al.
2021), and more particularly in New Zealand, which has
experienced one of the highest rates of agricultural land inten-
sification over recent decades (OECD/FAO 2015). Several
studies have demonstrated that land use changes have pro-
foundly impacted New Zealand's stream communities for sev-
eral taxonomic groups (Clapcott et al. 2012, Foote et al. 2015,
<figure><img src="image_12.png" /><figcaption>Figure 3. Population trends, species range shifts and relationships with taxonomy. Changes in population size (a), range size (b) and latitudinal range shifts (c) for each taxon. The name of the five most and/or least changing taxa are given on each plot. Relationships between taxonomy and changes in population size (d), range size (e) and latitudinal range shifts (f). The numbers around the taxonomic trees delineate major taxonomic orders: 1: Other, 2: Diptera, 3: Plecoptera, 4: Coleoptera, 5: Hemiptera, 6: Trichoptera, 7: Megaloptera, 8: Ephemeroptera.</figcaption></figure>

Joy et al. 2019). Such effects of human land-use could render freshwater ecosystems more susceptible to climate change (Verberk et al. 2016, Tickner et al. 2020), and our results indicate that their effects are currently operating in unison.

At the reach scale, changes in flow regimes were found to be the most important drivers of changes in taxonomic and functional β diversity. Riverine ecosystems are governed by patterns of temporal variation in flow regimes (Tonkin et al. 2018a). As climate and land-use change modifies the natural flow regime in many river systems, components of flow regimes are expected to shift, even under the most conservative climate change scenarios (Rood et al. 2008, Ficklin et al. 2018). The results of our hierarchical models indicate how influential flow driven habitat changes are to the organization of river macroinvertebrate communities. Nevertheless, while flow change was important for almost all facets of beta diversity, it had only a small role as a driver of taxonomic extinctions, where land-use change was the predominant driver. This may reflect the relative flexibility of most New Zealand invertebrates to unpredictable flow regimes, given its oceanic climate (Winterbourn et al. 1981, Tonkin et al. 2018a).

Contrary to our second expectation (E₂) our results identified greater amounts of increases in population and range size (i.e. winners) than decreases (i.e. losers) across the entire species pool. This suggests that the overall taxonomic homogenisation of these river macroinvertebrate communities recently observed (Mouton et al. 2020) is characterised by a greater amount of increase in population size and range size of winners than widespread declines of losers. Nonetheless, we found increasing rates of taxonomic extirpations at sites located at the north-eastern boundary of the North Island of New Zealand. Conversely, the greatest rates of taxonomic colonisations were located at the southern boundaries of each island. More importantly, almost a quarter of the taxa examined here also tracked the shifting isotherms by shifting their ranges towards the south pole. Given the geographic isolation of New Zealand and the observed latitudinal patterns of species distribution shifts here, the biodiversity of New Zealand's rivers could experience a 'cul-de-sac' effect, if it is to be exposed to more intense climate change (Sauer et al. 2011, Albouy et al. 2012). This trend may further be exacerbated for species dispersing exclusively along river networks
<figure><img src="image_13.png" /><figcaption>Figure 4. Differences in functional space occupancy between the species that experienced (a) population trends, and (b) latitudinal range shifts. Points size is proportional to the degree of change presented on each plot. Colour gradients also represent the degree of change presented on each plot. The black line delimits the convex hull (light grey polygon) occupied by the species pool within each bi-dimensional trait space.</figcaption></figure>

(Hylander and Ehrlén 2013, Bush and Hoskins 2017, Tonkin et al. 2018b).

We found that population trends were taxonomically structured, however no relationship was found with species range shifts (thus, partly accepting E₃₃). We also found that changes in population size and latitudinal range shifts of taxa were associated with distinct functional traits, confirming our expectation of a relationship between the functional characteristics of species and their vulnerability to climate and land-use change (E₃₃). Several mayflies (Ephemeroptera) and caddisflies (Trichoptera) exhibited decreasing abundances over the 26-year period, many species from these orders are widely considered as pollution-sensitive (Stark 1985, Usseglio-Polatera and Bournaud 1989). By contrast, we observed a rise in crustaceans (Crustacea), snails (Gastropoda) and scavenger beetles (Coleoptera), these groups being mainly composed of eurythermal taxa, typical from slow-flowing waters, and tolerant to a wide range of water quality conditions, including low oxygen but high nutrient concentrations (Stark 1995). Increasing population sizes and latitudinal range shifts were strongly related to high numbers of reproductive cycles per year (two or more) and descendants per reproductive cycle, respectively. Such r-selected strategies have been found to correlate with other increasing population sizes and/or species range shifts in freshwater and marine organisms' facing climate change. For example, McLean et al. (2018b) found that rapid warming drove marine pelagic fishes with r-selected life history traits to shift abruptly poleward. Similarly, Comte et al. (2014) found that species with high propagule pressure (i.e. r-strategists) and greater mobility, displayed the greatest range shifts in stream fishes facing climate change. Here, increasing population sizes and latitudinal range shifts were

also related to traits like long-life duration of adults and large body sizes, which are more typical of K-strategists' species. This result is in agreement with previous studies highlighting that such strategies can be promoted under climate change (del Cacho et al. 2012), especially in running waters experiencing climate-driven flow reduction (Floury et al. 2017).

## Conclusion

We found that rates of changes in macroinvertebrate taxonomic composition in New Zealand's River systems over a 25-year period are latitudinally structured but not changes in their trait composition. We observed increasing rates of taxonomic extirpations at the rear edge of mainland New Zealand but increasing taxonomic colonisations at the leading edges. Further, we found that most taxa tended to increase in population and range size rather than decrease. Changes in community composition correlated with recent climate and land-use change. Macroinvertebrate functional traits related to life-history and morphology explained changes in population size and species latitudinal range shifts. Our findings highlight the critical need to move beyond observation of emergent state-level variables to understanding the mechanisms underpinning taxonomic and functional reorganisation of biodiversity under ongoing environmental changes (Tonkin et al. 2019).

*Acknowledgements* – We thank Kevin Collier, Sylvain Dolédec, David Mouillot, the associate editor and four anonymous reviewers for their helpful comments on earlier versions of the manuscript.

*Funding* – JDT is supported by a Rutherford Discovery Fellowship administered by the Royal Society Te Apārangi (RDF-18-UOC-007).
**Transparent Peer Review**

The peer review history for this article is available at <https://publons.com/publon/10.1111/ecog.06148>.

**Data availability statement**

The data that support the findings of this study are available from the National Inst. of Water and Atmospheric research (NIWA) upon request. Restrictions apply to the availability of these data, which were used under license for this study.

**Supporting information**

The supporting information associated with this article is available from the online version.

**References**

Albouy, C. et al. 2012. Combining projected changes in species richness and composition reveals climate change impacts on coastal Mediterranean fish assemblages. – Global Change Biol. 18: 2995–3003.

Allan, D. et al. 1997. The influence of catchment land use on stream integrity across multiple spatial scales. – Freshwater Biol. 37: 149–161.

Allan, J. D. 2004. Landscapes and riverscapes: the influence of land use on stream ecosystems. – Annu. Rev. Ecol. Evol. Syst. 35: 257–284.

Anderson, M. J. et al. 2011. Navigating the multiple meanings of β diversity: a roadmap for the practicing ecologist. – Ecol. Lett. 14: 19–28.

Antão, L. H. et al. 2020. Temperature-related biodiversity change across temperate marine and terrestrial systems. – Nat. Ecol. Evol. 4: 927–933.

Blowes, S. A. et al. 2019. The geography of biodiversity change in marine and terrestrial assemblages. – Science 366: 339–345.

Booker, D. J. 2015. Hydrological indices for national environmental reporting. – Prepared for Ministry for the Environment, NIWA, p. 39.

Box, G. E. and Cox, D. R. 1964. An analysis of transformations. – J. R. Stat. Soc. B 26: 211–243.

Brondizio, E. S. et al. 2019. Global assessment report on biodiversity and ecosystem services of the Intergovernmental Science-Policy Platform on Biodiversity and Ecosystem Services. – IPBES Secretariat.

Brown, C. J. et al. 2013. Managing for interactions between local and global stressors of ecosystems. – PLoS One 8: e65765.

Bush, A. and Hoskins, A. J. 2017. Does dispersal capacity matter for freshwater biodiversity under climate change? – Freshwater Biol. 62: 382–396.

Cañedo-Argüelles, M. et al. 2020. As time goes by: 20 years of changes in the aquatic macroinvertebrate metacommunity of Mediterranean river networks. – J. Biogeogr. 47: 1861–1874.

Chessman, B. C. 2009. Climatic changes and 13-year trends in stream macroinvertebrate assemblages in New South Wales, Australia. – Global Change Biol. 15: 2791–2802.

Chevenet, F. et al. 1994. A fuzzy coding approach for the analysis of long-term ecological data. – Freshwater Biol. 31: 295–309.

Chin, A. et al. 2010. An integrated risk assessment for climate change: analysing the vulnerability of sharks and rays on Australia's Great Barrier Reef. – Global Change Biol. 16: 1936–1953.

Clapcott, J. E. et al. 2012. Quantifying relationships between land-use gradients and structural and functional indicators of stream ecological integrity. – Freshwater Biol. 57: 74–90.

Comte, L. et al. 2014. Species traits and phylogenetic conservatism of climate-induced range shifts in stream fishes. – Nat. Commun. 5: 5023.

Conti, L. et al. 2014. A trait-based approach to assess the vulnerability of European aquatic insects to climate change. – Hydrobiologia 721: 297–315.

Crabot, J. et al. 2020. Drying determines the temporal dynamics of stream invertebrate structural and functional beta diversity. – Ecography 43: 620–635.

Daskalova, G. N. et al. 2020. Rare and common vertebrates span a wide spectrum of population trends. – Nat. Commun. 11: 4394.

Davies-Colley, R. J. et al. 2011. Twenty years of New Zealand's National Rivers Water Quality Network: benefits of careful design and consistent operation 1. – J. Am. Water Resour. Assoc. 47: 750–771.

Dawson, T. P. et al. 2011. Beyond predictions: biodiversity conservation in a changing climate. – Science 332: 53–58.

del Cacho, M. et al. 2012. Effect of experimentally induced climate change on the seed bank of a Mediterranean shrubland. – J. Veg. Sci. 23: 280–291.

Dolédec, S. et al. 2006. Comparison of structural and functional approaches to determining landuse effects on grassland stream invertebrate communities. – J. North Am. Benthol. Soc. 25: 44–60.

Doledec, S. et al. 2011. Invertebrate community responses to land use at a broad spatial scale: trait and taxonomic measures compared in New Zealand rivers. – Freshwater Biol. 56: 1670–1688.

Dormann, C. et al. 2007. Methods to account for spatial autocorrelation in the analysis of species distributional data: a review. – Ecography 30: 609–628.

Dornelas, M. et al. 2014. Assemblage time series reveal biodiversity change but not systematic loss. – Science 344: 296–299.

Dornelas, M. et al. 2019. A balance of winners and losers in the anthropocene. – Ecol. Lett. 22: 847–854.

Dray, S. et al. 2020. R package 'ade4', analysis of ecological data: exploratory and euclidean methods in environmental sciences. – <https://CRAN.R-project.org/package=ade4>.

Drost, H.-G. et al. 2018. myTAI: evolutionary transcriptomics with R. – Bioinformatics 34: 1589–1590.

Fick, S. E. and Hijmans, R. J. 2017. WorldClim 2: new 1-km spatial resolution climate surfaces for global land areas. – Int. J. Climatol. 37: 4302–4315.

Ficklin, D. L. et al. 2018. Natural and managed watersheds show similar responses to recent climate change. – Proc. Natl Acad. Sci. USA 115: 8553–8557.

Floury, M. et al. 2017. Assessing long-term effects of multiple, potentially confounded drivers in ecosystems from species traits. – Global Change Biol. 23: 2297–2307.

Floury, M. et al. 2018. Climatic and trophic processes drive long-term changes in functional diversity of freshwater invertebrate communities. – Ecography 41: 209–218.

Foote, K. J. et al. 2015. New Zealand dairy farming: milking our environment for all its worth. – Environ. Manage. 56: 709–720.

Gotelli, N. J. et al. 2017. Community-level regulation of temporal trends in biodiversity. – Sci. Adv. 3: e1700315.
Haase, P. et al. 2019. Moderate warming over the past 25 years has already reorganized stream invertebrate communities. – Sci. Total Environ. 658: 1531–1538.

Heino, J. et al. 2009. Climate change and freshwater biodiversity: detected patterns, future trends and adaptations in northern regions. – Biol. Rev. 84: 39–54.

Hylander, K. and Ehrlén, J. 2013. The mechanisms causing extinction debts. – Trends Ecol. Evol. 28: 341–346.

IPCC 2014. Climate change 2014: synthesis report. – In: IPCC Core Writing Team et al. (eds), Contribution of working groups I, II and III to the fifth assessment report of the intergovernmental panel on climate change. IPCC, p. 151.

Johnson, C. N. et al. 2017. Biodiversity losses and conservation responses in the Anthropocene. – Science 356: 270–275.

Jowett, I. G. et al. 1991. Microhabitat preferences of benthic invertebrates and the development of generalised *Deleatidium* spp. habitat suitability curves, applied to four New Zealand rivers. – N. Z. J. Mar. Freshwater Res. 25: 187–199.

Joy, M. K. et al. 2019. Decline in New Zealand's freshwater fish fauna: effect of land use. – Mar. Freshwater Res. 70: 114–124.

Julian, J. P. et al. 2017. River water quality changes in New Zealand over 26 years: response to land use intensity. – Hydrol. Earth Syst. Sci. 21: 1149.

Kuczynski, L. et al. 2018. Concomitant impacts of climate change, fragmentation and non-native species have led to reorganization of fish communities since the 1980s. – Global Ecol. Biogeogr. 27: 213–222.

Kuhn, M. 2008. Variable importance using the caret package. – J. Stat. Softw. 28: 1–26.

Kuhn, M. et al. 2020. Package 'caret'. – <https://cran.r-project.org/web/packages/caret/index.html>

Laliberté, E. et al. 2014. Package 'FD': measuring functional diversity from multiple traits. – <https://CRAN.R-project.org/package=FD>.

Landcare Research 2015. Land Cover Database ver. 4.1, Mainland New Zealand, online.

Lavorel, S. et al. 2008. Assessing functional diversity in the field—methodology matters! – Funct. Ecol. 22: 134–147.

Leathwick, J. et al. 2010. Freshwater ecosystems of New Zealand (FENZ) geodatabase. – Dept of Conservation, Wellington, New Zealand, p. 51.

Legendre, P. 2019. A temporal beta-diversity index to identify sites that have changed in exceptional ways in space–time surveys. – Ecol. Evol. 9: 3500–3514.

Legendre, P. and Gauthier, O. 2014. Statistical methods for temporal and space–time analysis of community composition data. – Proc. R. Soc. B 281: 20132728.

Lenoir, J. et al. 2019. Species better track the shifting isotherms in the oceans than on lands. – Nat. Ecol. Evol. 4: 1044–1059.

Li, D. et al. 2020. Changes in taxonomic and phylogenetic diversity in the Anthropocene. – Proc. R. Soc. B 287: 20200777.

Loiseau, N. et al. 2020. Global distribution and conservation status of ecologically rare mammal and bird species. – Nat. Commun. 11: 5071.

Maclean, I. M. and Wilson, R. J. 2011. Recent ecological responses to climate change support predictions of high extinction risk. – Proc. Natl Acad. Sci. USA 108: 12337–12342.

MacLean, S. A. and Beissinger, S. R. 2017. Species' traits as predictors of range shifts under contemporary climate change: a review and meta-analysis. – Global Change Biol. 23: 4094–4105.

Magurran, A. E. and Henderson, P. A. 2010. Temporal turnover and the maintenance of diversity in ecological assemblages. – Phil. Trans. R. Soc. B 365: 3611–3620.

Magurran, A. E. et al. 2019. Temporal β diversity – a macroecological perspective. – Global Ecol. Biogeogr. 28: 1949–1960.

Mantel, N. 1967. The detection of disease clustering and a generalized regression approach. – Cancer Res. 27: 209–220.

McLean, M. et al. 2018a. A climate-driven functional inversion of connected marine ecosystems. – Curr. Biol. 28: 3654–3660.

McLean, M. et al. 2018b. Ecological and life history traits explain a climate-induced shift in a temperate marine fish community. – Mar. Ecol. Prog. Ser. 606: 175–186.

McLean, M. et al. 2019. Trait structure and redundancy determine sensitivity to disturbance in marine fish communities. – Global Change Biol. 25: 3424–3437.

Mouillot, D. et al. 2013. A functional approach reveals community responses to disturbances. – Trends Ecol. Evol. 28: 167–177.

Mouton, T. L. et al. 2020. Increasing climate-driven taxonomic homogenization but functional differentiation among river macroinvertebrate assemblages. – Global Change Biol. 26: 6904–6915.

Odum, E. P. 1950. Bird populations of the highlands (North Carolina) Plateau in relation to plant succession and avian invasion. – Ecology 31: 587–605.

OECD/FAO 2015. OECD-FAO Agricultural Outlook 2015–2024. – OECD/Food and Agriculture Organization of the United Nations, Paris, p. 148.

Oksanen, J. et al. 2014. Package 'vegan'. Community ecology package, Rpackagever.2. – <https://CRAN.R-project.org/package=vegan>.

Olden, J. D. et al. 2018. The Homogocene: a research prospectus for the study of biotic homogenisation. – NeoBiota 37: 23.

Outhwaite, C. L. et al. 2020. Complex long-term biodiversity change among invertebrates, bryophytes and lichens. – Nat. Ecol. Evol. 4: 384–392.

Pacifici, M. et al. 2017. Species' traits influenced their response to recent climate change. – Nat. Clim. Change 7: 205–208.

Paradis, E. et al. 2019. Package 'ape'. – Analyses of phylogenetics and evolution, ver. 2. – <https://CRAN.R-project.org/package=ape>.

Petsch, D. K. et al. 2021. A synthesis of land use impacts on stream biodiversity across metrics and scales. – Ecology 102: e03498.

Pinheiro, J. et al. 2017. Package 'nlme'. – Linear and nonlinear mixed effects models, ver. 3. – <https://CRAN.R-project.org/package=nlme>.

Poff, N. L. 1997. Landscape filters and species traits: towards mechanistic understanding and prediction in stream ecology. – J. North Am. Benthol. Soc. 16: 391–409.

Pyne, M. I. and Poff, N. L. 2017. Vulnerability of stream community composition and function to projected thermal warming and hydrologic change across ecoregions in the western United States. – Global Change Biol. 23: 77–93.

Quinn, J. M. and Hickey, C. W. 1990. Characterisation and classification of benthic invertebrate communities in 88 New Zealand rivers in relation to environmental factors. – N. Z. J. Mar. Freshwater Res. 24: 387–409.

Rood, S. B. et al. 2008. Declining summer flows of Rocky Mountain rivers: changing seasonal hydrology and probable impacts on floodplain forests. – J. Hydrol. 349: 397–410.

Ruhí, A. et al. 2017. Interpreting beta-diversity components over time to conserve metacommunities in highly dynamic ecosystems. – Conserv. Biol. 31: 1459–1468.

Sarremejane, R. et al. 2020. Stochastic processes and ecological connectivity drive stream invertebrate community responses to short-term drought. – J. Anim. Ecol. 90: 886–898.

Sauer, J. et al. 2011. Low mountain ranges: summit traps for montane freshwater species under climate change. – Biodivers. Conserv. 20: 3133–3146.
Scarsbrook, M. R. et al. 2000. New Zealand's National River Water Quality Network: long-term trends in macroinvertebrate communities. – N Z. J. Mar. Freshwater Res. 34: 289–302.

Shimadzu, H. et al. 2015. Measuring temporal turnover in eco-
logical communities. – Methods Ecol. Evol. 6: 1384–1394.

Smith, D. and McBride, G. 1990. New Zealand's National Network Water Quality Monitoring Network design and first year's operation. – J. Am. Water Resour. Assoc. 26: 767–775.

Stark, J. D. 1985. A macroinvertebrate community index of water quality for stony streams. – Ministry of Works and Development, Wellington, p. 52.

Tickner, D. et al. 2020. Bending the curve of global freshwater biodiversity loss: an emergency recovery plan. – BioScience 70: 332–342.

Tilman, D. et al. 1994. Habitat destruction and the extinction debt. – Nature 371: 65.

Tonkin, J. D. et al. 2017. Seasonality and predictability shape tem-
poral species diversity. – Ecology 98: 1201–1216.

Tonkin, J. D. et al. 2018a. Flow regime alteration degrades eco-
logical networks in riparian ecosystems. – Nat. Ecol. Evol. 2:
86.

Tonkin, J. D. et al. 2018b. The role of dispersal in river network metacommunities: patterns, processes and pathways. – Freshwater Biol. 63: 141–163.

Tonkin, J. D. et al. 2019. Prepare river ecosystems for an uncertain future. – Nature 570: 301–303.

Urban, M. C. 2015. Accelerating extinction risk from climate change. – Science 348: 571–573.

Usseglio-Polatera, P. and Bournaud, M. 1989. Trichoptera and ephemeroptera as indicators of environmental changes of the Rhone river at lyons over the last twenty-five years. – Regul. Rivers Res. Manage. 4: 249–262.

Van Looy, K. et al. 2017. At what scale and extent environmental gradients and climatic changes influence stream invertebrate communities? – Sci. Total Environ. 580: 34–42.

Vaughan, I. P. and Gotelli, N. J. 2019. Water quality improvements offset the climatic debt for stream macroinvertebrates over twenty years. – Nat. Commun. 10: 1956.

Verberk, W. C. et al. 2016. Field and laboratory studies reveal interacting effects of stream oxygenation and warming on aquatic ectotherms. – Global Change Biol. 22: 1769–1778.

Vié, J.-C. et al. 2009. Wildlife in a changing world: an analysis of the 2008 IUCN Red List of threatened species. – IUCN.

Winterbourn, M. J. et al. 1981. Are New Zealand stream ecosystems really different? – N. Z. J. Mar. Freshwater Res. 15: 321–328.

Wood, S. and Wood, M. S. 2015. Package ‘mgcv’. – R package ver. 1-29. – <https://CRAN.R-project.org/package=mgcv>.