Disentangling genetic from environmental effects on phenotypic variability of southern rock lobster (Jasus edwardsii) postlarvae

Abstract Environmental conditions experienced during the larval dispersal of marine organisms can determine the size‐at‐settlement of recruits. It is, therefore, not uncommon that larvae undergoing different dispersal histories would exhibit phenotypic variability at recruitment. Here, we investigated morphological differences in recently settled southern rock lobster (Jasus edwardsii) recruits, known as pueruli, along a latitudinal and temporal gradient on the east coast of Tasmania, Australia. We further explored whether natural selection could be driving morphological variation. We used double digest restriction site‐associated DNA sequencing (ddRADseq) to assess differences in the genetic structure of recently settled recruits on the east coast of Tasmania over 3 months of peak settlement during 2012 (August–October). Phenotypic differences in pueruli between sites and months of settlement were observed, with significantly smaller individuals found at the northernmost site. Also, there was a lack of overall genetic divergence; however, significant differences in pairwise FST values between settlement months were observed at the southernmost study site, located at an area of confluence of ocean currents. Specifically, individuals settling into the southernmost earlier in the season were genetically different from those settling later. The lack of overall genetic divergence in the presence of phenotypic variation indicates that larval environmental history during the dispersal of J. edwardsii could be a possible driver of the resulting phenotype of settlers.


| INTRODUC TI ON
Genetic connectivity in many marine invertebrates is contingent on the dispersal of pelagic larvae (Hedgecock, 1986). Highly structured populations typically result from larvae with short pelagic durations or nonplanktonic larval stages, whereas species with long pelagic larval durations (PLDs) have shown little genetic structure across large geographical scales (Sandoval-Castillo et al., 2018;Villacorta-Rath et al., 2016;Watts & Thorpe, 2006). Long PLDs, however, could lead to reduced recruitment success due to the larger window of time where predation and starvation can act upon individuals (Morgan, 2020).
Although it has been proven that fast growth promotes larval survival (Crean et al., 2011;Marshall et al., 2010), the driver of phenotypic variability is a complex interaction of maternal fecundity, maternal energy investment in the offspring and environmental factors (Cobb et al., 1997;González-Ortegón & Giménez, 2014;Green et al., 2014).
Indeed, the environment has profound effects on the metabolism and growth of larvae, which can be carried over throughout their life history (Marshall et al., 2003). Cohorts of marine invertebrate recruits exhibit high phenotypic divergence due to the wide range of environmental conditions experienced during the pelagic larval stage (PLD) (Rey et al., 2016). For example, if egg-hatching coincides with low food availability, larvae may experience poor survival or high postsettlement mortality (Giménez, 2006). Also, the offspring of females that become reproductive at different times of the reproductive season will have slightly different phenotypes at settlement, probably because of the environmental history experienced during dispersal (Giménez, 2006;Kunisch & Anger, 1984). Although evidence suggests the existence of phenotypic differences in invertebrate recruits across geographic and temporal scales (reviewed by Sanford & Kelly, 2011), it is unknown whether these differences are purely driven by dispersal history or if local adaptation also plays an important role in shaping recruitment (Hedgecock, 1986).
The present study explored the link between phenotype and genotype using the southern rock lobster, Jasus edwardsii, as a model species given its large reproductive output and protracted pelagic larval duration. This species is highly fecund, producing approximately 500,000 eggs per brood . Adults experience seasonal reproduction that extends from May to September, but egg-hatching periods vary slightly across populations through a latitudinal gradient (Powell, 2016). Jasus edwardsii exhibits the largest PLD of all species within the Palinuridae family, spending between 12 and 24 months before metamorphosing into a puerulus postlarva (Booth, 1994), giving a large window for environmental factors to shape larval phenotypes. Recruitment monitoring of this species in southeast Australia has evidenced large variability in abundance (Hinojosa et al., 2017;Linnane et al., 2014) and in size-at-settlement (Booth, 1979(Booth, , 1994, as it is expected for highly fecund species (Cobb et al., 1997). Additionally, sweepstakes reproductive success and natural selection have been proposed as possible mechanisms driving chaotic genetic patchiness in this species (Villacorta-Rath et al., 2018).
Based on previous evidence of phenotypic and genetic differences in recruits, we hypothesized that dispersal history drives phenotypic variability of recently settled individuals. The present study explored patterns of morphological and genetic differences in recently settled J. edwardsii along a latitudinal and temporal gradient in Tasmania, Australia. More specifically, the aims of this study were: (1) to compare size-at-settlement of J. edwardsii pueruli at three sampling sites during three consecutive settlement months, (2) to determine genetic differences in recently settled J. edwardsii across a geographical and temporal gradient, and (3) to investigate if the phenotypic variation was related to neutral and non-neutral genotypic variation.

| Sample collection
Lobster pueruli were caught in crevice collectors (Booth & Tarring, 1986) deployed at three sites on the east coast of Tasmania, southeast Australia ( Figure 1) (Gardner et al., 2000). Collectors were attached to permanent structures fixed on sandy substratum at three to nine meters of depth (Gardner et al., 2000). Puerulus collectors were serviced monthly by a diver who placed a mesh bag around each collector and attached a rope with a buoy to each bag. Collectors were then hauled into the boat, where they were cleaned (Gardner et al., 2000). All pueruli were removed from the collectors and immediately stored in 90% ethanol. For the purpose of this study, only recently settled pueruli (stages 1 and 2) during the months of August, September, and October 2012 were considered (Table 1).

| Comparing pueruli size-at-settlement
Pueruli individuals were removed from ethanol, blot dried to eliminate the excess ethanol, and weighed (g ± 0.01). To determine whether pueruli weight varied with settlement month, site, and the interaction of both, a two-way analysis of variance (ANOVA) was performed using the function "aov" implemented within the R package stats v4.2.1 (R Core Team, 2022). A post hoc examination of differences between individual weights was carried out using the Tukey's "honest significance difference" method using the function "TukeyHSD" of the R package stats v4.2.1.

| Analyses of raw sequencing data and reference catalog building
An initial quality check of raw indexed data was performed using FastQC v0.10.1. Data were then demultiplexed using the "process_ radtags" protocol from Stacks v1.29 (Catchen et al., 2011) and hard trimmed to 75 bp to ensure that Phred Quality Score (Q Score) of all reads was above 30. Finally, demultiplexed libraries were filtered for bacterial and viral content using Kraken-gcc v0.10.4 (Wood & Salzberg, 2014).
An initial run of the rad-loci pipeline, read alignment, and SNP calling was performed and samples with more than 25% missing data were discarded in order to maximize the final number of polymorphic loci obtained. A total of 70 samples were discarded and 193 samples with less than 25% missing data remained. The subsequent rad-loci pipeline run was therefore performed with a subset of 193 samples using VSearch v1.1.3 (https://github.com/torog nes/vsearch). Initially, all reads were pooled and clusters with a depth of at least 193 were retained, assuming that there would be at least one read per sample represented in each cluster. Subsequently, all clusters that were comprised of reads with more than a 4% mismatch (3 bp) were discarded.
Assuming that each member of a cluster was an allele, only clusters that had between two and 16 members were retained. Finally, only samples with a minimum of 10,000 alleles and 30% missing data were retained for subsequent read alignment and SNP calling.

| Read alignment, variant calling, and SNP filtering
Individual filtered reads of all samples were aligned to the reference catalog using bwa-intel v0.7.12 (Li, 2013). Variant calling was performed using the Genome Analysis Toolkit (GATK) v3.3_0 (McKenna et al., 2010), with a further correction of the variant call format (vcf) file performed to ensure the accuracy of the reference and alternate allele calls, and to filter out false positives. In the absence of a reference genome, the correction was made by calculating the ratio between the highest quality score over depth ("QD" in vcf file) and the lowest QD. If one sample at a specific position had a ratio threshold of 10, which corresponds to a 10% error on a Phred scale, the SNP at that position was replaced with missing data.
SNP filtering was performed in vcftools-gcc v0.1.13 to ensure that only bi-allelic data were present (−-min-alleles 2, −-max-alleles 2), to remove SNPs that were potentially in linkage disequilibrium (−-min-r2 0.2), to discard SNPs with a minor allele frequency (MAF) of less than 5% (−-maf 0.05), and to ensure that the minimum SNP depth was 5 (−-minDP 5). The maximum amount of missing data for each locus was set to 25% (−-max-missing 0.75) and individuals with more than 25% missing data were removed from subsequent analyses. Finally, only one SNP per locus (−-thin 75) was retained.
A total of 22 technical replicates remained after the removal of individuals with more than the missing data threshold. A Principal Component Analysis (PCA) was performed to visualize the spatial distribution of replicated samples using the function "dudi.pca" of the R package ade4 v1.7-19 (Dray & Dufour, 2007).

| Outlier SNP identification
Two outlier identification methods were used in order to minimize using BLAST+ v2.2.29. Queries with statistically significant e-values (E < 10 −6 ) and more than 84% identity were considered as valid alignments. Transcriptome sequences that provided significant alignments were annotated using the Trinotate pipeline (https://trino tate.github.io/) to determine whether they aligned with any known protein domain.

| Analyses of genetic diversity
Global F ST values and confidence intervals for each SNP panel were estimated using the R package mmod v1.3.2 (Winter, 2012).
Additionally, pairwise F ST values between sampling sites and years as well as confidence intervals were calculated in the R package hierfstat v0.04-22 (Goudet, 2005). A false discovery rate correction (FDR) was applied to calculated p-values using the function "p.adjust" of the R package stats v4.2.1 (R Core Team, 2022).
Additionally, a discriminant analysis of principal components (DAPC) and a PCA were performed to determine the possible number of genetic clusters based on allele frequencies of all three sampling sites across three settlement months using the R packages adegenet v2.1.8 (Jombart, 2008) and ade4 v1.7-19 (Dray & Dufour, 2007).
Genetic diversity was quantified using the standardized individual heterozygosity (sh). This metric represents the proportion of heterozygous loci over the mean heterozygosity across all markers, so that heterozygosity of all individuals is measured on the same scale (Coltman et al., 1999). The standardized individual heterozygosity was calculated for both SNP panels at (1) each site with all three settlement months combined and (2) at each settlement month at each site using the R package Rhh v1.0.1 (Alho et al., 2010). Differences between the maximum and minimum values of standardized individual heterozygosity were examined using a Mann-Whitney test. A Bonferroni correction was applied in order to account for multiple comparisons, and α was set at the 0.05/4 level.

| Isolation by time
Isolation by time (IBT) was estimated using Mantel tests between temporal and genetic distance using 9999 permutations within the R package ade4 v1.7-19 (Dray & Dufour, 2007). The function "dist.
gene" was used to compute the genetic pairwise distance between samples using the R package ape v5.6-2 (Paradis et al., 2004). The pairwise distance between two individuals is the number of loci for which they differ, divided by the total number of loci (Paradis et al., 2004). The temporal distance matrix was calculated using the

| Differences in time and size-at-settlement
The weight of individual puerulus differed significantly between sites and settlement months throughout the 2012 winter settlement peak. However, the interaction between the settlement site and month did not differ significantly (Table 2). Pueruli were lighter in weight in the most northerly site, Bicheno, and heavier at the midsampling site, South Arm ( Figure 2; Table 3).

TA B L E 3 Differences between mean
individual weight per month and site of settlement and confidence intervals at a 95% confidence level based on the "Tukey's Honest Significant Difference" method.
samples across multiple lanes did not introduce a large bias in the catalog-building process.

| Analyses of genetic diversity
There was no evidence of genetic differentiation in the neutral SNP panel across all sites and settlement months (global F ST = −0.001, n.s., Figures 4 and 5). However, pueruli settling into Recherche Bay during August was genetically distinct from individuals settling into Bicheno and South Arm (Table 4), corresponding to the phenotypic differences between the sites described above. In addition, individuals settling into Recherche Bay during August were genetically distinct from those settling during October (Table 4).

Allele frequency variation was driven by the settlement site,
although it accounted for a small proportion of the genetic variation between pueruli when analyzing the neutral SNP panel (Table 5).
On the other hand, the SNPs under putative positive selection exhibited a stronger overall level of genetic differentiation (Global In support of the small F ST differentiation between settlement months in Bicheno and South Arm, the isolation by time analyses showed no significant correlation between genetic distance and distance between recruitment time (Table 6). Similarly, the significant F ST values found among the Recherche Bay sampling periods were confirmed by the significant Mantel test at this site; however, the r value was very low (r = 0.087), indicating the existence of a very low correlation. Pueruli settling into Recherche Bay during the same month was genetically more similar to each other than individuals settling one or two months apart, as shown by the slight but significant regression ( Figure 6).

| DISCUSS ION
The present study examined pueruli size-at-settlement and genetic signatures during a winter settlement peak at three sites on the east coast of Tasmania, Australia. Significant differences in weight be-  More specifically, variations in size-at-settlement, measured as puerulus carapace length, have been observed in J. edwardsii arriving during different years and seasons at the same site (Booth, 1979;Booth & Tarring, 1986). It has been observed that stage 1 pueruli settling during winter in the northwest of New Zealand were significantly bigger than those settling during other times of the year (Booth & Tarring, 1986). Also, in a three-year study, Booth (1979) found significant differences in size-at-settlement of J. edwardsii arriving at one site in northeast New Zealand between most years and months of settlement. The authors attributed the phenotypic differences to: (1) different cohorts of phyllosoma metamorphosing into pueruli and settling simultaneously, and (2) to pueruli originating from different parental stocks (Booth, 1979). Here, we aimed to explore J. edwardsii pueruli differences in size-at-settlement further and investigate whether these could have a genetic origin. The lack of neutral genetic structure found in the present study indicates that it is unlikely that the phenotypic variation is directly linked to different larval sources. Thus, based on our results, we hypothesize that the observed phenotypic differences between pueruli settling into three sites in Tasmania during the study period were caused by the environmental history of settling individuals throughout dispersal.

| Dispersal history driving phenotypic variation
The environment has profound effects on the metabolism and growth of larval lobsters, which can be carried over into subsequent stages . Jasus edwardsii settling into east Tasmania can spend up to 2 years as phyllosoma larvae, encountering a large gradient of environmental factors (Bruce et al., 2007). Temperature and food availability are the main factors affecting the growth rates of larvae and can determine settlement success (Pechenik, 2006).
Therefore, individuals that encountered temperature and food conditions that promoted growth during the late phyllosoma stage, when they are likely closer to settlement grounds, could have a higher accumulation of energy reserves and metamorphose into larger pueruli (Jeffs et al., 2001). If the significant differences in pueruli weight across sites were caused by differences in environmental conditions during development, the persistence of size differences across the study sites would require that larvae were traveling in cohorts and experiencing different temperature and nutrition regimes, differences in size-at-settlement in J. edwardsii found in the present study. A mixture of collective dispersal and convergence of larvae at later stages has been described in a reef fish cohort settling into the same site (Shima & Swearer, 2016). If this also occurs in J.
edwardsii, then the genetic signature of settling individuals will be more randomly distributed, with no genetic divergence, as seen in Bicheno and South Arm but with a large variation in phenotypes of settlers due to divergent dispersal histories. Our data show that the total weight range (difference between maximum and minimum weight) varies across months within each sampling site, suggesting

TA B L E 5
Analysis of molecular variance (AMOVA) results using genetic distance as a function of (1) settlement month and (2)  Finally, the temporal differences in pueruli weight could be attributed to differences in environmental conditions related to the time of hatching. A latitudinal gradient of the incubation period due to differences in temperature has been described in the Dungeness crab, Cancer magister (Shirley et al., 1987) where individuals hatching at lower temperatures were larger than those hatching in warmer waters. Also, if time and place of hatching can determine where larvae will be advected, this would result in different larval cohorts experiencing divergent environmental conditions. It is therefore possible that females spawning at different times of a reproductive season can produce offspring with slightly different phenotypes (Kunisch & Anger, 1984;Shirley et al., 1987). This could have been the case in the present study, where individuals recruiting during different months exhibited significant differences in weight-at-settlement. The same trend was observed in the green crab, Carcinus maenas, where larval size at metamorphosis and survival of the first crab instar under starvation varied significantly across four supply events (Rey et al., 2016).
Individuals sampled earlier in the season exhibited smaller sizes than those sampled later, probably due to less favorable environmental conditions encountered by the former group (Rey et al., 2016).
The question that arises here is whether fitness in J. edwardsii settlers is related to size-at-settlement and can lead to successful recruitment. Although there is a general hypothesis that larger offspring are fitter than smaller offspring, the link between size and fitness is not completely understood (Giménez, 2006;Marshall et al., 2006). Larger individuals generally exhibit improved antipredatory behavior, are bolder, can forage for a longer time, and better resist periods of starvation (Dingeldein & White, 2016;Johnson et al., 2017). However, it is often the environment (biotic and abiotic factors) that ultimately provides the selective pressure and determines survival (Marshall et al., 2006). If smaller size represents a disadvantage for postsettlement survival, it would then result in fewer pueruli recruiting to the fishery and lower productivity from Bicheno compared with the other sampling areas studied here. This is actually the case. Average counts of puerulus at Bicheno are substantially higher than the other sites (Gardner, unpublished data), yet the Bicheno region has the lowest fishery productivity in Tasmania (Gardner et al., 2015). Therefore, it is possible that in J. edwardsii, smaller pueruli are less fit and less likely to recruit into the fishery.

| Genetic differentiation at the southernmost site
Lack of neutral genetic differentiation between sampling events suggests that populations sourcing pueruli into the east coast of Tasmania during winter months are genetically homogenous; however, there is a slight indication of genetic differentiation at the southernmost site.
Lack of overall genetic divergence but significant differences at smaller spatial or temporal scale is typical of chaotic genetic patchiness (Eldon et al., 2016), which is inline with previous evidence of chaotic genetic patchiness in spiny lobsters (reviewed by Silva et al., 2019).

TA B L E 6
Results of the Mantel tests between temporal and genetic distance at each settlement site using the neutral SNP panel. pueruli during the winter months. The slight differences in individual weights between months found and size-at-settlement could be a signature of each of these populations. However, the small sample size of the present study precludes the possibility of proving that the observed genetic divergence at Recherche Bay is due to chaotic genetic patchiness or IBT.

| CON CLUS IONS
The observed variation in size-at-settlement of pueruli observed in the present study is likely due to phenotypic variability during dispersal, rather than the genetic signature of the parental stock. writing -review and editing (equal). Carla A. Souza: Data curation (equal); formal analysis (supporting); investigation (supporting). Jan Strugnell: Conceptualization (equal); formal analysis (equal); investigation (equal); methodology (equal); project administration (lead); resources (lead); writing -review and editing (equal).

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

DATA AVA I L A B I L I T Y S TAT E M E N T
Genotype (VCF) file and reference loci sequences are available through Dryad, doi:10.5061/dryad.5mkkwh792.