Niche overlap of mountain hare subspecies and the vulnerability of their ranges to invasion by the European hare; the (bad) luck of the Irish

Niche conservatism is the tendency of related species to retain ancestral tolerances after geographic separation. We used Ecological Niche Modelling and Principal Components Analysis of bioclimatic and habitat variables to describe the extent of the species niche, and degrees of bioclimatic–habitat niche conservatism within the mountain hare (L. timidus) clade. Mountain hare niche space was contrasted with that of the European hare (L. europaeus), to shed light on species interactions in contact zones throughout Europe. All five subspecies of mountain hare had quantifiably distinct niches. Fennoscandian (L.t. sylvaticus, L.t. timidus) and highland (L.t. scoticus, L.t. varronis) subspecies, however, were most similar, exhibiting greatest apparent niche conservatism. They inhabit tundra, boreal forest and uplands, and, hence are presumed most similar to the ancestral form. The Irish hare was distinct, being consistently distinguished from other mountain hares in both 2D and nth dimensional (4D) niche space. The ecological distinctiveness of the Irish hare provides further evidence that it is an Evolutionarily Significant Unit, particularly vulnerable to displacement by introduced European hares with which it competes and hybridises. Projections under global climate change suggest that, by 2070, bioclimatic space for invasive European hares in Ireland will expand (by 79%) but contract for endemic Irish hares (by 75%), further facilitating their replacement. The near complete species replacement of the heath hare (L.t. sylvaticus) in southern Sweden, where the European hare has also been introduced, may suggest a similar fate may be in store for the Irish hare.


Introduction
The controversial concept, 'niche conservatism', is the tendency of emergent species to retain their ancestral ecological traits such that closely related species may be more ecologically similar than would be expected based on their phylogenetic divergence (Wiens et al. 2010). A niche comprises a multivariate set of abiotic and biotic conditions which facilitate the persistence of a species, and to which it is suitably adapted (Hutchinson 1957). However, studies investigating the relationships between a species' distribution and niche frequently fail to appropriately define their terms (Soberòn 2007). The fundamental niche is unconstrained by limiting biotic factors such as ecological competition, predation, dispersal ability, and environmental conditions (Hutchinson 1957;Wiens and Graham 2005). The realised niche is described as the fundamental niche constrained by limiting factors, i.e. the space occupied by, and the resources available to, an organism (Hutchinson 1957;Soberòn 2007). There are, however, issues inherent in the utilisation of these terms due to considerations of spatial resolution and biotic interactions (Araújo and Guisan 2006). Population viability depends on a degree of environmental stability or adaptive predictability, and hence, climate and habitat are key factors (Soberón and Peterson 2005;Jäkäläniemi 2011). In the absence of gene flow, populations may diverge genetically while occupying similar habitat to ancestral species, and hence, the species niche is conserved (Peterson et al. 1999;Wiens 2004). Niche adaptation may, therefore, be spatiotemporally stable and can be conserved during allopatric speciation via geographic isolation.
Environmental Niche Modelling (ENM) is increasingly used to estimate environmental suitability as a function of geospatial species occurrence relative to environmental variables (Phillips et al. 2006), thus capturing a species' niche space. Here, we investigate species occurrences and overlap in potentia, with dispersal constrained by climatic and habitat variables (hereafter, referred to simply as the species 'niche'). We aimed to describe the degree of niche conservatism within the mountain hare clade in Europe, capturing the niche of each subspecies. These are contrasted against the niche of the European hare, thus shedding light on observed differences in species interactions post-contact. We pay particular attention to the Irish hare due to its ecological equivalency to the European hare (i.e. temperate, lowland, grazing habits). We hypothesise that, given prolonged, postglacial isolation and ecological expansion: (1) the niche of the Irish hare is more ecologically distinct from other mountain hare subspecies than those subspecies are from one another and, (2) the niches of the Irish and European hare are more similar to each other than are those of other mountain hare subspecies, relative to the European hare (making the Irish hare more vulnerable to the impact of European hare invasion). We also model the predicted shift in the bioclimatic space suitable for Irish and European hares in Ireland, under projected global climate change. We hypothesised that (3) the bioclimatic space available for the Irish hare is likely to become increasingly unsuitable, and, (4) the bioclimatic space available for the European hare is likely to become increasingly suitable under warming temperatures due to their contrasting origins and differentially-adapted physiology (the former having Arctic, and the latter Middle Eastern, ancestry). Thus, we expect that the trajectory of global climate change is likely to be beneficial to the invader and detrimental to the native species.
In the Order Lagomorpha, the genus Lepus (hares and jackrabbits) is represented in Europe by five extant species, two of which, the mountain hare (Lepus timidus, Linnaeus, 1758) and the European hare (L. europaeus, Pallas 1837), are widely distributed. These species are readily distinguished phenotypically, by ear and limb length (all shorter in mountain hares, apart from the hind feet), head shape (convex in European hares), strongly contrasting black ear-tips (present in European hares), white-striped muzzle (present in European hares), ventral tail surface (black in European hares), body mass (lower in mountain hares) and pelage colour (darker and more uniform in mountain hares; Flux and Angerman 1990;Caravaggi et al. 2016). The European hare is a highly successful invasive species that has been introduced to a large number of countries worldwide (sensu Flux and Angerman, 1990). It is typically parapatric with the mountain hare, being separated by elevation or habitat, with narrow contact zones suggesting each species has a distinct niche separated by, for example, differences in habitat or climate (Amori et al. 2008). The invasion dynamics of the European hare and its interaction with native mountain hare populations are poorly understood (Thulin 2003;Reid 2011). Some contact zones between the species are largely stable (e.g. in the Alps and Scottish Highlands), though it is predicted that such elevationally-defined contact zones will shift upwards due to the effects of global climate change (Leach et al. 2015a). Other, more recently established contact zones are highly unstable with the larger European hare outcompeting and displacing the smaller mountain hare. Indeed, European hares have displaced mountain hares over much of southern Sweden (Jansson and Pehrson 2007) and part of southern Finland (Levänen et al. 2015) during the twentith century, and part of Ireland in the last few decades Reid 2011;Caravaggi et al. 2015Caravaggi et al. , 2016. Thus, post-introduction sympatry is a typically transient phenomenon (Thulin 2003).
The mountain hare is a circumpolar, arcto-alpine species complex, distributed from Ireland in the west, to Japan and Kamchatka in the east, and from the Alps in the south, to 75°N (Flux and Angerman 1990;Smith and Johnson 2008b). There are five extant European mountain hare subspecies differentiated by morphophysiological characteristics, behaviour, and ecology (Angerbjörn and Flux 1995). There is generally low genetic differentiation between subspecies, indicative of a post-glacial panmictic European population, which subsequently underwent fragmentation, isolation, and divergence (Hamill et al. 2006). The most widespread subspecies, recognised as the typical form (and thus presumed similar to the ancestral type), is the northern hare (L. timidus timidus, Linnaeus 1758) which inhabits tundra (in the north) and boreal forest (further south) in the Arctic and Fennoscandia (Angerbjörn and Flux 1995). Its diet varies seasonally, with hard woody material being consumed in winter, and grasses, sedges, and herbs, in late summer and autumn (Flux and Angerman 1990;Helle 1995). The heath hare (L.t. sylvaticus, Nilsson 1831) occurs in southern Sweden and Gotland (Winiger 2014). The subspecific status of this taxon is debated, with some regarding it as a synonym of L.t. timidus. However, many others recognise it as a distinct subspecies based on winter pelage, which is blue-grey rather than white (Lindström 1980;Suchentrunk et al. 1999;Thulin et al. 2003;Winiger 2014). Both the northern and heath hares are thus Fennoscandian mountain hare subspecies and geographically distinct from three isolated mountain hare populations, two of which are true highland mountain hares: the Scottish hare (L.t. scoticus, Hizheimer 1906) and the Alpine hare (L.t. varronis, Miller 1901). The former is widespread throughout montane habitats in Scotland, occurring up to 1300 m asl (Newey et al. 2011). The latter is generally found on forested slopes (Bisi et al. 2013;Rehnus et al. 2013) up to 3500 m asl (Thulin 2003;Rehnus et al. 2013), throughout the Alps (Angerbjörn and Flux 1995). The highland subspecies browse hard, woody plant material e.g. heather Calluna vulgaris (Flux and Angerman 1990). The Irish hare (L.t. hibernicus, Bell 1837) is endemic to the island of Ireland, where it has been isolated for 30,000-60,000 years (Hughes et al. 2006). One estimate placed the divergence of Irish hares from other mountain hares (specifically, Russian L.t. timidus) at ca. 360,000 years before present (Hughes et al. 2006). This subspecies possesses a comparatively high number of unique genetic forms (mitochondrial haplotypes) not shared by any other subspecies outside Ireland (Hughes et al. 2006). It exhibits considerable ecological plasticity, being found at all altitudes in Ireland, but is most common in the lowlands (Whelan 1985;. In contrast to other mountain hares it feeds predominantly on soft, mostly agricultural grasses, e.g. ryegrass Lolium perenne (Strevens and Rochford 2004). Nearly all mountain hare populations exhibit winter whitening as camouflage during winter snow cover (e.g. Hewson 1958), with one exception. The Irish hare has largely lost the trait, save for minimal whitening of the ear margins and feet (Flux and Angerman 1990). Such is the genetic, phenotypic, behavioural and ecological distinctiveness Niche overlap of mountain hare subspecies and the vulnerability of the Irish hare, that some contend it may warrant full species status (Hughes et al. 2006). It is as divergent from other mountain hare subspecies as the mountain hare is from other species such as the Arctic (L. arcticus, Ross 1819) or Alaskan (L. othus, Merriam 1900) hares (Paulo Prodöhl pers. comm.), whose taxonomic status and phylogenetic relationships with the mountain hare have been the subject of debate (e.g. Wu et al. 2005;MacDonald and Cook 2010).

Data sources and preparation
A total of 238,813 records of mountain hare subspecies and European hare found in Europe were obtained from a large number of sources, principally biodiversity data record centres, academics and ecologists (Tables S1, S2 in Supporting Information). Data were collected via a variety of methods, combinations of which differed between and within regions, countries and organisations, e.g. scientific surveys, hunting bags, opportunistic sightings by the public, ecological surveys, road casualties. Hereafter, we adopt the terms ''(sub-)species'' to refer to the mountain hare (including all subspecies) and the European hare, or ''subspecies'' when referring to the mountain hare only. Records were extracted during 2013-2014 and were sub-sampled by date (post-1950, to ensure consistency with current bioclimatic datasets), and geospatial accuracy (B1 km resolution). Furthermore, while there may be difficulties inherent in discriminating between sympatric species, we were unable to quantify observer bias. Duplicate records were removed, as were those considered erroneous based on known distributions of each (sub-)species (i.e. falling beyond the boundary of the International Union for Conservation of Nature range polygon; Smith and Johnston 2008a, b). Species-specific records that occurred within the known range of that species were, therefore, considered 'true' and retained, while those that occurred outside the known range were considered 'false' and removed. The range polygons for each mountain hare subspecies were extracted from the parent IUCN range polygon and sub-divided into geographically isolated populations i.e. Ireland, Scotland, and the Alps, whilst the Fennoscandian mountain hare subspecies ranges were delineated according to Bergengren (1969). Due to a lack of sufficiently precise data in northern Fennoscandia and much of central Europe, Northern and European hares appeared erroneously 'absent' from parts of their known range. Furthermore, data exhibited considerable sample bias (Yackulic et al. 2013), with large numbers of records occurring around urban centres, particularly in the UK and Sweden. There are a number of methods available for accounting for sample bias (see Fourcade et al. 2014), including the utilisation of target background points (Phillips et al. 2009) or bias grids (Elith et al. 2010). Targetbackgrounds are defined as background points drawn from occurrences of a focal class (e.g. lagomorphs, herbivorous mammals). Thus, background data will exhibit similar spatial bias to that of the modelled species (Phillips et al. 2009). Similarly, a bias grid is a surface scaled to represent survey effort (Elith et al. 2010), a quantity unknown for almost all ([99.9%) of our data. However, data manipulation (i.e. removing data in over-sampled regions) may be effective in reducing or removing bias (Phillips et al. 2009). Thus, in order to reduce sample selection bias, presence records were thinned using OccurrenceThinner version 1.04 downloaded from www.phycoweb.net/ software. OccurrenceThinner uses probability algorithms to remove occurrence records based on an associated kernel density grid. The probability that an occurrence will be removed is proportional to occurrence density described by the kernel density grid (Verbruggen et al. 2013). Due to the extremely high density of occurrences in some regions (e.g. urban areas in the UK and southern Sweden), data were sequentially thinned to appropriate densities which were informed a priori by densities of records elsewhere in the species range. A priori thinning aimed to equalise the densities of occurrence records on a landscape scale, and, hence, produce ecologically relevant models. A total of 9075 records were used in modelling (see Table S2 for species specific pre-and post-thinning occurrence counts, Fig S1 for occurrence distribution maps). 10,000 background data points (i.e. pseudo-absences) were generated randomly within the range of each individual (sub-)species, analogous to the Restricted Background approach detailed in Fourcade et al. (2014).
Climate data were downloaded from WorldClim (www.worldclim.org) at 30 arc-second (ca. 1 km 2 ) resolution. Species records were associated with mean data from 1950 to 2000 for current models only. Three raw-format (mean temperature, precipitation seasonality and temperature seasonality) and three composite (Hilliness Index, Normalised Difference Vegetation Index (NDVI), and water balance) environmental variables were used (Table S3). Eight land cover variables (coniferous forest, crops, mixed forest, moorland and heathland, pasture, peat bog, scrub and sparse vegetation; see Table S3 for vector filenames) were obtained from the CORINE Land Cover 2006 (EEA 2010).
Shapefile and raster creation and manipulation were carried out using ArcGIS 10.2.2 (ESRI 2011).

Environmental Niche Modelling
MAXENT is a popular presence-only modelling tool (Phillips et al. 2006, which uses a maximum entropy approach, i.e. the probability distribution which best represents the data is the one with the largest entropy. Despite its widespread use, MAXENT has been criticised due to its vulnerability to overfitting and the use of logistic output to estimate absolute occurrence probabilities (e.g. Royle et al. 2012). Such limitations may be mitigated against by careful a priori data manipulation, to closer approximate the assumptions of the model, e.g. that occurrence data represent unbiased independent samples, constant probability of detection, and that detectability is independent of model variables (Yackulic et al. 2013). Indeed, MAXENT has been shown to consistently outperform other comparable modelling techniques (e.g. Wisz et al. 2008;Tarkesh and Jetschke 2012). While MAXENT is relatively robust against collinear variables (Rodriguez-Robles et al. 2010;Kuemmerle et al. 2012), several climatic variables exhibited strong collinearity; exploratory models suggested a strong cumulative influence. Variables with the greatest permutation importance, i.e. mean temperature (collinear with minimum temperature, maximum temperature) and annual water balance (collinear with minimum precipitation, maximum precipitation, mean precipitation), were retained. Climatic and environmental variables with a mean permutation importance of \ 2 (complex cultivation, human influence index, inland marsh, natural grassland, radiation, snow, urban, number of months with positive water balance) were also removed. ENMs were run using linear, quadratic, product and threshold features with clamping and extrapolation disabled, for 50 replicates. Presence records were split randomly into a 75% training set and a 25% test set, with cross-validation.
Models for the Irish and European hare were projected under global climate change at time-slices for the current period (2010-2014), 2050s and 2070s. IPCC Fifth Assessment Report Coupled Model Intercomparison Project Phase 5 (CMIP5) future climatic data for the Representative Concentration Pathway (RCP) 8.5 for 2050 (averaged across 2041-2060) and 2070 (average for 2061-2080) were downloaded from WorldClim at 1 km 2 grid cell resolution. RCP 8.5 indicates a mean average global temperature increase of 2°C by the 2050s and 3.7°C by the 2070s. All variables were averaged across five Global Circulation Models (GCMs), CNRM-CM5, GFDL-CM3, GISS-E2-R, Had-GEM-ES and MIROC-ESM-CHEM, thus reducing model error (sensu Pierce et al. 2009). Originally described as ''extreme climate change'', this climate scenario now appears to best fit observed climatological trends (following Leach et al. 2015a; Table S3). Changes in predicted range extent were calculated using Max SSS, i.e. the sum of test specificity plus sensitivity, which is effective when using presence only data, and is not affected by pseudo-absences (Liu et al. 2013). A major caveat of this approach is that CORINE habitat variables were kept constant when projecting into future time-slices as no robust predictions are available for how land cover will respond under future climatic conditions. However, this approach is consistent with most studies that project species ranges into future conditions (e.g. Acevedo et al. 2012).

Model evaluation
Models were evaluated using the Area Under the Curve (AUC; Fielding and Bell 1997) of the Receiver Operating Characteristic (ROC) curve, a model-accuracy assessment measure that is independent of prevalence (McPherson et al. 2004). The classification of AUC values follows a commonly-used, yet arbitrary ranking system based on suggestions by Swets (1988), Greiner et al. (2000). Values between 0.9 and 1.0 are considered excellent, 0.9-0.8 good, 0.7 and 0.8 average and \0.7 poor. However, where ROC curves are constructed from presence-only data, the maximum possible AUC is\1 (Wiley et al. 2003), and it is not possible to determine optimal performance (Phillips et al. 2006). Nevertheless, relative performance may still be inferred, given that an AUC of 0.5 describes random prediction (Phillips et al. 2006). We also tested the omission rate (proportion of true occurrences misidentified), sensitivity (proportion of presences which are correctly predicted), specificity (proportion of absences which are correctly predicted), proportion correct (proportion of the presence and absence records correctly identified), and True Skill Statistic (TSS), calculated using SDMTools package (Van der Wal et al. 2012) in R (version 3.2.2). TSS is a prevalence-independent metric derived from threshold sensitivity and specificity. Values range from -1 to ?1 and test the agreement between the expected and observed distribution, and whether the outcome could be predicted due to chance (Allouche et al. 2006

Niche overlap and equivalency
The similarity of (sub-)species continuous-surface probability models (i.e. geographic niche overlap) were evaluated using the niche overlap metric, I (Warren et al. 2008). This method calculates pairwise overlap between models, producing values between 0 (no overlap between niche models) and 1 (identical niche models). Warren's I is based on the probability (p x,i , p y,i ) of a species (X or Y) occurring in a given cell (i); defined by the ENM. In contrast to Schoener's D (Schoener 1968), another commonly-used metric, Warren's I treats p x and p y as probability distributions with no biological assumptions, and, hence, is more appropriate for presence-only analyses (Warren et al. 2008). Niche overlap metrics were calculated for contemporary and future climate-projected models using the R package fuzzysim (Barbosa 2015).
Niche equivalency tests were used to assess whether paired-species ENM overlap values (I) were significantly different from a one-tailed normalized null distribution of comparative overlap values. Null distributions were generated by comparing ENMs of two focal species to random subsets drawn from pooled presences, where the number of extracted (i.e. 'null') presences were equal to the number of observed presences for each species. This was repeated 100 times for each species pair (Warren et al. 2008). Ecological niches were said to be non-equivalent if paired-species overlap values were significantly lower than those of the null distribution (P B 0.05). Niche equivalency tests were carried out using ENMTools (Warren et al. 2010) and using only contemporary data.

Ecological distance
Principal Component Analysis of occurrence records and associated data was used to reduce bioclimatic and habitat variables associated with all species records to four hypothetical axes with eigenvalues [1, describing ecological niche space, using core R functions. A multifactorial General Linear Model (GLM) was used to establish differences in Principal Components (PC1 through PC4) between each (sub-)species with Bonferroni pairwise post hoc test for multiple comparisons used to identify niche space differences. Biplots of paired Principal Component Axes were used to plot the proximity of each (sub-)species in 2D space. For each pairwise plot the mean Mahalanobis distance (De Maesschalck et al. 2000) was calculated between: (1) all pairwise comparisons of mountain hare subspecies excluding the focal subspecies (i.e. the Irish hare); and (2) all pairwise combinations including the focal subspecies. Mahalanobis distances were calculated using the R package StatMatch (D'Orazio 2015). The n-dimensional Euclidean distance between each pair of (sub-)species was also calculated across all four Principal Components simultaneously, thus deriving a single measure of distance between (sub-)species in multidimensional (4D) niche space. Euclidean distances were calculated using the R package pdist (Wong 2013). A t test was used to test for significance of differences between the two groups (mountain hares including and excluding the Irish hare).
The predicted probabilities of mountain hare (sub-)species presence closely approximated the actual range extent of each (sub-)species (Fig. 1). Niche space for the European hare was predicted northward beyond its northern (invasive) range edge in Sweden extending west into southern Norway, southward beyond its southerly (natural) range edge in northeastern Iberia and in all directions around its current invasive range in Northern Ireland (Fig. 1f).

Ecological (dis)similarities
Geographic niche overlap measures derived from continuous-surface probability models described potential overlap between six (sub-)species pairs (I C 0.4; Table 3). Almost all pairwise comparisons between hare (sub-)species and ENMs generated using randomly selected background points did not differ from null distributions, and, hence, their niches can be said to be similar. Only four pairwise comparisons between were found to be significantly different (i.e. less similar than expected by chance; P B 0.05), though the relationship was unidirectional rather than reciprocal: the Alpine hare was distinct from the Scottish hare and the European hare; the Irish hare was distinct from the European hare; and the Northern hare was distinct from the Irish hare (Table 3). Thus, the ecological niche of the Alpine hare, for example, was more distinct from that of the European hare than would have been expected by chance, but not vice versa. Our results suggest that while the ecological niches of hare (sub-)species in Europe are similar, they are not identical.
All Principal Component values varied significantly between (sub-)species (Table S4). The biplot of PC1 and PC2 (accounting for 41% of cumulative variation) suggested that the niches of Fennoscandian mountain hare subspecies were more similar to one another than they were to any other hare (sub-)species (Fig. 2). Both highland mountain hare subspecies were also more similar to one another than they were to any other hare (sub-)species. The Scottish hare occupied a similar, yet slightly more productive environment, suggested by a more positive value on Niche overlap of mountain hare subspecies and the vulnerability PC1, indicating higher NDVI. Variation in Irish hare niche space not only did not overlap with any other mountain hare subspecies, but its centroid was further away from other mountain hare subspecies than it was from the European hare, which was associated with agricultural crops. The Irish hare was associated with temperate, highly productive pastures (Fig. 2). Other pairwise comparisons between remaining Principal Components showed fewer distinct differences (Fig. S2), as they accounted for less variation, Niche overlap of mountain hare subspecies and the vulnerability yet almost all pairwise post hoc tests between (sub-)species were significant (Fig. S3). The Irish hare was most similar to the heath hare and European hare on PC1, the European hare on PC2 (Fig. 2) and the Scottish hare on PC3 and PC4 (Fig. S2). The European hare did not differ significantly from the Alpine hare on PC1, the heath hare on PC2 (Fig. 2) and the northern hare on PC4 (Fig. S2).

Ecological distance
All pairwise mountain hare subspecies comparisons including the Irish hare had significantly longer Mahalanobis distances on each Principal Component biplot, than all pairwise comparisons excluding the Irish hare, except between PC2 and PC3 (Fig. 3c). The Irish hare was on average 13% further away from other mountain hares than they were to one another across all four Principal Component Axes. Pairwise Mahalanobis distances were greatest on those axes explaining the greatest percentage variation in bioclimatic and habitat variables. Mahalanobis distances on Principal Component biplots did not differ significantly for other mountain hare subspecies (Fig. 3b, d, e) with the exception of the Alpine hare, which had significantly shorter distances than all pairwise combinations excluding itself on PC1:PC2, and significantly longer distances on PC2:PC3 and PC3:PC4 (Fig. 3a). The Alpine hare was on average 7% further away from other mountain hares than they were to one another across all four Principal Component axes. Consequently, when using a single metric for the mean nthdimensional (4D) Euclidean distance between species pairs, those mountain hare pairs involving the Irish hare were significantly further away from other mountain hares than pairs excluding the Irish hare were from each other (t df=5 = 3.66, p = 0.015; Fig. 4).

Predicted range shifts
Geographic niche overlap metrics derived from continuous-surface probability models suggest that niche overlap between the European hare and several mountain hare subspecies will increase in coming decades (Table 5). Most notably, overlap with the heath hare is predicted to be almost complete by 2070 (I = 0.91). Geographic niche overlap between the European hare with the Scottish (I = 0.81) and Northern (I = 0.73) mountain hares increases by 2050 but is largely stable thereafter ( Table 5).
Projections of predicted probability of occurrence for Irish and European hares in Ireland under future climate change suggest major changes in the likely distribution of their respective envelopes over the next half century. ENMs for the Irish hare predicted the current range to cover the whole of Ireland at 83,497 km 2 (Fig. 5a). By 2050, the suitable bioclimatic envelope was predicted to contract westward with the most suitable habitat remaining in the northwest, the total range extent declining to 35,461 km 2 (Fig. 5b), and declining further to 21,107 km 2 by 2070 (Fig. 5c). Thus, ENMs predicted a 75% contraction in the species suitable bioclimatic envelope over the next 50 years. Conversely, the current bioclimatic envelope of the European hare was predicted to be restricted mostly to Northern Ireland at 12,417 km 2 (Fig. 5d). By 2050, it expanded south and westward to 53,874 km 2 (Fig. 5e), expanding further to 66,312 km 2 or 79% of the island by 2070 (Fig. 5f).

Discussion
All five European subspecies of mountain hare were found to have quantifiably distinct niches. Fennoscandian subspecies (the northern and heath hares) were more similar to each other than any other subspecies while the highland subspecies (the Scottish and Alpine hares) were also more similar to each other than any other subspecies. The Irish hare was not only distinct, but had zero overlap with other subspecies and was consistently distinguished from other mountain hares in both 2D and nth dimensional (4D) space. Moreover, the niche space of the Irish hare was more similar to that of the European hare than any other mountain Niche overlap of mountain hare subspecies and the vulnerability hare. Thus, all subspecies had separate niches, but the Irish hare stands out as being uniquely different, with little or no commonality to other mountain hares. During the Last Glacial Maximum (LGM), ancestral mountain hares likely maintained a panmictic European population which inhabited highlands, boreal forest, and tundra, in snowy, cold conditions (Angerbjörn and Flux 1995). Northern, Scottish and Alpine hares appear, therefore, to exhibit niche conservatism, retaining their ancestral ecological traits and environmental distributions, such that they are ecologically similar despite their geographic isolation. The extent of their distributions may have been constrained by interspecific interactions with the European hare that invaded lowland areas after the LGM. Moreover, agricultural intensification and changes in land management are likely to have further reinforced elevational and latitudinal separation.
Contact zones between European hares and mountain hare subspecies in the Scottish highlands and the Alps are largely stable. This stability can be largely attributed to differences in habitat and dietary preferences. Scottish hares are most abundant on moorlands, where up to 90% of their diet can be comprised of heather (Hewson 1962). The Alpine hare is found on forested slopes (Rehnus et al. 2013) at altitudes of up to 3,000 m asl (Bisi et al. 2015) throughout the Alps where it exhibits a flexible foraging strategy, but positively selects forest habitats for protection and cover (Rehnus et al. 2013). These preferences stand in contrast to the lowland habitat and soft agricultural grasses and herbs preferred by the European hare  Schai-Braun et al. 2015), thus limiting the potential for interspecific competition. It has been suggested, however, that rapid climatic changes will facilitate the expansion of European hare populations at the expense of mountain hares (Leach et al. 2015a). This is supported by our models and niche overlap metrics which describe an increasing probability of geographic overlap and, hence, interspecific interaction between European hares and neighbouring mountain hare subspecies in coming decades. It is important to remember, however, that our models do not account for land use change and, indeed, the habitat vacated by mountain hare subspecies may not be favourable for the expansion of the European hare (sensu Bisi et al. 2015).
In contrast, the heath hare and, most notably the Irish hare, are adapted to a lowland ecology in the absence of contact with the European hare. The heath hare was isolated at the most southerly tip of the Scandinavian Peninsula, far from Russian European hares, and buffered by an expansive northern hare population to the north, while the Irish hare was isolated on an island. Thus, in contrast to other mountain hare subspecies, the niches of the heath and Irish hares were, and indeed, are particularly  vulnerable to invasion by the ecologically similar European hare, with their populations susceptible to collapse upon contact with the latter after anthropogenic introductions (Thulin 2003;Caravaggi et al. 2015). After the introduction of the European hare into southern Sweden, it expanded northwards, rapidly outcompeting and hybridising with the heath hare (Suchentrunk et al. 1999;Thulin et al. 2003;Winiger 2014) which is now said to be virtually extinct (Carl-Gustaf Thulin pers. comm.), and, to a lesser extent, the Northern hare . The continued northward expansion of the European hare into Sweden has only been slowed by the southern extent of the snowline and lower temperatures (here, the Northern hare was associated with a high degree of seasonality in both temperature and precipitation, and coniferous forest in our ENMs); conditions to which the European hare is not well adapted. However, our models and niche overlap metrics suggest that a rapidly changing climate will facilitate the continued expansion of the European hare northward in Sweden. It is highly likely that this will be to the detriment of the native Northern hare. A similar pattern of species replacement and northern range restriction may be underway in Finland, following the natural expansion of the European hare from Russia (sensu Tiainen and Pankakoski 1995;Syrjälä et al. 2005;Levänen et al. 2015).
The Irish hare occupies temperate lowlands, where it is most abundant in agricultural pastures (Whelan 1985;Reid et al. 2006 and feeds predominately on grasses (Strevens and Rochford 2004). Our ENMs predicted Irish hare range extent almost exclusively using temperature seasonality (a measure of climatic stability). Seasonal temperatures in Ireland are relatively stable due the maritime influence of the warm Atlantic Conveyor, from a mean daily temperature of 19°C in the summer to 2.5°C in winter (www.met. ie). In contrast, Great Britain and central or northern Europe experience greater seasonal variation, the latter experiencing freezing winters with long periods of snow cover and hot summers, which are rare in Ireland. Thus, when interpreting the results of the Irish hare model, one must acknowledge that while the modelled niche space is necessarily restricted to conditions within Ireland, the potential niche of the Irish hare must be broader given its observed ecological plasticity and adaptability. It is impossible to quantify, however, the degree to which its true niche has been truncated, if indeed that is the case. Moreover, while the Irish hare has evolved to exploit a wide range of habitats and food items, adapted over thousands of years, this does not necessarily equate to short-term adaptive potential. Our analyses provide additional well-quantified ecological evidence, to add to observed genetic, phenotypic and behavioural differences, by which the Irish hare might be judged to warrant full species status (e.g. Arribas and Carranza 2004). No single line of evidence is sufficient to resolve the issue of specific-status, but taken together, the evidence for the Irish hare being at the very least an  (Reid 2011), with most failing to become established . Our ENMs predict that most of Ireland is (and presumably was) unsuitable for the European hare, providing a potential explanation as to why most introductions failed. At present, there is a relatively range-restricted population of introduced European hares in Northern Ireland (Caravaggi et al. 2015), the only region of Ireland currently predicted by our ENMs as being suitable for the species. Their range expanded three-fold between 2005 and 2012/2013 (Caravaggi et al. 2015) with a core range populated solely by the invader being established recently (Caravaggi et al. 2016). Moreover, there is a high degree of multi-generational, bidirectional hybridisation in areas of sympatry (Hughes et al. 2009;Prodöhl et al. 2013). Certainly, the Irish hare shares much in common with the European hare, though the latter has a greater association with arable land, typical of its more easterly and central European distribution (Smith et al. 2005). Indeed, previous studies have demonstrated that both species exhibit comparable niche breadths and almost complete niche overlap in Ireland Caravaggi et al. 2015). Niche overlap metrics in the present study are, on the face of it, contradictory. However, it must be remembered that the European hare niche described herein is reflective of its entire range, whereas niche similarities described by previous studies were confined to Ireland. Niche overlap estimates derived from ENMs are associated with geographic overlap of generated probability surfaces (Warren et al. 2008). Our estimates, therefore, must be placed in context of the small range of the European hare in Ireland at present, and future projections of European hare range increase and Irish hare range decrease. Climatological projections under future climate change suggest Ireland is likely to get warmer and drier, though with heavier winter rainfall , favouring an increase in arable agriculture . ENMs suggest such changes will result in a dramatic reduction in the suitable bioclimatic envelope for the Irish hare. As such, while the Irish hare exhibits considerable ecological plasticity and adaptability, the species may struggle under increased climate instability. Remaining suitable areas are likely to be in the cool, wet west of Ireland in habitats suboptimal for the European hare, such as peat bogs or the north-west uplands. Conversely, the bioclimatic envelope suitable for the European hare is likely to expand in Ireland in future, with conditions and any increase in arable cropland becoming increasingly favourable. Indeed, it has been suggested that invasive interactions may be more likely in areas with greater than average climatic instability (Leach et al. 2015b). The European hare may, therefore pose a direct threat to the ecological integrity of the Irish hare in the short-term (i.e. over the next 30 years or so) only, after which areas of sympatry are likely to be temporally transient. Our ENMs predict that the west of Ireland will remain suitable for the Irish hare and unsuitable for the European hare by 2070. It should be cautioned, however, that no mountain in Ireland is high enough to have a permanent snowline, nor any habitat suboptimal for the European hare expansive enough, in terms of individual patch size, to provide refuge for the Irish hare, should post-European hare introduction population dynamics mirror those of Sweden. Consequently, only occupied offshore islands in the northwest are likely to provide refuge for the Irish hare in the long-term.
Control of invasive species is frequently recommended and increasingly common (e.g. Courchamp et al. 2003). Most successful eradications have occurred on islands, where recolonisation is less likely (e.g. Imber et al. 2000;Howald et al. 2007). Invaders may be removed via biological (e.g. the introduction of a predator, competitor or pathogen), physical (e.g. shooting, trapping) or chemical (poisoning) methods (Courchamp et al. 2003). The use of any one or combination of control methods requires careful evaluation due to the potential for significant unintended consequences on non-target species and the wider ecosystem. Control or eradication of range-restricted non-native species is eminently feasible given the application of appropriate techniques, and observation of control/removal criteria (e.g. Bomford and O'Brien 1995). Once invasive species become widespread, however, eradication may be impractical and applied management becomes increasingly difficult and expensive. This is certainly the case with regards to the European hare populations in Sweden. The European hare was introduced to Sweden in the late 19 th century (Lonnberg 1905), and rapidly expanded across southern Sweden, completely displacing the heath hare from [9,000 km 2 by 1999, and establishing a considerable zone of sympatry which also extends into the southerly extent of the northern hare (sensu Bergengren 1969;Jansson and Pehrson 2007). The situation in Scotland is more nuanced, as the European hare is considered a priority species in Great Britain, despite its alien origin, as is the native mountain hare. Our niche overlap projections suggest that interspecific competition between European and Scottish hares will become increasingly common. We may expect, therefore, the Scottish hare to come under increasing pressure from what we can reasonably assume to be an ecologically dominant competitor, leading to spatial displacement. It must be remembered, however, that our projections did not account for changes in land use. Given the different habitat preferences of Scottish and European hares (e.g. Hewson 1962, Schai-Braun et al. 2015, competition and displacement may be less severe than our models suggest and may be mediated by habitat management to the benefit of the Scottish hare (e.g. the maintenance of heather moorland and other upland habitats).
Given that the current European hare population in Ireland may have been introduced as recently as the 1970s (Caravaggi et al. 2015) with its extent and numbers expanding rapidly, policy makers and conservationists in Ireland would do well to take heed of the Irish hare's likely future prospects by using Sweden as a case study example. The authorities of both political jurisdictions of Ireland are signatories to the Convention on Biological Diversity (UNEP 1992), the Bern Convention (1979), the European Habitats Directive (EEC 43/1992), and the EU Regulation 1143/2014 on Invasive Alien Species (Official Journal of the European Union [OJ] 2014), and hence, are obliged to address invasive species issues. Lessons from the heath and Northern hare in Sweden highlight the precarious, unstable nature of mountain and European hare interspecific population dynamics, and suggest that continued inaction from authorities in Northern Ireland will only serve to facilitate the continued expansion of the invader, to the detriment of the endemic.