Seascape genetics of the spiny lobster Panulirus homarus in the Western Indian Ocean: Understanding how oceanographic features shape the genetic structure of species with high larval dispersal potential

Abstract This study examines the fine‐scale population genetic structure and phylogeography of the spiny lobster Panulirus homarus in the Western Indian Ocean. A seascape genetics approach was used to relate the observed genetic structure based on 21 microsatellite loci to ocean circulation patterns, and to determine the influence of latitude, sea surface temperature (SST), and ocean turbidity (KD490) on population‐level processes. At a geospatial level, the genetic clusters recovered corresponded to three putative subspecies, P. h. rubellus from the SW Indian Ocean, P. h. megasculptus from the NW Indian Ocean, and P. h. homarus from the tropical region in‐between. Virtual passive Lagrangian particles advected using satellite‐derived ocean surface currents were used to simulate larval dispersal. In the SW Indian Ocean, the dispersion of particles tracked over a 4‐month period provided insight into a steep genetic gradient observed at the Delagoa Bight, which separates P. h. rubellus and P. h. homarus. South of the contact zone, particles were advected southwestwards by prevailing boundary currents or were retained in nearshore eddies close to release locations. Some particles released in southeast Madagascar dispersed across the Mozambique Channel and reached the African shelf. Dispersal was characterized by high seasonal and inter‐annual variability, and a large proportion of particles were dispersed far offshore and presumably lost. In the NW Indian Ocean, particles were retained within the Arabian Sea. Larval retention and self‐recruitment in the Arabian Sea could explain the recent genetic divergence between P. h. megasculptus and P. h. homarus. Geographic distance and minimum SST were significantly associated with genetic differentiation in multivariate analysis, suggesting that larval tolerance to SST plays a role in shaping the population structure of P. homarus.


| INTRODUC TI ON
Seascape genetics, where population genetics draws on ecology, oceanography, and geography to explain connectivity and spatial genetic structure, was recently reviewed by Selkoe, D'Aloia, et al. (2016).
The scalloped spiny lobster Panulirus homarus is widespread in the Indo-West Pacific, extending northwards from southeastern Africa to the Arabian Sea, and eastwards to Japan, Indonesia, and northern Australia (Holthuis, 1991). Panulirus homarus comprises several putative subspecies, and its evolutionary history and phylogeography have been popular subjects in the recent literature (Farhadi, Farahmand, Nematollahi, Jeffs, & Lavery, 2013;Farhadi et al., 2017, Lavery et al., 2014Reddy, MacDonald, Groeneveld, & Schleyer, 2014;Singh, Groeneveld, Al-Marzouqi, & Willows-Munro, 2017). Analyses of mitochondrial and nuclear sequences from several genetic markers and microsatellite loci have come to broadly similar conclusions, that at least three or four separate evolutionary lineages exist, occupying distinct geographic ranges.
The genetic information largely agrees with morphological groupings (Berry, 1974a). Panulirus homarus rubellus (red megasculpta form) from the Southwest (SW) Indian Ocean is the most divergent and possibly a separate species based on phylogenetic analyses using mitochondrial and nuclear markers (Farhadi et al., 2017;Lavery et al., 2014;Singh et al., 2017). Panulirus homarus "brown" is restricted to the Marquesas Islands in the Central Pacific, and P. homarus homarus (green microsculpta form) is widespread in intervening areas (Farhadi et al., 2017). The status of P. homarus megasculptus (green megasculpta form) from the Northwest (NW) Indian Ocean (mainly Arabian Sea) remains unclear. Neither Singh et al. (2017) nor Lavery et al. (2014) could find a phylogenetic basis for distinguishing it from P. h. homarus from the SW Indian Ocean based on mtDNA and nuclear loci, but Farhadi et al. (2013), Farhadi et al. (2017) found evidence for reduced gene flow between the NW Indian Ocean and populations in Kenya and Tanzania, further to the south.
The Western Indian Ocean extends from eastern South Africa to the Arabian Gulf, including waters around Madagascar and several small island states, such as Mauritius, Seychelles, and Comoros ( Figure 1). It is a tropical/subtropical region with high biodiversity and endemism, brought about by diverse habitats along the continental shelf, and around isolated islands or subsurface plateaus (van der Elst et al. 2009;Halo, Malauene, & Ostrowski, 2017). Two western boundary currents dominate the coastal oceanography: the Agulhas Current to the south and the Somali Current to the north. The Somali Current reverses its flow direction seasonally, under the influence of the SW monsoon (flows northwards) and the northeast (NE) monsoon (flows southwards) seasons (Schott & McCreary, 2001). It also receives waters from the East African Coastal Current along coastal Tanzania and Kenya. Water from the South Equatorial Current enters the northern Mozambique Channel, forming a series of eddies which propagate southwards through the channel (Halo et al., 2017;Otwoma & Kochzius, 2016). These waters are augmented by eddies originating from the East Madagascar Current at the southern end of the channel (de Ruijter et al., 2004;Ridderinkhof, Bars, Heydt, & Ruijter, 2013), and their confluence forms the upper Agulhas Current off southern Mozambique and eastern South Africa. Spalding et al. (2007) subdivided the Western Indian Ocean into 12 marine ecoregions, or "areas of relatively homogenous species composition, clearly distinct from adjacent systems." Biogeographic forcing agents characteristic of ecoregions include upwelling cells (e.g., the Central Somali Coast), nutrient inputs from freshwater influx (Sofala Bight/Swamp Coast), the influence of ocean currents (Northern Monsoon Coastal Current), bathymetric or coastal complexity (East African Coral Coast), or differences in temperature regimes or sediments. Toward the south, the subtropical Natal ecoregion is characterized by high species diversity, resulting from an influx of species from the Indo-Pacific, transported by the Agulhas Current (Bustamante & Branch 1996;Spalding et al., 2007). Madagascar comprises two separate ecoregions, Southeast Madagascar with cooler water and upwelling cells, and Western and Northern Madagascar where this upwelling is negligible or absent (Spalding et al., 2007).
We used biological information, habitats, and ocean surface currents in a seascape genetics approach to explore the genetic structure observed within the P. homarus subspecies complex in the Western Indian Ocean. The hypothesis that isolation by distance is the main driver of differentiation among subspecies was tested.
Alternatively, biogeographic breaks determined by ocean circulation patterns and environmental gradients may have constrained gene flow, resulting in genetic differentiation between the three groups.
Specific aims were to: (a) evaluate potential phylogeographic boundary zones; (b) use virtual particle tracking simulations to infer larval dispersal patterns based on surface currents; and (c) determine the extent to which genetic structure correlates with environmental variables and simulated larval movement.  Dao, Todd, and Jerry (2013) and 13 loci (Pho G01, Pho G03, Pho G21, Pho G22, Pho G25, Pho G27, Pho G30, Pho G32, Pho G35, Pho G36, Pho G42, Pho G53, and Pho G58) described by Delghandi et al. (2015), was chosen for PCR amplification. Forward primers were labeled with a fluorescent dye on the 5′ end. The microsatellite panel was partitioned into six multiplex reactions and a single reaction for Orn 17 (Supporting Information Table S1).

| Microsatellite sampling
PCR products were carried out in 10 µl volume reactions, which contained 5 µl of 1× Kapa 2G Fast Multiplex Mastermix (Kapa Biosystems), 0.2 µmol/L of each primer, and 1 µl of 10-50 ng of F I G U R E 1 Map of Panulirus homarus sampling locations and the number of samples collected per locality. The major current systems (adapted from Lutjeharms, 2006) are also shown genomic DNA. For multiplexes A, B, C, and Orn 17, cycling conditions were as follows: an initial denaturation at 95°C for 5 min,  (Chapuis & Estoup, 2007) was used to estimate the null allele frequencies for each marker and population. Bootstrap resampling over loci was conducted with 100,000 replicates. The excluding null alleles (ENA) method was used to compare null allele corrected and uncorrected global and population F ST values using paired t tests.

| Genetic structure and gene flow
Population genetic structure between the three subspecies and between different localities was evaluated using STRUCTURE v. 2.3.4 (Pritchard, Stephens, & Donnelly, 2000). K was set from 1 to 12. Each K was run for 10 iterations. The admixture model was chosen and run with and without sampling location as a prior. A total of 1,000,000 MCMC iterations were performed after a burn-in period of 250,000 iterations. Structure Harvester (Earl & von Holdt, 2012; https://tay-lor0.biology.ucla.edu/structureHarvester/) and Pophelper (Francis, 2017; https://pophelper.com/) were used to identify the number of genetic clusters (K) that fit the data best (Evanno, Regnaut, & Goudet, 2005). An analysis of molecular variance (AMOVA) was performed using Arlequin, to assess the magnitude of genetic differentiation between the groups identified by the STRUCTURE analysis.
The role of seascape features in population genetic structure, and potential contact zones, was inferred using the Geneland v.4.0.3 R-package (Guillot, Mortier, & Estoup, 2005). This method couples multilocus microsatellite genotype data and geospatial information by using geographic sampling locality co-ordinates as priors, enabling the detection of subtle population structure, and also to infer the locations of genetic discontinuities or transition zones (Guillot et al., 2005). The analysis consisted of 10 independent runs of 1,000,000 MCMC iterations, sampling every 100 generations. The number of predefined genetic clusters (K) ranged from 1 to 12. The spatial model was employed, with and without accounting for null alleles. The correlated allele frequencies model was chosen to detect subtle genetic structure (Guillot, 2008).
The null hypothesis of isolation by distance was examined by testing the relationship between genetic and geographic distance using a Mantel test. The analysis was carried out using the program IBDWS (Jensen, Bohonak, & Kelley, 2005, https://ibdws.sdsu. edu/~ibdws/). The genetic distances (F ST /(1 -F ST ) were regressed against the geographic distances (straight-line distances on Google Earth) with 1,000 randomizations between populations.

| Environmental factors
The relationship between genetic differentiation of the P. homarus groups (pairwise F ST ), geography, and environment was tested using the multivariate multiple regression approach of a distance-based redundancy analysis (dbRDA) (Legendre & Anderson, 1999). Monthly composite oceanic variable data of ocean productivity chlorophyll a (chl a, mg/m 3 ), sea surface temperature (mean, maximum and minimum SST, °C), and turbidity (diffuse attenuation coefficient KD490; m −1 ) were obtained from the GMIS Marine Geodatabase (mcc.jrc. ec.europa.eu), based on the NASA MODIS terra satellite at a 4-km spatial resolution. Data were collected from 2012 to 2015, coinciding with the genetic sampling regime. A matrix of genetic differentiation between sampling sites (F ST ) was used as the response variable, and mean measurements of environmental variables (chl a, mean, max and min SST, turbidity), and larval recruits (based on estimates from ocean modeling simulations as the number of larvae that reach settlement at each site on the coast) were used as the explanatory variables. The dbRDA analysis was carried out using the "capscale" function in the Vegan package (Oksanen et al., 2015) in R. The correlation among variables was measured using Pearson's correlation coefficient in Vegan. Where variable pairs were highly correlated (Rellstab, Gugerli, Eckert, Hancock, & Holderegger, 2015), the biologically more relevant one to P. homarus was used. Minimum SST was considered an important correlate because it decreases along a gradient from tropical (P. h. homarus preference) to subtropical conditions (P. h. rubellus preference) (Berry, 1971a(Berry, , 1974a. SST is expected to affect phyllosoma larvae that disperse in the upper layers of the water column. Turbidity was selected because P. homarus juveniles and adults prefer turbid conditions in the shallow subtidal zone (Berry, 1971a;Holthuis, 1991). The Euclidean geographic distance matrix was transformed into a continuous rectangular vector by principal co-ordinates analysis (PCoA) using the "pcnm" function in Vegan. The "ordistep" function was used to identify the most significant variables to include in the model (Supporting Information   Table S2). The "varpart" function in Vegan was used to estimate the amount of genetic differentiation due to each predictor variable. Significance of the full model was assessed by comparing it to the null model containing only the intercept, and significance of each predictor variable was tested by comparing models which test each predictor after removing the effects of the other predictors.
Significance tests were carried out using the "ANOVA.cca" function in Vegan.

| Larval dispersal simulations
We explored the effects of ocean circulation on larval dispersal and settlement, and hence genetic differentiation, in the virtual particle tracking tool, Parcels v.0.9 (Probably A Really Computationally Efficient Lagrangian Simulator; Lange & van Sebille, 2017).
Lagrangian particle trajectories were simulated using Parcels with the equation: where X is the three-dimensional position of a particle and v (x, ) is the three-dimensional velocity field at that location from an ocean general circulation model.
Ocean velocity data were obtained from the GlobCurrent Project (Johannessen et al., 2016; https://www.globcurrent.org/) which provides daily average surface ocean current estimates for the world's oceans at 25-km spatial resolution by merging data from multiple satellite altimeter missions with gravity models and combining these with Ekman currents (at surface and 15 m depth) estimated from drogued surface drifter data, Argo floats, and near-surface winds.
The u-and v-component geostrophic velocities are defined by the following equations (Saraceno, Strub, & Kosro, 2008): where υ represents the zonal component, ν is the meridional component, is the variation in the sea surface height, x∕ y is the distance between grid points, g is the gravitational acceleration, and f is the Coriolis force.
These products are combined with the Ekman currents esti- Known life history parameters of P. homarus (Table 1)

| Genetic diversity
Summary statistics of 271 P. homarus individuals genotyped across 21 microsatellite loci are shown in Table 2 and in Supporting   Information Tables S3A and B. Microchecker did not detect evidence for scoring errors due to stuttering or large allelic dropout, but indicated that null alleles may be present at ten loci (Orn 11, Orn 12, Orn 16, G21, G25, G27, G30, G32, G53, and G58). All loci were included in subsequent analyses because no significant difference was found between ENA-corrected versus non-ENA-corrected F ST values, overall and for each population (p > 0.05).
However, out of the 12 populations, Orn 5 deviated significantly from HWE in 5 populations, Orn 12 in 6, Orn 16 in 3, Orn 17 in 5, Orn 32 in 1, G21 in 5, G22 in 2, G32 in 6, G36 in 2, G42 in 9, and G53 in 1 population. The analyses were also conducted after removing loci that deviated from HWE, which resulted in a smaller dataset with 8 loci from 257 individuals. The results from analyses of the smaller dataset did not contradict those obtained from using the full dataset, and hence, all loci were retained (see Supporting Information Appendix S1

| Genetic structure
Genetic differentiation between the subspecies could clearly be distinguished using the clustering analyses implemented in STRUCTURE and Geneland. STRUCTURE analysis indicated that K = 3 was optimal for the P. homarus subspecies dataset (Figure 2a,b). The AMOVA indicated that there was 4.5% difference between each of the groups identified by the STRUCTURE analysis (

| Environmental factors
The "ordistep" function identified geography, minimum SST, and larval recruits of 2009 as the most significant variables affecting the response matrix of genetic variation (Supporting Information Table   S2). This model was significant when compared to the null model (F = 6.20, p = 0.001). The first axis accounted for 55% of the fitted variation and 35% of the total variation, and the second axis explained 42% of the fitted variation and 27% of the total variation ( Figure 4). Variance partitioning indicated that minimum SST explained 18% of the total variation (F = 6.98, p = 0.004), and that larval recruits of June 2009 could explain 8% of the variation, but was not significant (F = 1.82, p = 0.15). Geography explained 13% of the total variation (F = 4.67, p = 0.019) and was significant (Figure 4).

| NW Indian Ocean and Kenya
At the NW Indian Ocean sites (Oman, Yemen) and Kenya, the peak breeding season of P. homarus is in June. Particle density plots are shown for June ( Figure 5) and January (Supporting Information Figure S1) along with dispersal trajectories for January and June

| SW Indian Ocean
Particle density plots for SW Indian Ocean sites (Mozambique, Madagascar, and South Africa) for the peak breeding season in January are shown in Figure 6. Density plots are also shown for June (Supporting Information Figure S4) along with particle trajectories for January (Supporting Information Figure S5) and June (Supporting Information Figure S6)

| Genetic structure and boundary zones
We examined the fine-scale population genetic structure and phy-  Table S4). Farhadi et al. (2013Farhadi et al. ( , 2017 ) also found evidence of genetic divergence between Arabian Sea and more equatorial P. homarus populations, based on hypervariable control region sequences and microsatellites. The two groups are morphologically distinguishable (microsculpta vs. macrosculpta forms) and inhabit distinct geographic ranges, with some overlap in Oman (J. Groeneveld, personal communication). Nevertheless, more conserved mitochondrial and nuclear DNA markers could not differentiate the two groups (Lavery et al., 2014;Singh et al., 2017), thus suggesting a recent divergence.
The simulated particle trajectories strongly support the retention of P. h. megasculptus larvae within the NW Indian Ocean.
but the eddies weaken seasonally (Bower & Furey, 2012), giving rise to eastward dispersal of particles released in January, into the Arabian Sea. Panulirus homarus megasculptus in Oman has a prolonged breeding season with multiple broods (Al-Marzouqi, Al-Nahdi, Jayabalan, & Groeneveld, 2007;Al-Marzouqi, Groeneveld, Al-Nahdi, & Al-Hosni, 2008), and drifting larvae are expected to be present in the Arabian Sea throughout the year. Ocean circulation is influenced by the seasonally reversing Somali Current and eddies which develop along the coasts of Oman, Yemen, and Somalia (Schott & McCreary, 2001). These eddies likely retain larvae in the NW Indian Ocean, thereby facilitating the incipient genetic structure observed between P. h. megasculptus and P. h. homarus (our results and see also Farhadi et al., 2013). The oceanographic boundary is not restricted to dispersal of lobster larvae and also prevents northwards dispersal of starfish Acanthaster planci larvae from east Africa into the northern Arabian Sea (Vogler et al., 2012).

| Environmental drivers of genetic discontinuities
The multivariate analysis dbRDA was used to investigate the correlation between genetic variation, environmental, and spatial factors across the seascape. Minimum SST was correlated with genetic variability (Figure 4) (Spalding et al., 2007). We propose the Delagoa Bight ecoregion as a contact zone between P. h. homarus and P. h. rubellus. The Bight is a coastal indentation at the southern end of the Mozambique Channel, and both its location and unique oceanographic features may have contributed to the genetic discontinuity seen here. The

| A biogeographic perspective
Bight is affected by eddies moving gradually southwards through the Mozambique Channel, and also by cross-channel eddies originating from the East Madagascar Current, on its path around the southern tip of Madagascar (Cossa, Pous, Penven, Capet, & Reason, 2016;Lutjeharms, 2006). Additionally, the Bight is located near the upper reaches of the Agulhas Current (Lutjeharms, 2006), and particles released here were frequently advected southwestwards.
Particles released at Mozambican sites in the peak January breeding season of P. h. rubellus were often entrained in the Agulhas Current or its inshore filaments, and dispersed southwestwards.
Some particles retroflected eastwards into the SW Indian Ocean (presumably lost) or leaked into the SE Atlantic basin through ring shedding, eddies, and filaments (Gordon, 2003 to swim and position themselves in the water column to benefit from water movements at different depths (Bradford et al., 2005;Chiswell & Booth, 1999, 2008. Vertical migrations to deeper habitats during the day would reduce predation, with phyllosomas rising to the surface at night to feed (Bradford et al., 2005;Butler et al., 2011).
Directional swimming is also not included in the particle model. Jeffs, Montgomery, and Tindle (2005) demonstrated that late-stage larvae receive sensory cues directing them toward the shore; cues may include acoustic reef sounds (Hinojosa et al., 2016), chemical (Hadfield & Paul, ), celestial, magnetic field, or electrosense signals . The postlarval puerulus stage is also a strong swimmer, able to cross the shelf using stored energy reserves (Phillips & Olsen, 1975;Wilkin & Jeffs, 2011). It is therefore clear that behavior will affect the dispersal of larvae (Butler et al., 2011). Particles released in southeast Madagascar were propagated across the Mozambique Channel on several occasions, reaching the African shelf. Directed swimming by larvae dispersed to here may have assisted them in reaching coastal settlement habitats in Mozambique or eastern South Africa.
Migrations and spawning behavior of adult spiny lobsters may mitigate the influence of the Agulhas Current on the downstream dispersal of larvae (Groeneveld & Branch, 2002;Santos, Rouillard, & Groeneveld, 2014), through determining the location and timing of larval release. Vertical migration of larvae in the water column has been shown to favor larval retention close to local shores (Phelps, Polton, Souza, & Robinson, 2015;Rivera et al., 2013 in the Agulhas Current, we suggest that none, or very few of them would have been returned to the coast. The particle model used here may overstate the magnitude of larval dispersal. The spatial resolution of the GlobCurrent model is 25 km, and coastal submesoscale processes are therefore not resolved. Fine-scale hydrographic processes that may "push" or "pull" particles toward or away from the coast are most likely crucial in determining larval settlement patterns, but could not be resolved with the model used. While global operational ocean models (e.g., https://hycom.org/or https://marine.copernicus.eu/) are approaching spatio-temporal resolutions to resolve submesoscale processes, they are unable to represent the correct levels of variability, particularly in the highly dynamic Agulhas Current (Meyer, Braby, Krug, & Backeberg, 2017).
Particle tracking simulations were constrained to January and June to provide seasonal contrast to dispersal patterns. The simulations included the peak breeding seasons of P. h. rubellus in the SW Indian Ocean (January) and of P. h. megasculptus and P. h. homarus in the NW Indian Ocean and Kenya (June). Simulations were further restricted to 2009 and 2010 and showed considerable inter-annual variability-a factor that has been implicated in recruitment patterns of species with drifting larvae (Botsford, 2001;Chelton, Bernal, & McGowan, 1982;Griffin, Wilkin, Chubb, Pearce, & Caputi, 2001

ACK N OWLED G M ENTS
We thank two anonymous reviewers whose comments on earlier drafts helped to substantially improve this work. We thank

DATA ACCE SS I B I LIT Y
Data are available from the Dryad Digital Repository: https://doi. org/10.5061/dryad.hf068m7.