Genetic structure and life history are key factors in species distribution models of spiny lobsters

Abstract Aim We incorporated genetic structure and life history phase in species distribution models (SDMs) constructed for a widespread spiny lobster, to reveal local adaptations specific to individual subspecies and predict future range shifts under the RCP 8.5 climate change scenario. Location Indo‐West Pacific. Methods MaxEnt was used to construct present‐day SDMs for the spiny lobster Panulirus homarus and individually for the three genetically distinct subspecies of which it comprises. SDMs incorporated both sea surface and benthic (seafloor) climate layers to recreate discrete influences of these habitats during the drifting larval and benthic juvenile and adult life history phases. Principle component analysis (PCA) was used to infer environmental variables to which individual subspecies were adapted. SDM projections of present‐day habitat suitability were compared with predictions for the year 2,100, under the RCP 8.5 climate change scenario. Results In the PCA, salinity best explained P. h. megasculptus habitat suitability, compared with current velocity in P. h. rubellus and sea surface temperature in P. h. homarus. Drifting and benthic life history phases were adapted to different combinations of sea surface and benthic environmental variables considered. Highly suitable habitats for benthic phases were spatially enveloped within more extensive sea surface habitats suitable for drifting larvae. SDMs predicted that present‐day highly suitable habitats for P. homarus will decrease by the year 2,100. Main conclusions Incorporating genetic structure in SDMs showed that individual spiny lobster subspecies had unique adaptations, which could not be resolved in species‐level models. The use of sea surface and benthic climate layers revealed the relative importance of environmental variables during drifting and benthic life history phases. SDMs that included genetic structure and life history were more informative in predictive models of climate change effects.


| INTRODUC TI ON
Species distribution models (SDMs) are popular numerical methods which correlate geo-located species occurrence data with environmental datasets to infer habitat suitability (Elith & Leathwick, 2009).
SDMs have many applications in applied ecology and biogeography, such as predicting suitable sites for species, range shifts, and spread of invasive species and predicting future distribution patterns as a result of climate change (Melo-Merino et al., 2020). SDMs have been used extensively in terrestrial ecosystems. Their use in marine ecosystems has increased significantly in recent years, enabled by the rapid growth of remote sensing technology and increasing availability of high-resolution marine environmental data available on global databases (reviewed by Melo-Merino et al., 2020;Robinson et al., 2017). Even so, using SDMs to model distributions in marine environments remains complex, because of the dynamic nature of physical and biological variables in these systems (Palumbi et al., 2019;Robinson et al., 2011). In marine environments, speciesenvironment relationships occur at multiple spatial and temporal scales, but in most cases environmental data are restricted to sea surface and benthic layers, with limited data for midwater habitats (Assis et al., 2017). Life history phases of marine invertebrates often inhabit different discrete habitats, as drifting larvae dispersing in sea surface or midwater layers, and as juveniles and adults in benthic environments (Robinson et al., 2011). Water movements (ocean currents, eddies, fronts) are highly variable and together with larval duration and behavior determine the spatial reach of dispersal processes (Butler et al., 2011). The combination of dispersal processes, location of benthic habitats, and life history traits regulates recruitment, survival, and reproduction, and ultimately determines geographic distribution patterns (Soberon & Peterson, 2005).
Exposure to habitat changes can strongly influence the phenotypic plasticity of a species, directly contributing to population fitness, and adaptive genetic variation and evolution (Reed et al., 2011).
Local adaptation typically occurs when populations near the edge of a species range develop tolerance to more extreme environmental conditions, sometimes accompanied by a loss of tolerance to formerly suitable conditions (Huey & Hertz, 1984). Over time, locally adapted populations may become genetically differentiated, giving rise to genetically structured populations across their geographic range. Local adaptation and genetic structure within widespread marine invertebrates have been documented for planktonic taxa (reviewed in Sanford & Kelly, 2011), spiny lobster (Gaeta et al., 2020;Lavery et al., 2014;Singh et al., 2018;Truelove et al., 2017), crabs (van Tienderen & van der Meij, 2017), oysters (Burford et al., 2014), and sea urchins (Carreras et al., 2020). Recent studies have shown that incorporating local adaptation and genetic structure in SDMs can be informative in predictive models of climate change effects (Cacciapaglia & van Woesik, 2018;Hällfors et al., 2016).
Spiny lobsters (Palinuridae) support economically valuable fisheries wherever they are found. Lobster species are often regional icons, and they are among the most researched animals on earth (Phillips, 2013). Spiny lobsters occupy a trophic level as benthic consumers and are essential to the function and maintenance of reef ecosystems where they occur at high densities (Briones-Fourzán & Lozano-Álvarez, 2015). Spiny lobsters are highly fecund, and after hatching, larvae drift in ocean currents for several months to years before settling on the seafloor (George, 2005). Larvae can be dispersed over vast distances by prevailing ocean currents (Chiswell et al., 2003;Groeneveld et al., 2012) or retained near their origin by submesoscale processes and larval behavior (Butler et al., 2011;Chiswell & Booth, 1999;Kough et al., 2013;Singh et al., 2019). The effects of environmental change on spiny lobster distribution and abundance have been well-documented for several species, particularly in relation to fisheries management (Caputi et al., 2010(Caputi et al., , 2013Chávez & García-Córdova, 2011;Cockcroft et al., 2008).
Here, we chose the scalloped spiny lobster Panulirus homarus (Linnaeus, 1758) as a model organism for testing SDM habitat suitability predictions at an ocean-wide scale. P. homarus is widespread in the tropical and subtropical Western Indo-Pacific (Chan, 2010;Holthuis, 1991) and inhabits shallow coastal reef patches exposed to waves and currents. Females release larvae with a prolonged (4-6 months) drifting dispersal phase (George, 2005;Phillips & Booth, 1994). The long drifting larval phase is susceptible to the effects of climate change, especially variability in the strength and direction of ocean currents and water temperature (Singh et al., 2018(Singh et al., , 2019. P. homarus comprises three subspecies based on geography and morphological characteristics (Berry, 1974 (Farhadi et al., 2017;Lavery et al., 2014;Singh et al., 2017). Farhadi et al. (2013Farhadi et al. ( , 2017, have also shown molecular evidence that P. h. megasculptus is diverging from P. h. homarus. The aims of this study were to investigate the effects of genetic partitioning on predictions of current and future areas of habitat suitability for the commercially important spiny lobster Panulirus homarus, using SDMs. Genetic and oceanographic barriers derived from Singh et al. (2018) were incorporated in the SDMs, and the results were contrasted with SDMs that treated P. homarus as a single homogenous unit (i.e., all subspecies combined). We took both larval dispersal and benthic life history phases into account by modeling sea surface and benthic environmental data layers.
SDMs were used to forecast future changes in suitable habitats for P. homarus (to the year 2,100) under a maximum greenhouse gas emission scenario: representative concentration pathways (RCP 8.5), which depicts a worst-case scenario for the rise in global emissions (IPCC 2013). RCP 8.5 assumes that by the year 2,100, the radiative forcing level will reach 8.5 W/m 2 , sea level will have increased by >90 cm, and seawater temperature will be approximately 3°C warmer. We predicted that SDMs that incorporated genetic partitioning would reveal local adaptations that could not be discerned from models that treated P. homarus as a homogenous unit and that predicted shifts in the geographic distribution of the three subspecies under the RCP 8.5 climate change scenario would be influenced by local adaptations.

| Study region
The Western Indo-Pacific biogeographic realm (Spalding et al., 2007) covers the western and central parts of the Indian Ocean, including the east coast of Africa, Red Sea, Gulf of Aden, Persian Gulf, Arabian Sea, Bay of Bengal, and Andaman Sea, as well as the coastal waters surrounding Madagascar, the Seychelles, Comoros, Mascarene Islands, Maldives, and Chagos Archipelago. P. homarus also occurs along the western edge of the Central Indo-Pacific, including Indonesia, Thailand, Taiwan, northern Australia, and as far to the east as Japan (Holthuis, 1991). The study region included the entire known distribution range of P. homarus, extending across tropical and subtropical latitudes.

| Environmental and geographic occurrence data
Both larval dispersal and benthic life history phases were taken into account by modeling sea surface and benthic environmental data layers. Monthly averaged climatic variables (current velocity, water temperature, and salinity) for the present sea surface (to model larval habitat suitability) and benthic (to model adult habitat suitability) layers (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014) as well as forecasted data layers for the RCP 8.5 climate change scenarios for the year 2,100 were downloaded from Bio-ORACLE (Ocean Rasters for Analysis of Climate and Environment, http://www.oracle.ugent.be; Assis et al., 2017;Tyberghein et al., 2012) at 5 arcmin spatial resolution (9.2 km). Predictor variables were chosen if present and future projections were available for the study area and for both surface and benthic layers (Table 1). A variance inflation factor analysis (VIF, Marquardt, 1970) was carried out by regressing each predictor variable against the others to assess potential collinearity. The VIF analysis was performed in the R package usdm (Naimi, 2012), and variables with a VIF > 10 (an indication of collinearity) were excluded from subsequent analyses (Chatterjee & Hadi, 2006).
Georeferenced occurrence data for P. homarus were obtained from the Global Biodiversity Information Facility (www.gbif.org) and the Ocean Biogeographic Information System (www.iobis.org). Occurrence data obtained from the published literature were also added (Lavery et al., 2014;Reddy et al., 2014;Singh et al., 2018;Singh et al., 2019).
Lobster fishing locations (described by Al-Breiki et al., 2018;Amali & Wulan Sari, 2020;Fielding & Mann, 1999;Hariri et al., 2000;Jong, 1993;Mashaei & Rajabipour, 2002Priyambodo et al., 2017;Sanders & Bouhlel, 1984;Senevirathna & Munasinghe, 2014;Sultana et al., 2009;Thangaraja et al., 2015) were used to derive additional occurrence data after validation by plotting representative locations on Google Earth (earth.google.com/web/). The occurrence data consisted of 114 unique occurrence points throughout the Western Indo-Pacific and the western part of the Central Indo-Pacific ( Figure 1). Grouping of the populations by subspecies was informed by morphological identification of specimens (where available), and genetic and geographic information (Lavery et al., 2014;Singh et al., 2018Singh et al., , 2019. Principal component analyses (PCAs) implemented in the R package FactoMineR (Lê et al., 2008) using all present surface and benthic environmental variables were conducted to examine whether the occurrence data of each subspecies are linked to environmental variables.

| SDM construction and forecasting
The machine-learning method MaxEnt 3.4.1 (Phillips & Dudík, 2008) was used to develop present-day and forecasting SDMs (based on the RCP 8.5 climate change scenario) for P. homarus (subspecies combined) and for each subspecies separately. SDMs were based on marine sea surface and benthic climate layers ( Table 1). The MaxEnt models were generated using 10 bootstrap replicate runs for each model and with 10,000 background points. During optimization 25% of the data were used as test data and 75% as training data. The final models were run using all the data (100%). A logistic output was chosen, and the average of the 10 replicates was used for further analyses. The performance of each model was evaluated using the receiver operating characteristic area under the curve (AUC, Fielding & Bell, 1997). A permutation approach and jackknife tests were used to assess the importance of the predictor variables in predicting the species distribution. Initial calibration models were run with all 18 climate variable layers, but subsequent analyses were only conducted with layers which were not excluded in the VIF analysis (see above). GRASS GIS v. 7.0.4 (Neteler et al., 2012) was used to import raster maps of habitat suitability from MaxEnt and to calculate the area (km 2 ) of habitat suitability. Pixels on the raster maps were classified into probabilities: not suitable (0-0.39), low suitability (0.4-0.59), medium suitability (0.6-0.79), and high suitability (0.80-1.0). The inferred distribution range was constrained to the eastern hemisphere (0-180°E) and to latitudes between 50°N and 50°S, to exclude ocean regions considered unlikely to be naturally colonized by P. homarus. The Mediterranean Sea was also excluded from the present-day SDMs because P. homarus has not been recorded there.

| PCA outputs
For surface environmental variables, the first PCA axis explained 40%, and the second axis 25% of subspecies niche suitability TA B L E 1 The receiver operating characteristic area under the curve (AUC) and standard deviation (SD), percentage contribution, and permutation importance for surface and benthic predictor variables for Panulirus homarus subspecies combined and for individual subspecies

| Sea surface variables
The present-day and RCP 8.5 SDMs using sea surface predictor variables for subspecies combined and for individual subspecies were robust, with all AUC values above 0.9 (

| Benthic variables
The present-day and RCP 8.5 SDMs based on benthic climate layers were robust for the three subspecies combined and for individual subspecies, with all AUC values above 0.9 (Table 1)

| Species distribution probabilities for presentday and the RCP 8.5 climate change scenario
Present-day estimates of highly suitable sea surface area for larval dispersal phases were substantially larger than the size of benthic areas for juvenile and adult P. homarus (subspecies combined), and for individual subspecies (Table 2). Benthic habitat areas were spatially contained within estimated surface habitat areas. Suitable benthic habitat areas were closely aligned with shallow coastal waters (Figures 3-6), matching known seafloor habitats of P. homarus (Holthuis, 1991

| D ISCUSS I ON
At the onset of this study, we predicted that SDMs based on genetically partitioned subspecies of P. homarus would reveal local adaptations that cannot be discerned when they are modeled as a single species and that shifts in geographic distribution under the RCP 8.5 climate change scenario would be influenced by unique adaptations. Both predictions were strongly supported by the analyses presented. Sea surface SDMs for the three subspecies combined attributed variance roughly equally to current velocity (40%), water temperature (30%), and salinity variables (30%), but individual models could identify the relative importance of variables unique to each subspecies. A similar pattern was observed for benthic SDMs: Roughly equal variance attributed to current velocity (32%), water temperature (34%), and salinity variables (34%) for combined subspecies, but specific variables became relatively more important in models of individual subspecies. Overall, our findings aligned well with the theory that natural selection drives populations exposed to more extreme environmental conditions to adapt to become better suited to its local conditions than other members of the same species (Bell & Collins, 2008;Huey & Hertz, 1984).
Panulirus homarus rubellus inhabits a narrow continental shelf flanked by the upper reaches of the strong western boundary Agulhas Current. Drifting larvae released in this environment can either be retained over the narrow shelf by submesoscale processes, dispersed downstream along the coast, or become entrained in the Agulhas Current and presumably lost offshore (Singh et al., 2019).
Recent gene flow and population genetic structure analyses supported a hypothesis that some larvae are retained nearshore by lee eddies and counter currents, and that net gene flow direction was toward the southwest (Singh et al., 2018(Singh et al., , 2019. In the present study, the sea surface SDM for P. h. rubellus attributed nearly 80% of variance to current velocity variables, and the benthic SDM, some 55%. The SDM results therefore agree with gene flow and oceanographic information, confirming that P. h. rubellus is adapted primarily to survival in a strong ocean current regime. Two other spiny lobster species occurring in the same region have similarly adapted to survival in strong ocean currents, by undertaking benthic migrations to redress downstream dispersal of larvae (Groeneveld, 2002;Groeneveld & Branch, 2002;Santos et al., 2014).

Panulirus homarus megasculptus inhabits the NW Indian Ocean
but is excluded from the Red Sea and Persian Gulf, potentially because larvae are unable to tolerate high salinity water which forms in enclosed water bodies where evaporation exceeds precipitation (Jha et al., 2010;Kumar & Prasad, 1999). The hypersaline water then exits the Red Sea into the Gulf of Aden and the Persian Gulf into the Sea of Oman (Kumar & Prasad, 1999). We suggest that these areas form biogeographic transition zones to the successful dispersal of influences of cold-water upwelling cells and solar heating (Marra & Barber, 2005). Tolerance to sudden changes in water temperature would therefore be an adaptive advantage to benthic species unable to move away. P. h. megasculptus is therefore adapted to tolerate a wide range of water temperature in benthic and salinity in drifting larval phases.
F I G U R E 2 PCA plot of Panulirus homarus subspecies using all (a) surface and (b) benthic predictor variables. The ellipses represent 95% confidence intervals of subspecies groupings. The percentages by axes indicate how much variation is explained by the principal components. CV, current velocity; LT, long-term  (Araujo & Guisan, 2006); that geo-located occurrence data were accurate; and that sea surface layers could explain pelagic habitats without including midwater layers (Bentlage et al., 2013).
For benthic life history phases, occurrence records obtained from regional databases (GBIF and OBIS), the published literature and locational data from commercial fisheries were verified by plotting on Google Earth (earth.google.com/web/), and the final occurrence dataset represented the known distribution of P. homarus well, including geographic information on subspecies distribution patterns (see Figure 1). are therefore likely to influence spatial dispersal patterns, larval survival, and benthic settlement success of postlarvae (Caputi et al., 2013). An important caveat of the SDMs used in this study is that they do not consider biotic factors or species interactions in predictions of suitable habitats (Hutchinson, 1957;Soberon & Peterson, 2005) but rely only on abiotic environmental conditions, strengthened by information on dispersal capabilities (or movements) obtained from a previous process-oriented study (Singh et al., 2018).

F I G U R E 4 Surface and benthic habitat suitability maps for the
In conclusion, genetic partitioning of a widespread marine spiny lobster revealed unique local adaptations, which were not apparent in SDMs that excluded the genetic information. Individual subspecies were adapted to distinct combinations of ocean current, water temperature, and salinity regimes. Drifting larvae and benthic life history phases were adapted to different combinations of sea surface and benthic environmental variables considered. Highly suitable habitats for benthic phases were spatially enveloped within sea surface habitats suitable for drifting larvae. SDMs predicted that highly suitable habitats for P. homarus subspecies will decrease by the year 2,100, based on the RCP 8.5 climate change scenario.

ACK N OWLED G M ENTS
We would like to thank the two anonymous reviewers for their comments which greatly enhanced the quality of this work.

CO N FLI C T O F I NTE R E S T S
The authors declare that there are no competing interests.

DATA AVA I L A B I L I T Y S TAT E M E N T
Data used in this study are publicly available at https://www.biooracle.org/, www.gbif.org, and www.iobis.org.