Increasing Cervidae populations have variable impacts on habitat suitability for threatened forest plant and lichen species

Large herbivores play a key role in temperate and boreal forest ecosystems. Cervidae (deer) population densities and community structure have undergone drastic changes in many parts of the world over the past decades, often with deer populations increasing. Many studies show impacts of Cervidae on multiple ecosystem properties, including vegetation and biodiversity, at local spatial scales. At larger spatial scales, however, impacts of changing Cervidae populations on forest ecosystems are less known. Although both abiotic and biotic dimensions contribute to shaping species’ niches, abiotic variables are generally given prominence when modelling species habitats and ranges. This is despite biotic changes, including changes in trophic structure, being an important component of global environmental change. In this study, we examined the potential contribution of Cervidae densities to the habitat suitability for rare plant and lichen species across the temperate and boreal forests of Norway, where cervid densities have increased over the past 60 years. We also examined how these changes in herbivore communities may have shaped habitat suitability for rare lichens and plants and discuss the results in light of continuing shifts in herbivore assemblages. We ran habitat suitability models for 47 species of rare plants and lichens, which were selected based on herbivory reported as a criterion for placement on the national red list for species. Climate (temperature and precipitation), forest (forest type and productivity), soil pH and Cervidae densities (moose Alces alces, red deer Cervus elaphus and roe deer Capreolus capreolus) were used as independent variables. Densities of one or more of the three Cervidae species were inferred to be associated with the distribution of 14 (ten lichen, one bryophyte and three vascular plant species) of these 47 species. We found a range of habitat suitability associations with Cervidae densities, including positive, negative and hump-backed responses. Increases in Cervidae densities over the past 60 years may have led to different spatial trends in habitat suitability across the 14 species. Our results suggest that Cervidae densities are associated with the distribution of rare forest plant and lichen species differently at large spatial scales; experimental studies should test the causality of these associations. If causal, this implies that Cervidae management should find a balance between high and low densities to conserve several plant and lichen species. The preponderance of epiphytic lichens species, for which habitat suitability was associated with Cervidae densities, calls for field studies to focus on Cervidae impacts on forest lichens.


Introduction
Large herbivores are important drivers of ecological and ecosystem dynamics (Danell et al., 2006). Since large herbivores often have ecosystem engineering roles, recent losses of herbivore species and the biases in species loss towards large-bodied herbivores, are assumed to have consequences for ecosystems (Estes et al., 2011;Ripple et al., 2015). In many parts of the world, herbivore communities today are dominated by homogenous livestock populations and are functionally poor compared to native herbivore assemblages (Gill, 2015;Svenning et al., 2016). In parts of Europe, however, agricultural abandonment of outlying land, previously used for extensive agricultural purposes, has led to reductions in livestock density and replacement by wild herbivores; low predator densities, increases in woody plant cover and climatic changes have attributed to this shift (Navarro and Pereira, 2012;Espunyes et al., 2019;Speed et al., 2019).
In temperate and boreal forest ecosystems, densities of Cervidae species (i.e. deer species; referred to as cervids hereafter) are increasing T in response to changing management of both wildlife populations and agricultural and silvicultural systems (Côté et al., 2004). Increasing cervid populations affect forest structure, biodiversity and ecosystems (Côté et al., 2004;Kolstad et al., 2017;Boulanger et al., 2018). Variable effects of cervid browsing on forest plant species and communities have been observed (Fuller and Gill, 2001;Rooney and Waller, 2003;Speed et al., 2014), and there is no clear consensus of how plant species alpha diversity responds to cervid browsing (Bernes et al., 2018). The lack of a clear effect of browsing cervids on plant diversity is likely due to the cooccurrence of both direct (e.g. browsing, trampling) and indirect (e.g. altered competitive interactions, dispersal) impacts of cervids on plant communities (Bernes et al., 2018;Kolstad et al., 2018;Kolstad et al., 2019). The meta-analyses undertaken by Bernes et al. (2018) showed that the species richness of vascular plants (herbs) and bryophytes were greater at higher levels of herbivory. The path analyses by Kolstad et al. (2019) demonstrated that cervid impacts on species richness are likely to be mediated through other biotic variables, including densities of dominant species and shrub biomass.
Cervids have been shown to affect the abundance or presence of plant species at a local scale (Côté et al., 2004;Speed et al., 2014;Bernes et al., 2018;Kolstad et al., 2019). Yet, cervid population trends and impacts tend to operate at relatively large spatial scales (Herfindal et al., 2015;Angelstam et al., 2017). Therefore, the large-scale habitat suitability of plant species may potentially be affected by changes in cervid populations. Such a question can be addressed through habitat suitability modelling. This approach is usually undertaken using abiotic variables alone (e.g. climate). However, niche theory acknowledges the importance of biotic variables to species' distributions. Hence, there is a need to incorporate biotic variables in models of habitat suitability in order to model species' Eltonian niches (i.e. also including biotic dimensions of the niche; Wisz et al., 2013;Trainor and Schmitz, 2014). Previous studies have shown that spatial maps of large herbivore densities can be used to model habitat suitability of plant species influenced by large herbivores Speed et al., 2019). Speed and Austrheim (2017) showed that large herbivore densities (sheep and reindeer) were associated with the distribution of rare vascular plants in mountain ecosystems in Norway. However, the greatest changes in Norwegian herbivore communities have occurred in forested ecosystems. While grazing livestock dominated lowland forests in the mid-20th century, browsing cervids became dominant in the late 20th century and into the 21st century . Norwegian forests are habitat for almost half of the national red listed species, far more than any other ecosystem (Henriksen and Hilmo, 2015). Although fungi and invertebrates are the dominant forest taxa on the red list for species, there are also many threatened species of lichens, vascular plants and bryophytes in the forest ecosystems of Norway (Henriksen and Hilmo, 2015). In this study, we investigate whether forest cervids are associated with the habitat suitability of rare plants and lichens across the boreal forest ecosystems of Norway. We also address how recent changes (1940s vs present) in cervid densities across Norway may have altered habitat suitability for these species.

Species occurrence data
The study area was defined as areas of forest within Norway for which forest productivity has been classified, as defined using the AR50 land cover map (Norwegian Forest and Landscape Institute, 2007). The potential study species used here were limited to species of vascular plants, bryophytes and lichens on the Norwegian Red List for species in 2015 (Henriksen and Hilmo, 2015) that (i) were red listed due to landuse change with herbivory specified in description of impact, and (ii) had forest identified as the main (or one of the main) habitats for the species. For these 64 species (Table S1), all georeferenced species occurrence records located in Norway were downloaded from the Global  Table S3.
Biodiversity Information Facility: lichens (GBIF, 2018a); bryophytes (GBIF, 2018b); vascular plants (GBIF, 2018c). Only records with coordinate uncertainty less than or equal to 1414 m (the diagonal of a 1 km grid cell) and species with at least 20 records within forests were used in further analyses. This resulted in a list of 47 species (bold text in Table S1), of which two are bryophytes, 20 lichens and 25 vascular plants.

Environmental data
To characterise Norwegian climate, we downloaded WorldClim bioclimatic variables for Norway at 30 arc second resolution (Fick and Hijmans, 2017). This dataset was resampled to a 1 km grid using the nearest neighbour approach. We selected two variables to represent the principle climate axes of this dataset, namely annual precipitation and mean temperature of the warmest quarter. See Speed and Austrheim (2017) for further information on the climatic variable selection. Soil pH was used to characterise soil and bedrock across Norway. We used soil pH at 5 cm depth, from Soil Grids at a 1 km resolution (Hengl et al., 2014).
We used the AR50 land cover map of Norway to quantify forest type (coniferous, deciduous or mixed) and forest productivity (low, medium or high). The metabolic biomass (biomass 0.75 to account for allometric scaling of body size) of forest cervid species (roe deer Capreolus capreolus L., red deer Cervus elaphus L. and moose Alces alces L.) in Norwegian municipalities in contemporary (2015) and historic (1949) times (estimated from hunting and livestock statstics; Speed et al., 2019) were used to characterise cervid densities. While reindeer (Rangifer tarandus L.) are present and forage within limited parts of the study region, they are predominantly a tundra species in Norway, and were therefore not included in this study. Livestock species were not included due to very low contemporary densities in forest regions . All environmental (abiotic and biotic) variables are shown across Norway in Fig. S1. The temporal change (difference between 2015 and 1949 densities) in cervid biomass densities across Norway are shown in Fig. S2. We checked all environmental variables for potential collinearity using visual plots, correlation coefficients and variance inflation factors (VIF). Correlation coefficients between all  Table S2 for variable importance of all variables (i.e. including also climate and land cover).
pairs of environmental variables were all below 0.55 (Fig. S3) and all VIFs were below 2.7 (Table S2), substantially below the critical thresholds of 0.7 and 10 respectively (Dormann et al., 2013). We therefore proceeded with these eight environmental variables.

Modelling
There are a wealth of methods available to model species' habitat suitability, each with their own advantages and disadvantages (Guisan et al., 2017). The best choice of habitat suitability model is likely to vary between species. Thus, when modelling multiple species, the preferable option is often to use an ensemble modelling approach. This involves fitting models with a range of methods and averaging the output across the modelling methods (Guisan et al., 2017). Since the current study sets out to model habitat suitability of many (47) species, we selected the ensemble modelling approach here. Multiple models were developed for each species using six commonly applied methods. These were general linear models (GLM), general additive models (GAM), random forest (RF), boosted regression trees (BRT), mixture discriminant analysis (MDA) and flexible discriminant anlaysis (FDA). We selected this group of models as a broad and robust range of approaches with good performance and lack of overcomplexity (Breiner Fig. 4. Response curves plotting habitat suitability against cervid biomass densities, for species where one of the cervid biomass variables had a variable importance of > 0.25. Mean and standard errors of predictions shown with solid and dashed lines respectively. Note that y axis scales differ across plots while the x axis scale is the same within each column (cervid species). Lichen and plant species are ordered down each column, alphabetically by species within higher taxa (lichens followed by vascular plants). Some species had variable importance values of > 0.25 for both moose and red deer so response curves for both cervids are shown for these. Cervid biomass densities are metabolic biomass. Guisan et al., 2017).
The ensemble of models were fitted and evaluated using the package sdm (Naimi and Araújo, 2016) in the R environment (R Core Team, 2018). VIFs were estimated using the usdm package (Naimi et al., 2014). Since the species data used is presence-only (i.e. we do not have true species-absence data), we also created a background dataset. To account for spatial bias in species observations in Norway, we used weighted background data. This was a sample of 1000 points in the forested areas of Norway weighted by the spatial distribution of plant species occurrence records (Phillips et al., 2009). This ensures that model predictons contrast the environment at locations with presences of the selected species against regions of the study area where plant occurrence data has been sampled, rather than against unsampled areas. The following iight environmental variables were fitted: mean summer temperature, annual precipitation, soil pH, forest type, forest producvtivity, moose biomass density, red deer biomass density and roe deer biomass density (2015 densities). We used densities of the forest cervids individually, since they differ in diet (moose and roe deer are browsers, and red deer is an intermediate feeder, Hofmann, 1989) and distribution within Norway .
All models were replicated using 5-fold cross validation (Guisan et al., 2017) and evaluated in terms of the area under curve (AUC) statistic. Cross validation partitions the presence dataset into a test and training datasets. Models are fit with the training dataset and evaluated against the test dataset in order to evalute the model against 'independent' data. This is repeated k times (here k = 5) with 1/k presence observations as the training dataset in each cross validation (Hijmans, 2012).
The variable importances of the environmental variables were estimated as mean and standard errors from the five fold cross validation runs. These summarise the relative contribution of each variable to the models. Cervid biomass density response curves were plotted for each red listed species if a cervid biomass density variable importance was > 0.25 averaged across models and replicates (with 8 independent variables, the null expectation is a variable importance of 0.125). Model predictions (averaged across methods and replicates) were plotted over the whole of Norway (assuming unlimiting dispersal of every species). Model predictions were made using contemporary (2015) since that is temporally closest to the time of sampling of the majority of the plant and lichen occurrence records. Model predictions were also projected against historic (1949) cervid densities to examine the changing habitat suitability with altered cervid densities.

Results
Across the 47 species for which habitat suitability models were run (Table S1), model fits were generally very good (median AUC 0.92, range 0.66 to 0.99, Fig. 1, Table S3). Model fit was higher for bryophytes (mean AUC = 0.98, sd = 0.002) than for either lichens (0.88 ± 0.07) or vascular plants (0.90 ± 0.06). The GAM models  Table 1 Overview of habitat requirements and potential pathways for cervids to influence the distribution of the 14 species modelled in this study. were the best fitting (0.952 ± 0.053), whereas the FDA models had the lowest fit (0.885 ± 0.075, Table S3).
Across species, the climate variables of mean summer temperature and annual precipitation contributed most to the habitat suitability models (Fig. 2a, Table S3). These were followed by soil pH. The variable importances of cervid biomass densities were, on average across species, comparable to those of forest productivity and forest type (Fig. 2a). The variable importances of temperature and precipitation were negatively correlated (Fig. 2b). The variable importance of roe deer biomass was positively correlated with the variable importance of temperature, while the importance of red deer biomass was negatively correlated with the importance of temperature (Fig. 2b).
Of the three cervid species biomasses, the highest variable importance was found for moose biomass in the case of the vascular herb Campanula barbata (variable importance mean across methods and replicates ± SEM of 0.65 ± 0.01, Fig. 3, Table S4). The species for which the variable importance of red deer was the greatest was the foliose epiphytic lichen Pectenia cyanoloma (0.47 ± 0.01) and for roe deer it was the vascular herb Vicia cassubica (0.35 ± 0.01).
For the 14 species for which densities of one of the three cervids was an important (> 0.25) driver of species habitat suitability, prediction maps extrapolated across the whole of Norwegian productive forests were made using both recent (2015, Fig. S4) and historic cervid densities (1949, Fig. S5). For some species (Gyalecta flotowii, Herbertus stramineus, Menegazzia subsimilis, Opegrapha vermicellifera, Pectenia cyanoloma, Staurolemma omphalarioides and Thelotrema macrosporum), changes in cervid biomass density have led to increases in model predictions in the west of Norway and decreases in the south east of the country (Fig. 5). For other species, increases in habitat suitability with changing cervid biomass densities are seen in the south east (Scorzonera humilis and Vicia cassubica). Changes in herbivore densities lead to increases in predicted habitat suitability of Campanula barbata in central Norwegian forests, and for Rinodina distuncta along coastal forests of central Norway. Model predictions of Galium sterneri and Phaeophyscia kairamoi decreased following the changes in herbivore densities, particularly in south east Norway (Fig. 5).

Discussion
In this study, we set out to (i) establish whether cervid densities were associated with the habitat suitability of rare forest plants and lichens in Norway, and (ii) investigate whether historical changes in cervid densities may have increased or decreased habitat suitability for these species. Our findings show that cervid densities are indeed associated with habitat suitability of some forest plant and lichen species (Fig. 3). Positive, negative and hump-backed associations between cervid densities and habitat suitability were all observed (Fig. 4), and recent changes in cervid density may have had differential impacts on habitat suitability between species and in space (Fig. 5). While field experiments and meta-analyses have demonstrated impacts of cervid browsing on diversity and specific plant taxa at the local scale (Speed et al., 2014;Bernes et al., 2018;Boulanger et al., 2018;Kolstad et al., 2019), our study indicates associations between cervid densities and plant and lichen species' habitat suitability at far larger spatial scales.
Our study species were plants and lichens on the national Norwegian red list in forest habitats with herbivory specified in the description of the threat status. Our results indicate that the habitat suitability of 14 of the 47 study species was associated with the biomass density of one or more cervid species. Given that the study species were selected from the red list on the basis that changing herbivory lay in part behind their threatened status, this number is perhaps lower than might be expected. This low proportion may relate to the spatial scale and the relative roles of wild cervids and domestic livestock.
Potential pathways by which cervids may influence the distribution of plant and lichens species may be direct (feeding, trampling) or indirect (e.g. changed vegetation structure or composition): these potential pathways are outlined in Table 1. The majority (10 of 14) of the study species for which cervid densities contributed to habitat suitability were lichens. This may be surprising since the forest cervids (i.e. not Rangifer tarandus) included in this study rarely feed upon lichens (Tremblay et al., 2005;Latham and Boutin, 2008). However, just as the impacts of cervid browsing on plant diversity are often indirect (Beguin et al., 2011;Bernes et al., 2018;Kolstad et al., 2019), the impact of cervid densities on lichens is also likely to be indirect (Moore and Crawley, 2014). Trees are the substrate for many epiphytic species, including most of the lichens in this study (Table 1, Table S1). Therefore, the impacts of browsing cervids on the recruitment of trees (particularly deciduous species, Myking et al., 2013;Kolstad et al., 2018) is likely to be the dominant mechanism behind these findings (Table 1). Lichen species' habitat suitability was influenced more often by red deer than roe deer with moose being intermediate. Epiphytic lichens responded negatively or with hump-backed responses to moose density, but positively (or hump-backed) to red deer density. This may be due to the negative effect of moose browsing on deciduous tree abundance (Kolstad et al., 2018) and the frequency of bark-stripping by moose (e.g. Bergqvist et al., 2001); see Table 1 for further information. Conversely, red deer have been shown to prevent saplings from shading lichens on mature deciduous trees (Moore and Crawley, 2014).
The distribution of epiphytic lichens relates to the distribution of host tree species (Hedenås et al., 2003). Epiphytic lichens are generally more abundant and species rich in old-growth forests, and are particularly widespread on deciduous species (Table S1). Therefore, disturbances to host tree species mediated through intense browsing or forestry management will impact the distribution of epiphytic lichens (Table 1). Cervid densities are strongly tied to forest management through forage availability (Milner et al., 2013;Wam et al., 2016), and habitat suitability models are correlative. We cannot, therefore, ascertain a causal link between cervid densities and habitat suitability of threatened plant and lichen species. The positive associations (and the increasing side of the hump backed relationships) may also be caused by the plant or lichen species sharing habitat preferences with the herbivore species, and not a direct impact of herbivory. This pattern may appear independent of any direct interaction between the cervids and the plants or lichens. Still, shared habitat preferences do indicate the potential for a cervid to influence (in any direction) the habitat suitability for a plant or lichen species. Nevertheless, our findings show that field-based research into cervid impacts on lichen species in boreal and temperate forests is clearly warranted and should be a priority for species of conservation concern.
The four vascular plant species for which habitat suitability were inferred to be affected by cervid densities were all herbaceous perennials: Campanula barbata (positive response to moose density) Galium sterneri and Scorzonera humilis (positive responses to roe deer density) and Vicia cassubica (unimodal response to roe deer density (Fig. 3). It should be noted that these four species have narrow geographic distributions (Fig. 5). All four are found in both grasslands and forests. The roe deer and moose are both browsers, so the increases in habitat suitability with increasing cervid densities within forests may be due to increasing canopy openness with increased cervid densities (Mathisen et al., 2010;Eichhorn et al., 2017;Kolstad et al., 2018). The increase in habitat suitability for rare non-forest specialists with increasing cervid densities mirrors the findings of Boulanger et al. (2018) who found lower richness of herbaceous species, particularly of non-forest species, inside ungulate exclosures.
Our study demonstrates the importance of including biotic factors, in addition to abiotic factors, when assessing habitat suitability or species distributions (Wisz et al., 2013;Trainor and Schmitz, 2014;Speed and Austrheim, 2017). This allows species' Eltonian niches to be more fully quantified, but also to infer how habitat suitability may be affected by a changing environment. The habitat suitability of thermosensitive species (species for which temperature was a strong determinant of habitat suitability) was also affected by roe deer biomass densities (Fig. 2b). This may be due to the high variation in roe deer biomass densities in warmer parts of Norway (Fig. S3). In contrast, species for which habitat suitability was affected by red deer biomass tended to be thermo-insensitive. The implication of this is that the management of threatened plant and lichen species will require a greater focus on roe deer densities, as the climate warms. However, the predicted range expansion of red deer in Norway is also likely to influence future habitat suitability (Milner et al., 2006;Jarvie and Svenning, 2018).

Implications for management in a changing world
The increase in large browsers and intermediate feeders represent a strong biotic change for forest ecosystems, particularly set against intensifying forestry and climate change. Changing herbivory levels have been inferred to be important causes of threat for red-listed forest species in Norway, but studies at relevant scales are lacking. Our study shows that red listed plant and lichen species analysed here demonstrate different associations with cervid biomass densities (Fig. 4). Our results show that the effects of six decades of herbivore change seem to have had different impacts on species across Norway (Fig. 5). This suggests that there is no single cervid management option that will increase habitat suitability across all threatened plant and lichen species. When positive, negative and humped-backed responses of rare plants and lichens to cervid density are prevalent, potential compromise approaches to managing cervids could be (i) maintaining intermediate densities of cervids or (ii) diversifying the range of cervid intensities in space. A similar impact of ungulate herbivores was found in the case of threatened vascular plants and ungulate densities in mountain regions of Norway . However, livestock (sheep and semi-domesticated reindeer in mountain regions) can generally be managed at smaller spatial scales (farms, grazing and herding districts) than cervid species. The cervids are generally wider ranging and management processes in our study region operate hierarchically (Danielsen, 2001).
In this study we used cervid densities as independent variables. We did not include livestock densities, although the overarching change in forest large-herbivore communities over the last six decades can be characterised as a replacement of grazing livestock by browsing forest cervids (a form of rewilding; Speed et al., 2019). However, since future changes in herbivore assemblages in forests in Norway are likely to involve further changes in wild herbivore communities (Jarvie and Svenning, 2018) and a continued near-absence of livestock, the focus on wild species is most relevant also in the future.
We have considered changes in herbivore assemblages over a sixdecade period in this study. This is relatively long-term for ecological studies. However, it has been argued that the role of herbivores in forest communities should be contextualised with a far longer perspective, accounting for the evolutionary history of plant-herbivore interactions developed when herbivore assemblages were more functionally diverse (Fløjgaard et al., 2018, i.e. from the Pleistocence and before). In Norway, there has been a gradual expansion of red deer distribution and a decline in moose abundance throughout the Holocene, with increasing anthropogenic limitation of populations in more recent times (Rosvold et al., 2013). This is consistent with widespread loss of browsing herbivores globally (Janis et al., 2000;Bakker et al., 2016;Malhi et al., 2016). The dominance of positive responses of habitat suitability to cervid density found in our study (Fig. 4) suggests that low densities of wild browsing herbivores in forests over past decades may indeed have been a threat to certain species. However, current high densities of cervids in regions with a functional absence of predators, can limit the recruitment of keystone deciduous tree species, and this is also of biodiversity concern (Kolstad et al., 2018), particularly given the occurrence of multiple threatened epiphytic species. There is thus a need for cervid management to find a middle line (Milner et al., 2013) as well as a requirement for field based research to increase understanding of cervid impacts on forest lichens.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.