Quantifying the distribution and potential biotic interactions between deer and flora using species distribution modelling

ABSTRACT

Invasive species are ranked as one of the leading drivers of global biodiversity loss. To mitigate their impact, we must understand the future risks caused by invasive species, particularly to flora of conservation concern. Here we used species distribution modelling (SDM) to project the current and future (RCP45 and RCP85 2050) distributions of four deer species and 13 plant species of conservation concern for the island of Ireland, quantifying changes in distributions and overlap. Large areas of suitable habitat for the deer species were predicted with high accuracy across all counties, with future climate scenarios identifying an expansion in sika deer distributions and a decrease in muntjac and fallow deer distributions. Red deer declined under the moderate climate change scenario but increased under the worst-case projection. Future projections predicted the (local) extinction of six (out of 13) endangered and vulnerable plant species. An expansion in distributions was observed for four plant species; however, these areas had large overlap with the future predictions of deer, placing further pressures on these plant species. These findings suggest that targeted conservation and management measures are required to alleviate the pressures on ‘at-risk’ plant species due to grazing from native and non-native deer.

KEYWORDS:

• climate change
• conservation
• forestry
• invasive species
• species distribution modelling

Introduction

Climate change and concomitant human activities are placing ecosystems under unparalleled pressures (Chen et al. 2011; Linders et al. 2019; Pejchar and Mooney 2009). The recent Intergovernmental Science Policy Platform on Biodiversity and Ecosystem Services (IPBES) report (Díaz et al. 2019) identified that over 1 million species are threatened with accelerated extinction. This deterioration is directly linked to human activity, with the spread of invasive species ranked as one of the primary drivers of this loss. Invasions by non-native species directly and indirectly affect biodiversity and, as a result, ecosystem functions (Linders et al. 2019). In order to mitigate the potentially devastating scenarios for biodiversity and ecosystems, it is imperative to better understand the future extinction risk caused by invasive species (Bellard, Bernery, and Leclerc 2021).

Deer species, including sika (Cervus nippon), muntjac (Muntiacus reevesi) and fallow (Dama dama), are preponderant examples of invasive species in Europe (Saggiomo, Esattore, and Picone 2020). Following their introduction throughout Europe, sika deer have undergone rapid population increase (Carden et al. 2010; Homolka and Heroldová 2003), becoming naturalized in some countries with extensive and invasive populations found in Germany, the Czech Republic, United Kingdom and Ireland (Carden et al. 2010; Homolka and Heroldová 2003; Pérez-Espona, Pemberton, and Putman 2009). Muntjac deer are also expanding rapidly across Europe, with predictions that this species will expand its range to encompass every European country (Ward, Richardson, and Mergeay 2021). Large-scale browsing and grazing by non-native deer species has cascading effects on biodiversity due to direct interactions with other species, including plants, invertebrates, smaller species and other herbivores (Akashi and Nakashizuka 1999; Homolka and Heroldová 2003; Rooney and Waller 2003). For example, Pollard and Cooke (1994) identified that muntjac deer preferentially grazed on lower leaves of honeysuckle (Lonicera periclymenum) which are the preferred egg laying sites of the white admiral butterfly (Ladoga Camilla), leading to their population decline. This subsequently raises numerous management and control issues (Carden et al. 2010).

Overabundance of deer, and the additional effects that are associated therewith, such as effects on biodiversity or ecology, is extensively discussed in their native ranges (Iida, Soga, and Koike 2016; Kato and Okuyama 2004; Ohashi et al. 2014; Uno et al. 2019). For example, Iida, Soga, and Koike (2016) studied the effects of sika deer overpopulation on flora and invertebrates in Western Hokkaido in Northern Japan. Two sites were investigated, the shore of Lake Troya (a control site) and Nakanoshima Island (experimental site). During the 1950’s and 60’s, three sika deer were introduced on the island for tourism (Kaji, Koizumi, and Ohtaishi 1988). Over this period there was a significant increase in sika deer density on the island (Ikeda et al. 2013), which resulted in a dramatic decrease in both understory vegetation and beetle diversity. This study quantifies how populations of sika deer can alter native species richness, abundance and likely the ecological functions of an area.

Saggiomo, Esattore, and Picone (2020) conducted a recent bibliometric analysis which presented all published scientific work on sika deer that has been accomplished to date. It was found that most studies came from Asia, Europe, Africa and North America. The literature on native populations focused mainly on the effects of sika overpopulation and its effects on vegetation such as the dwarf bamboo, whereas the literature on sika deer in their introduced range mainly explored hybridization (Bailey 2013; Gélin et al. 2019; Huxel 1999). This is in contrast to the impact of muntjac and fallow deer on native and non-native vegetation that has been widely researched (Pollard and Cooke 1994; Cooke and Farrell 2001, Morecroft et al. 2003; Husheer and Frampton 2005; Marozas et al. 2009; Fuller and Gill 2010; Cooke 2021; Cunningham et al. 2022). The consensus so far is that deer species, particularly non-native species, negatively impact vegetation diversity in regions, with cascading effects on insect and bird species. However, studies have not been conclusive, with differences across species and geographic regions. Studying one species in isolation can lead to species-specific results that are not generalizable to the wider and connected ecosystem (Meyer et al. 2019), particularly for organizations and practitioners managing these ecosystems. Subsequently, extensive studies exploring the impact of multiple deer species distributions and populations on the flora of non-native ecosystems are lacking (Saggiomo, Esattore, and Picone 2020).

In Ireland, between 1978 and 2008 there was a range expansion of 565%, 353% and 174%, for red, sika and fallow deer, respectively (Carden et al. 2010). The increase in both human and deer populations has heightened the potential for human-deer conflict in recent years, with research predominantly focusing on their impact to the agriculture and forestry sectors (Hawthorne, Colombaroli, and Mitchell 2021; Liu and Nieuwenhuis 2014; Wilson, Wirz, and Byrne 2021), deer-vehicle collisions (Liu, Nieuwenhuis, and McCullagh 2018), their role as disease vectors and/or hosts (Crispell et al. 2020; Kelly, Mullen, and Good 2021; McDonough and Holloway 2020), hybridization with the native red deer (Carden et al. 2010; Gélin et al. 2019; Huxel 1999; McDevitt et al. 2009; Smith et al. 2014), and their impact on the native vegetation (Perrin, Mitchell, and Kelly 2011), although the latter is the most depauperate research area at a national scale.

The diet of deer species is varied, with that of the non-native sika and native red deer overlapping in south-west Ireland, causing competition (Burkitt 2009). Although the bark from saplings and grasses are the preferred diet of sika, they have been known to browse and graze unconditionally on forbs, berries, heather and other vegetation (Mann and Putman 1989; Yokoyama, Kaji, and Suzuki 2000). Muntjac deer have also been found to graze plant species that are avoided by other deer species, such as bluebells (Hyacinthoides non-scripta), which can shift the equilibrium of native ecosystems (Dolman et al. 2010). In Ireland, the Red Data List of Irish Vascular Plants (National Botanical Gardens 2017) identifies that there are currently seven critically endangered plant species, 52 endangered species and 69 vulnerable species. Subsequently, any expansion of non-native (or even native) deer distributions in the future may place further pressures on these plant species that are already at risk due to climate change and habitat degradation. Despite the potential impact of invasive species, many of their effects can be mediated through vegetation abundance, species richness and diversity (Linders et al. 2019). To that end, non-native deer must be focused on as they are a key species for conservation and management efforts, particularly in supporting the Irish Biodiversity Action Plan (2019). To achieve sustainability in ecosystems and biodiversity, the prevention of any further spread, or establishment of an invasive species has been recognized as a more efficient and a less costly method of control than large-scale culling or barrier control (Barbet-Massin et al. 2018). Moreover, if it is known where future spreads may occur, then adaptation efforts can be implemented to ensure minimal impact on any native and/or flora of conservation concern. Therefore, the initial step for effective management is quantification of how deer and flora distributions will change in the future (Carden et al. 2010).

Species distribution models (SDMs) are one of the most important GIScience application areas (Burch et al. 2020; Laffan, Skidmore, and Franklin 2014; Miller et al. 2020), and a powerful spatial analytical tool for studying the geographic distribution of a range of taxa (Franklin 2010; Peterson et al. 2011), providing a methodological framework for researchers and practitioners to quantitatively assess the relationship between species distributions and environmental factors. SDMs project relationships in both environmental and geographic space using a variety of statistical methods and machine learning algorithms (Miller and Holloway 2017) and have been widely used for various applications, including projecting suitable aquaculture sites (O’Mahony, de la Torre Cerro, and Holloway 2021), forecasting global change impacts on biodiversity (Ehrlén, Morris, and Buckley 2015), locating nature reserves (Araújo et al. 2004), and predicting future ranges of invasive species (Barbet-Massin et al. 2018; McMahon et al. 2021). For example, McMahon et al. (2021) modelled the distribution of 15 non-native invasive grass and forb species of conservation concern in western North America. The aim of the study was to characterize the species environmental niches and identify potential invasion risks. The results showed that the species are likely to expand their ranges and coupled with continued increases in aridity, temperature and burned areas the risk of future invasions is high. Through the use of GIScience and SDMs, this study was able to aid in identifying possible new invasion sites and prioritize the prevention and control of the invasive species.

Despite the fact that non-native deer species are of ecological and economic concern in Ireland, coupled with the ability of SDMs to project current and future distributions, there is a paucity of research that focuses on predictions of deer distributions in Ireland, particularly related to native vegetation. Such information is vital in determining which locations require regulation, management, or conservation. Here we explore four main research questions with the overall aim of identifying the factors that drive deer distributions and to reveal areas at greatest risk to biodiversity loss for the island of Ireland: 1) What is the current suitable distribution of the four deer species currently established in Ireland? 2) What are the current suitable distributions of endangered/vulnerable plant species? 3) How will the distributions of the deer species and plant species differ under future climate change scenarios? 4) Where will the deer species and plant species future distributions overlap?

Methods

Data collection

Presence data for the four deer species and 13 endangered/vulnerable plant species were obtained from the Global Biodiversity Information Facility (GBIFa-q) using the rgbif package (Chamberlain et al. 2021) in R version 4.1.2 (R Core Team 2021). Following data cleaning as of December 2021, 309 records for C. nippon, 299 records for C. elaphus, 835 records for D. dama and 281 for M. reevesi were returned, while plant species were retained if they had at least 10 occurrence records. All data observations were returned for the period after 2000 only, meaning any observations prior to then were disregarded. Plant species were identified from The Red Data List of Irish Vascular Plants (National Botanical Gardens 2017) and are outlined in Table 1. In total, 13 plant species were used including 5 endangered and 8 vulnerable plant species. Most plant species in The Red Data List of Irish Vascular Plants had little or no occurrence records, meaning they were omitted from the analysis.

Table 1. Information on the species used in the modelling framework along with IUCN status.

Name	Scientific Name	Group	IUCN Status
Red deer	Cervus elaphus	Mammal	Least concern
Sika deer	Cervus nippon	Mammal	Least concern
European fallow deer	Dama dama	Mammal	Least concern
Reeves muntjac	Muntiacus reevesi	Mammal	Least concern
Chives	Allium schoenoprasum	Plant	Vulnerable
Nettle-leaved Bellflower	Campanula trachelium	Plant	ENDANGERED
Cornflower	Centaurea cyanus	Plant	ENDANGERED
St Dabeocs heath	Daboecia cantabrica	Plant	Vulnerable
Mackay’s Heath	Erica mackaiana	Plant	Vulnerable
Alder Buckthorn	Frangula alnus	Plant	Vulnerable
Wood Crane’s-bill	Geranium sylvaticum	Plant	ENDANGERED
Meadow Barley	Hordeum secalinum	Plant	ENDANGERED
Hairy Bird’s-foot-trefoil	Lotus subbiflorus	Plant	Vulnerable
Alpine Meadow-grass	Poa alpina	Plant	Vulnerable
Northern Yellow-cress	Rorippa islandica	Plant	Vulnerable
Subterranean Clover	Trifolium subterraneum	Plant	ENDANGERED
Spring Vetch	Vicia lathyroides	Plant	Vulnerable

All observations that were utilized within this study came from human observations, predominantly from citizen science programmes such as the National Biodiversity Data Centre and iNaturalist. For all but four of the species, coordinate uncertainty was reported for over 90% of the observations for each species, with the majority of the species reporting uncertainty at less than 1000 m, which is the approximate spatial resolution of the SDM. P. alpina, R. islandica and V. lathyroides had over 50% of their observations with a high uncertainty; however, these represented some of the smaller sample sizes. Recent research has identified that data quantity is more important than spatial and/or other bias in the predictive performance of SDMs (Gaul et al. 2020). Due to the low number of observations reported for these species, we opted for data quantity over spatial bias and uncertainty introduced by using these locations where uncertainty in spatial positioning was higher than the resolution of the SDM. Please see the Supplementary Material for a breakdown of this information of the GBIF observations.

Environmental data included climate and land cover data. Climate data were obtained from the WorldClim database (Fick and Hijmans 2017) using the sdmpredictors package (Bosch et al. 2017) in R version 4.1.2 (R Core team 2021). Corine land cover (EEA, 2018) data was used to calculate the percentage of arable, forest, grass, pasture, urban and water cover. All variables had a spatial resolution of 30 arc seconds, approximately 1 km. To represent the future climatic predictors, the Community Climate System Model version 4 (CCSM 4.0) was selected under the Representative Concentration Pathways (RCP) 4.5 and 8.5. RCP was chosen over the Shared Socio-economic Pathways (SSP) due to the unavailability of SSP at the required scale in the WorldClim database, which would hinder reproducible research. This model is a coupled climate model for simulating the earth’s climate maintained by the Climate and Global Dynamics laboratory at the National Centre for Atmospheric Research, USA (Gent et al. 2011).

In this study, RCP4.5 and RCP8.5 were selected to model the future distributions of four deer species and the 13 endangered/vulnerable plant species in the years 2050 (average predicted conditions of 2041–2060). RCP 4.5 is an ‘intermediate’ scenario with an emissions peak around 2040 and a steady decline thereafter. RCP 8.5 is the basis for the ‘worst-case’ scenario were emissions to continue to increase until 2100 (Meinshausen et al. 2011). We opted to assume that land cover will remain consistent in the future, which is most likely a simplification of reality. While there have been recent global products of future land use, such as ESRI’s Living Atlas Land Use 2050 data, this currently is not projecting in our opinion realistic covers. For example, the product predicts the vast removal of cropland across Ireland, being replaced with grassland, the conversion of Killarney forests into swampland and the conversion of large parts of Dublin into cropland. A recent study by Kang et al. (2022) found that ESRI’s current land cover data had the lowest overall accuracy when land cover was assessed. Therefore, in the absence of a realistic future projection, we feel that the assumption that land cover remains unchanged is justified and is standard practice when future land cover projections are missing in SDM research (e.g. Stanton et al. 2012).

Data analysis

SDMs generate a prediction of habitat suitability in both environmental and geographic space. The most common response variables used are data representing species presence-absence or presence-only (with corresponding pseudo-absences), which are coupled with a set of corresponding and georeferenced environmental variables (Franklin 2010). Species-environment relationships are then estimated using a range of statistical approaches that can then be extrapolated in space and time, generating a habitat suitability map that identifies the probability that the environmental conditions at that location are favourable for the study species. A categorical representation of habitat suitability is created by using a threshold that delineates between suitable-unsuitable, with habitat suitability values above and below this predicted as suitable and unsuitable, respectively. Figure 1 provides a methodological workflow of the SDM steps used in this study.

Figure 1. Conceptual workflow of the species distribution modelling framework implemented in this study.

GBIF provides information on presence-only, meaning we did not have a true consideration of absence. As such, 10000 pseudo-absences (or background points) were randomly drawn from the study area, following the established standards set out by Barbet-Massin et al. (2012) to ensure a binomial function could be applied to the statistical models. According to Barbet-Massin et al (2012, 2018), the results of SDM projections can be impacted by the location of pseudo-absences. Therefore, we decided to replicate data analysis three times with different pseudo-absence selections to account for such model uncertainty.

Three different statistical models were implemented within the biomod2 package (Thuiller et al. 2016) in R version 4.1.2 (R Core team 2021). The methods we selected for initial analysis included Generalized Linear Model (GLM), Maximum Entropy (Maxent) and Random Forest (RF). Studies have shown that different modelling techniques can produce different results for the same species and datasets, with ensembles across different statistical models suggested (Araujo and Guisan 2006; Elith and Graham 2009; Thuiller 2004). However, the results for GLM and Maxent were poor in comparison to RF (results not shown), and when coupled with recent research that has identified there is no particular benefit to using an ensemble approach across different statistical models over individually tuned statistical models (Hao et al. 2020), we opted to progress this research using only RF. RF constructs many decision trees using bootstrap samples (Liaw and Wiener, 2002) averaging these results into one projection. By default, RF can be modelled in the form of classification and/or regression trees (Zhang et al. 2019) and Artificial Neural Networks (ANN) (Lek and Guégan 2000; Ripley 1996).

Models were evaluated after parameterization to justify the acceptance of projections for their intended purpose (Araujo and Guisan 2006). We performed 3-fold cross-validation of the data by splitting the occurrence records and pseudo-absences using an 80:20 split into training and testing data following best practice in the SDM literature when independent test data are not available (Hijmans 2012; O’Mahony, de la Torre Cerro, and Holloway 2021; Peterson et al. 2011). The model’s discriminatory power between suitable and unsuitable was assessed using three different metrics which can be used to provide an accurate assessment of the models. The area under the curve (AUC) statistic of a receiver operating characteristic (ROC) (Hanley and McNeil 1982) is a threshold-independent metric and widely used in SDM research (Franklin 2010). It has a value range between 0 and 1 with a value below 0.5 deemed no better than a random selection and a value of 1 representing a highly accurate model. True Skill Statistic (TSS) measures the performance of models generating presence-absence predictions with values of 0 indicating no agreement and 1 indicating perfect agreement (McHugh 2012). Cohen’s kappa corrects the overall accuracy of model predictions by the accuracy expected to occur by chance (Allouche, Tsoar, and Kadmon 2006). It ranges from values −1 to + 1, where + 1 indicates perfect agreement and values of 0 or less indicate a performance no better than random (Cohen 1960).

For each pseudo-absence configuration, we therefore ran three different RF models to obtain a consensus distribution, calculating an ensemble forecast distribution. This meant that for each species, we had nine possible distribution projections. The projections from these individual models were ensembled using the median habitat suitability value. The accuracy values were also reported based on the mean of the nine individual models. The ensemble model uses a continuous probability distribution which was then discretized into a binary suitable-unsuitable map using the probability threshold that maximized the TSS value (Liu et al. 2005).

The same SDM procedure was undertaken for the 13 plant species. To determine plant species richness under current conditions and future scenarios RCP45 and RCP85 2050, the binary maps produced for each plant species were stacked. For example, if all 13 species were projected in the same 1 km grid cell, then the value would be 13.

The Biomod.RangeSize function within the biomod2 package was used to calculate the current and future distribution in square kilometres (km²) for each species. It was also used to calculate the overlap of the four deer future distributions with each plant species future distribution.

Results

We projected the four deer species distributions with high accuracy values across the island of Ireland for all accuracy metrics (Figure2, Table 2). The current habitat suitability of C. nippon is approximately 1257 km² (Figure 2(D), Table. 3) with suitable areas identified in all counties with the exception of Galway and Mayo. There are four large geographic areas of continuous suitable climate and habitat, these are the Wicklow mountains, the upland areas of County Kerry, the upland areas of West Cork, the upland area between Carlow and Wexford and some suitable areas in the Burren region of Clare. The current habitat suitability of M. reevesi is approximately 231 km² (Figure 2g, Table. 3) with a suitable area mainly confined to the eastern Peninsula of Strangford Lough, Co. Down. The current habitat suitability of C. elaphus is approximately 571 km² (Figure 2j, Table. 3) with suitable areas identified in all counties, although there are two large areas of continuous suitable climate and habitat, these are Killarney, Co. Kerry and Co. Donegal. The current habitat suitability of D. dama is approximately 382 km² with suitable areas in all counties (Figure 2m, Table. 3). However, there are six large geographic areas of continuous suitable climate and habitat, these are Phoenix Park, Co. Dublin, the Mourne mountains, Co. Down, the Tyrone-Monaghan border region, Co. Sligo, the uplands of the Laois/Offaly border and Co. Waterford with bordering areas of southern Tipperary and Kilkenny.

Figure 2. Species richness of flora of conservation concern projected into current and future scenarios (a-c). Distributions of sika deer (d-f), muntjac deer (g-i), Red deer (j-l) and Fallow deer (m-o).

Table 2. Mean accuracy metrics for the nine iterations of the random forest models run per species. AUC=area under the curve. TSS=true skill statistics.

Species	TSS	KAPPA	AUC
Allium schoenoprasum	0.738	.470	0.876
Campanula trachelium	0.777	.759	0.884
Centaurea cyanus	0.713	.253	0.850
Cervus elaphus	0.537	.434	0.853
Cervus nippon	0.758	.646	0.937
Daboecia cantabrica	0.850	.367	0.938
Dama dama	0.753	.731	0.940
Erica mackaiana	0.993	.963	0.999
Frangula alnus	0.970	.647	0.985
Geranium sylvaticum	0.501	.611	0.724
Hordeum secalinum	0.774	.682	0.943
Lotus subbiflorus	0.969	.875	0.987
Muntiacus reevesi	0.911	.912	0.979
Poa alpina	0.926	.944	0.961
Rorippa islandica	0.598	.316	0.791
Trifolium subterraneum	0.999	.911	0.999
Vicia lathyroides	0.888	.734	0.948

Table 3. The current, future RCP45 2050 and future RCP85 2050 range size for the four deer species and 13 flora species.

Name	Current range size (km^2)	Future range size RCP45 2050 (km^2)	Future range size RCP85 2050 (km^2)
Cervus elaphus	571	392	2835
Cervus nippon	1257	6392	14446
Dama dama	1459	1088	543
Muntiacus reevesi	231	0	0
Allium schoenoprasum	382	22213	50842
Campanula trachelium	22	0	0
Centaurea cyanus	13	4	9
Daboecia cantabrica	508	24471	10557
Erica mackaiana	4	0	0
Frangula alnus	16	4	14
Geranium sylvaticum	14	0	0
Hordeum secalinum	354	5226	25992
Lotus subbiflorus	357	84625	84625
Poa alpina	4	0	0
Rorippa islandica	43	11	6
Trifolium subterraneum	5	0	0
Vicia lathyroides	16	0	0

We also projected the distributions of the 13 flora species of conservation concern with medium to high accuracy (Table 2), with AUC scores ranging from 0.724 to 0.999, kappa scores ranging from 0.253 to 0.963 and TSS scores ranging from 0.501 to 0.999, suggesting good predictive models. While kappa had the lowest individual scores, this is not out of line with criticism of its predictive power (Pontius and Millones 2011). When distribution outputs were stacked to represent species richness of species of conservation concern (Figure 2a-c), the areas of high flora species richness did not overlap with the current distribution of three deer species, with D. dama being the only exception as it is the only deer species to have a current area of continuous suitable climate and habitat in Co. Dublin where current flora species richness (of conservation concern) is highest. For flora of conservation concern in general, high levels of species richness mainly occur in and/or around regions with a mixture of urban and agricultural environments such as Dublin, Waterford, Cork and Dundalk (Figure 2a). Other areas that have lower levels of flora species richness occurring outside these regions included the more rural regions of West Galway/Connemara, the Burren, the Shannon Estuary, Killarney, small areas of West Cork, the Three Sister Estuary (Waterford Harbour) and the East coast of Co. Wexford.

The area of high ‘at-risk’ floral species richness in the climate scenario RCP45 2050 is projected to increase (Figure 2b), with large areas of the east reporting a projected higher richness. These are areas where the models project a suitable climate and habitat for multiple endangered/vulnerable plant species, which occur in the south, south-east and mid-west of Ireland (Figure 2e). The climate scenario RCP85 2050 predicts a substantial increase in plant species richness across the island compared to current plant species richness (Figure 2c). Areas with the highest levels of species richness include large parts of the east and south of the island. Also, areas of high flora species richness are not just confined to the urban areas.

The area projected as suitable for C. nippon substantially increases from the current 1,257 km² when projecting into future climate scenarios, with suitable areas of 6,392 km² and 14,446 km² for RCP45 and RCP85, respectively (Table 3). The climate scenario RCP45 2050 reveals that the habitat of C. nippon will increase, expanding its current primary locations, such as the Wicklow mountains, the upland areas of County Kerry, the upland areas of west Cork, the upland areas of Carlow/Wexford border and the Co. Clare region. Also indicating more large areas of continuous suitable climate and habitat in the mid-west, midlands, southern, eastern and north-eastern parts of the island (Figure 2e). The climate scenario RCP85 2050 reveals that the habitat of C. nippon will substantially increase, expanding its current primary locations and with an increase in areas of suitable climate and habitat in all regions of the island with the exception of the northwest (Figure 2f).

For M. reevesi, there are no suitable areas projected for both future climate scenarios (Figure 2h-i). The area projected as suitable for C. elaphus decreases from the current 571km² to 392km² for the future climate scenario RCP45, but has a large increase for the future climate scenario RCP85 with an area of 2835km² (Table 3). The climate scenario RCP85 2050 reveals that the habitat of C. elaphus will be substantially increased, with a large increase in Kerry/West Cork region and a slight decrease in current primary locations (Figure 2k, i). The area projected as suitable for D. dama decreases from the current 1,459km² when projecting into future climate scenarios, with suitable areas of 1,088km² and 543km² for RCP45 and RCP85, respectively (Table 3). The climate scenario RCP45 2050 reveals that the habitat of D. dama will decrease, losing some of its current primary locations, such as Co. Waterford with bordering areas of southern Tipperary and Kilkenny, and Phoenix Park, Co. Dublin. However, other primary locations have a smaller decrease in size (Figure 2n). The climate scenario RCP85 2050 reveals that the habitat of D. dama will be substantially decreased, losing some of its current primary locations and a new suitable area in the Kerry/West Cork region. With the other suitable areas being in small areas of the midlands and North-East of the island (Figure 2o).

Out of the 13 plant species, four plant species have a range increase from their current range for both future climate scenarios. These species are A. schoenoprasum, D. cantabrica, H. secalinum and L. subbiflorus, with all species having a large increase compared to their predicted current range size (Table 3). There are six plant species, C. trachelium, E. mackaiana, G. sylvaticum, P. alpina, T. subterraneum and V. lathyroides which have no suitable geographic locations for both future climate scenarios with a projection of extinction reported by our models. Three plant species have a decrease in both future climate scenarios. These species are C. cyanus, F. alnus and R. islandica (Table 3).

When the distributions of C. nippon and the 13 plant species are considered together, all six of the flora species with a predicted future area of suitability are predicted to have some overlap with the future range of C. nippon (see Table 4, 5); however, the areas of overlap vary in size. For example, F. alnus is the only flora species whose future range size is predicted to have a 100% overlap with C. nippon for both future climate scenarios. For the future scenario RCP45 2050, R. islandica had a predicted overlap of 18.18%, D. cantabrica 11.59%, A. schoenoprasum 9.99%, H. secalinum 9.20% and L. subbiflorus 7.55% (Table 4). For future scenario RCP85 2050, D. cantabrica had an predicted overlap of 32.99%, H. secalinum 25.04%, A. schoenoprasum 21.39%, L. subbiflorus 17.07% and R.islandica 16.67% (Table 5).

Table 4. The future RCP45 2050 overlap of the four Deer species in km² and per cent with the 13 flora species.

Overlap RCP45 2050	Cervus elaphus (km^2)	Cervus elaphus (%)	Cervus nippon (km^2)	Cervus nippon (%)	Dama dama (km^2)	Dama dama (%)
Allium schoenoprasum	19	.09	2219	9.99	33	0.15
Campanula trachelium	0	.00	0	.00	0	0.00
Centaurea cyanus	0	.00	0	.00	0	0.00
Daboecia cantabrica	259	1.06	2835	.00	735	3.00
Erica mackaiana	0	.00	0	.00	0	0.00
Frangula alnus	4	100	4	100	4	100
Geranium sylvaticum	0	.00	0	.00	0	0.00
Hordeum secalinum	0	.00	481	9.20	3	0.06
Lotus subbiflorus	392	.46	6392	7.55	1088	1.29
Poa alpina	0	.00	0	.00	0	0.00
Rorippa islandica	0	.00	2	18.18	0	0.00
Trifolium subterraneum	0	.00	0	.00	0	0.00
Vicia lathyroides	0	.00	0	.00	0	0.00

Table 5. The future RCP85 2050 overlap of the four Deer species in km² and per cent with the 13 flora species.

Overlap RCP85 2050	Cervus elaphus (km^2)	Cervus elaphus (%)	Cervus nippon (km^2)	Cervus nippon (%)	Dama dama (km^2)	Dama dama (%)
Allium schoenoprasum	840	1.65	10875	21.39	145	0.29
Campanula trachelium	0	.00	0	.00	0	0.00
Centaurea cyanus	0	.00	0	.00	0	0.00
Daboecia cantabrica	312	2.96	3483	32.99	286	2.71
Erica mackaiana	0	.00	0	.00	0	0.00
Frangula alnus	13	92.86	14	100	13	92.86
Geranium sylvaticum	0	.00	0	.00	0	0.00
Hordeum secalinum	76	.29	6509	25.04	46	0.18
Lotus subbiflorus	2835	3.35	14446	17.07	543	0.64
Poa alpina	0	.00	0	.00	0	0.00
Rorippa islandica	0	.00	1	16.67	0	0.00
Trifolium subterraneum	0	.00	0	.00	0	0.00
Vicia lathyroides	0	.00	0	.00	0	0.00

When the distributions of C. elaphus and the 13 plant species are considered together, five of the flora species with a predicted future area of suitability are predicted to have some overlap with the future range of C. elaphus (see Table 4, 5); however, the areas of overlap are not extensive. F. alnus is the only flora species whose future range size is predicted to have an extensive overlap with C. elaphus for both future climate scenarios with a 100% predicted overlap in scenario RCP45 2050 and a predicted overlap of 92.86% in scenario RCP85 2050 (Table 4). The future RCP45 2050 distribution for H. secalinum has no predicted overlap with C. elaphus and the lowest predicted overlap of 0.29% for the future RCP85 2050 distribution. All the other flora species are predicted to have an overlap of less than 3.5% with C. elaphus for both future scenarios, largely due to the restricted south-west distribution of C. elaphus under future climate scenarios.

When the distributions of D. dama and the 13 plant species are considered together, four of the flora species with a predicted future area of suitability are predicted to have some overlap with the future range of D. dama (see Table 4, 5); however, the areas of overlap are not extensive, with all the flora species predicted to have an overlap of less 3% for both future climate scenarios (Tables 4, 5). However, there are two species, C. cyanus and R. islandica who are both predicted to have no future overlap with D. dama for both future climate scenarios. Although, the future RCP45 2050 distribution for F. alnus has the highest predicted overlap of 100% with D. dama; the RCP85 2050 distribution also has the highest predicted overlap for all flora species with a predicted overlap of 92.86%.

Discussion

The main aim of this research was to identify the geographic distribution of the four deer species currently found on the island of Ireland (three of which are non-native) and to reveal areas that may be at risk of biodiversity loss, focusing on endangered and vulnerable plant species. In doing so, we aimed to support recent calls in the literature to support the decision-making spectrum of related management and conservation dilemmas on the flora in the non-native ranges of deer species (Carden et al. 2010; Saggiomo, Esattore, and Picone 2020). This research is to our knowledge the first attempt to use a species distribution modelling framework to estimate the current and future distributions of deer and plant species of conservation concern at an all-Ireland level. We projected species distributions to high accuracy levels (Table 1), with large areas of currently suitable habitat for deer species across counties. Future climate scenarios identified an increase in sika deer, while muntjac and fallow deer are projected to reduce their distributions under climate change. Red deer appeared to decline under the moderate climate change scenario but increase under the worst-case projection. However, climate change projections appeared to project the local extinction of six (out of 13) endangered and vulnerable plant species (Table 3). Moreover, while an increase in area was identified for four plant species, several of these areas overlap with the future projections of deer species (Tables 4, 5), which could cause further pressures on these ‘at-risk’ plant species.

Our results are incongruent with the current thinking that all deer species will continue to expand their habitat in Ireland (Carden et al. 2010). While we projected an expansion of sika under both RCP scenarios, which is congruent with predictions (Carden et al. 2010), we identified a decrease in range size for fallow and muntjac deer, and a decrease under RCP45 for red deer, albeit with an increase under RCP85. This does not imply that all deer species will not be present in currently occupied areas at the time of this prediction, it simply forecasts the most suitable habitat for the modelled species in relation to the environmental variables of that period. A central assumption of SDM assumes that a species is at equilibrium with its environment, which is where it occurs in all suitable sites, while being absent from all unsuitable sites (Araújo, Pearson, and Rahbek 2005). However, introduced species often violate this assumption as occurrence records of the species do not always inform its potential distribution (Gallien et al. 2012). This could have been improved by using global distribution records from the native ranges (O’Mahony, de la Torre Cerro, and Holloway 2021); however, given the length of time the deer species have been present on the island of Ireland (e.g. sika were introduced in the 1860s – Purser et al. 2009), these species are most likely in equilibrium with the environment. As such they are no longer expanding their ranges, which may account for the reduction in range size projected by our models. Muntjac are the species that has been present on the island for the shortest period, meaning the low habitat suitability projected for this species may be an artefact of the assumed equilibrium, and as such should be treated with caution.

Moreover, we opted to threshold our habitat suitability projections from a continuous probability to a binary output. This was to provide the framework to test the areas of overlap with ‘at-risk’ plant species. While there is a low probability of occurrence in these regions based on climate and land cover, it is feasible for the species to be found there. This is particularly relevant for invasive species, which are often released and dispersed through less-than-ideal environmental conditions (Barbet-Massin et al. 2018).

Despite this, expansion of invasive deer species in other countries has slowed, with muntjac range expansion dropping to 2–3% in England from 2011 to 2016 compared to an expansion of almost 12% in 2007 (Ward, Richardson, and Mergeay 2021). When this concept is coupled with the lower-than-average forest cover in Ireland compared to other European countries (<1% at the turn of the century – Mitchell 2000), it is likely that deer species in Ireland are inhabiting areas of lower suitability, or at the very least climatically different areas than it may favour given complete forest coverage. We opted to keep land cover consistent in the models, in part due to a lack of reliable alternatives, which is in line with the SDM literature (Stanton et al. 2012). However, with active afforestation policies in Ireland (i.e. Forest Action Plan 2022) that are aiming to increase the area of forests as well as the native components of these woodlands, it is likely that our models are an underestimation of future geographic distributions of all deer species. Despite this, our models provide a baseline of the minimum potential impact to ‘at-risk’ plant species, and as such will support and inform future management strategies.

We projected a local extinction of six ‘at-risk’ plant species under future climate scenarios. The local extinction of any species could have cascading effects on the ecosystem, but our models have projected that almost 50% of the plant species considered here may go extinct (Figure 2). Given that the 13 species included in this analysis were the only ones to reach a minimum number of observations in GBIF, the reality is that there may be a much higher rate of extinction for flora species of conservation concern, as these are in theory the most common of the endangered and vulnerable plant species. As such, there is a pressing concern for research to be directed towards understanding how these species will respond to climate change. It should be noted here that the drivers of vascular plant extinctions are nuanced, with numerous pressures globally and within Ireland (Nic Lughadha et al. 2020; O’Rourke, Lughadha, and Bacon 2022). Therefore, the potential for extinction based on these climate models must be considered alongside many other complex and interlinked ecosystem functions and pressures, such as habitat loss, exploitation, intensification of agriculture and non-native herbivory.

Contrary to the extinction risk, we did identify that some ‘at-risk’ species are projected to increase their geographically suitable area. ‘At-risk’ plant species richness appeared to increase nationally, with some of the highest projected richness occurring in urban and agricultural areas (Figure 2e, f). We posit here the importance of the landscape matrix (i.e. multiple habitats) and urban wildflower meadows that have become increasingly popular in city centres in recent years to allow citizens to observe and enjoy the pleasures of nature (Bretzel et al. 2016) and more recently support the biodiversity crisis (Mody et al. 2020). These green spaces help to mitigate the negative impacts of urbanization on biodiversity such as habitat loss, fragmentation, decreased population size and loss of genetic diversity across a multitude of species (Lambert, Cawkwell, and Holloway 2021; Vega and Küffer 2021). Moreover, Ireland has been reported as a location whereby species have expanded their geographic ranges in recent years in response to climate change (Bellard et al. 2018; Donnelly et al. 2014), and while such research warrants further investigation, our findings suggest that certain plant species that are currently restricted locally on the island may benefit from changes to the existing climate.

For example, the Biodiversity Audit of the Áras an Uachtaraín in Dublin in 2019–2020 identified 297 species of plant in the 130-acre site (Gaughran and Stout 2020). While significantly larger than the average urban greenspace in Ireland, the importance of these refuges to native species is clear, with the dispersal of plants from the plantation sites to support urban range expansion well documented (Gilbert and Anderson 1998). Moreover, given the avoidance of urban centres by deer species in Ireland, urban areas could provide a refuge for many plant species, including endangered and vulnerable species to not only persist, but thrive under future climates. Almost 50% of the plant species included in this analysis are available to purchase from horticulturalists, suggesting this urban refuge may support these species.

Conversely, it is also important to consider the areas predicted to have a low richness of ‘at-risk’ species. Linders et al. (2019) state that one of the most effective mediation strategies for invasive species is to support vegetation abundance, species richness and diversity, so that any grazing pressure is ‘shared’ among species. It can be seen from Figure 2b-c that areas of lower species richness under the two RCP scenarios are found in the south-west, including Co. Kerry and Cork, the north-west and large areas of Northern Ireland. When these outputs are compared with the future outputs of C. nippon and C. elaphus distributions (Figures 2e-f and k-l), the area in the south-west of Ireland can be identified as an area that may need to be explicitly monitored due to the presence of one or two ‘at-risk’ plant species and a large geographic distribution of deer. Therefore, any management strategy should look to ensure that there is a diverse range of vegetation available, such that certain species are not subject to increased and unsustainable grazing. This is particularly pertinent given the multifaceted and nuanced implications of forest management, with Hawthorne, Colombaroli, and Mitchell (2021) suggesting that the complete exclusion of grazing by deer in the region would be undesirable due to the potential encroachment of woody species leading to the decline of woodland specialists.

There was one plant species predicted to have a 100% overlap with three deer species under climate change. F.alnus is predicted to occur in the west of Ireland in and around current suitable areas in the Burren, Killarney and West Cork Region. This plant is mainly found in damp, peaty soil and the wooded areas on limestone soils (Goodwin 1943) some of which are land covers deer species do not frequent. Therefore, it may be possible that despite high overlap with deer species, F.alnus may occur in land cover less suitable; however, they should be monitored regularly for expansion of their ranges.

Based on the discrepancy of the findings, it is difficult to ascertain the exact role of climate or non-native pressure on the varying fortunes of the ‘at-risk’ plant species. While our models suggest extinction for 50% of the plant species studied here based on climate, which is congruent with other global studies (Román-Palacios and Wiens 2020), it is clear that the presence and overlap of non-native species could further propagate such considerations. Moreover, the importance of restoring functional ecosystems or trophic rewilding has become increasingly important within the conservation debate (Svenning et al. 2016), with studies suggesting that restoring native animal populations to their functional roles can enhance natural carbon capture and storage, and potentially negate climate warming (Schmitz et al. 2023). Extinct fauna in Ireland includes the grey wolf, the great elk and the brown bear, to name just a few (Hickey 2013), meaning that the native and non-native deer within the (largely restricted) forest ecoystems are not the top-herbivore and are currently uncontrolled by any natural predator, which further propogates their risk to endangered flora. Therefore, a diversification of plant flora in a region, coupled with a ‘rewilding’ of the ecosystems is currently considered the most effective management practice to support the impact of non-native deer species on native plants.

We only considered abiotic variables in model parameterization, including climate and habitat. However, deer are herbaceous species whose distribution is dependent on food sources. We did not explicitly include biotic factors, such as the distribution of food sources or competition, in the model as has been suggested (Peterson et al. 2020; Soberon and Peterson 2005). In part this was due to the fact that we wanted to compare the distributions of endangered and vulnerable plant species with those of the deer species, and any inclusion of vegetation in parameterization of these distributions would have introduced circularity into the model. However, given the importance of climatic variables in determining deer distributions (results not shown), it is possible that such considerations were incorporated in the model as latent variables (Schymanski et al. 2013). While this is important to note, so too is the fact that perfect statistical independence between environmental variables is perhaps unrealistic of real-world variables, with propagation of environmental variables most likely present in all SDM research (Holloway and Miller 2017). It is infeasible to omit all variables based on the assumption that they may already be included elsewhere in the model, and when coupled with the high accuracy of the models (Table 1), our overarching aim to accurately predict future distributions of these species has been achieved.

Further, we also assumed unlimited dispersal capacity of all species, which is in line with species who have high dispersal capacities or the potential for translocation by humans (Holloway, Miller, and Gillings 2016). We also note the importance of the low number of occurrences. While this has been found to impact distribution projections, research has found that models can be reliably and accurately predicted using as little as three observations (van Proosdij et al. 2016). To ensure a minimum standard of reliability, we opted to use a minimum of ten observations for the plant species. However, of the 128 species listed on the Red Data Plant List, only 13 had a minimum of ten sightings on GBIF. Future research collecting verified plant species locations that are made available through open data repositories is needed to improve model calibration and support biodiversity conservation and supports recent calls for more comprehensive surveys and biodiversity audits in Ireland (Gaughran and Stout 2020).

Due to these limitations in data availability, some of our species reported high uncertainty in the spatial positioning (see Supplementary Material). As stated, recent research has identified that data quantity is more important than spatial bias and uncertainty in the predictive performance of SDMs (Gaul et al. 2020), meaning we opted to use these data points to ensure we could undertake this research. Despite these caveats, we do report high accuracy values (Table 1), allowing confidence in the discussion of the important environmental drivers of the plant species distributions. Moreover, of the three species that reported high uncertainty, only P. islandica has a lower than average accuracy metric, meaning the uncertainty has not appeared to overly impact model predictions. However, both P. aplina and V. lathyroides are two of the species projected to go extinct, but represent two of the smallest current distributions of plants. More research and data are therefore needed to investigate whether these species are projected to become extinct, or whether this is an artefact of data availability and their cryptic nature. Given that the 13 species included in this analysis were the only ones to reach a minimum number of observations, the reality is that there may be a much higher rate of extinction risk for flora species of conservation concern.

Finally, there persists a need to explore habitat suitability of deer and plant species in conjunction with other measures of biodiversity, such as density and abundance. New methods, such as Poisson point process models, are emerging as suitable tools to predict the distribution and abundance of species from presence-only data (Renner and Warton 2013; Schank et al. 2017; Warton and Shepherd 2010). Similary, new data fusion techniques are being used to combine presence-only data and presence-absence across multiple spatial scales that are resolving some of the issues related to data sparsity (Morera-Pujol et al. 2022). Given the uncertain relationships identified, future studies should focus on developing predictions of both the distribution and density of deer to inform forest management and biodiversity conservation.

Conclusion

It has been identified that over 1 million species are threatened with accelerated extinction (Díaz et al. 2019), with the spread of invasive species ranked as one of the primary drivers of this loss. Here we utilized a species distribution modelling framework, which is a prominent GIScience technique, to assess the current and future habitat suitability of four deer species (three non-native and one native) and 13 endangered/vulnerable plant species. To our knowledge, this is the first time this applied research has been undertaken at an all-Ireland scale. We then used the predicted future distributions to identify how much the deer and plant species distribution will overlap. We predicted large areas of current suitable habitat for the deer species across all counties. Future climate scenarios identified an increase in sika deer, while fallow deer are projected to reduce their distribution under climate change and our results predict a local extinction for muntjac. Also, our results predict the local extinction of six (out of 13) endangered and vulnerable plant species. Moreover, while an increase in area was identified for four plant species, several of these areas overlap with the future projections of deer species, which could cause further pressures on these ‘at-risk’ plant species. However, ‘at-risk’ plant species richness appeared to increase nationally, with some of the highest projected richness occurring in urban and agricultural areas. Finally, our results should support the conservation of these endangered/vulnerable plant species from the negative impact of native/non-native deer species and can be used to mediate the impact on the ‘at-risk’ plant species through different conservation and management measures. Not only has this research provided practitioners with a product with which to identify areas to prioritize conservation activities, saving their time, resources and research efforts, it has also resulted in a better understanding of the effects of climate change on projected distributions and extinction rates.

Acknowledgement

We would like to thank the editors, reviewers, and Astrid Wingler for their comments and suggestions.

Disclosure statement

No potential conflict of interest was reported by the author(s).