Genetic variation in the small bivalve Nuculana inaequisculpta along a retreating glacier fjord, King George Island, Antarctica

. - Climate change is strongly influencing regions of Antarctica but the consequences on microevolutionary processes have been little studied. Patterns of population genetic diversity were analysed in the Antarctic bivalve Nuculana inaequisculpta (Protobranchia: Nuculanidae) from a fjord with 70 years of documented climate-forced glacier retreat. Thirty-nine individuals from five sites at different distances from the glacier terminus were collected, and the COI gene was sequenced from each individual. No statistically significant genetic differentiation was found between sites nor a significant correlation between the proximity of glaciers and genetic diversity, suggesting a high dispersal capability and therefore, a planktonic larval stage for this species. Nevertheless, we encourage increasing the sample size and number of loci in future studies to confirm our findings.


Introduction
West Antarctica is a hotspot for physical Climate Change, and ~90% of all glaciers along the West Antarctic Peninsula (WAP) are rapidly retreating (Cook et al. 2016). Marine species, especially those in the shallows, are subject to drastic environmental changes, including those which may negatively impact species survival (e.g., increased ice scour) or favor colonization of new available habitats exposed by retreating glaciers. Several studies have documented ecological impacts of glacial retreat and melting on marine communities (e.g., Moon et al. 2015, Sahade et al. 2015, Barnes et al. 2020, Bascur et al. 2020. However, it is less well known how such perturbations may drive microevolutionary responses of marine species at population genetic scales. There is a lack of studies examining this question in Antarctic ecosystems. Studies testing patterns of genetic diversity in relation to climate-change forced distributional shifts in other areas have found mixed responses. For instance, simulation studies have found that the edge of a shifting range, under climate change, has low levels of neutral genetic diversity (Cobben et al. 2011, Arenas et al. 2012. In contrast, some empirical data, for example of Larix decidua, have shown high genetic diversity at the leading edge of an expanding population following a retreating glacier in the Swiss Alps (Pluess 2011).
The genetic population response of a species after perturbation is also related to aspects of its biology, such as its dispersal capability. For instance, a species with low dispersal capability might be expected to present a stronger pattern of genetic diversity (e.g., decreasing genetic diversity along an expansion transect) than a species with high dispersal capability. This is because migration between demes is not fast enough to dilute the pattern (Klopfstein et al. 2006). In contrast, high dispersal capability would be likely to drive high gene flow between populations and therefore, rapid homogenization of any incipient genetic differentiation (Slatkin 1987).
In Antarctica, it has been argued that climatic conditions have filtered for life-history strategies that reduce or lack larval stages in benthic organisms, especially those which are planktotrophic (Pearse et al. 2009). However, planktonic larvae have been described for many Antarctic benthic species (e.g., Stanwell-Smith et al. 1997, 1999, and their actual number may have been greatly underestimated due to the lack of detailed morphological and molecular studies on planktonic larvae (Thatje 2012). Studies around Antarctica have shown that developmental type influences patterns of genetic differentiation across barriers. For example, Poulin et al. (2014) found that species lacking planktonic larvae were more divergent across the Antarctic Circumpolar Current than species with planktonic larvae. Similarly, Muñoz-Ramírez et al. (2020) found that Aequiyoldia eightsii, a species with a lecitotrophic larvae, showed genetic patterns that were congruent with ocean current dispersal. Dispersal capability of species likely influences colonization rates and genetic patterns on ecological time scales. However, no studies have examined the genetic consequences at the temporal scale of current climate change in Antarctica.
Nuculana inaequisculpta (Lamy, 1906) (Fig. 1D) is a small bivalve species that belongs to the Protobranchia subclass, a taxonomic group that contains the most ancient bivalves in the world (Yonge 1939, Villarroel & Stuardo 1998. It is an abundant infaunal species, patchy distributed in the South Shetland Islands and along the Antarctic Peninsula (Cattaneo-Vietti et al. 2000, Lovell & Trego 2003, Gordillo et al. 2017. This species lives in benthic, muddy marine habitats and feeds mainly on organic deposits in sediments (Zardus 2002). It occurs as deep as ca., 800 m, and range in size from 2.5-16.0 mm in shell length (Cattaneo-Vietti et al. 2000, Gordillo et al. 2017. In general, protobranchia species have two separate sexes, with a similar (1:1) sex ratio, high fecundity, and lecithotrophic larval development with a pericalymma larva (Zardus 2002). Recent bioenergetic analyses have found that Nuculana inaequisculpta are sensitive to food availability and exhibit lower fitness in recently colonized areas, perhaps due to lower food supply and higher perturbations near the glacier fronts (Bascur et al. 2020).
In this study, sequences of the Cytochrome C Oxidase Subunit I (COI) gene obtained from Nuculana inaequisculpta were analysed to investigate patterns of genetic diversity along a fjord with a retreating glacier. A strong pattern of decreasing genetic diversity towards the glacier front could be considered indicative of low dispersal capability as expected under a direct development life-history strategy or benthic larvae, and would suggest impacts of current climate change on population diversity. In contrast, a lack of correlation would indicate high dispersal capability, as expected for planktonic larval stages, and highlight the importance of high dispersion as a mechanism of resilience against climate-change induced perturbations.

Sampling sites and molecular data collection
Thirty-nine individuals of N. inaequisculpta were collected from five sites (6-8 individuals per site) at Marian Cove, King George Island, Antarctica, a fjord with a significant glacier retreat documented since 1950 (Park et al. 1998, Cook et al. 2016, Barnes et al. 2020 (Fig. 1

Genetic differentiation and diversity
To estimate levels of population differentiation, F st values between sites were calculated using the R-packages ADEGENET (Jombart 2008), and HIERFSTAT (Goudet 2005). Genetic diversity indexes, haplotype richness (S), haplotype diversity (Hd) and nucleotide diversity (π) for each site and for the entire fjord were estimated using PEGAS (Paradis 2010). To visually display relationships between haplotypes, a haplotype network was estimated using the median-joining network algorithm of Bandelt et al. (1999) implemented in the program POPART (Leigh & Bryant 2015).
Based on results of low genetic diversity for the closest site to the glacier terminus, a test was conducted to estimate the probability of obtaining the lowest π value for that site under a model of panmixia and chance. For this, a panmictic population was assumed and samples were randomly sampled to build a frequency distribution of 1,000 simulated π values, selecting six samples at each iteration to calculate π. Subsequently, the empirical π value was compared to the frequency distribution of random π values. Support for a significantly low value is given if the empirical π value is lower than the 95% of the simulated π values.

Results and discussion
A total of eighteen haplotypes were present among the thirty-nine individuals of N. inaequisculpta at Marian Cove. The average number of nucleotide differences between samples was 0.402. Genetic structure between populations was low. Genetic differentiation was greater between sites MC2 and MC5 (F st = 0.1057), whilst the lowest differentiation values were between MC3 versus MC1 and MC3 versus MC4 (both F st = 0.0279). However, none of the F st values between sites were statistically significant at the 0.05 level. A similar lack of structure was observed in the haplotype network ( Fig. 2A). A central, high-frequency haplotype was present at all sites at relatively high frequency. Most of the remaining haplotypes were derived from that, with low frequency and separated by only few mutational steps. This star-like topology is typical of populations that have undergone a recent bottleneck and subsequent population expansion (Slatkin & Hudson 1991). Living space on the Antarctic shelf has undergone cyclic expansions and contractions due to the advance and retreat of the Antarctic Ice Shield during glacial cycles (Augustin et al. 2004). These major historical events can explain the signal of demographic expansion, which could have initiated after the Last Glacial Maximum. However, this could also be explained, at least in part, due to recent colonization of the fjord. Data are not available to distinguish these hypotheses, but larger geographic sampling, including populations from the entire distribution of the species could help resolve this question in a future study.
Haplotype diversity (Hd) across sites varied from 0.643 (in MC1) to 1 (in MC2) (Fig. 2B). There was no clear tendency for Hd to increase or decrease with glacier proximity, and the lowest values were found at the opposites sites of the transect (Hd MC1 = 0.643; Hd MC5 = 0.8). Nucleotide diversity tended to increase with proximity to the glacier from sites MC1 (π MC1 = 0.0032) to MC4 (π MC4 = 0.0053), but it strongly decreased at site MC5, with MC5 presenting the lowest value among all sites (π MC5 = 0.0021) (Fig. 2C). No significant correlation was found between geographic proximity to the glacier terminus and genetic diversity for both Hd (r= -0.307; P-value= 0.6157) and π (r= 0.06; P-value= 0.9216), suggesting high dispersal capability for N. inaequisculpta, and a planktonic larval stage. However, the lowest π value found at MC5 (the closest site to the glacier terminus) suggests that strong perturbations (e.g., ice scouring, glacier meltwater discharge, high sedimentation) near the glacier terminus may impact genetic diversity. To test if this low value of genetic diversity could imply a deterministic force (e.g., impact of perturbations on genetic diversity), the probability of having the lowest π value at MC5 was estimated under population panmixia and chance. The probability of obtaining a π value equal or lower than the observed value (π= 0.00201) only by chance was P= 0.1897, and therefore, not significant. This indicates that

Figure 2. A) Haplotype network, B) Haplotype diversity (Hd), and C) Nucleotide diversity (π) for Nuculana inaequisculpta based on COI sequences.
Black dots and dashes in the haplotype network represent nodes and mutational steps, respectively / A) Red de haplotipos, B) Diversidad haplotípica (Hd) y C) Diversidad nucleotídica (π) para Nuculana inaequisculpta basada en secuencias COI. Los círculos negros y las líneas cortas perpendiculares sobre las líneas que conectan haplotipos corresponden a nodos y pasos mutacionales, respectivamente there is no enough evidence that the lowest π value at MC5 is due to deterministic forces (e.g., recent colonization of the area or perturbations). However, evidence of low fitness in Nuculana specimens from near the glacier terminus has been recently provided by Bascur et al. (2020), who found, by conducting bioenergetic analyses, that individuals from MC5 had the lowest energy content in their tissues. This congruence between bioenergetic and genetic results strongly suggests an impact of glacierderived perturbations on the species's performance.
The pattern of low genetic diversity towards MC1 (the farthest to the glacier terminus) is puzzling, and may suggest a more complex scenario including low competitive abilities. If N. inaequisculpta is a species that exploits its dispersal capability to reach recently emptied habitats, but is a generally weak competitor, populations that are present in longer-established communities (e.g., those which are further away from the glacier terminus) may be at disadvantage, and therefore maintain lower population sizes. This could explain low genetic diversity at MC1 and the progressively higher diversity towards MC4 (although not significant). Indeed, an additional site (MC6) furthest from the glacier terminus (not shown on the map) was also sampled, but it yielded no individuals of N. inaequisculpta, so it was not included in this study. This MC6 site was rich and diverse in other species though (data not shown), suggesting that the community was old and long established. This seems to support the idea that N. inaequisculpta could be a weak competitor given its absence in this site. However, it is recognized that this is speculative, and as a tentative hypothesis it will require further study. Therefore, it is proposed that the comparative genetic study of multiple species with different life history strategies (e.g., with and without planktonic larval stage) and competitive abilities will allow testing this hypothesis in future studies.
There are clear limitations to the inferential power of the marker used here (one single mitochondrial gene) and the low number of individuals used, so these results should be taken with caution. Nevertheless, the results provided in this study are a first step towards a more comprehensive understanding of the patterns and processes governing the genetic diversity of this bivalve species, which can be largely improved in future studies by increasing the sample size and the adoption of more powerful genomic markers such as RADSeq-derived single-nucleotide polymorphisms.