Genomic analyses indicate resilience of a commercially and culturally important marine gastropod snail to climate change

Sampling

Eight rocky shore locations (0 to 5 m depth) were selected for sampling (Table 1) spanning the known range of T. militaris from Hastings Point (northern NSW) to Jervis Bay (southern NSW) (Fig. 1), representing a seven-degree latitudinal and a 3.9 °C annual mean sea surface temperature gradient. Additionally, to sample a known cross-shelf gradient of temperature (less confounded by distance and unconfounded by latitude), three locations at varying distances from shore and under different temperature regimes were selected: Nambucca Heads (mainland), Split Solitary Island (3 km from the mainland) and South Solitary Island (7 km from the mainland). These represent inshore, mid-shelf, and offshore positions, respectively, with an annual mean sea temperature varying by about 1 °C from inshore to offshore (Malcolm, Jordan & Smith, 2010). Notably, South and Split Solitary islands are remote from shore and are also likely to have very low to no harvesting pressure. Between 25 and 30 T. militaris individuals were collected from each location (Table 1). As internal foot muscle tissue was found to yield clean, high molecular-weight DNA, it was necessary to narcotise specimens. Field collections were carried out under a scientific collection permit granted by the NSW DPI Fisheries under section 37 of the NSW Fisheries Management Act 1994 (#P01/0059(A)-4.0).

DNA extraction and genotyping

Approximately 25 mg of foot muscle tissue was sampled from each specimen using sterilised scalpel and forceps, avoiding the inclusion of mucous-rich epidermal tissue. Tissue samples were immediately placed in 2.0 mL Eppendorf snap-lock microcentrifuge tubes containing 500 µL of hexadecyltrimethylammonium bromide (CTAB) lysis buffer (100 mM TrisHCl, 20 mM EDTA, CTAB 2% w/v, NaCl 1.5 M) and refrigerated at 4 °C for 2 weeks. Proteinase K (30 µL at 20 mg/mL) was added to the samples which were incubated in a Allsheng shaking incubator at 60 °C, 200 rpm overnight. Samples were cooled to room temperature and purified by addition of an equal volume of 24:1 chloroform isoamyl. After centrifugation, the aqueous phase was retained and DNA was precipitated by the addition of 800 µL of dilute CTAB buffer (100 mM TrisHCl, 20mM EDTA, CTAB 2% w/v) to each tube which were incubated in an Aosheng MSC-100 shaking incubator at 60 °C at 400 rpm until DNA/CTAB complexes were visible. After centrifugation, the DNA formed a pellet which was twice cleaned with 70% ethanol by repetitive inversion. The DNA pellet was air-dried to remove residual ethanol and subsequently resuspended in sterile lab-grade water. Extracted DNA was quality checked using NanoDrop, Qubit assay and gel electrophoresis.

For single nucleotide polymorphism (SNP) genotyping, 20 μL of extracted DNA was sent to Diversity Arrays Technology Pty Ltd (Canberra, Australia) (DArT). The DArT organisation provides a process pipeline of whole-genome profiling, without the need for a reference genome (Jaccoud et al., 2001). High-throughput DArTseq technology was used to genotype Turbo militaris DNA. Here, the PstI-based complexity reduction method (Wenzl et al., 2004) was applied for the enrichment of genomic representation with single copy sequences. This method involved the digestion of DNA samples with a cutting enzyme PstI, paired with a set of secondary frequently cutting restriction endonucleases, ligation with site-specific adapters, and amplification of adapter-ligated fragments. Post digestion with a restriction enzyme pair, a PstI-overhang-compatible oligonucleotide adapter was ligated, and the adapter-ligated fragments were amplified in adherence to standard protocol (Wenzl et al., 2004). To develop SNPs, the DArTseq technology was optimized using two PstI-compatible adapters corresponding to two different restriction enzyme overhangs. The genomic representations were generated following the procedures described by Kilian et al. (2012). Next-generation sequencing technology was implemented using HiSeq2000 (Illumina, USA) to detect SNP markers. Sequence data was analysed using DarTsoft14 and DArTdb (Kilian et al., 2012).

SNP calling

In total, 208 individual T. militaris were initially genotyped with the DArTseq^™ platform yielding a total of 19,837 SNP loci with a mean read depth of 15.31 and 10.73% missing data. To improve SNP quality, while optimising the number of loci available for population genomic analyses, quality control filters, based on the descriptive statistics from the DArTseq^™ pipeline, were applied to data using the R package dartR v.2.7.2 (Gruber et al., 2019; Mijangos et al., 2022). Prior to SNP quality control, we checked for the presence of genetically related individuals, as their inclusion can lead to biased genetic estimations of downstream analyses. We calculated a similarity genetic distance matrix for individuals on the proportion of shared alleles per pairs of individuals with the function gl.propShared in dartR. Two individuals, one from Nambucca Heads (NAM07) and one from South Solitary Island (SSI16) were characterised as closely related and were removed from downstream analyses.

Two SNP datasets were generated for analyses: dataset ‘1’ for the analysis of overall genetic structure, and dataset ‘2’ for tests of GEAs and the identification of candidate loci. Both datasets were generated by retaining a single SNP per tag, removing secondaries, applying a locus and individual call rate of 80%, reproducibility threshold of 80%, and a Hamming distance threshold of 0.2 to control for the influence of linkage disequilibrium between loci. SNPs were called for dataset ‘1’ by applying a minor allele frequency (MAF) threshold of 0.03 and removing all loci departing from Hardy–Weinberg expectations. In contrast, SNPs were called for dataset ‘2’ by setting MAF to 0.01, and not filtering out SNP loci deviating from by Hardy–Weinberg expectations (all loci included). Finally, we also used poppr (Kamvar, Tabima & Grunwald, 2014) in the R package to calculate the number of private alleles found in each population, and remove these SNPs using gl.drop.loc in package dartR. After filtering, a total of 3,527 and 6,852 SNP loci for 206 individuals (Table 2) were retained for data sets ‘1’ and ‘2’, respectively.

Tests for population differentiation

Several estimates of genetic diversity were generated using the poppr package in R, including observed (HO) and expected (HE) heterozygosity and allelic richness (AR) (Joop Ouborg, Angeloni & Vergeer, 2010). For AR, allele counts were rarefied by the minimum number of individuals genotyped using the allelic.richness command in the R package hierfstat (Weir & Goudet, 2017). Statistical differences in genetic diversity measures among sites was estimated using the Hs.test function in the R package adegenet (Jombart, 2008). Departures from random mating were calculated using F_IS (inbreeding coefficient) for the overall dataset and for each sample location using the basic.stats function in hierfstat.

Tests for population genetic structure were subsequently performed, by calculating global population differentiation (F_ST) with 95% confidence limits (Weir & Cockerham, 1984), and population pairwise measures of F_ST with significance determined using permutation (999) in the dartR. Multiple testing effects were corrected for using the Benjamini–Hochberg FDR procedure (Benjamini & Hochberg, 1995). An analysis of molecular variation (AMOVA) was performed in the R package poppr (Kamvar, Tabima & Grunwald, 2014), using a model that partitioned variation among sample sites and within sample sites, with significance based on a randomization test with 999 permutations. A mantel test of isolation by distance (IBD) was carried out using the gl.ibd (Rousset, 1997) function in the R package dartR (Gruber et al., 2019) with results visualised in a scatterplot. A discriminant analysis of principal components (DAPC) was performed using adegenet package in R (Jombart, 2008). The find.clusters function was used to detect the number of clusters in the population. The best number of subpopulations has the lowest associated Bayesian Information Criterion (BIC). A cross validation function (Xval.dapc) was used to confirm the correct number of PCs to be retained. Finally, we used sparse non-negative matrix factorisation (sNMF) implemented in the R package LEA (Frichot & François, 2015) in R. This algorithm estimates the genetic ancestry components for each sample. For this study, 15 runs were performed with α = 100 for each K value (one to eight). The selection of the best number of putative ancestral populations was guided by the cross-entropy criterion (wherein, for K, a plot of the cross-entropy curve formed a ‘knee’) and the results from the best run were visualised using the barplots function.

Environmental variables

Physical oceanographic data were downloaded from the Copernicus Marine Environment Monitoring Service (https://marine.copernicus.eu), using a 20-year historical time-series encompassing 2001–2020 (daily temporal resolution) matched to each sampling location. These variables included temperature, water flow and productivity: sea surface temperature (SST) at 0.05° spatial resolution; remotely-sensed eddy kinetic energy (EKE) at 0.25° spatial resolution; and remotely-sensed chlorophyll a concentration (CHLa) at 0.04° spatial resolution (Table 2). The native spatial resolutions of oceanographic variables were used when matching daily data to sampling locations. Mean values were calculated for each variable by averaging the daily data. In addition, the lowest and highest SST for each location was extracted, and an absolute temperature range variable (temp.range) was calculated by subtracting the minimum from the maximum. A pairs.panels scatter plot of matrices was generated in the R package psych to confirm a lack of (Pearson’s) correlation between variables (r² ≤ 0.8), which were subsequently used in GEA analyses (described below).

Genome–environment association (GEA) analyses

To detect putative genomic signatures of selection, tests for GEAs were conducted using two complementary models: latent factor mixed models (LFMM2) (Caye et al., 2019), and the Bayesian method available in BayPass V1.01 (Gautier, 2015). BayPass and LFMM accounted for potentially confounding allele frequency differences due to population structure in a mixed linear model framework, but in different ways: LFMM estimates GEAs when simultaneously correcting for population structure with latent factors, while BayPass uses a neutral covariance matrix constructed from population allele frequencies. The subset of environmental variables described above were used as predictor variables for both analyses.

Genotype-environment associations were explored with Baypass V1.01 (Gautier, 2015) under the auxiliary (AUX) covariate model (-covmcmc and -auxmode flags). The first core model (without the environmental data) was run to estimate a covariance matrix (Ω) of population allele frequencies, which is an approximation of genomic differentiation between populations caused by demographic history. In order to reach convergence and reproducibility of the MCMC estimates, five independent runs, each with a randomly chosen seed were performed using default parameters, except for: pilot runs length of 1,000 iterations, number of sampled parameter values of 1,000, and a burn-in period length of 2,500 iterations. Secondly, the average of the five covariance matrices were used as input for the auxiliary covariate model to detect evidence of an association corrected for population structure. Environmental variables were scaled using the “-scalecov” option and the same running parameters as the core model were applied. Finally, the strength of association between genotype and the covariates was assessed by calculating the average of the log-transformed Bayes Factor (BF) in deciban units (dB) for each locus and environmental predictor. Significance was determined following Jeffrey’s criterion for decisive associations (BFis ≥ 20) (Jeffreys, 1939).

Latent factor mixed models (LFMM2) tested for linear relationships between environmental variables and genetic variants with random latent factors using a least-square method. Population structure was inferred by estimating individual ancestry coefficients based on sparse non-negative matrix factorisation (SNMF) method implemented in the snmf function in the R package LEA v3.10.2 (Frichot & François, 2015). Ancestry coefficients were determined for 1–8 ancestral populations (K) by generating an entropy criterion that evaluates the fit of the statistical model to the data using a cross-validation technique (Frichot & François, 2015). The K with the lowest cross-entropy value using 100 repetitions for each K value was selected. Subsequently, the optimal factor, K = 1, was used to inform the LFMM to identify whether allele frequencies were correlated with any of the environmental variables. To increase the statistical power of associations, missing genotype data were imputed via the ‘impute’ function in the LEA package, using the most common allele frequency observed in each K with the method ‘mode’. Next, we used the function lfmm_ridge to compute a regularised least-squares estimate using a ridge penalty. Individual associations between each SNP frequency and each environmental variable were assessed using statistics test calibrated using genomic inflation factor (function lfmm_test). Corrections for multiple comparisons were applied with the Benjamini-Hochberg algorithm with a false discovery rate (FDR) threshold of 5% (Benjamini & Hochberg, 1995). Significance associations were determined using a threshold of 0.001, since the probability of finding a false positive result increases with lower thresholds (Ahrens, Byrne & Rymer, 2019).

The gradient forest (GF) algorithm was subsequently used to describe the strength of the associations of spatial, environmental variables and candidate loci and to map spatial patterns of allelic turnover in climate space (Ellis, Smith & Pitcher, 2012; Fitzpatrick & Keller, 2015). Gradient forest is a machine learning method initially developed to model the turnover of ecological community assemblages in relation to environmental gradients (Ellis, Smith & Pitcher, 2012). Recently, this method has been adapted as a landscape genomics toolbox, substituting allele frequencies at genetic loci for species to model allelic turnover in climatic space (Fitzpatrick & Keller, 2015). The turnover functions in gradient forest allow for inference of the environmental predictors driving observed changes in allele frequency (Fitzpatrick & Keller, 2015). Analyses were implemented in the R package gradientForest (Ellis, Smith & Pitcher, 2012), using a regression tree-based approach to fit a model of responses between genomic data and environmental variables (Capblancq et al., 2020). Specifically, adaptive genetic variation turnover were modelled on the seascape using the candidate SNPs (derived from LFMM, and BayPass) set as the response variables. The machine learning algorithm partitioned allele frequencies at numerous splits values along each environmental gradient and calculated the change in allele frequencies for each split. The split importance (i.e., the amount of genomic variation explained by each split value) was cumulatively summed along the environmental gradient and aggregated across alleles to build a non-linear turnover function to identify loci that are significantly influenced by the predictor variable (Ellis, Smith & Pitcher, 2012). The analysis was run over 500 regression trees for each of the four environmental variables with all other parameters at default settings. The cumulative goodness-of-fit among SNPs was represented as an R² value indicating how well a predictor explained changes in allele frequency and which predictors were most important in predicting genomic changes. The resulting multidimensional genomic patterns were summarised using principal component analysis (PCA), allowing the relative importance of predictor variables on allelic turnover to be visualised. Finally, using the top gradient forest model, we interpolated genetic composition and allelic turnover across the sampling range of eastern Australia.

Overall population genetic structure

Patterns of genetic diversity did not differ greatly across the eight sample locations (H_O = 0.186–0.212; AR = 1.630–1.640; Table 3). Most sites showed a weak excess of heterozygotes (F_IS = 0.048–0.118), however these estimates did not differ significantly from zero (p > 0.01; Table 3). Overall genetic differentiation was found to be significant, but weak and close to zero (global F_ST = 0.002, p < 0.001) indicating a lack of genetic structure among sampling locations (Table 3). These findings are further supported by weak, yet significant, estimates of genetic differentiation among all pairs of sampling locations (F_ST = 0.000–0.004, Table 4). AMOVA also indicated a lack of overall genetic structure indicating genetic variance attributed to differences among sites to be non-significant (0.212, p > 0.05) while the majority of variance was explained by genetic variation between individuals within sites (99.86, p > 0.05). Similarly, DAPC and LEA analyses indicated a lack of genetic structure, both identifying a single population cluster (K = 1). Finally, Mantel tests revealed no significant relationship between genetic differentiation and distance between sampled locations (R2 = 0.03356, p = 0.278) providing further evidence of panmixia.

Genotype x environment associations

BayPass and LFMM each identified a number of candidate SNP loci exhibiting significant genotype-by-environment associations for each of the environmental predictors tested (Fig. 2). BayPass detected between 0 and 2 SNPs with significant correlations (log₁₀(BF) > 20) for each of the environmental predictors, with zero overlap in candidates between predictors. In contrast, LFMM detected between 5 and 12 SNPs that were significantly correlated with each of the environmental predictors (Fig. 2), but with only two loci overlapping between the EKE and SST predictor variables. Concordance of candidate SNP loci between methods was low (three loci only; Fig. 2) but expected given these methods have varying sensitivities to detecting loci under selection, use different methods for controlling for demography, and adopt different association algorithms. Given that T. militaris exhibits panmixia, its lack of population structure was likely to influence the inference of omega matrix and K clusters on BayPass and LFMM respectively. Furthermore, LFMM tests for relationships between individual-based allele frequencies whereas BayPass is at population level. Additionally, sampling design is likely to have had an influence, wherein a strategic sampling design accommodated to environmental heterogeneity and spatial variation on a landscape is essential to potentially identify and validate patterns of local adaptations across natural populations. Overall, GEAs indicated that a larger proportion of SNPs were significantly associated with temp range (13) and CHLa (11), followed by EKE (7) and SST (5).

Gradient Forest modelling used the unique candidate adaptive loci detected by both LFMM and BayPass (Luo et al., 2021), with 39 in total, four in BayPass and 35 in LFMM, but three were shared between both analyses and one was shared within LFMM (variables EKE and SST) (thus 39—3—1 = 35 unique candidate loci). Gradient forest analyses showed 14 of the 35 candidate SNP loci to be significantly correlated with environment (R² values > 0; mean = 0.08, range 0.003–0.3). Overall, EKE, CHLa and SST were found to be the most important predictors of genomic variation, while temp range had less of an effect (Fig. 3). Turnover functions from the GF model show the weighted cumulative importance values and sharp turnovers for all environmental predictors, but again with temp.range having lower importance relative to all other predictors (Fig. 3A). Biplots based on the first two principal components captured approximately 99% of the total variation and point to EKE and SST as the most prominent drivers of genomic variation (Fig. 3B). A spatial depiction of genomic composition in multi-dimensional climatic space based on PCA is provided in the allelic turnover map (Fig. 3C). The map indicates that the turnover of putatively adaptive allelic variation tracks closely with latitude, with the genomic composition of northern and southern most sampling locations being distinct from those from geographically intermediate locations (Fig. 3D).

Evidence of panmixia

Our analyses point to a lack of genetic structure across the entire distribution of T. militaris indicating widespread gene flow along the eastern seaboard of Australia. Such genetic patterns are also found amongst other trochoid taxa with long pelagic larval phases which are expected to facilitate long distance dispersal (Berry et al., 2019; Díaz-Ferguson et al., 2010; Nikula, Spencer & Waters, 2011; Silliman, Grosholz & Bertness, 2009). Several other eastern Australian marine invertebrates also exhibit high gene flow facilitated by larval traits including the Crown-of-thorns sea star (Acanthaster spp.) (Pratchett et al., 2015), the surf bivalve Donax deltoides Lamarck, 1818 (Murray-Jones & Ayre, 1997; Miller et al., 2013) and the black sea-cucumber Holothuria (Mertensiothuria) leucospilota (Brandt, 1835) (Chieu et al., 2023). Here gene flow and population structure has been linked to the long-distance dispersal of pelagic larvae facilitated by a fast-flowing East Australian current (EAC). While knowledge of the reproductive biology or larval competency of T. militaris is poor, our results suggest that the species may also generate long-lived, planktotrophic larvae contributing to high levels of biological connectivity across its distribution (Cowen & Sponaugle, 2009). Overall, these findings are consistent with previous genetic studies on trochoid taxa indicating population admixture over vast geographical areas. In this case we have provided evidence of gene flow among potentially locally adapted populations spanning major environmental gradients.

Evidence of local adaptation

Despite a lack of overall genetic structure across the sampling distribution, significant genotype associations with heterogeneous habitat features were observed across the sampling distribution, including associations with annual mean sea surface temperature, sea surface temperature range, EKE and productivity (CHLa). Drift processes leading to neutral genetic structure are often suppressed in broadcast spawning marine organisms with large population sizes (Gélin et al., 2017; Palumbi, 2003; Pinsky & Palumbi, 2014), but numerous studies have shown that adaptive genetic divergences can still be established and maintained under strong selection pressure (Hendry, 2017; Nosil, 2012; Schluter, 2000). In fact, many studies have demonstrated that adaptive variation can be maintained despite high levels of gene flow specifically in marine invertebrates, including gastropod snails (Miller et al., 2019; Sandoval‐Castillo et al., 2018). However, these findings are based on correlative tests only and should be interpreted with caution, as controlled mechanistic experiments are needed to validate these patterns and drivers of putative adaptive differentiation (Savolainen, Lascoux & Merilä, 2013; Stinchcombe & Hoekstra, 2008). Also, while our sampling regime was designed to correct for geographical distance, we cannot rule out the possible influence of artefactual associations (i.e., false positives) and SNP associations with other environmental factors varying with latitude (such as light intensity, dissolved oxygen levels and rainfall-driven variation in salinity). Nevertheless, our findings are consistent with many studies of marine species which demonstrated genetically determined clines related to climatic variables in Australia and overseas (Poloczanska et al., 2013, 2016; Wernberg et al., 2016).

Resilience to future climatic challenges

The coastal waters of south-eastern Australia are a recognised climate change hotspot with warming occurring above the global average with changes enhanced by strengthening of the EAC (Cresswell, Peterson & Pender, 2016; Malcolm et al., 2011; Oliver et al., 2017). Predictions anticipate warming to continue to increase, generating deteriorating marine conditions for sessile taxa that are already at their thermal limits (e.g., Davis, Champion & Coleman, 2021). Within the EAC, the effect of warming and changes to circulation may also counter one another, with warming enhancing larval survival, but a strengthened current reducing larval supply to the coast by restricting cross-current larval dispersal (e.g., in lobster, Cetina-Heredia et al., 2015). Our findings suggest that T. militaris may have the capacity to adapt to future climatic challenges, assisted by both widespread gene flow across environmental gradients and the availability of standing genomic variation for selection to act upon. Findings of putative adaptive genetic differentiation associated with temperature, nutrients and ocean currents is a particularly important finding, suggesting that standing genetic variation may be available for selection to act on to counter future environmental change, assisted by widespread gene flow, high fecundity (Romolo & Trijoko, 2021) and highly probable short generation time in this species (Kimani, 1996; Romolo & Trijoko, 2021).

While our findings suggest that T. militaris is likely to be generally resilient to shifts in the physical ocean climate, aggregations in some areas may still be vulnerable to risks of maladaptation. Those most at risk are likely to include locally adapted populations, where projected local changes in climate are high, and connections to non-local aggregations are relatively weak (Hoffmann & Sgrò, 2011). Recent climate projections indicate that many low energy embayment habitats are likely to experience greater increases in SST than open coastal habitats (Guyondet et al., 2015; Scanes, Scanes & Ross, 2020; Vila-Concejo et al., 2007); Also, biophysical models suggest that the biological connections between low energy embayment and high energy open coastal habitats can be weak in some marine invertebrates from south-eastern Australia (Riginos et al., 2016; Treml et al., 2015). Consequently, it is possible that locally adapted aggregations from low energy embayment habitats may be most vulnerable to climate change effects, where gene flow is unlikely to assist local aggregations in adapting to warming sea surface temperatures via the migration of thermally adapted genotypes. In such cases adaptive management strategies might be needed, including the assisted migration of thermally adapted genotypes to populations showing signs of climate stress. Such approaches are being widely advocated as a tool for “climate proofing” threatened marine and terrestrial animal and plant communities (Aitken & Whitlock, 2013; Prober et al., 2015; Layton et al., 2020; Hoffmann, Miller & Weeks, 2021). Although the northernmost Hastings Point population is not genetically isolated, its position at the northern trailing edge may render it vulnerable to stochastic events, such as heatwaves (Ab Lah et al., 2018; Mamo et al., 2019). Vulnerability in trailing edge populations can be amplified by genetic impoverishment through loss of individuals, without replenishment through immigration (gene-flow), consequently exposing these populations to the risk of localised extinction (Clark et al., 2020; Coleman et al., 2020). Indeed, this species has undergone prior range shifts indicating that the trailing (warm) edge is likely to be vulnerable to ongoing warming.

Implications for fisheries management

Sustainable fisheries management requires information on factors that influence the resilience of individual fishing stocks to fishing pressure and environmental disturbance (Astles et al., 2006; Kenny et al., 2018). This includes understanding the geographic boundaries of biological populations and the recruitment potential of individual stocks persisting within and across these populations (Binks et al., 2019; Roughgarden, Iwasa & Baxter, 1985) and how this might be altered under climate change. In the case of T. militaris, the presence of a single panmictic population unit across its distributional range, generally means that the opportunity for recolonisation following depletion events (overharvesting or environmental disturbance) are enhanced for central and southern sub-populations. However, for the northern trailing edge, the opportunity for repopulation with genotypes from pools of genetic diversity further south may be hampered by the dominant poleward flow of the EAC. Furthermore, increasing human population and harvesting pressure may reduce local abundance (Cooling & Smith, 2015) and the opportunity for thermal adaption and repopulation.