Specific niche requirements underpin multidecadal range edge stability, but may introduce barriers for climate change adaptation

To investigate some of the environmental variables underpinning the past and present distribution of an ecosystem engineer near its poleward range edge.


| INTRODUC TI ON
Climate change poses a major threat to biodiversity through driving the global redistribution of species (Parmesan et al., 1999). At leading edges, warming temperatures can drive range extensions of organisms (Burrows et al., 2011), although cool summer or winter temperatures and extreme weather events can affect survival or reproductive success (Firth et al., 2011;Hutchins, 1947), which can limit range extension. If the leading edge adjoins a stretch of unfavourable habitat or a physical barrier to movement, individuals of a species may have 'nowhere to go' (Loarie et al., 2009). If the trailing edge 'catches up' with the leading edge, a species may be susceptible to local extirpation or even be at risk of global extinction (Marzloff et al., 2018), unless individuals can adapt to changing environmental conditions.
In the marine environment, mesoscale nearshore circulation patterns are modified by shoreline configuration, bathymetry and tides creating eddies or fronts, sharp discontinuities of physical and chemical variables characterized by different temperatures and/or salinities and hence density regimes (e.g. Pineda, 1994;Simpson & Hunter, 1974). Seasonal tidal mixing fronts, such as those found in northwest France (Le Boyer et al., 2009) or in the Irish Sea (Simpson et al., 2009), form within a few days and persist for only a few months before their dissipation. Other fronts, such as the southern Antarctic Circumpolar Current front (Orsi et al., 1995) can persist for longer and are quasi-permanent features in the ocean. Tidal mixing fronts separate vertically well-mixed and stratified water columns. They are often enriched in nutrients on the stratified side as the surface layer is usually shallower than the photic zone, leading to enhanced primary production (Sharples, 2008). When they occur, the spatially stratified environmental conditions they represent can act as barriers to dispersal for a range of taxa including fish and invertebrates (Galarza et al., 2009) and can adversely affect individual survivorship, thus making fronts potentially strong structural determinants of species distributions signified by clusters of range edges around them (e.g. Gaylord & Gaines, 2000;Sagarin et al., 1999).
Given that dispersal is a key process shaping the distribution of marine species by mediating connectivity and recruitment within and among patches (James et al., 2019), the presence and duration of fronts could affect the dispersal and population dynamics of a species, thereby affecting long-term trends in abundance and distribution (Ayata et al., 2010;Pringle et al. 2017). For instance, Banks et al. (2007) showed how a front created a persistent cold-water barrier that has prevented further advancement and colonization of the south coast of Tasmania by the invasive urchin Centrostephanus rodgersii. Consequently, range extension (or lack of) and species responses around fronts and other barriers to movement (e.g. habitat continuity, Fenberg & Rivadeneira, 2019) may well be idiosyncratic to more general predictions of range extensions under climate change scenarios. Hence, there is a need for more work in these areas if we are to better understand patterns and trends in biogeography and inform conservation.
A number of fronts are found around the islands of Ireland and Britain which straddle a biogeographic boundary between cold 'Boreal' and warmer 'Lusitanian' waters. Forbes (1858) produced the first map of the broad distributional limits of marine species of Britain and Ireland with a line depicting the 'general limit of southern types' (herein "Forbes' Line," Figure 1). In response to recent warming, many species have exhibited range extensions beyond previously known limits that coincided with this line in Britain; both northwards through the Irish Sea and eastwards along the English Channel (Mieszkowska et al., , 2014Southward et al., 1995).
In contrast to the surrounding warmer waters of the Atlantic Ocean and Celtic Sea, the Irish Sea is dominated by cooler, mixed water in summer (Simpson & Hunter, 1974, Figure 2a Devoy, 2008;Pingree & Griffiths, 1978). For more detail on spatio-temporal variation in sea surface temperature see Figure S1. (b) Coastal configuration, spatial changes in wave climate and variation in spring tidal ranges around Ireland (map redrawn from Devoy, 2008 with wave heights from Bricheno & Wolf, 2018). Maps created by Tim Absalom, University of Plymouth GeoMapping Unit Fronts separate the cold water of the Irish Sea from the warmer Atlantic Ocean and Celtic Sea, respectively ( Figure S1). Within the Irish Sea, a patch of warm water develops in the summer that is bounded by the Western Irish Sea Front and the Western Irish Sea gyre. This creates a mosaic of patches of warmer and cooler waters around the island. Crisp (1989) was one of the first to describe the influence of tidal fronts on the distribution of intertidal taxa. He described a positive association between the reef-forming worm Sabellaria alveolata (Linnaeus, 1767), the topshell Phorcus lineatus (da Costa, 1778) and chthamalid barnacles (all warm-affinity taxa) with the patch of warm summer-stratified water bounded by the Western Irish Sea Front ( Figure 2a). Since Crisp (1989), surprisingly little research has been done, on either investigating patterns of distribution of intertidal taxa in relation to tidal fronts, or indeed the mechanisms underpinning these patterns (but see Ayata et al., 2010). The presence of so many tidal fronts combined with a complex coastline ( Figure 2b) and a mosaic of mixed and stratified waters in a biogeographic boundary region presents an ideal study system for testing hypotheses about proximate environmental factors determining the range edges of climate change indicators at their poleward edge of distribution.
The honeycomb worm, S. alveolata is protected under the EU Habitats Directive because it is considered an ecosystem engineer that creates a range of biogenic structures (sheets, veneers, hummocks to the largest reefs in Europe, Curd, et al., 2019), which provide habitat for a multitude of other species throughout its geographic range (Cole and Chapman, 2007;Plicanti et al., 2017;Jones et al., 2018). Its planktotrophic larvae settle and metamorphose preferentially on the cemented sand tubes of conspecific adults (Wilson, 1968), typically in areas where rocky reefs abut sandy beaches that can supply coarse sand for tube-building (Firth et al., 2015). It is often described to be distributed from North Africa to southwest Scotland (Dubois et al., 2002), but it is also present in Ireland (Culloty et al., 2010). Despite very little being known about its distribution in Ireland (only 40 records in Oceanographic Biodiversity Information System (OBIS) database, representing 1.6% of total records), Crisp (1989) identified a leading edge on the east coast in the northern Irish Sea. It is possible that another leading edge exists on the northwest coast where it was also known to occur in the past (Duerden, 1895). In Britain, where there is a strong heritage of sustained observations (Hawkins et al., 2013;2,357 records in OBIS database, representing 93.6% of records), it is considered a climate change indicator that is expected to benefit from climate change (Mieszkowska et al., 2013) due to exhibiting population increases in response to warming temperatures  and decreases in response to extreme cold events (Crisp, 1964;Firth et al., 2015). Bush (2015) identified the leading edge in Britain to be in the northern Irish Sea at Auchenmalg, southwest Scotland.
Here, we investigated some of the environmental variables underpinning the distribution of S. alveolata around Ireland. We combined historical data from a range of sources with contemporary data from broad-scale field surveys to describe past and present distributions to build up a detailed picture of the fine-scale distribution.
Building on existing comprehensive unpublished datasets from the 1950s and 2000s, we revisited 60 target locations, alongside many more (283) in the 2010s to quantify spatio-temporal changes in distribution and abundance during a period of rapid warming. Then, using a habitat suitability modelling approach, we assessed some of the environmental drivers which we considered to be important in underpinning the distribution of S. alveolata. The outputs provide insights into some of the environmental variables underpinning the biogeography of this key ecosystem engineer and will help to inform conservation and management plans for a species that may be susceptible to climate change and other anthropogenic stressors.  Table S1 for descriptions), some estimated abundances per unit area, whilst others simply noted presences. Where a SACFOR classification was not specified, we assigned 'present'. All records that were not collected by the authors were carefully scrutinized and any dubious records were omitted.

| ME
Comprehensive island-wide surveys were conducted at the same 60 locations in the 1950s (Southward & Crisp, 1954), 2000s (Simkanin et al., 2005) and in the 2010s. The 2010s surveys were done by authors of this paper, who also revisited an additional 75 locations sampled in either the 1950s or 2000s, along with 208 new locations to chart the broad-scale distribution of S. alveolata on different shore types. For these surveys, S. alveolata abundance was estimated following a 30-min search at each location (typically by two people, covering up to 1000 m 2 depending on the topography) using the modified semi-quantitative SACFOR abundance scale (Table S1).
Following the methodology of Firth et al. (2015), 'Abundant' and 'Superabundant' were combined as 'Abundant' to ensure comparability with Crisp and Southward, as 'Superabundant' was adopted in the 1970s. Changes in abundance were considered in two ways-a conservative two-category change (e.g. Frequent to Abundant) and a less conservative one-category change. Intertidal habitat substrata were categorized as rock and boulders, coarse substrate, sand, muddy sand, sandy mud, mud, and mixed sediments. Due to lack of habitat information, five data points were excluded from analysis from the Connemara coastline (central-west region). To assess sensitivity to alternative distribution modelling techniques, in the first instance we fitted generalized linear models (GLM), generalized additive models (GAM), multivariate adaptive regression splines (MARS), random forests (RF) and boosted regression trees (BRT), as well as support vector machines (SVM) models to presence/absence data using the sdm R package (Naimi & Araújo, 2016). Using default model settings, all methods yielded similar results both in terms of variable influences on predicted probability of presence and cross-validation accuracy levels (i.e. area under the curve (AUC) varied from 0.81-0.91 for GAM and SVM, respectively). Following guidance in Elith et al., (2008), BRT models underwent further tuning leading to being adopted after achieving the highest AUC (~0.92 on 10-fold random cross-validation versus ~0.88 on 10-fold cross-validation with spatial blocks) of all models. Briefly, this AUC score was achieved by fitting a Bernoulli distribution, implemented in the gbm.step function of the dismo package v. 1.1-4 (Hijmans et al., 2017)

| Distribution modelling
with a learning rate of 0.001, tree complexity of 5, and a bag fraction of 70% (meaning that at each iteration of the model, 70% of the data are drawn at random, without replacement, from the training set).
Cross-validation accuracy levels were based on 10-fold cross-validation and spatial blocks assuming a ~67 km autocorrelation range as estimated on the 11 quantitative environmental layers (excluding substrate type) using the blockCV R package (Valavi et al., 2019).
Folds were stratified by prevalence so that each fold roughly contained the same proportion of presences and absences. Relative importance of each explanatory variable (determined by the reduction in root mean square error (RMSE) following introduction of a variable) and partial dependencies were produced using the ggbrt package (Jouffray et al., 2019).

F I G U R E 3
Map illustrating the locations around Ireland where Sabellaria alveolata was ever recorded as present. Full circles = known presences; empty circles = known absences; Boxes illustrate discretely bounded regional sub-populations: the northeast; the south; Dingle Peninsula; Galway Bay; northwest; Lough Swilly.

| Mapping of past and present distributions
A total of 981 records spanning 182 years (1836-2018) and the whole coastline of Ireland were collated during this study, comprising 319 presences (33%) and 662 absences (not seen) (67%) (data available in Curd, Cordier, et al., 2020). 776 (79%) records were collected by either the authors of this paper or by Crisp and Southward. The remaining 205 (21%) were from a combination of online databases, unpublished theses, reports, published papers and museum specimens. The vast majority of the records (954; 97%) were from intertidal surveys (Figure 3a).
Given the paucity of subtidal records, we focus on the intertidal dataset.
Sabellaria alveolata exhibited a discontinuous distribution with six discretely bounded regional sub-populations (Figure 3a). It was found at a large number of locations along the south coast (Cullenstown, Co. Wexford to Galley Head, Co. Cork), in Galway Bay, and in the northwest (Killala Bay, Co. Mayo to Fintra Beach, Co. Donegal). It was also found in smaller pockets on the Dingle Peninsula in the southwest, in Lough Swilly on the north coast, and along the northeast coast (Coney Island, Co. Down to Howth, Co. Dublin, where it was typically less abundant (i.e. exhibited lower SACFOR estimates)).
The Lough Swilly sub-population abuts the Islay Front; the southern sub-population abuts the Celtic Sea Front to the east, and the northeast sub-population is bounded to the north and south by the

| Distribution modelling
Cross-validation of BRT predictions across independent spatial blocks yielded an area under the curve (AUC) of 0.88, suggesting a high degree of reliability in presence/absence classification.
Moreover, BRT spatial predictions of S. alveolata distribution also qualitatively match regional patterns of observed occurrences Note: Colonizations denote when Sabellaria alveolata is classed as present after being classed as Not Seen (absent) in the preceding comparable sampling period. 1-step and 2-step changes indicate either a one or two step SACFOR change (e.g. Common to Abundant = 1-step, Frequent to Abundant = 2-step). Extirpations denote when S. alveolata was classed as Not Seen after being classed as present in the preceding comparable sampling period. a Denotes three locations in the northeast where S. alveolata was present in 1950s, not seen in 2000s but then recolonized the same locations in the 2010s. See Figure 4 for graphical representation.

TA B L E 1
The number of locations representing each category of abundance change between sampling periods where S. alveolata was present on at least one of the three sampling periods substrate type were the four most influential variables (in decreasing order) (Figure 6b). Probability of presence marginally increases with exposure to waves, up to a threshold of 1.8 m in significant wave height, beyond which S. alveolata occurrence is highly un-

likely. Sabellaria alveolata predicted occurrence is low in areas
where tidal amplitude is <1.8 m, and reaches an optimum where tidal amplitude is between 1.8-2.7 m. Its occurrence is also most likely in regions of higher stratification (stratification index > 3.2, Figure 6b). It typically occurs on rock and boulder, sand, muddy sand and mixed sediment but is typically absent from sandy mud or coarse sediment-dominated substrates (Figure 6b). Minimum sea surface temperature was the fifth most important environmental variable (with a ~ 7% contribution to model predictions). However, its influence is limited to probability of occurrence marginally increasing when minimum temperatures are above ~ 6°C (See Appendix S2).
We also considered the interactions between the most influential variables ( Figure S3). There was a significant interaction between mean wave height and tidal amplitude, stratification index, substrate type, minimum SST and wind speed respectively. In general, when wave height is >1.8 m (i.e. in exposed areas), probability of S. alveolata occurrence is marginal and other variables have little influence on model predictions. Sabellaria alveolata is predicted to thrive in areas where exposure to waves is moderate (wave height between 1.3-1.8 m), in particular where tidal amplitude is moderately high (between 2-2.5 m), in stratified waters (stratification index > 3.5), when minimum SST is above 6°C and when wind speed is between 11-12 m/s.  Clougherhead: 12 records 1952Clougherhead: 12 records -2018Clougherhead: 12 records , present in 17% of records (O in 1958Clougherhead: 12 records , 1984. Balbriggan 1958Balbriggan -2016Balbriggan :4 records, present at 75% of records (O in 1958Balbriggan , C in 2013Balbriggan , F in 2016

| D ISCUSS I ON
We provide a comprehensive account of broad-scale (>7,400 km) and long-term (182 years) changes in the distribution and abundance of a widespread ecosystem engineer towards its poleward distributional edge around Ireland. We collated 981 records of presences (33% of records) and not found-presumed absences (67%). Collating absences over consecutive surveys not only provides strong evidence of true regional absences, but also enables much greater power in statistical model predictions. Sabellaria alveolata exhibited a discontinuous distribution with multiple discretely bounded regions, three of which persistently aligned with tidal fronts. Two of these fronts Lough Swilly, with a persistent distribution gap in between. Whilst it was restricted to six geographic regions, the sub-populations within these regions and the previously recognized edges appeared to be relatively stable through time, with evidence of local increases within the regions but no evidence of range extensions beyond the previously identified poleward range edge in the northeast (Crisp, 1989).
Similarly, Bush (2015) found the poleward range edge in Britain to be persistent between the 1980s and 2012.
Despite mounting evidence that S. alveolata responds positively to warming temperatures and negatively to cold winters (Crisp, 1964;Firth et al., 2015;Mieszkowska et al., 2005), its two leading edges in Ireland appear to be stable. As such, it is perhaps more appropriate to refer to these as poleward-most range edges, rather than leading edges which may imply some degree of movement. In the face of continued warming, if the trailing edges catch up with the polewardmost edges, this important ecosystem engineer may be susceptible to extirpation and even extinction, unless it can adapt or disperse through other means (e.g. intentional or unintentional humanassisted migration (Firth, Knights, et al., 2016;Woodin et al., 2014)).
The northeast appeared to be highly 'dynamic', with evidence of declines, potential local extirpations and recolonizations. In this region, it is considered that a combination of reduced wave weights and tidal amplitude, coupled with being surrounded by cooler summer waters outside the Western Irish Sea Front in addition to being bounded by unsuitable substrate to the south make these populations potentially vulnerable to extinction risk as discussed below.
Our modelling suggests a combination of suitable wave height, tidal amplitude, local water stratification and substrate type are pri- where tidal amplitude is moderately high (between 2-2.5 m) and where water is stratified in the summer (stratification index > 3.5).
Whilst minimum sea surface temperature (SST) was the fifth most important variable, it only had a marginal influence on predicted probability (~7%) of S. alveolata occurrence. In this region, SST variability was relatively narrow (i.e. 10.5-12°C and 4-7°C for mean and minimum SST, respectively), and thus, there is relatively little variance in temperature to explain. This is further exacerbated by temperature data being averaged over multiple years (Bates et al., 2018).
Future work should consider the influence of extreme and prolonged periods of low temperatures, particularly during the reproductive season, when gametes and larvae may be negatively affected by low water temperatures .
Through comparing 60 locations across three time periods (1950s, 2000s, 2010s), we found that population increases were observed at one time or another at all but three locations. We also discovered that S. alveolata underwent declines in the 2000s at 42% of locations where it was recorded as present in the 1950s. Of these declines, 80% were extirpations, 75% of which were in the northeast.
The fact that it was subsequently recorded at all of these locations in the 2010s suggests eventual recovery following the potential occurrence of a damaging extreme event or a sustained period of colder weather following the 1950s surveys (Éireann, 2019). Sabellaria alveolata and other taxa (Firth et al., 2015, see Mieszkowska et al., 2007 for work on Phorcus; Southward, 1991  Abundant on the SACFOR scale, meaning that they were unlikely to be forming high-density reefs which would buffer individuals from thermal extremes and other environmental perturbations (Bertness & Leonard, 1997). Southward and Crisp (1954) identified that rocky shore communities in the northeast were susceptible to excessive winter cooling of sea surface temperature due to the easterly and south-easterly aspect of the shoreline. Furthermore, some of the locations in this region are considered vulnerable to physical damage from snail (winkle) harvesters through trampling and breaking off chunks of reef (Preston & Portig, 2001). The dynamic nature of the northern part of S. alveolata's range is likely due to a combination of factors including cooler temperatures (exacerbated by the easterly aspect, Firth, White, et al., 2016), damage from trampling and harvesting and limited dispersal or immigration in the presence of multiple tidal fronts, rendering these populations highly vulnerable.
These populations are perhaps even more vulnerable because they are likely to be genetically isolated. Population genetic analysis, whilst based on a single population from Galway Bay, shows that S. alveolata from Ireland are genetically differentiated from other European and North African populations (Muir et al., 2020).
Genetic structure at the regional scale around Ireland has not yet been evaluated for this species, but hydrodynamic modelling studies suggest that barriers to dispersal could exist, most likely due to tidal fronts (Robins et al., 2013), restricting connectivity among discrete populations. For instance, the Western Irish Sea Front and gyre can promote local larval retention, self-recruitment and genetic isolation of mussel populations within the Irish Sea (Gosling et al., 2008), and the Islay and Celtic Sea Fronts can further prevent larvae entering the Irish Sea (Lynch et al., 2004;Wilmes et al., 2019).
Genetic diversity of S. alveolata is greatest in its northern part of the distribution in the Irish Sea and English Channel due to the preservation of genetic diversity in refugia over Pleistocene glacial cycles (Rigal, 2005), making this hotspot of genetic diversity vulnerable to population instability. For example, population collapse following the extremely cold winter of 1962/63 (Firth et al., 2015) is evoked as the likely cause of reduced genetic variation in populations in the eastern Irish Sea (Nunes et al., accepted). Despite stochasticity being evident in local and regional populations (i.e. within the discrete regions shown in Figures 3 and 4), our results suggest a complex population structure in S. alveolata and 'portfolio effect' (Schindler et al., 2015) at an island-scale, arising from an interplay between connectivity, gene flow with local adaptation, environmental pertubations and recruitment. Hotspots of genetic diversity are likely in a fragmented network of metapopulations with only occasional genetic interchange with different source populations. Such occasional recruitment may lead to Allee effects (decreased population growth or fitness at small population sizes), but any exchanges may buffer these and prevent local adaptation. Ongoing investigations into the dispersal and population genetics of S. alveolata around Ireland will cast further light on this.

Sabellaria alveolata reefs and their surrounding sedimentary
habitats are known to support species of commercial importance (Dubois et al., 2002;Plicanti et al., 2016;Schimmenti et al., 2015) and the worms themselves are also collected as fishing bait, particularly in the Mediterranean (Gambi et al., 1992). In France, S. alveolata reefs have long been exploited for bivalves, causing the authorities to restrict harvesting activities in 1970 after a marked decline of the Mont-Saint Michel reef (Dubois et al., 2002). Recently, the non-native oyster, Magallana (formerly Crassostrea) gigas has become increasingly abundant on S. alveolata reefs in Mont-Saint Michel (Dubois et al., 2006), exacerbating physical damage to the reefs directly through overgrowth on worm tubes and indirectly through physical damage by harvesters (Dubois et al., 2002). Further information on damage caused by harvesting and associated invasive species would be beneficial for informing on how multiple stressors affect the vulnerability of this species.
The historical reach of this study coupled with the fine-scale mapping revealed persistent regional sub-population boundaries that align with tidal fronts, two of which are northernmost range edges. Crisp (1989) was the first to point out that intertidal fauna showed discontinuities around fronts in Ireland due to the impact of mixing through tidal energy. In that paper, he also suggested that long stretches of sandy coast that were subjected to severe scouring contribute to the biogeographic distribution patterns of both warmadapted and cold-adapted species. The absence of S. alveolata in the southern Irish Sea between Dublin and Carnsore Point can largely be explained by the existence of long stretches of sandy coast. Here, we provide quantitative evidence that a combination of suitable substrata and a suite of hydrographic variables (tidal amplitude, stratification index and in particular, wave height) may act as environmental filters (Webb et al., 2002) and are likely to be key explanatory variables that underpin the distribution and persistence of S. alveolata.
These can be used to predict the presence of this ecologically im-

| CON CLUS IONS
In 1858, Edward Forbes drew his "general limit of southern types" has a discontinuous distribution with two poleward-most range edges and six discretely bounded regional sub-populations. Using a combination of statistical model outputs and expert judgement, we identified a potential suite of specific mesoscale hydrographic niche requirements and thresholds required for population success.
We also showed that the sub-population boundary edges that align with tidal fronts were largely stable over the past 62 years; this is despite (a) global and regional increases in sea surface temperatures that would predict range extension by this species, and (b) local and regional population fluctuations within the range that are associated with climate change and extreme weather events. This suggests that a combination of specific niche requirements associated with frontal regions have facilitated the long-term regional stability of this species, but may in the future act as barriers to range extension. Thus, these environmental boundaries may limit the capacity of this important ecosystem engineering species to redistribute in response to future climate change.

ACK N OWLED G M ENTS
Dedication: This paper is dedicated to Edward Forbes FRS FGS (1815-1854) for his pioneering work on biogeography. Alfred Russel Wallace is often referred to as the "Father of biogeography" and Wallace's Line between Borneo and Sulawesi is globally recognized.
The line drawn by Forbes representing the "general limit of southern types" is less well known. Perhaps it was the coining of the term "Wallace's Line" and advocating of the concept by the influential Thomas Henry Huxley that facilitated this line becoming the most famous biogeographic demarcation in the world. Since then, a large number of major biogeographic boundary zones have been identified globally, but few have actually been named. We propose the naming of "Forbes' Line" (Figure 1)

PEER R E V I E W
The peer review history for this article is available at https://publo ns.com/publo n/10.1111/ddi.13224.

DATA AVA I L A B I L I T Y S TAT E M E N T
All data used in this paper are available in: Curd, Cordier, et al. (2020).