Climate and environmental changes driving idiosyncratic shifts in the distribution of tropical and temperate worm reefs

Abstract

An increasing number of studies have forecast the potential responses of marine life to future climate change. This study predicts how the distributional range of temperate and tropical worm reefs (WRs) might respond to climate and environmental changes (CECs). Compared with current distributions, the tested hypotheses were: (i) under a low CO₂ concentration and active atmospheric carbon capturing scenario (RCP2.6), both tropical and temperate WRs will maintain their current distributions and face only slight multi-directional biogeographic changes along the century; and (ii) under a high CO₂ concentration scenario (RCP8.5) WRs will shift toward higher latitudes, with marked changes for tropical species and slight changes for temperate species, specifically at the end of the 21st century. The hypotheses were tested using species distribution modelling, and exploratory statistical analyses were performed to tune model settings. Under scenario RCP2.6, in the middle of the century, areas of suitable habitat are predicted to slightly increase for the temperate WRs and conversely contract for tropical WRs. At the end of the century, multi-directional shifts without range retraction were predicted for both species, but tropical WRs showed major changes in their distribution. Under scenario RCP8.5 and throughout the century, multi-directional shifts increased the areas of suitable habitat for temperate WRs, whereas tropical WRs experienced shifts toward high latitudes and significant retraction at low latitudes. Results indicate that biogeographic range shifts are idiosyncratic for temperate and tropical WRs depending on the CECs scenarios considered.

INTRODUCTION

Climate and environmental changes (CECs) have synergistic impacts on marine ecosystems (e.g. altering temperature, ocean circulation, nutrient concentrations and seawater pH) and can drive spatio-temporal changes in the distribution and abundance of species, with cascading effects on ecosystem functioning (Harley et al., 2006; Hoegh-Guldberg & Bruno, 2010; Burrows et al., 2011; Gattuso et al., 2015). The expected sign of CECs on marine life reaches from genetic and physiological adaptations to distributional shifts (Harley et al., 2006; Poloczanska et al., 2013; Gotelli & Stanton-Geddes, 2015). Under climate and environmental changes, the general prediction for marine and terrestrial species distribution, based on observations and theories, is a shift towards high latitudes, where more favourable conditions might exist, or alternatively, multi-directional distribution shifts are expected to occur (Burrows et al., 2011; VanDerWal et al., 2013; García Molinos et al., 2015; Lenoir & Svenning, 2015). Marine organisms are valuable for assessing the possible impacts and magnitude of CECs on species distribution (Monaco & Helmuth, 2011; Riebesell & Gattuso, 2015); one such example is ecosystem engineers, organisms that create and provide persistent habitats for the associated biodiversity (Jones et al., 1997; Coleman & Williams, 2002). Understanding the responses of such species under CEC scenarios has become a challenge for marine ecologists. Moreover, it is crucial in forecasting possible shifts in the distributional range of ecosystem engineers under CEC (Przeslawski et al., 2008).

A number of approaches (e.g. laboratory and field studies, long-term multigenerational experimental research, and modelling) have been used to assess the impacts of CECs on marine biodiversity (Wernberg et al., 2012; Andersson et al., 2015; Riebesell & Gattuso, 2015). Numerical model simulations of species distribution generate crucial predictive models that can be useful for conservation planning (Andersson et al., 2015). There is strong evidence that recent observations of sea urchin and fish distribution shifts support modelling forecasts (Poloczanska et al., 2016). Although there are some examples of modelling studies using multiple drivers and multiple scenarios showing shifts in marine biodiversity in the global scale (Jones & Cheung, 2015; García Molinos et al., 2015) or considering ecosystem engineer organisms (Couce et al., 2013), there have been few attempts to forecast climate-driven impacts on tropical and temperate species on a large geographic scale (Saupe et al., 2014). Previous modelling studies have indicated that the expected effects on the distribution of marine organisms are as follows: (i) species will exhibit poleward shifts in their distribution, (ii) tropical species will be highly vulnerable due to the contraction of suitable habitats, and (iii) the geographic distribution of suitable habitats will decrease more quickly and markedly under a scenario with high CO₂ concentrations (García Molinos et al., 2015; Jones & Cheung, 2015). Alternatively, within groups of solely tropical or temperate species, no consistent pattern among species would occur due to varied responses in the direction of shifts or an increase/decrease in the suitable area (Saupe et al., 2014). Although in recent years a large number of modelling studies have shown distributional shifts under CECs, this approach is not yet common to assess the responses in closely related ecosystem engineer species.

In temperate and tropical coastal areas around the world, some species of sabellariids build so-called worm reefs (WRs) (Zamorano et al., 1995; Dubois et al., 2002; Bremec et al., 2014; Faroni-Perez, 2014; Firth et al., 2015; Nishi et al., 2015; Faroni-Perez et al., 2016), which are important biogenic structures, following hermatypic corals and stromatolites (Goldberg, 2013). WRs represent habitats of high biodiversity on which hundreds of species of invertebrates and fishes rely for feeding, nursery, refuge and migration (Gherardi & Cassidy, 1994; Dubois et al., 2002; Ataide et al., 2014). There is strong evidence, based on historical and recent surveys that the distribution and abundance of WRs are already changing in response to recent CECs (Firth et al., 2015). In the present study, multiple drivers and scenarios were used to target future biogeographic changes in the distribution of two species of polychaetes of the Sabellariidae family that are ecosystem engineers (Coleman & Williams, 2002). Among the sabellariids, the genus Phragmatopoma is endemic to the American coastline, contains four valid species and is mainly registered in intertidal and shallow subtidal zones (Capa & Hutchings, 2014; Faroni-Perez et al., 2016). Here, the focus was on two reef-building species, Phragmatopoma caudata, Krøyer in Mörch, 1863, and P. virgini, Kinberg, 1866, as they are genetically close related and occur in distinctive areas (Nunes et al., 2016). In terms of marine biogeographic provinces (Briggs & Bowen, 2012), P. caudata occurs in the tropical regions of the Western Atlantic, ranging from the Florida coast to the southern coast of Brazil, whereas P. virgini occurs in the temperate regions of the South-east Pacific Ocean ranging from Ecuador to Patagonia and northwards up to Uruguay in the South-west Atlantic Ocean.

Therefore, in order to assess if CECs will cause distribution shifts in tropical and temperate WRs, present-day habitat suitability distributions were modelled, as well as under IPCC-CMIP 5 projections for two different time intervals (mid-century: 2040–2059 and late-century: 2080–2099). The tested hypotheses was as follows: (i) under a low CO₂ concentration and active atmospheric carbon capturing scenario (RCP2.6), both tropical and temperate WR species will maintain their current distributions and face only slight multi-directional biogeographic changes across the century; and (ii) under a high CO₂ concentration scenario (RCP8.5), both species will shift toward high latitudes, with marked changes for tropical WRs and slight changes for temperate WRs, specifically at the end of the century.

MATERIALS AND METHODS

Datasets: species occurrences and environmental predictors

For the dataset of species occurrences, an extensive online search was performed using the keyword ‘Phragmatopoma’ in various reference databases (i.e. Web of Science, SCOPUS, SciElo, Google Scholar), the catalogue of the Smithsonian National Museum of Natural History, and the NONATOBase (Pagliosa et al., 2014). Additionally, original data on 57 new records was added to the dataset. Each occurrence data was georeferenced (WGS84) and record sites with low accuracy were discarded to reduce noise in the analysis. Overall, 97 literature sources were included in the analysis, resulting in 284 and 94 occurrence records for P. caudata and P. virgini, respectively.

The primary datasets of environmental drivers used are from the Coupled Model Intercomparison Project Phase 5 (IPCC, 2013; http://cmip-pcmdi.llnl.gov/cmip5/). The model selected was the HadGEM2-ES, from Met Office Hadley Centre, UK (Jones et al., 2011; Martin et al., 2011). The global climate model provided by HadGEM2-ES includes all components of the earth system (i.e. physical atmosphere and ocean components with the addition of schemes to characterize aspects of the Earth system) (Collins et al., 2011; Martin et al., 2011; Zelinka et al., 2013), and have been used to forecast shifts in distribution of marine and terrestrial species (Good et al., 2013; Saupe et al., 2014; García Molinos et al., 2015; McQuillan & Rice, 2015; Xavier et al., 2016). The monthly output data of the concentration of carbon, chlorophyll, dissolved nitrate, ammonium and silicate, hydrogen ion, alkalinity, air and sea surface temperature, precipitation, salinity and water runoff was acquired from HadGEM2-ES (Table 1). The dataset was divided into a ‘current’ period, which included the recent past (1986–2005), whereas future projections for the Representative Concentration Pathways (RCPs) 2.6 and 8.5 were subdivided into the ‘middle of the century’ (i.e. mid-century 2040–2059) and the ‘end of the century’ (i.e. late-century 2080–2099). These future projections represent contrasting scenarios: the RCP2.6 characterizes the optimistic and environmentally friendly future with low CO₂ concentrations, and the RCP8.5 represents the worst-case scenario with the highest CO₂ concentrations. Each environmental predictor was rescaled to a resolution of 0.08 × 0.08°, using bi-linear interpolation. Additionally, for each predictor, the interannual monthly mean was calculated from the primary data. The annual values (average, minimum and maximum) were then calculated totalling 36 environmental drivers layers (see online supplementary material).

Table 1. Niches and their environmental drivers with code and units used for modelling habitat suitability.

Selection of predictors

The predictors, here classified according to the use of the source by the worms (i.e. biogenic, osmotic, thermal and trophic; see Table 1), were selected based on environmental variables representative of the key dimensions of the niche for the studied species (Austin, 2007). This multi-dimensional niche approach was used as a step selection of variables that can influence the distribution of Phragmatopoma species. Following a hierarchical framework to reduce multi-collinearity, an array of a peer-to-peer correlation matrices using the set of layers from the ‘current’ time period was created that excluded highly correlated predictors (Pearson, r > 0.60). After this step, only nine out of initially 36 predictors remained, two for each niche dimension, except the trophic, which had three predictors per species. Subsequently, to select for each species the most representative predictor for each pre-established niche dimension, a set of generalized linear models (GLMs) was built, exploring the possible combinations of variables and considering one predictor per niche. The final sets of niche predictors per species were assessed by the most influential predictive capacity across models. The predictive capacity was ranked based on the number of times a given predictor was significant across GLM models (see online supplementary material). Final layers of predictors selected are as follows: pHmin (biogenic) and tosmin (thermal) for both the species, and frivermax (osmotic), simax (trophic) selected for P. caudata, and prmax (osmotic) and detocmax (trophic), for P. virgini. Predictors’ names and units are described in Table 1. Statistical analysis of these layers can be found in the online supplementary material.

Modelling

Maximum entropy modelling (MaxEnt) was used to model the current distributions of P. caudata and P. virgini and then project their future habitat suitability according to the future IPCC-CMIP 5 projections. MaxEnt is a correlative algorithm that uses the occurrence data of species coupled to environmental data for modelling species distributions (Phillips et al., 2006; Phillips & Dudik, 2008; Elith et al., 2011) and can be used to predict the future impacts of CECs on marine organisms (Couce et al., 2013; Verbruggen et al., 2013; Saupe et al., 2014). Generally, the models produced by MaxEnt are influenced by the extent of the study area and species occurrence, which are data input premises. Thus, the MaxEnt algorithm performs better when the study area for the calibration model does not include areas outside the species range occurrence (Anderson & Raza, 2010; Elith et al., 2011). For this reason, the model is calibrated with the immediately surrounding region for each species known occurrence. Subsequently, the background area of each species was referred strictly to the occurring latitudinal span and longitudinally delimited by the closest cell grid representing the continental shoreline, as the two species inhabit the intertidal or shallow subtidal zones. Additionally, the OccurrenceThinner (Verbruggen et al., 2013) was used to filter the data and reduce the geographic sampling bias since MaxEnt better estimates the contribution of the niche environmental variables for each species, when not overestimating the model performance (Elith et al., 2011; Veloz et al., 2013). OccurrenceThinner uses the occurrence records and a raster of a kernel density function of the study area to reduce occurrence records, omitting more records from densely sampled regions. This filtering resulted in 42 and 21 occurrence records for P. caudata and P. virgini, respectively.

For the species-specific model tuning, various potential sources of complexity from the MaxEnt features were investigated, considering methods proposed by Elith et al. (2010, 2011), Shcheglovitova & Anderson (2013) and Radosavljevic & Anderson (2014) (see online supplementary material). Then, the fitted model was assessed based on the performance of explanatory predictors and MaxEnt features. The maps and results reported herein are based on the MaxEnt average (5 runs), cross-validated, logistic output, and with the clamping setting as true. The model predictions were tested by partial receiver operating characteristic analyses (p-ROC AUC) with E = 0.5 of the maximum error rate, 50% random points and 500 iterations for the bootstrap (Peterson et al., 2008; Peterson, 2012).

Post modelling

To quantify the shifts in the distribution of tropical and temperate WRs, a metric based on total grid cell values (i.e. each pixel suitability) was applied. The total numbers of grid cells were 4515 for P. caudata and 3322 for P. virgini. For each projection, grid-cell values were extracted and frequency distribution and the delta value were computed. The delta value is the relative difference in total grid-cell frequencies between each future projection and the present day. Although the threshold of minimum training presence values was 0.11 for P. caudata and 0.09 for P. virgini, the threshold value of 0.2 was set to define unsuitable areas, pixel values ≤0.2 are included in the non-suitability category (see online supplementary material). Consequently, changes in delta values correspond to a gain or a loss of suitable area. The remaining grid cells (i.e. pixel values: comprised between 0.21 and 1) correspond to the suitable area, and for which changes in the delta value were categorized as follows: 0.21 ≤ low ≤ 0.4; 0.41 ≤ medium ≤ 0.6; 0.61 ≤ high ≤ 0.8; 0.81 ≤ very high ≤ 1. This analysis attempted solely to characterize and quantify ‘potential’ shifts in the distribution of each species based on their areas of environmental suitability.

Processing all modelling

All software used in this study is listed in the online supplementary material.

Abbreviations

CECs:

Climate and environmental changes.

WRs:

Worm reefs built by Sabellariidae.

RCPs:

Representative concentration pathways; radiative forcing values for scenario set, containing emission, concentration and land-use trajectories (van Vuuren et al., 2011).

RCP2.6:

Mitigation scenario leading to a very low forcing level.

RCP8.5:

Very high baseline emission scenario.

RESULTS

Modelling analysis and contribution of multiple drivers

MaxEnt models showed good predictive performances for the models (Table 2). For the present-day period, the predicted areas of suitable habitat were in agreement with the current distribution of both species. The current areas of suitability cover 44.7% and 60.8% of the total extent of Phragmatopoma caudata and P. virgini, respectively. For both species, the predictors regarding the trophic and osmotic niches showed an intermediary influence in the model, although they had an inverse influential order (Table 3). The most influential predictor for both species was sea surface temperature (tos), related to the thermal niche. In contrast, pH (related to the biogenic niche) showed the lowest contribution to the model for both species (Table 3).

Table 2. Summary of model statistics showing the receiver operating characteristic curves for MaxEnt with E = 0.5 of acceptable omission error according to p-ROC AUC values.

Table 3. Influence of environmental drivers used in modelling habitat suitability. Codes used for predictors are defined in Table 1.

Future trends of habitat suitability for worm reefs

All models predicted future differences in the distribution of habitat suitability for tropical and temperate WRs (Figure 1). For the ‘optimistic scenario’ (RCP2.6), in contrast with the current distributions of both WRs, a gain of suitable areas was predicted for the end of the century. However, idiosyncratic responses were forecast; in the mid-century, tropical WRs would lose area, while temperate WRs would gain area. At the end of the century, tropical WRs would markedly increase their distribution, while temperate WRs would slightly lose area compared with the mid-century situation (Figure 1). Under the ‘business as usual’ scenario (RCP8.5), which predicts the highest temperature increase, models also showed idiosyncratic responses of tropical and temperate WRs (Figure 1). Tropical WRs tend to lose suitable areas across the century, and a majority of changes would occur up to the mid-century. In the 50s, it is expected that non-suitable areas increase by 13.84% for a total of 16.48% in the 90s (Table 4). The temperate WRs tend to gain suitable areas throughout the century and by the end of the century, all biogeographic areas would be suitable (Table 4). However, the major gain would occur up to the mid-century.

Fig. 1. Predicted changes in total suitable areas of the tropical Phragmatopoma caudata and the temperate Phragmatopoma virgini worm reefs considering current and future predictions for RCP2.6 (low emissions) and RCP8.5 (high emissions) scenarios, and the time-scales of the 50s (mid-century: 2040–2059) and the 90s (late-century: 2080–2099).

Table 4. Projected changes in the potentially suitable areas for topical (Phragmatopoma caudata) and temperate (P. virgini) worm reefs with contrast to their current habitat suitability. Δ represents the difference in frequencies of pixel suitability between current and scenarios; RCP2.6 (low emissions) and RCP8.5 (high emissions) and between the 50s (mid-century: 2040–2059) and the 90s (late-century: 2080–2099) time period. In parentheses, the total frequencies (F%) given in percentage of grid cells per projection.

BIOGEOGRAPHIC SHIFTS

Tropical worm reefs

In contrast to the current distribution, the RCP2.6 model shows opposite shifts during the century for P. caudata (Figures 1 & 2A–C). In the middle of the century, the model projected a reduction of the suitable area (~8%) in the Caribbean region, Venezuela and the northern Gulf of Mexico due to the loss of grid cells with low and medium classes of suitability. On the other hand, increments in high (~8%) and very high (~3%) classes of habitat suitability were found along the east coast of Brazil and the western part of the Gulf of Mexico (Figure 2B). Interestingly, for the end of the century, the model predicted more advantageous conditions with increasing areas and potential suitability, demonstrating the potential of species range expansion. With respect to the total area predicted as suitable, the models showed an increase from ~45% to ~69%. It is noteworthy that in the late century, 24% out of these 69% suitable areas are expected to reach the classes of high and very high suitability (Table 4). In general, areas largely in the Caribbean region, the West Indies, the western coast of the Gulf of Mexico and the coast of Brazil are herein projected as suitable (Figure 2C), with only small gaps of non-suitable areas in northern South America and the northern part of the Gulf of Mexico.

Fig. 2. Maps* show spatial variation in current (A) and future (B–E) potential suitable area distributions for tropical worm reefs (Phragmatopoma caudata) predicted given the RCP2.6 (low-emission) and RCP8.5 (high-emission) scenarios. Each projection (i.e. current, 50s and 90s) is based on a dataset covering two decades. Pixel colours correspond to the following: 0–0.2 = area of non-suitability (cold colours: medium to dark blue), 0.4 ≤ low ≤ 0.2 1(colours: light blue to green); 0.41 ≤ medium ≤ 0.4 (colours: green to light green); 0.61 ≤ high ≤ 0.8 (colours: light green to yellow); 0.81 ≤ very high ≤ 1 (warm colours: orange to red). *Continental shoreline of the tropical regions of the Western Atlantic, ranging from Florida coast (USA) to the southern coast of Brazil.

Under the RCP8.5 model, a considerable contraction was observed of the total area of potential distribution from the current (~45%) to the middle (~31%) and the end of the century (~28%; Figure 2A, D, E; Table 4). In general, an increase of non-suitable areas is projected for the entire northern part of South America, the Caribbean region, the West Indies, and the south-western and northern parts of the Gulf of Mexico. Furthermore, of the total remaining suitable area, a significant amount is classified into the low and medium classes of suitability (mid-century: ~21%, late-century: ~14%; Table 4), demonstrating vulnerability due to a massive loss of suitability. Finally, the remaining suitable areas with high and very high classes of suitability (mid-century: ~10%, late-century: ~14%) are projected for the north-western region of the Gulf of Mexico and the south-eastern part of Brazil (Figure 2D, E).

Temperate worm reefs

Under the RCP2.6 modelling outputs, Phragmatopoma virgini showed a gradual trend of a spreading distribution with increasing values of the potential suitability of habitat over the century (Figure 3A–C; Table 4). It is noteworthy that in contrast to present-day conditions, the suitability increases by 17.25% by the mid-century, reaching 78% of potential habitat suitability over the biogeographic province. Overall, projected areas that gain suitability are found towards the north of Peru and South Patagonia (Argentina and Chile), demonstrating multi-directional distribution shifts. In contrast to the mid-century, a slight increment (~2%) of the area of suitability for the end of the century was found. Furthermore, increasing values in the potential suitability for the low (~13–10%) and high (~5%) classes was observed throughout the century, regardless of a slight loss (1.6–0.75%) of the medium, indicating that P. virgini might retain its distributional range under this optimistic scenario. The biogeographic areas with increases in potential suitability were Ecuador, along the Peru coastline, and in Chile, mainly for the central part at Regions III–VII (i.e. Atacama, Coquimbo, Valparaíso, O'Higgins and Maule).

Fig. 3. Maps* show spatial variation in current (A) and future (B–E) potential suitable area distributions for temperate worm reefs (Phragmatopoma virgini) predicted given the RCP2.6 (low emissions) and RCP8.5 (high emissions) scenarios. Each projection (i.e. current, 50s and 90s) is based on a dataset covering two decades. Pixel colours correspond to the following: 0–0.2 = area of non-suitability (cold colours: medium to dark blue), 0.4 ≤ low ≤ 0.2 1(colours: light blue to green); 0.41 ≤ medium ≤ 0.4 (colours: green to light green); 0.61 ≤ high ≤ 0.8 (colours: light green to yellow); 0.81 ≤ very high ≤ 1 (warm colours: orange to red). *Continental shoreline of the temperate regions ranging from Ecuador in South-east Pacific Ocean and Uruguay in South-west Atlantic Ocean southward to Patagonia.

Under the RCP8.5 model, a gain of suitable areas during the century was also predicted (Figure 3A, D, E; Table 4). The projections show that by mid-century, 83.2% of the biogeographic province will be suitable for P. virgini. In this period, increases in the classes of low (~16%) and high (~10%) suitability were projected, notably demonstrating species gain in suitable distribution. The potential suitability is expected to increase along the coast of Ecuador, Peru and Chile (i.e. north and central coast) in the Pacific Ocean and in Uruguay and Argentina (i.e. Buenos Aires towards south of Chubut) in the Atlantic Ocean. In general, only small gaps of remaining non-suitable area in Southern Chile (i.e. Regions XI Aisén and XII Magallanes and Antártica Chilena) and Argentina (i.e. Tierra del Fuego and Santa Cruz) were observed. Finally, regarding the projection for the end of century, no grid cell is projected to be a non-suitable area, indicating that all of the P. virgini biogeographic provinces will be potentially suitable. Coupled to this evident gain of suitable area, the increments were projected into the medium (~10%) and high (~20%) classes of suitability.

DISCUSSION

Biogeographic shifts of temperate and tropical worm reefs are idiosyncratic to the climate and environmental change scenarios considered. Overall, models showed an increase in suitable habitats for temperate WRs, with a varying magnitude depending on scenarios and time periods. For the tropical WRs, responses were generally more antagonistic, also depending on the scenario and period considered. This study was the first attempt at modelling the future distribution of polychaetes using multiple drivers under IPCC-CMIP 5 scenarios. Even though the values from different modelling approaches should not be comparable (Guillera-Arroita et al., 2015), projections from a previous study using multiple drivers under IPCC-CMIP 5 scenarios have also shown that among marine invertebrates, idiosyncratic responses to CECs are expected to occur (Saupe et al., 2014).

The assumptions for temperate and tropical WRs under RCP2.6 and RCP8.5 scenarios did apply only in part. To clarify, under RCP2.6, multi-directional changes in habitat suitability would occur, but projections showed that at the end of the century, the magnitude of changes is either slight or marked for the temperate and tropical WRs, respectively. Additionally, in the case of scenario RCP8.5, marked changes in the distribution pattern were observed across the century for both species, but a consistent southward shift with suitable retraction at low latitudes was found for tropical WRs only, whereas temperate WRs exhibited continuous multi-directional shifts. Therefore, the present results suggest that the common assumptions of higher latitude range shifts cannot be globalized, as the level of CEC conditions (i.e. the different scenarios and regional impacts) and species tipping points (i.e. tropical, temperate or pandemic) also need to be considered.

Environmental variables and insights into scales

In the present models, drivers’ contribution rates differed between congeners, corroborating previous studies, which showed that CECs would not equally impact distinct biogeographic regions (Brierley & Kingsford, 2009), and that temperate and tropical species would have different responses. Tropical species are expected to be more vulnerable than temperate ones (Rosa et al., 2014). Although demonstrating different patterns for either reef builder, sea surface temperature was the most important predictor, whereas pH played a minor role in explaining the distribution of suitable areas. This result is congruent with previous modelling studies conducted for bio-engineering corals (Couce et al., 2013). However, it is possible that the model could have underestimated the role of pH for WRs habitat suitability. In the absence of actions to reduce CO₂ emissions, within the next few decades, seawater pH values will represent a negative impact for marine invertebrates (Orr et al., 2005; Shirayama & Thornton, 2005). Furthermore, the underwater solidification of cement used by sabellariids to build their tubes, and, consequently, the reefs, appears to be initially dependent on pH differences between the worms’ secretory glands and seawater (Stevens et al., 2007; Wang & Stewart, 2012). Further experiments are required to better assess the vulnerability of WRs to pH under future ocean acidification scenarios and shed light on the accuracy of models to predict the impact of low pH on shifts of species distribution. The projected increase of CO₂ into the atmosphere and subsequent environmental changes throughout this century will enhance either risk: the vulnerability of tropical WRs at low latitudes and the spread of temperate WRs at high latitudes, both of which will have cascading effects on ecosystem functioning. Although warming seawater by itself would not be the unique driver to constraining biodiversity (Gattuso et al., 2015), the thermal predictor had the highest contribution to models.

Species distribution models and the method herein used consider that occupied grid cells estimate suitable environmental conditions, and the set of abiotic traits in these occupied cells is more favourable than adjacent unoccupied cells and non-similar abiotic traits. The geographic scale dealt with in this study encompasses a large part of the American continent. Both worm species modelled have a wide-range distribution at the continental scale, indicating that differences in local and micro-habitat variables have a relatively low importance in understanding and predicting species distribution at larger spatial scales. Spatially, the predictor layers used in this study provide sensible results at continental scales that are not necessarily accurate to local scale. Care should be applied to considerations about suitable pixel and grid cell size due to some peculiarities at fine scale, which might not be favourable for these species considering biotic (e.g. tipping points and intrinsic settlement ability) and abiotic factors (e.g. substratum availability). For example, in a grid cell a river runoff can drop coastal marine salinity to unfavourable values. In fact several WRs occur in coastal areas adjacent to estuaries (see in Faroni-Perez, 2014; Aviz et al., 2017; Eeo et al., 2017). During low tides WRs can cope with rain. Therefore, WRs often undergo salinity fluctuation (i.e. euryhaline to mesohaline), but no empirical data indicate that continually low salinity values (i.e. mesohaline to oligohaline) are favourable. The species modelled have wide-range distributions, however, at local spatial scales WRs do not have a homogeneous distribution. For example, low recruitment ratios in P. caudata were associated with the passage of cold fronts that dispersed larvae along coastline (Faroni-Perez, 2014). Although larvae can drift to different locations, they have preferences to settle gregariously onto specific substrata, a behaviour associated with sensory organs (Faroni-Perez et al., 2016). Larvae settle preferentially on existing reef or hard substrata, and substratum trait (type, density, over area, heterogeneity) might exert significant influence on the models. However, marine species distribution modelling lacks high resolution intertidal geological layers encompassing the American continents and hence this predictor was not considered in the present study. In some places around the world sabellariid reefs can be found directly on soft bottoms (Faroni-Perez et al., 2016), but that is not common among species of the genus Phragmatopoma.

Present results indicated that distribution shifts in WRs would occur. Studies have demonstrated that cost-effective methodological approaches, such as scaling up (Ferrari et al., 2016; González-Rivero et al., 2016), can be applied to target regional areas (e.g. 10–10,000 km) toward acquiring empirical data on distribution shifts. Further studies are encouraged to use the scaling up approaches for monitoring WRs and their associated biodiversity by a well-established temporal assessment as a key measure to detect distribution shifts as response to short-term CEC events, such as heatwaves or cold spells. Global impacts of CECs on marine biodiversity are inevitable (Harley et al., 2006; Poloczanska et al., 2014; Gattuso et al., 2015), and identifying them as a matter of urgency would promote better strategies for management actions attempting to reduce the adverse effects for biodiversity conservation.

Tropical and temperate worm reefs under the RCP2.6

Under the optimistic scenario, projections of suitable areas will persist for congeners at the end of the century, indicating the benefits for WRs to guarantee range maintenance. Similar results of suitable area persistence were found for the gastropods Bulla occidentalis and Conus spurius (see online Table S2.1, MaxEnt, all variables in Saupe et al., 2014), and interestingly, the range of these molluscs co-occurs with tropical P. caudata. Across the century, temperate reefs will continually gain suitable areas, whereas tropical reefs will face an initial reduction up to the middle of the century followed by recovery and the gain of newly suitable areas at the end of the century. A decrease in habitat suitability up to the middle of the century followed by a recovery in suitability were also found for several tropical and temperate molluscs, such as Anomia simplex, Dinocardium robustum, Lucina pensylvanica, Conus spurius, Crepidula fornicata, Melongena corona, Oliva sayana, Strombus alatus and Terebra dislocata (see online Table S2.1, MaxEnt, all variables, 2041–2060 and 2081–2100 in Saupe et al., 2014). These projections appear very likely when taking into account that the RCP2.6 scenario anticipates that increasing greenhouse gas emissions will peak in the middle of the century, with a subsequent drop to zero due to low emissions and the active removal of atmospheric greenhouse gases (IPCC, 2013). Then, under this peak-and-decline scenario, biodiversity will face alterations in the environment, which are expected to result in mild changes. Moreover, the present results are also consistent with a previous study, which showed that shifts by the middle of the century (2040–2059) will be much slower under the RCP2.6 scenario compared with the RCP8.5 scenario (Jones & Cheung, 2015). This highlights the importance of promoting effective actions on carbon emission reduction and minimizing the impacts on biodiversity.

Multi-directional shifts will be expected at the end of the century for both temperate and tropical WRs, and potential cascading effects are not likely to occur. Therefore, the majority of the new area will have a low class of potential habitat suitability. For instance, for tropical WRs, most of the projected gain in suitability will occur in the Caribbean region, but it is not evident that WRs will spread toward that area. Firstly, these areas are hotspots of biodiversity well known for their tropical coral reefs and fish communities (Floeter et al., 2008), and sessile organisms undergo high rates of space competition (Lopez-Victoria et al., 2006). Secondly, as experiments have demonstrated that polyps can survive under this optimistic scenario (Fine & Tchernov, 2007), coral reefs may not face high mortality, and space competition rates would remain high, minimizing chances for the larval settlement of sabellariids. Thirdly, in the Caribbean region, many constraints prevent the establishment of WRs, a reason why they commonly do not occur where coral reefs exist (Goldberg, 2013). Overall, under this environmentally friendly scenario, results indicate some crucial perspective: temperate and tropical WRs will not have their existence threatened, which benefits the associated fauna as well. There is not much evidence of negative cascading effects for the established community.

Tropical and temperate worm reefs under RCP8.5

Under the business-as-usual scenario, projections of suitable areas will markedly decrease for tropical WRs and increase for temperate WRs, either case indicating crucial issues for ecosystem functioning. The singular pattern of changes between congeners can be related to species-specific responses under future variability in the conditions of multiple drivers (Poloczanska et al., 2013). Some marine species can have plastic responses, such as thermal acclimation (Somero, 2002; Helmuth et al., 2006), a phenomenon particularly common in invertebrates with planktonic dispersal (Sanford & Kelly, 2011). Temperate (i.e. cold-water adapted) species are supposed to be more tolerant to global warming compared with their tropical (i.e. warm-water adapted) counterparts, mainly because the latter live close to their thermal tolerance limits (Tewksbury et al., 2008; Somero, 2010; Nguyen et al., 2011). This could support the habitat persistence and amplification of temperate WRs and the range contraction at the tropics found for tropical WRs. The response found for temperate WRs, however, is not uni-directional and may promote spatial homogenization, since habitat suitability amplifies across the whole biogeographic province. CECs promote a trend towards spatial homogenization of marine communities due to species distributional shifts or even decreases in survivorship rates (Garrabou et al., 2009; García Molinos et al., 2015). For instance, several marine organisms were affected by mass mortality events or have retracted their distribution ranges, whereas others have extended their distribution, mostly followed by recent climatic anomalies, such as heatwaves and winter storms (Garrabou et al., 2009; Wernberg et al., 2012; Firth et al., 2015).

At the end of the century, the majority of areas at low latitudes will become unsuitable for tropical WRs, most likely due to CEC-mediated shifts of distribution towards higher latitudes. This result matches predictions that the impacts of CECs would be severe in the tropics (Tewksbury et al., 2008; Poloczanska et al., 2013). Poleward expansion and equatorial contraction in distribution were also found by modelling coral reefs (Couce et al., 2013). Currently, CECs are already causing distribution shifts towards higher latitudes in several tropical marine species, from phytoplankton to fishes (Wernberg et al., 2011; Poloczanska et al., 2014). The loss of suitability in the tropical Western Atlantic Ocean predicted for tropical WRs will probably mediate cascade effects on local communities. Consequent to a WR range contraction at low latitudes, associated species are also expected to decline due to a loss of habitat. Generally, present results corroborate several modelling forecasts for the distribution of fishes, bio-engineering corals and other invertebrates, demonstrating a marked decline in habitat suitability in tropical regions that would promote local extinction (Cheung et al., 2009; Couce et al., 2013). Then, since hotspots of biodiversity are mainly located in tropical and subtropical areas (Brown, 2014), model predictions indicate grave ecological issues for ecosystem functioning, and marine biodiversity at low latitudes will be threatened.

The projected range extension of temperate WRs would incorporate ecological issues for local communities. For example, heat-tolerant species can display plastic responses to CECs or may even be able to shift their ranges over time for new latitudes tracking thermal niches (Somero, 2002). Within the genus Phragmatopoma, life-history traits indicate they are r-strategists. For instance, WRs can have an average of 65,000 individuals m⁻² (Faroni-Perez, 2014) with each female able to spawn 1500–2000 oocytes (McCarthy et al., 2003), and the dispersal and substrate colonization rates by juveniles are high (Zamorano et al., 1995; Faroni-Perez, 2014). Moreover, WRs can dominate the substratum and exclude less competitive species (Taylor & Littler, 1982). Therefore, it is likely that P. virgini will have the ability to persist and colonize new suitable sites along the South-American coasts, which would promote negative cascading effects. Since WRs provide ecosystem services for many organisms, it is very likely that the associated species may also extend their ranges by tracking the existence of reefs. Thus, the multi-directional persistence of temperate WRs might facilitate the dispersion and biological invasions of the associated fauna into new territories. Some potential biological invasions in marine ecosystems promoted by CECs have already been reported, for example to solitary and colonial ascidians and encrusting bryozoans (Stachowicz et al., 2002), and it is suggested that they will be more intense at high latitudes (Cheung et al., 2009; de Rivera et al., 2011). Another plausible issue is that if dense aggregations of WRs can exert physical pressure resulting in space competition (Taylor & Littler, 1982), then the population dynamics of pre-established adjacent organisms such as algae, anemones and molluscs will face disturbances. Thus, the enlargement of habitat suitability found for temperate WRs may drive key changes in community composition and productivity by outcompeting sensitive species.

In the Atlantic Ocean, WRs formed by the tropically distributed P. caudata are flanked in the south by Sabellaria nanella followed by P. virgini. End of century projections forecast a southward increase in habitat suitability for P. caudata and a northward increase in habitat suitability for P. virgini, making all three taxa potentially overlap their distributions around latitudes 25° and 40°S. The future trend towards sympatry among these three western Atlantic sabellariids represents an interesting ecological scenario. There is a potential for reef forming organisms to coexist through spatial and temporal niche partitioning (Knowlton & Jackson, 1994). Sympatry can drive novel patterns of resource utilization and habitat occupation. Some sympatric gregarious sabellariid species are known to exist in the same reef habitats, e.g. Neosabellaria cementarium and Idanthyrsus ornamentatus (Posey et al., 1984; Pawlik & Chia, 1991), but mostly WRs are biogenic structures of a single species. Among WRs, spatial niche tuning can be mediated by pelagic-benthic processes through the larvae sensory organs and the intrinsic ability to select specific reef microhabitats to settle (Pawlik & Chia, 1991; Faroni-Perez et al., 2016). Moreover, zonation preferences can lead to local niche partitioning. For example, the majority of S. nanella records come from subtidal areas, although recently some records have been noted in the intertidal (Bremec et al., 2014), while the majority of P. caudata and P. virgini records are from the intertidal. Niche overlap between species has been described to WRs such as between N. cementarium and I. ornamentatus (Posey et al., 1984) and Sabellaria alveolata and S. spinulosa (Holt et al., 1997). Temporal niche partitioning can also occur through physiological and reproductive strategies. For example, biological rhythms of gametogenesis, spawning and recruitment can be different in phylogenetically closely related sabellariids (see in Eckelbarger, 1979; McCarthy et al., 2003; Culloty et al., 2010; Faroni-Perez, 2014; Faroni-Perez & Zara, 2014; Aviz et al., 2017). Another possible outcome of the potential future distributional overlap is the formation of a hybrid zone. In the laboratory, some sabellariids are able to interbreed, for example Phragmatopoma caudata, P. californica, Idanthyrsus ornamentatus and Neosabellaria cementarium (Smith & Chia, 1985; Pawlik & Chia, 1991), however developed larvae are abnormal. To date there is no evidence that P. caudata, P. virgini and S. nanella have similar spawning periods and would interbreed in situ. Synchronous spawning, interbreeding and hybrid larvae deficiency would reduce fecundity, leading to local WRs decline of one or more species. These results encourage further studies to investigate whether a more finely tuned niche partitioning can support sympatric WR species under future CEC scenarios.

Overall, the present results highlight that the distribution of two important ecosystem-engineers species will shift, and depending on the future CEC scenarios, the species will have idiosyncratic, poleward or multi-directional responses. Although under the more optimistic scenario (RCP2.6), model predictions indicate the minimization of biodiversity loss and abrupt community changes, the predictions for the worst-case RCP scenario (RCP8.5) indicate range retraction/expansion and shifts for the distribution of worm reefs. Finally, it should be considered that in order to make more exact model predictions about the magnitude of multiple-drivers (e.g. salinity, nutrients and pH) on future species distribution, experimental approaches are needed to resolve potential physiological thresholds and authenticate the modelling calibration.