Combined mechanistic modelling predicts changes in species distribution and increased co‐occurrence of a tropical urchin herbivore and a habitat‐forming temperate kelp

Identify climate change impacts on spawning and settlement of a tropical herbivore (Tripneustes), optimal habitat of a temperate kelp (Ecklonia) and implications for these species regions of interaction under future climate.


| INTRODUC TI ON
Globally, a large proportion of terrestrial and marine species are responding to climate change by undergoing range redistributions (Chen, Hill, Ohlemüller, Roy, & Thomas, 2011;Poloczanska et al., 2013), with effects on biodiversity and ecosystem structure (Pecl et al., 2017). Warming is one of the primary causes leading to species range shifts (Chen et al., 2011;Poloczanska et al., 2013). All species have an optimal thermal range where they best function; outside of this range, performance decreases and extreme temperatures lead to demise (Angilletta, 2009;Pörtner & Knust, 2007). Ocean isotherms are shifting at an equal or faster rate than terrestrial isotherms, resulting in a greater range shift for marine species, as they move to new areas that encompass their thermal range (Burrows et al., 2011;Poloczanska et al., 2013;Sen Gupta et al., 2015).
Dispersal is a central process in distributional shifts (García Molinos, Burrows, & Poloczanska, 2017;Travis et al., 2013); it is particularly important for marine organisms as most species rely on the release of reproductive propagules into the water column where they are then carried by ocean currents (Kinland & Gaines, 2003) to replenish populations across their biogeographical range. Oceanic conditions strongly influence dispersal, as temperature often triggers reproduction and spawning events (Fincham, Rijinsdorp, & Engelhard, 2013;Yang & Rudolf, 2010) and ocean currents determine propagule trajectories influencing settlement patterns (Werner, Cowen, & Paris, 2007). Thus, diagnosing the impact of climate-driven circulation changes on larval dispersal is important to anticipate changes in species ranges (e.g. Cetina-Heredia, Roughan, van Sebille, Feng, & Coleman, 2015).
The Biotic-Abiotic-Mobility (BAM) framework highlights the three key factors that should be considered when predicting species range shifts in response to climate change (Peterson, Papeş, & Soberón, 2015). However, changes in species distributions are commonly predicted using correlative models such as species distribution models (SDMs; Peterson et al., 2015), which focus mostly on abiotic factors, by statistically relating known occurrence data of a species to environmental variables in order to predict a species range (Guisan & Zimmermann, 2000). Differently, mechanistic models use a biophysical approach to associate population processes with abiotic conditions. For instance, Rodríguez, García, Carreño, and Martínez (2019) forecasted Millepora alcicornis future habitat constructing a physio-climatic predictor which combined the coral maximal quantum yield of photosynthesis with sea surface temperatures (SSTs) into a raster reflecting coral optimal conditions. Nonetheless, there are other evolutionary and ecological processes which determine species distribution, such as reproduction and dispersal (i.e. mobility; Leroux et al., 2013), which are often lacking or simplified in correlative or mechanistic approaches (Robinson et al., 2011). Leroux et al. (2013) modelled the impact of climate change for North American butterflies, through a reaction-diffusion framework, which considered the population growth rate (which depends on both abiotic conditions and biological interactions), climate envelope (i.e. abiotic conditions) and diffusion rate (i.e. rate of movement) for the butterflies. Though Leroux et al. (2013) integrated all three factors of the BAM framework in the model, the research highlighted the need to also consider changes to species interactions which are often not considered but can be fundamental in determining ecosystem function and structure (Blois, Zarnetske, Fitzpatrick, & Finnegan, 2013).
Tropicalization offers a striking example of how distribution shifts of tropical species into temperate habitats can create novel biotic interactions, with consequences cascading throughout entire ecosystems (Vergés, Steinberg, et al., 2014). Tropicalization has been seen along Western Boundary Currents (WBC) such as the East Australian Current (EAC), whose climate-driven intensification (Wu et al., 2012) contributes to the greater transport of warm tropical water and larvae being carried into temperate ecosystems (Vergés, Steinberg, et al., 2014). Globally, this has resulted in new herbivore-algae interactions that have led to phase shifts, as areas once inhabited by dense kelp forests turn into turf algae dominated reefs, with tropical consumers preventing the re-establishment of temperate kelp (Bennett, Wernberg, Harvey, Santana-Garcon, & Saunders, 2015;Vergés et al., 2016Vergés et al., , 2019Vergés, Tomas, et al., 2014). This example highlights the importance of understanding and identifying climate-mediated changes in species interactions to understand future ecosystem composition and enable adaptive management decisions (Koehn, Hobday, Prachett, & Gillanders, 2011).
We develop a new combined mechanistic modelling approach to predict distributional range shifts and changes in regions of co-occurrence that allow for biotic interactions under climate change.
We use a case study that forecasts future range shifts of a tropical herbivore, the tropical sea urchin Tripneustes gratilla (hereafter Tripneustes) and a temperate habitat-forming kelp Ecklonia radiata (hereafter "Ecklonia"). Ecklonia forms ecologically and economically important kelp forests that support a high diversity of organisms  tropicalization et al., 2019). Tripneustes is an important herbivore that preferentially consumes Ecklonia when a range of seaweeds are on offer, as demonstrated through both aquaria assays (Dworjanyn, Pirozzi, & Liu, 2007) and field experiments (Vergés et al., 2016). Further, Tripneustes is able to overgraze seagrass and macroalgae forests when present in high densities, as was seen in east Africa and Lord Howe Island (LHI), respectively (Eklöf et al., 2008;Valentine & Edgar, 2010). To develop our model, first, we characterize the fundamental niche and the realized niche (sensu Hutchinson (Vandermeer, 1972)) of thermal tolerance of Ecklonia. We then use temperature predictions under future Representative Concentration Pathway 8.5 (RCP 8.5) carbon emissions scenario to predict its distribution in the future. RCP 8.5 is a "business as usual" scenario where emissions and the population continue to grow throughout the 21st century (van Vurren et al., 2011). Second, we model the impact of a changing climate on the dispersal and geographical patterns of settlement of Tripneustes using the Connectivity Modelling System (CMS; Paris, Helgers, van Sebille, & Srinivasan, 2013). CMS is a Lagrangian-tracking algorithm that allows implementing organism's life traits to simulate the advection of particles such as larvae by ocean currents. Finally, we combine these results to determine the new overlapping areas where biotic interactions are likely to occur under this future climate scenario.
Through this case study approach, we seek to answer the following two questions: (a) how will spawning and settlement of Tripneustes larvae change in the future along the east coast of Australia? And (b) could Tripneustes pose an increased threat to Ecklonia in the future by increasing the overlap, and therefore the potential to interact, between these two species?

| Model details
The oceanographic data used to forecast the impact of climate change on the dispersal of Tripneusteus and distribution of Ecklonia came from the Ocean Forecasting Australian Model version 3 (OFAM3; Oke et al., 2013) configured to downscale and project future climate as detailed in Feng, Zhang, Sloyan, and Chamberlain (2017) and Zhang et al. (2016). The climate projection includes changes in temperature and circulation. OFAM3 is an eddy-resolving ocean model (1/10° spatial resolution) that has been thoroughly validated in the study region of eastern Australia (Oke et al., 2013).
Due to the complexity of the systems represented in climate models, some physical processes are parameterized; such parameterizations can lead to biases or spurious trends denominated climate drift (Sen Gupta et al., 2012). An advantage of the OFAM3 configuration is that it reduces model drift by using non-adaptive relaxation (Zhang et al., 2016). Relaxation is a forcing term that prevents the modelled fields to diverge from a known climatology; in the OFAM3 configuration such forcing does not depend on differences between the climatology and the model state; thus, climate change signals are not masked (Zhang et al., 2016). Additionally, the climate signal used to force this OFAM3 configuration is an ensemble of climatologies from 17 CMIP5 climate models , averaging out non-systematic biases of individual climate models, and therefore avoiding the introduction of errors associated with a specific climate model.

| Optimal temperature range distribution of Ecklonia
Ecklonia is distributed throughout temperate Australia and its northern range edge lies within the temperate-tropical transition zone, where further expansion equatorward is mostly limited by temperature (Lüning, 1990;Martínez et al., 2018). We characterized the realized thermal niche of Ecklonia: its actual distribution with respect to temperature by using the mean temperature range of its present habitat in eastern Australia. Specifically, we use depth-averaged (surface to the 75 m isobath) temperatures from the contemporary scenario (2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015) to calculate the temperature mean and standard deviation (SD) at Brunswick Heads (28.6°S) and at the southernmost point of Tasmania (43.6°S) over this 10-year period. The mean temperature + SD (23.9°C + 1.7°C) at Brunswick Heads was identified as the upper temperature threshold where presently Ecklonia survives at all depths (i.e. surface to 75 m). Similarly, the mean temperature − SD (14.4-1.3°C) in southern Tasmania was identified as the lower temperature threshold for Ecklonia's present habitat range. We then mapped the proportion of time (%) that daily depth-averaged temperature fell within this range in both contemporary and future scenarios. We found that at the northernmost location where Ecklonia is still observed at depth (>25 m in Henderson), these temperature conditions were met 70% of the time. We then used the 70% contour in the future scenario to infer the equatorial limit of Ecklonia's distribution in the future. We did not focus on changes in the poleward limit of their distribution because future temperatures in those regions are not projected to shift outside of Ecklonia's optimal thermal range (Martínez et al., 2018), and poleward expansion is indeed limited by the continental shelf.
We then characterized the fundamental thermal niche of Ecklonia by compiling reported upper survival threshold from the literature.
This climate envelope determines the potential niche (sensu Grinnell; Vandermeer, 1972), that is, the potential geographical distribution considering temperature as the only limiting factor. Wernberg, de Bettignies, Joy, and Finnegan (2016) found that Ecklonia could tolerate short-term temperatures of up to 26.5°C for ~45 min, after which net photosynthetic rate, an indicator of stress (Hurd, Harrison, Bischof, & Lobban, 2014), significantly decreased. We constructed maps of the proportion of time that the daily depth-averaged (surface −75 m isobath) temperature was >26.5°C for both the contemporary and future scenarios. At Ecklonia's northernmost location (i.e. Henderson), the daily depth-averaged temperature was >26.5°C in a contemporary scenario, <20% of the time. Therefore, the 20% contour in the future scenario was used to predict the equatorial limit of Ecklonia's distribution based on tolerance to maximum temperatures. We then superimposed the realized thermal niche onto these maps and compared the output of distribution models based on the two different niche indicators.
In order to explore changes to overlap between Ecklonia and

| Spawning and settlement behaviour of Tripneustes
A Lagrangian framework, which uses as reference the position of a particle over time, is ideal to study the transport of larvae by ocean currents. We used the CMS (Paris et al., 2013) to simulate the dispersal of virtual Tripneustes sea urchin larvae. CMS is a Lagrangian algorithm that allows simulating the advection of particles by ocean currents while incorporating species-specific life traits. In this study, larvae are advected with 3-dimensional velocity fields produced by the ocean model (OFAM3), that is with ocean current vectors defined by north-south, east-west and up-down components, which are used to advect larvae both horizontally and vertically. In addition to advection by the velocity fields, diffusion is implemented following Cetina-Heredia et al. (2015. Spawning and settlement criteria are dictated by Tripneustes reproductive behaviour and larval development. Tripneustes is a fast-growing sea urchin commonly found in tropical/subtropical regions globally (Lawrence & Agatsuma, 2013).
To simulate contemporary and predict future dispersal of Tripneustes larvae, we based our particle releases on known biological spawning behaviours of adult Tripneustes and simulated advection over a time period that corresponds to its pelagic larval duration (PLD).
The conditions that enable spawning of Tripneustes are temperaturedependent (Chang-Po & Kun-Hsiung, 1981). Mos, Cowden, Nielsen, and Dworjanyn (2011) describe that Tripneustes could be consistently induced to spawn every 4-6 weeks when kept at a constant 25°C and fed an abundance of macroalgae. Rahman, Tsuchiya, and Uehara (2009) found that healthy embryonic development of Tripneustes occurred at temperatures between 22 and 29°C.
Tripneustes has been found up to a depth of 75 m (Lane, Marsh, VandenSpiegel, & Rowe, 2000). We combined these conditions and allowed the daily release of particles, from the surface up to the 75 m isobath at all model grid-cell locations, if the mean temperature fell within 22-29°C over the past 6 weeks. A maximum of one particle could be released per day at any location (given by latitude, longitude and depth). As we had no larval supply or adult population data for specific locations (Everett et al., 2017), we could not include this information. Hence, spawning was based on empirical data of temperature suitability for spawning, which also appeared to reflect present distribution and abundance of Tripneustes. For instance, in our model Tripneustes was able to spawn all year round in the tropics as has been suggested by Malay, Juinio-Menez, and Villanoy (2000), with less spawning in temperate regions where abundance is very low (Williamson, 2015). Hence, spawning in our model reflected only "potential spawning" of larvae.
Tripneustes has a PLD of 15-52 days (Scholtz, Bolton, & Macey, 2013). Therefore, we considered virtual larvae could settle after 15 days and up to 52 days if they came close to the coast (i.e. anywhere between the coastline and the 75 m isobath). No other settlement conditions were considered due to lack of information or model resolution, and hence, we only consider "potential settlement."

| Data analysis of spawning and settlement across regions and seasons
To understand changes to potential spawning and settlement in both contemporary and future scenarios, we looked at differences in the total yearly number of virtual larvae released or settled within subtropical, tropical and temperate regions, as well as into downstream regions (Tasmania and LHI; Figure S1).
Understanding shifts between these regions can help determine whether Tripneustes is likely to contribute to the tropicalization of temperate Australia in the future (i.e. an increase in the proportion of Tripneustes inhabiting temperate ecosystems and a poleward shift in Tripneustes distribution). It also enabled us to target specific areas of interest, such as LHI, where a population outbreak of Tripneustes in the past caused significant loss of canopy seaweeds (Valentine & Edgar, 2010). In addition, we investigated monthly differences in total number of larvae spawned and settled within regions. We also recorded the PLD and distance (km) that larvae travelled before settling.
Tropical and subtropical regions were classified according to the IMCRA 4.0: Provincial Bioregions (Commonwealth of Australia, 2006). Tropical-temperate transition zone was determined as subtropical ( Figure S1), with areas north (>24.5°S) considered tropical and areas south (<30.5°S), temperate. The temperate region originally included a section that combined mainland south-eastern Australia and Tasmania, since the maximum depth within the Bass Strait (between mainland Australia and Tasmania) is mostly shallower than 75 m. These two regions were subsequently split into two using shallower 50 m isobaths in this region, with each settlement region extending roughly halfway into Bass Strait ( Figure S1). These regions remained the same for contemporary and future scenarios.
Larval trajectories were analysed using MATLAB (version R2018a) to quantify spawning, settlement and determine changes in connectivity between sources and sinks.
Changes to potential spawning and settlement within specific regions were analysed using the statistical platform R version 3.4.4.
(R Core Team, 2018). To determine statistical differences between the total number of particles spawned and settled within a specific region, an unpaired two-sample Wilcoxon test was used with contemporary and future scenarios as the arguments and the number of years within each scenario as the replicates. To ensure years were F I G U R E 1 Maps showing proportion of time that temperature falls within Ecklonia's current temperature range on the east coast of Australia at two spatial scales. (a) and (b) represent a contemporary scenario (2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015). (c) and (d) represent a future RCP 8.5 scenario (2090-2100). The subtropical region lies between the dotted lines, north is the tropical, and south is the temperate region. Thin blue line indicates the 75 m isobath and the thick black line encompasses Ecklonia's current and predicted distribution. White indicates 0% of time. RCP, Representative Concentration Pathway independent samples, we corroborated autocorrelations of total spawning and settlement for each region across years. These were not strongly autocorrelated, and therefore, assumptions of independent sample units were met.

| Temperature as a predictor of Ecklonia's future distribution
In the future scenario, Ecklonia's northern distribution experiences a poleward contraction of ~530 km (from 27°S to 31.8°S) based on shifts in its realized thermal niche (Figures 1 and 2). This range contraction is supported by forecasts based on changes in its fundamental thermal niche, that is high temperatures (>26.5°C) above Ecklonia's maximum thermal tolerance threshold ( Figure S2). In most northern areas of Ecklonia's predicted new range, high temperatures are expected to be experienced 0%-20% of the time ( Figure S2c,d), as is presently experienced in Ecklonia's northern distribution ( Figure S2a

| Changes in potential spawning and settlement of Tripneustes
Modelled spawning potential in a contemporary scenario adequately reflected present Tripneustes occurrence. In the contemporary scenario, greatest spawning of Tripneustes occurred in its usual habitat range (the tropics and subtropics). Spawning then gradually decreased towards Sydney, where Tripneustes is still found but in low abundances (Williamson, 2015). No spawning occurred south of 35.8°S, which is close to Merimbula, the furthest point where Tripneustes has been recorded (Williamson, 2015). Hence, temperature-based spawning in the model appears to appropriately reflect Tripneustes' broad present distribution (Figure 3a), suggesting this is a suitable abiotic driver to predict future spawning grounds and range shifts.
Temperate regions experienced a 3.6-fold increase in potential spawning of Tripneustes (Figure 4a; Table S1). In a contemporary scenario, potential spawning within temperate regions typically occurs less than 40% of the time, never exceeds 70%, and no spawning occurs poleward of 35.8°S (Figure 3). In a future scenario, spawning potential in temperate regions increased dramatically, with areas suitable for spawning 50%-100% of the time stretching from 30.5-35.8°S and with the southernmost possible spawning location predicted in Tasmania at 42.5°S (Figure 3). This equates to a ~650 km poleward range expansion in suitable spawning grounds (Figure 2).
Temperatures within the subtropics also became more suited to enable spawning (Figure 3), resulting in an 8.1% increase in future spawning ( Figure 4a; Table S1). Conversely, the model predicts the tropics will become less suitable for the spawning of Tripneustes, which decreases by 43% in the future (Figure 4a; Table S1). In a contemporary scenario, temperatures are suitable for spawning 40%-100% of the time in the tropics, while temperatures are only suitable for spawning 10%-70% of the time in the future (Figure 3).

F I G U R E 2
Maps of change in the distribution of Ecklonia and Tripneustes. The black region indicates the area where Ecklonia will be lost according to changes in temperature. The blue outlines areas where Ecklonia and Tripneustes presently co-exist and will continue to co-occur in the future based on model simulations for (a) potential spawning and (b) potential settlement of Tripneustes larvae. Merimbula denotes the present southernmost recorded location of Tripneustes (Williamson, 2015). The orange indicates areas where Ecklonia and Tripneustes may co-exist in the future only based on model simulations for (a) potential spawning and (b) potential settlement of Tripneustes larvae. The green indicates areas where Ecklonia will exist in the future but will likely not interact with Tripneustes Contemporary settlement extends further south than spawning grounds reaching regions within the continental shelf between mainland Australia and Tasmania ( Figure 5); however, south of 35.2°S, the southernmost spawning latitude, settlement density (i.e. within grid cells of ~10 km 2 ) accounts for less than 0.01% relative to total settlement. In the future scenario, settlement occurs as far south as 43.6°S along the coast off Tasmania and settlement densities larger than 0.01% extend ~400 km poleward to 39.2°S, relative to those in the contemporary scenario. Maximum settlement densities within temperate latitudes off Tasmania increase two orders of magnitude, F I G U R E 3 Maps showing proportion of time that temperature is suitable for spawning of Tripneustes in a contemporary scenario (2006-2015; a) and future RCP 8.5 scenario (2090-2100); b. Blue line indicates the 75 m isobath. The subtropical region lies between the dotted lines, north is the tropical, and south is the temperate region. White indicates 0%. Thick blue line indicates southernmost point where Tripneustes has been recorded to date (Williamson, 2015). RCP, Representative Concentration Pathway Temperate regions experienced the greatest changes to potential settlement, with a 1.5-fold increase in temperate zones (not including Tasmania) and a 677-fold increase in Tasmania, where previously settlement was close to zero ( Figure 4b; Table S1). Settlement in the tropics and subtropics experienced a decrease of 43% and 13.7% with respect to contemporary settlement, respectively (Figure 4b; Table S1). However, potential settlement was still greater in the tropics and subtropics than within temperate regions. LHI experienced no significant change in settlement (Figure 4b; Table S1).

| Main Tripneustes larval sources
Climate change leads to substantial changes on the connectivity between urchin populations in different regions. Under the contemporary scenario, only 21.7% of the larvae that settle in temperate regions are sourced from within temperate latitudes, whereas in the future nearly half of the larvae (45.1%) settling in temperate reefs will be sourced from within that same region (Figure 6a,b).
Additionally, the proportion of larvae released in the subtropics that are predicted to settle in temperate regions also changes, increasing from 14.6% to 22.6% (contemporary and future climate, respectively; Figure 6c,d). LHI experienced the greatest changes to larval sources. In a contemporary scenario, both the tropics and subtropics are significant sources of larvae for LHI urchin populations, with the subtropics contributing most larvae (Figure 6a). In a future scenario, self-recruitment becomes the largest source of larvae to LHI, and almost no larvae are sourced from the tropics (Figure 6b).

| Phenology changes
The most obvious changes to phenology occur in the tropics and subtropics for both potential spawning and settlement (Figures 7   and 8). In the tropics, spawning and settlement in a contemporary scenario are highest across eight months of the year from May to December (late autumn-early summer). In a future scenario, this is halved to only four months of the year and the timing shifts from July to October (mid-winter-early spring). In the subtropics contem-

| Larvae dispersal distance
In both, a contemporary and future scenario most larvae settled within 15 days ( Figure S3a,b) and travelled <400 km ( Figure S3c,d) while potential for dispersal was >3,500 km.

| D ISCUSS I ON
In this study, we have used a combined mechanistic modelling approach to create a predictive map that investigates climate-driven changes in the habitat suitability of a habitat-forming kelp species (Ecklonia), climate-driven changes in the distribution of an ecologically important tropical herbivore (Tripneustes) that account for dispersal processes, and lastly potential changes in the future interaction between these two species caused by new areas of cooccurrence. Overall, we predict that Ecklonia's range will contract poleward substantially due to warming; additionally, we find that climate-driven changes in the oceanic environment shift the spawning and settlement range of Tripneustes, causing its poleward expansion, and an increased area of interaction between this herbivore and Ecklonia with potential consequences on ecosystem health if kelp is overgrazed (Figures 2, 3 and 5). Our case study explores a new method to model species range shifts that incorporates all three aspects of the BAM framework including the anticipation of species interactions in new regions of co-occurrence, a concept that is often overlooked in SDMs. As climate change strengthens WBCs and creates global warming hotspots alongside the eastern coasts of many continents (Cetina-Heredia et al., 2015;Hobday & Pecl, 2014), this approach is relevant and can be applied to predict future shifts and interactions between foundation species and consumers in other temperate ecosystems globally aiding adaptive management.

F I G U R E 6
Connectivity matrices showing source of Tripneustes larvae that settle/sink in a region. Panels on the left represent a contemporary scenario (2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015). Panels on the right represent a future RCP 8.5 scenario (2090-2100). (a) and (b) indicate where the % of larvae that settle/sink in a region were originally sourced from, with values totalling 100% in the vertical direction. (c) and (d) indicate where the % of larvae released from a specific source/region then settle/sink, with values totalling 100% in the horizontal direction. LHI, Lord Howe Island; NaN, no spawning; SUB, subtropics; TAS, Tasmania; TMP, temperate region; TRP, tropics; RCP, Representative Concentration Pathway

| Future loss of kelp caused by warming and a range-shifting tropical herbivore
Our findings predict a smaller range contraction for Ecklonia based on temperature changes alone (Figures 1 and 2) than a previous study by Martínez et al. (2018), which used an SDM approach through MAXENT to also model the future distribution of Ecklonia. Their Furthermore, future conditions used in MAXENT were obtained from Bio-ORACLE, which provides ocean data produced from averaging three CMIP5 climate models (Assis et al., 2017). In contrast, OFAM3 future projection (Zhang et al., 2016) uses an ensemble of 17 CMIP5 models as forcing. Our approach provides a more optimistic outlook for kelp distribution in the future where temperature alone may not cause the severe retraction of habitat that was estimated by Martínez et al., (2018). However, further losses to those we predicted may be expected with the increased frequency of marine heatwaves (Frolicher, Fischer, & Gruber, 2018). Our predictions of Ecklonia's distribution are based on mean temperatures, possibly underestimating potential severe decimation of Ecklonia populations due to isolated marine heatwave events, which were not explored in this study . zones where this kelp is declining (Vergés et al., 2016). Hence, an increased abundance of Tripneustes could facilitate a greater contraction of Ecklonia's forecasted range. Further, Tripneustes is also able to tolerate higher temperatures (22-29°C) (Rahman et al., 2009) than Ecklonia. Therefore, as extreme heat events are expected to increase in frequency (Frolicher et al., 2018), Tripneustes may pose an even greater grazing risk to Ecklonia during these periods of thermal stress. Evidence for a range expansion of Tripneustes is already emerging in Western Australia, where increases in abundance have been observed following a marine heatwave event (Smale, Wernberg, & Vanderklift, 2017). Furthermore, phenological changes in the timing of spawning and settlement peaks for Tripneustes (Figures 7 and 8) may have unknown consequences for other species that interact with Tripneustes, such as creating a mismatch between predators and the availability of larvae as food (Durant, Hjermann, Ottersen, & Stenseth, 2007). If a mismatch occurred in favour of the larvae, this could cause Tripneustes outbreaks. Tripneustes is forecasted to spread into Tasmania where kelp forests are already being decimated by range-shifting Centrostephanus rodgersii (Ling, Johnson, Frusher, & Ridgway, 2009), novel interactions between kelp and Tripneustes in these areas where previously they did not co-occur (Figure 2, orange area) could lead to a significantly greater loss of kelp than was predicted by using optimal temperature range in our model (Figure 2, black area). Future management techniques, such as ensuring Tripneustes predators are not overfished (Eklöf et al., 2008;Ling et al., 2009) or potentially harvesting urchins as a new commercial fishery (Scheibling, Hennigar, & Balch, 1999) would aid in decreasing the risk of population outbreaks in the future.

| Drivers of change in Tripneustes settlement patterns
Within the aspects considered in the BAM framework, warming (abiotic factor) rather than changes in larval pathways (movement) appeared to be a greater driver of poleward range expansion for Tripneustes. This may also be the case for species with short PLD for This was the case in our simulations, where maximum settlement occurred in the first 15 days ( Figure S3a,b). Future research may consider separating the effects of a strengthened EAC from warming, to determine precisely the dispersal mechanisms for increased poleward range shift of Tripneustes into temperate ecosystems.
Conversely, in the tropics potential spawning and settlement is predicted to decrease due to warming temperatures reducing suitability as a year-round spawning location (Figures 3-5 and 7). This could contribute to reduced health of tropical coral reef ecosystems in the future, as sea urchins are considered to have an important ecological role in maintaining the balance of coral-algae competition on coral reefs (Coyer, Ambrose, Engle, & Carroll, 1993). Removal of important herbivores such as sea urchins can cause overgrowth of algae and facilitate ecosystem phase shifts from coral to algae dominated reefs (Coyer et al., 1993;Hughes, Reed, & Boyle, 1987). Consequently, a reduction in Tripneustes urchins may have negative consequences to tropical coral reefs, particularly if other important herbivores like fish also shift their ranges poleward (Vergés, Steinberg, et al., 2014).

| Connectivity implications for LHI
The results from LHI have particular conservation relevance due to the unique status of this island as a UNESCO World Heritage site (Environment Australia, 2002) and given records of past outbreaks of Tripneustes in this region (Valentine & Edgar, 2010). They also highlight the importance of considering dispersal when predicting species range shifts. Even though we did find that warming increased the suitability of LHI as a year-round spawning site (Figures 3 and 4a), we found no significant difference in settlement ( Figure 4b) and there appeared to be less variability in monthly settlement rates (Figure 8), suggesting no particular greater risk of outbreaks. However, we discovered a significant change in source regions supplying Tripneustes larvae to LHI (Figure 6a,b). In the future scenario, most larvae are self-recruited from LHI itself, the proportion of larvae exported from temperate regions remains unchanged and small, and there is a drastic reduction in supply from the subtropics and tropics (only significant sources in the contemporary scenario). These changes are most likely driven by alterations to ocean circulation; as the EAC strengthens with increased flow in a poleward direction, eastward flow towards New Zealand diminishes (Oliver & Holbrook, 2014), which may reduce larval transport from mainland Australia to South Pacific Islands such as LHI of all species with a planktonic phase. Furthermore, the EAC is projected to separate further south (Oliver & Holbrook, 2014); therefore, the EAC eastward flowing extension may also extend poleward and bypass LHI. These changes imply that LHI may become more isolated from tropical/subtropical Australia in the future, which could alter ecosystem composition and reduce population resilience through gene flow changes decreasing genetic diversity (Sgrò, Lowe, & Hoffmann, 2011), as species rely more on self-recruitment to maintain populations. The changes in connectivity patterns revealed here are likely to influence many other species and highlight the benefits of using a Lagrangian approach that considers dispersal mechanisms.

| Limitations to the dispersal model
We included the most current research on the biological properties of Tripneustes into our model, however, there are still knowledge gaps that limited our dispersal simulations. For example, in our model Tripneustes spawning was based solely on temperature determined by laboratory experiments in controlled environments (Mos et al., 2011;Rahman et al., 2009), as there is currently no certainty of spawning cues for wild populations. Furthermore, ocean acidification may have important impacts on future populations of calcifying organisms like Tripneustes, reducing reproductive output (Dworjanyn & Byrne, 2018) and larval success (Sheppard Brennand et al., 2010). Ocean pH was not projected with the downscaled ocean climate model; therefore, we could not incorporate these effects. Further research into these areas would improve representation of Tripneustes larvae and are factors to consider if using this approach for other species.
Additionally, the very recent discovery of a new species of sea urchin Tripneustes kermadecensis, which overlaps in range and has long been confused with T. gratilla, may mean that some of the biological properties used in our model may have potentially come from studies that unknowingly used T. kermadecensis or hybrids of the two species in their research. Bronstein et al. (2019) suggest that populations poleward of Sydney are likely to belong to T. kermadecensis. Given that frequency of suitable spawning temperature in our models is only high north of Sydney ( Figure 3) and this coincides with the southernmost distribution of T. gratilla inferred from genetic analysis (Bronstein et al., 2019), we are fairly confident that T. gratilla's distribution was well represented in our model. Nonetheless, findings should be interpreted with some caution until further genetic sampling is undertaken to determine with certainty the full range of the two species.

| Concluding remarks
This study provides a novel combined modelling approach that accounts for all three aspects of the BAM framework (Peterson et al., 2015) to address the limitations that other SDMs and mechanistic models encounter when modelling species distributions. The model uses output from OFAM3 to map the contemporary and future realized and fundamental thermal niches of a dominant temperate kelp species, and a Lagrangian particle-tracking framework to forecast and map range shifts of a tropical sea urchin due to climate-driven changes in dispersal. The output from these two approaches is then combined into one predictive map, which shows potential for strengthened interactions between the species in the future due to new regions of co-occurrence, possibly leading to greater kelp loss through increased grazing pressure. We were also able to inadvertently discover changes to connectivity of Tripneustes between source-sink locations. Our approach of niche characterization, dispersal modelling and distributional prediction is transferable to other marine ecosystems and can be used to determine range shifts and forecast changes to species interactions between other foundation species and consumers. It can assist managers worldwide to predict future ecosystem composition and important species interactions, and potentially prepare adaptive measures to better manage changing marine ecosystems.

ACK N OWLED G EM ENTS
CSIRO's Ocean Downscaling Strategic Project provided the OFAM3 model output. P.C.-H. and this work were partially supported by Australian Research Council Linkage grants LP 150100064 and LP16100162 to MR. This work was also supported by the NSW Environmental Trust and by Australian Research Council Discovery grants ARC DP170100023 and DP190102030 to AV. We thank Jason Everett for his feedback and two anonymous reviewers whose comments helped improve the manuscript.

DATA AVA I L A B I L I T Y S TAT E M E N T
The larval trajectories generated for this study are available through