Pollution control can help mitigate future climate change impact on European grayling in the UK

We compare the performance of habitat suitability models using climate data only or climate data together with water chemistry, land cover and predation pressure data to model the distribution of European grayling (Thymallus thymallus). From these models, we (a) investigate the relationship between habitat suitability and genetic diversity; (b) project the distribution of grayling under future climate change; and (c) model the effects of habitat mitigation on future distributions.

tions of future flow scenarios obtained from the Centre for Ecology and Hydrology.
To examine potential mitigation effects within habitats, models were run with manipulation of orthophosphate, nitrite and copper concentrations.

Results:
We mapped suitable habitat for grayling in the present and the future. The Full model achieved substantially higher discriminative power than the Simple model.
For low suitability habitat, higher levels of inbreeding were observed for adaptive, but not neutral, loci. Future projections predict a significant contraction of highly suitable areas. Under habitat mitigation, modelling suggests that recovery of suitable habitat of up to 10% is possible.

Main conclusions:
Extending the climate-only model improves estimates of habitat suitability. Significantly higher inbreeding coefficients were found at immune genes, but not neutral markers in low suitability habitat, indicating a possible impact of environmental stress on evolutionary potential. The potential for habitat mitigation to alleviate distributional changes under future climate change is demonstrated, and specific recommendations are made for habitat recovery on a regional basis.

| INTRODUC TI ON
Human ecosystem alteration, such as habitat loss and degradation, invasive species and overharvesting, can cause environmental stress (Brown, Saunders, Possingham, & Richardson, 2013;Crain, Kroeker, & Halpern, 2008), which can directly impair the adaptive potential of populations and increase vulnerability to extinction (Charmantier & Garant, 2005;Frankham, 2005;Hoffmann & Hercus, 2000). These stressors are thought to be impacting biodiversity in a way comparable to historic mass extinction events (Ceballos et al., 2015) with unprecedented declines in global biodiversity now occurring.
These impacts on biodiversity are now exacerbated by climate change through cumulative or synergistic effects (Brook, Sodhi, & Bradshaw, 2008). Climate change is expected to become a major threat to biodiversity over the next decades (Bellard, Bertelsmeier, Leadley, Thuiller, & Courchamp, 2012;Thomas et al., 2004), and the integration of climate change predictions into current conservation and biodiversity planning is therefore essential (Araújo, Cabeza, Thuiller, Hannah, & Williams, 2004;Heller & Zavaleta, 2009). Such integration requires knowledge of the sensitivities of species to various climatic parameters in order to assess vulnerability to climate change (Hulme, 2005). A synergistic approach incorporating habitat change and degradation with exposure and sensitivity to climate change and adaptive capacity will further improve vulnerability analysis (Dawson, Jackson, House, Prentice, & Mace, 2011;Williams, Shoo, Isaac, Hoffmann, & Langham, 2008). The investigation of differences in habitat quality may serve as a predictive framework to assess evolutionary potential across the range of a species because there is a close relationship between environmental conditions and genetic diversity (Charmantier & Garant, 2005;Frankham, 2005;Hoffmann & Hercus, 2000).
Species distribution modelling (SDM) can be used to identify key environmental parameters that affect the distribution of a species, by combining species occurrence with environmental data (Elith & Leathwick, 2009). These models are powerful tools to assess specific sensitivity to environmental change and to predict the influence of climate on species distribution (Thomas et al., 2004). The implementation of SDMs has been successful in conservation and reserve planning (Rodríguez, Brotons, Bustamante, & Seoane, 2007), invasive species management (Jiménez-Valverde et al., 2011), epidemiology research (Puschendorf et al., 2009) and predicting potential effects of climate change on biodiversity (Hijmans & Graham, 2006).
The classical approach of SDMs is to model "bioclimatic envelopes" (Pearson & Dawson, 2003) because climate is usually the dominant factor in determining species ranges (Araújo & Peterson, 2012) and climatic variables are therefore considered sufficient to describe changes in distribution (Bucklin et al., 2015). As a first approximation, this approach can describe or predict niche requirements and indicate where tolerance limits are exceeded under conditions of climate change (Pearson & Dawson, 2003). However, improvements in model performance can be made by the addition of non-climate parameters, particularly at small spatial scales, where other local factors can become dominant (Pearson & Dawson, 2003;Pearson, Dawson, & Liu, 2004;Stanton, Pearson, Horning, Ersts, & Reşit Akçakaya, 2012). This approach allows development of climate change mitigation strategies through the identification of non-climate-related drivers of biodiversity change (Pereira et al., 2010). Sutton and Soto (2012) show that predictive performance is highest when there are interactions between climatic and other variables. These authors highlight the importance of including variables such as land cover to study their effect in combination with climate (Stanton et al., 2012).
European grayling (Salmonidae, Thymallus thymallus) is a protected species (listed in appendix III of the Bern convention) of conservation and economic importance (Swatdipong, Primmer, & Vasemägi, 2010), that, even in comparison with other salmonids, shows high sensitivity to high temperature (Ibbotson et al., 2001;Jonsson & Jonsson, 2009) and exhibits narrow water quality requirements (Oberdorff, Pont, Hugueny, & Porcher, 2002;Uiblein, Jagsch, Honsig-Erlenburg, & Weiss, 2001). As such, grayling make an ideal indicator species of habitat quality and climate change. We aim to investigate the influence of climate and other habitat quality parameters on the distribution of this species. Specifically, we (a) compare a Simple climate-only model with one that incorporates the sensitivity of grayling to habitat quality parameters (i.e. current velocity, land use, water chemistry and predator density); (b) test associations of habitat quality and neutral and adaptive genetic diversity; (c) forecast grayling distributions under future climate change scenarios; and (d) investigate the effects of habitat improvement on future distributions. The results are discussed in the context of conservation management.  Table 1). The eastern part of the longitudinal range of European grayling (adjacent to the Ural Mountains) was omitted to avoid bias due to low numbers of records available in that area. Subsequently, we focussed on a subset study area (Figure 1 (Orr et al., 2010). These data sets were used in the Full model with the aim to assess the sensitivity of grayling to non-climatic parameters that impact on habitat suitability. The same subset study area was used in the Simple model to allow for comparison of model performance when using only climate data.

| Development of the Full model
Unique occurrence points (2,846) recorded since 1990 for European grayling were obtained from the GBIF database in 2014 (www.gbif. org) (Figure 1). This time frame was selected as there is a lack of recently updated records in GBIF, so that only using records made after 2000 would significantly underestimate the current distribution. To account for uneven sampling effort and to reduce spatial clustering, subsampling of occurrence records was done using the SDM toolbox (Brown, 2014)  The scale is the same as the resolution of environmental layers and was chosen due to the generally high site fidelity of grayling, aiming for a maximal resolution of small-scale differences in habitat preference (Nykänen, Huusko, & Lahti, 2004). Only occurrence points that had records for all environmental parameters were used for model tuning and evaluation. There were 441 occurrence points for the Simple model and 283 for the Full model.  Table 1). Land cover data for the year 2000 were retrieved from the European Environmental Agency at a resolution of 100 m. Predation by the great cormorant (Phalacrocorax carbo, hereafter cormorant) on inland fish populations has increased significantly over the past thirty years (Callaghan, Kirby, Bell, & Spray, 1998), having severe local impact in some cases (Vetemaa et al., 2010). To investigate the impact of this predator, GBIF records of cormorant since 1990 were used to estimate spatial patterns of abundance using a Gaussian kernel density F I G U R E 1 Right: Study area used for the selection of bioclimatic variables; Left: Subset study area; unique occurrence records of Thymallus thymallus are shown in red; these were spatially rarefied at a scale of 30 km 2 for the entire European distribution (right) and at a scale of 1 km 2 for the subset study area (left); populations genotyped at microsatellite and MH loci are shown in black, river lines are shown in grey Simple model study area for parameter selection function, which smooths individual occurrence counts per grid cell (Silverman, 2018), which was 10 km 2 here, using the SDM toolbox (Brown, 2014) for ArcGIS 10.3 (ESRI). As this estimate of density may be biased and confounded by sampling intensity, we additionally generated a categorical presence/absence layer based on the observation of occurrence records within a 10-km 2 grid cell and compared model predictions. Measurement of biochemical oxygen demand (BOD), concentration of calcium, chloride, copper, nitrate, nitrite and orthophosphates, pH, suspended solids and total ammonium were obtained from long-term surveys conducted by the UK Environment Agency. As it is necessary to acquire data from a large number of occurrence and background points, a good spatial coverage of measurements was ensured (minimum of 6,000 sites) after filtering using the following approach. To remove outliers for each site, all records since the year 2000 were averaged for the variables above and sites were removed if either the mean or the standard deviation exceeded a modified z-score of 3.5 (Iglewicz & Hoaglin, 1993). Water temperature measurements were derived from the River Surface Water temperature database as described in Orr et al. (2010). For each region, all records since 1985 were extracted and were again filtered if measurements exceeded a modified z-score of 3.5. Of the remaining records, averages were calculated per site, separately for each season. The time period was chosen to maximize the records per site in order to resolve local differences despite strong fluctuations. River flow data were obtained from the Future Flows and Groundwater Levels data set (Prudhomme et al., 2013, specifically the national maps of changes in river flow statistics (Prudhomme, 2012).  Table 1 lists variables tested for relative impact on habitat suitability for grayling.

| Variable selection, model calibration and evaluation
We used model selection to avoid overfitting the data. Briefly, we used a stepwise variable selection process including species-specific tuning of regularization and AICc as an evaluation metric (see Warren, Wright, Seifert, & Shaffer, 2014;Zeng, Low, & Yeo, 2016 transferability of the model, we used a "masked geographically structured" evaluation approach, where the study area was divided into k geographic blocks and evaluation was done on each block in turn (Radosavljevic & Anderson, 2014

| Testing association of habitat suitability and population inbreeding
To examine associations between habitat suitability and genetic  (Hochberg & Benjamini, 1990). When a significant relationship between habitat suitability and genetic diversity was identified (see RESULTS ), we further tested the relative importance of effective population size (from Dawnay et al., 2011) and management background (native non-stocked, native stocked or introduced) on genetic diversity using a Random Forest analysis in the R package randomForest (Liaw & Wiener, 2002).

| Future distribution under differing climate change scenarios
Data for future climate scenarios were downloaded from the WorldClim database (Hijmans et al., 2005)

| Evaluating the effects of habitat improvement
To evaluate the effect of habitat improvement as a mitigation strategy under climate change conditions, future projections were also done under scenarios where either orthophosphate, nitrite or copper concentration was artificially manipulated independently or targeting all three parameters in combination. Here, concentrations were capped at the maximum values that had been shown not to impact on habitat suitability estimates for grayling (0.15 mg/L, 0.02 mg/L and 2 μg/L for orthophosphate, nitrite or copper concentration respectively; see RESULTS). These variables were selected, because they were shown to affect habitat suitability for grayling in the Full model (see RESULTS) and could be targeted through management within a habitat improvement context.

| Comparison of Full and Simple model
The Full model (including non-climatic parameters) showed a better relative performance than the Simple model (which was restricted to just the bioclimatic variables; delta AICc = 4,374). The habitat suitability predicted by the Simple and Full model, respectively, across the study area is shown in Figure 2. The classified predictions of the two models agree in 66% of the study area ( Figure 2). In 20% of the study area, the Simple model predicts higher suitability than the Full model, and in 14% of the area, the Simple model predicts lower suitability than the Full model ( Figure 2). Tuning of regularization resulted in an optimal random multiplier of 2.2, which repre-

| Full model
When using the Full model, 11% of the area was classified as unsuitable (below MTP), 61% as low/medium suitability (below ESS) and 28% as high suitability (above ESS; Figure 2). Climatic environmental parameters contributed 51% of the final set of variables used to build the model. These were precipitation of the wettest quarter (Bio 16), isothermality (Bio 3), described as the ratio of the mean diurnal range (Bio 2) to the mean annual range (Bio 7), maximum temperature in the warmest month (Bio 5) and mean water temperature in summer and the 90th (Q10) percentile and 5th (Q95) of annual flow, representing highest and lowest flow rates, respectively. Water chemistry parameters had a total contribution of 37% to the final model (Table 2). These parameters were the concentrations of calcium, dissolved copper, nitrite and orthophosphates. Further, land cover had a contribution of 12% and cormorant density had a contribution of less than 1% (Table 2).
This was also the case using a presence/absence layer of cormorant occurrence observations within 10 km 2 as an alternative modelling input, which showed minor differences in model predictions (average per cell standard deviation was 0.007).
As a main model output, the sensitivity of grayling to different environmental parameters was quantified. A probability of grayling presence above 50% was observed for ( Figure S1). Highest suitability in regard to land use was observed for fruit trees and berry plantations, inland marshes and broadleaf forest. Land use classifications with a probability of presence below 50% were observed for urban fabric and industrial or commercial areas, natural grasslands and moors and heathlands ( Figure S1).

| Association of habitat suitability and population inbreeding
No differences were observed between high and low suitability habitats for observed or expected heterozygosity and allelic richness.
Inbreeding coefficients (F IS ) were significantly higher in low suitability habitats, defined by the Full model than in high suitability habitats for MH class II markers (clustered Mann-Whitney-Wilcoxon test, p < .001) but not for neutral markers (Figure 3). Random Forest analysis also showed high relative importance of habitat suitability estimates compared to effective population size and management background in explaining variation in MH II F IS values (Table 3).
Habitat suitability was the parameter that mostly decreased the mean square error (MSE) ( Table 3). When using the mean decrease in Gini index as measure of variable importance, both habitat suitability and effective population size were ranked most important (Table 3).
No significant differences in genetic diversity were observed between areas of low and high suitability defined by the Simple model.

| Future projections
Under conditions of climate change, the projections for 2050 predict predominantly a significant loss of high suitability habitat ( Figure 4; Table 4). This was estimated as a reduction in highly suitable areas of 21% for the RCP 2.6 and the RCP 8.5 scenario of change (Table 4)

TA B L E 2 Relative contribution and
permutation importance for all variables included in the Full model

| Evaluating the effects of habitat improvement
Decreasing orthophosphate, nitrite or copper concentration in areas where they exceed the estimated tolerance of grayling showed recovery of suitable habitat is possible under current and future conditions (Table 5, Figures 5 and 6). An overall recovery of suitable habitat of up to 7%, 6% and 4% relative to the total study area is predicted to be possible by making locally optimal improvements in either orthophosphate, nitrite or copper concentration under current conditions and predictions for 2050 for the RCP 2.6 and 8.5 scenarios, respectively (Table 5). If the target parameters are optimized in combination, the total improvements of suitability classifications are further increased and estimated 14%, 10% and 6% in total for current and future scenarios (RCP 2.6 and 8.5), respectively (Table 5). Possible improvements would greatly affect areas in which grayling is currently present (Figures 5 and 6). The highest net total gain in estimated suitability that can be achieved reveals locally optimal strategies for improvement, targeting either orthophosphate, nitrite or copper concentration ( Figure 6).
Highest local increases in habitat suitability of 26% were observed for reductions in orthophosphate concentrations under current conditions ( Figure 6). For other locations, a reduction in nitrite or metal pollution had higher effects on suitability achieving relative maximum improvements of 20% and 17% under current conditions ( Figure 6 ).

| D ISCUSS I ON
Here, we modelled habitat suitability for grayling across the UK and   (Alabaster & Lloyd, 2013;Foley et al., 2005).

| Habitat suitability and immune genetic diversity
We found that inbreeding coefficients (F IS ) were significantly higher in areas identified as low suitability habitats for adaptive,  (Wahlund, 1928) is unlikely to explain this result, as both types of markers would be equally affected. Inbreeding caused by a recent population decline or higher variation in family size as a consequence of the reduced habitat quality would also be expected to be reflected by neutral markers. A technical artefact, such as a higher probability for the presence of null alleles at the MH II in samples from low suitability habitats, seems unlikely, as samples were prepared in duplicates and random order (see Huml et al., 2018). Taking other potential drivers, such as effective population size and management background into account in the Random Forest analysis, habitat quality best explains the observed patterns of F IS at immune genes. There is evidence that the prevalence of infectious disease and the susceptibility of hosts can be increased under unfavourable environmental conditions (Austin, 2007;Schmidt-Posthaus & Wahli, 2015).
This includes metal pollution and eutrophication (Shirakashi & El-Matbouli, 2010;Wedekind, Gessner, Vazquez, Maerki, & Steiner, 2010), which were also identified to be among the main factors impacting on habitat suitability for grayling in this study. While our small sample size (N = 10) demands cautious interpretation of our results, we suggest directional selection pressure on MH II genes as a plausible explanation. This could be due to a higher prevalence of opportunistic pathogens under conditions of reduced habitat quality or environmental stress (Austin, 2007;Boutin, Bernatchez, Audet, & Derôme, 2013;Schmidt-Posthaus & Wahli, 2015), resulting in significantly higher F IS observed at the MH II for low-quality habitats, but with no observed effect of habitat quality for microsatellites. This is of high conservation relevance as periods of increased pressures of disease-mediated directional selection can lead to significant losses of immune genetic diversity and the future potential to face pathogens (Coughlan et al., 2006). However, a direct measure of selection pressures, for example through the assessment of shifts in the microbiome composition in grayling in F I G U R E 4 Projections of habitat suitability estimates for 2050 for RCP 2.6 (left) and 8.5 (right): Averages across 100 replicates are shown (the mean standard deviation was 0.06 for both the 2.6 scenario and the 8.5 scenario of change across the study area), with warmer colours indicating higher suitability; habitat classifications were done using the minimum training presence threshold (MTP) and equal training sensitivity and specificity (ESS) threshold

| Future climate predictions
Significant reductions in suitable range for grayling under future climate predictions were also demonstrated. As climate change occurs, species will need to respond either adaptively or by tracking range changes (Hoffmann & Sgro, 2011). That there is a scope for adaptive responses to climate change in salmonids has been shown by several studies. Eliason et al. (2011) showed that populations of sockeye salmon (Onchorhynchus nerka) exhibit different tolerance limits to temperatures, reflecting historic temperature ranges of their habitat. An adaptive population divergence in response to temperature has been shown by Kavanagh, Haugen, Gregersen, Jernvall, and Vøllestad (2010)  Maximum summer temperature was the parameter subject to most change under future projections ( Figure S2). While air temperatures generally show a linear relationship to water temperatures, the latter are also affected by flow, water volume, shading and wind shelters and deviations from linearity have been particularly shown, when maximum air temperatures exceed 25°C (Erickson & Stefan, 2000;Webb, Clack, & Walling, 2003). This was also evident from the river surface water temperature data set used here (Orr et al., 2010 (Sutton & Soto, 2012). Further investigations on local thermal refugia within grayling rivers that are expected to suffer from temperature increases, such as those in the South East, and management actions, such as the creation of riparian woodland zones are recommended (Malcolm et al., 2008). Also, the effect of groundwater abstraction should be evaluated in this context and in regard to increasing the risk of low summer flows, which can be critical for grayling (Riley et al., 2009).

| Mitigation of climate change impacts
One of the most important messages of this study is that the mod- Organic pollution and eutrophication associated with high concentrations of phosphate and nitrate are one of the major anthropogenic impacts on freshwater systems (Birk et al., 2012;Blabolil et al., 2016). Harmful effects for fish are caused by oxygen depletion, resulting from increased phytoplankton growth (Elshout, Dionisio Pires, Leuven, Wendelaar Bonga, & Hendriks, 2013).
Because of temperature dependencies, the impact of eutrophication is expected to further increase under climate change conditions (Moss et al., 2011). Metals have toxic effects on fish, mainly caused by a disruption of mechanisms important for ion regulation (Alsop & Wood, 2011). Metal pollution is considered a serious threat to freshwater ecosystems (Förstner & Wittmann, 2012) and is associated with historic mining activities within the UK (Macklin, Hudson-Edwards, & Dawson, 1997).

| CON CLUS IONS
Reducing non-climate-related environmental stress has been highlighted to be among the most important management action in the face of climate change (Heller & Zavaleta, 2009). This is stressed by the findings on low habitat suitability impacting on adaptive genetic variation and likely evolutionary potential, which is thought to be capable to promote climate change adaptation in salmonids (Eliason et al., 2011).
While this study does not try to give accurate predictions of the future distribution of grayling, its main goal was to show the potential of different habitat improvement strategies to increase habitat suitability for grayling under conditions of climate change and to give specific suggestions on local priority actions. One caveat of our study is that we used Maxent as a single modelling approach. While Maxent generally performs well in comparison with other techniques (Aguirre-Gutiérrez et al., 2013), the choice of modelling technique has been identified as a major source of uncertainty in SDMs (Thuiller, Guéguen, Renaud, Karger, & Zimmermann, 2019). Therefore, an ensemble SDM approach, which combines several techniques, may be a better alternative to comprehensively quantify uncertainty of predictions (Thuiller et al., 2019). Further, while in this study we successfully identify a number of non-climatic parameters that can be targeted to mitigate negative effects of climate change, we did not exhaustively evaluate all variables that may impact on habitat suitability of grayling. For example, here we focussed on extreme conditions of river flow (Q10 and Q95), given their particular relevance for the survival of juvenile salmonids (Riley et al., 2009;Warren et al., 2015).
However, seasonal patterns and variation in the flow regime, as well as the relevance of groundwater influxes, may be another important aspect to consider in future studies. Also, habitat fragmentation through hydropower development and dam construction has been shown to adversely affect grayling populations and maintenance of population connectivity is considered another important conservation priority (Junge, Museth, Hindar, Kraabøl, & Vøllestad, 2014).
Effective conservation of grayling across the higher latitudinal range of the species distribution, such as the UK and Scandinavia, is particularly warranted to safeguard the species, as the effects of climate change are expected to be less drastic than across its lower latitudinal range. The importance of evaluating priorities and invoking management actions on habitat improvement within the continental distribution of grayling is clear, given an expected stronger effect of climate change on habitat suitability for grayling at lower latitudes.
In summary, in the UK rivers should be managed to reduce levels of orthophosphates, copper and nitrate to ensure the best prospects for grayling under a changing climate. Evaluation of non-climate parameters should become routine in species distribution modelling for conservation management in a climate change context.

ACK N OWLED G EM ENTS
We thank the Grayling Research Trust and Manchester Metropolitan University for funding for this project. Thanks are also given to Prof.
Steven Weiss and Dr. Henri Persat for their input and guidance throughout the study and to two anonymous reviewers for insightful comments on an earlier version of the manuscript.

DATA AVA I L A B I L I T Y S TAT E M E N T
Occurrence data used in this study and an R-script used to fine-tune model parameters are available on Dryad (https ://doi.org/10.5061/ dryad.msbcc 2ftz). All other data used are available publicly or upon request from the sources given in the manuscript.

S U PP O RTI N G I N FO R M ATI O N
Additional supporting information may be found online in the Supporting Information section.