Phylogeography and genetic diversity of the microbivalve Kidderia subquadrata, reveals new data from West Antarctic Peninsula

It is well established that Antarctic biodiversity has been strongly influenced by rapid climatic fluctuations during the Quaternary. Marine invertebrates from Antarctica constitute an interesting lens through which to study the impacts of the last glacial periods as glaciation impacted the distribution and intraspecific genetic variation of these animals. However, the impact on the spatial genetic distribution and historical demography of local processes in areas adjacent to the West Antarctic Peninsula (WAP) is less clear. Here we present new genetic information on the bivalve Kidderia subquadrata, a small mollusk that inhabits intertidal rocky island ecosystems throughout the WAP. Using a phylogeographical approach, we examined the spatial patterns of genetic diversity in this brooder species to test the hypothesis of strong genetic structure in incubating organisms and the hypothesis of glacial refugia in organisms with limited dispersion. We found evidence of strong genetic structure among populations of the WAP and a recent expansion in the South Shetland Islands. Our findings are concordant with the predictions that incubating organisms, abundant in Antarctica, present a strong genetic structure among their populations and also support the hypothesis of glacial refugia in organisms with limited dispersion. The effect of the coastal current pattern in the WAP is suggested as a driver to the local spatial dynamics of the genetic diversity distribution. Although genetic information about this microbivalve is still scarce, the knowledge reported here has increased our understanding of the evolutionary patterns of this organism that is endemic to the Southern Ocean.

www.nature.com/scientificreports/ has been defined as an in situ refuge associated to volcanoes or areas of geothermal activity 18 . These refuges have been particularly informative in reconstructing the periglacial and postglacial history of marine organisms 19,20 . Finding a pronounced geographical structure may provide evidence of putative refugial areas [21][22][23] . Refugial populations are usually composed of subsets of the genetic diversity and long-term isolation of populations within geographically separate refugia will lead to genetic differentiation 14,24 . Despite the fact that several genera of bivalves have survived the glacial periods, a reduced number of these persist in Antarctic benthic systems and their biogeography patterns have yet to be exhaustively studied. The review of Seymour Island Paleocene molluscs 25 lists 57 species and 41 genera, of which only 12 genera are still represented in the SO; most of the remaining genera occur at present in the seas north of the Polar Front 26 and a few still live around Antarctica (12.5%) 25 . Moreover, the analysis of species richness in bivalve families revealed that the distribution has high richness north of the Antarctic peninsula and low richness to the south in, for example in families like Pectinidae, Nuculidae, Mytilidae, Gaimardiidae [26][27][28][29][30][31][32][33][34][35] . Antarctic bivalves have a varied range of reproductive strategies. In some the sexes are separate (gonochoric organisms) while others are simultaneous hermaphrodites 36,37 , some display indirect development and external fertilization (planktonic larvae) 38 while others exhibit direct development in which the embryos are retained by the female, who provides parental care to the offspring (brooder species) [39][40][41][42][43][44][45] .
The microbivalve K. subquadrata Pelseneer (1903) 46,47 is a brooder species with a 4.5 mm average shell length 47 . They incubate their offspring until the juvenile stage and have been reported inhabiting the rocky intertidal of islands adjacent to the West Antarctic Peninsula (WAP) 48 . No genetic information on this species has been reported prior to this work. In this paper, we aim to increase the available information about the Antarctic biodiversity and microbivalve species. To that end, we have developed a spatial analysis of the genetic diversity in populations of K. subquadrata from the WAP to advance the understanding of the evolutionary history of these microbivalves that inhabit Antarctica.
Based on the absence of the larval phase, we hypothesized that due to a low dispersal potential in this species we would find a high spatial genetic structure among populations with evidence of in situ refugia in the SSIs. To test this hypothesis, we used the mitochondrial molecular marker, Cytochrome oxidase I (CoxI) and a suite of phylogeographic analyses to estimate the spatial genetic diversity and the historical demography of this species.

Results
Genetic diversity and spatial genetic structure. The study included seven localities across the WAP ( Fig. 1). No insertion/deletions or stop codons were detected in the sequence dataset. In total 37 haplotypes were recovered (Table S1), with a whole dataset haplotype diversity of 0.722 ± 0.032 and 51 polymorphic sites ( Table 1). The pairwise FST distance revealed a high diferentiation among sites, with Doumer island being high and significantly differentiated from the other locations (Table S1). The IBD-Mantel test indicated that the relationship between the geographic distance and the linearized genetic distance was not significant among islands (r = 0.196; p-value = 0.25). The general AMOVA showed high genetic structure among populations (F ST = 0.87, p-value = 0.0001). The highest values of genetic differentiation measured with SAMOVA analysis were detected when the populations were separated into 4 groups (Group 1: Signy Island; Group 2: King George, Penguin, Greenwich, Deception Islands; Group 3: Livingston Island and Group 4: Doumer Island) with values of F CT = 0.92 (p-value = 0.029) ( Table 2). The model based on the Bayesian clustering algorithm and implemented in Geneland detected three clusters for the dataset (K = 3), with a high posterior probability of cluster membership (p-value = 0.9). Thus, demonstrating the strong genetic structure of K. subquadrata populations from the WAP. Cluster 1 (Fig. 2) includes samples exclusively from Signy Island; Cluster 2 includes samples from King George, Penguin, Greenwich, Deception and Livingston Islands; and, finally, Cluster 3 includes samples from Doumer Island.
The genealogical reconstruction of the haplotype network comprised 37 different haplotypes and showed a central haplotype (H1, Fig. 3) with the highest frequency (49.6%) that is distributed among five islands (Signy, King George, Penguin, Greenwich and Deception Islands) in a star-like topology. Additionally, most of the islands presented unique haplotypes with frequencies of 0.6% (Signy), 3.4% (Deception), 3.9% (Penguin), 5.6% (King George) and 7.8% (Greenwich). Ten mutation steps separated H1 from the samples of Livingston Island, where 4 unique haplotypes (3.4% of the total haplotypes) were detected. Another 10 mutational steps separated these groups from the second most frequent haplotype (H2, 19.5%); this group include exclusively samples from Doumer Island, where 2 unique haplotypes were found. A unique haplotype from Signy Island (H3) is separated from those of Doumer Island by one mutational step (Fig. 3).
Historical demography. The Tajima D test and the Fu Fs statistic (Table 1) were negative but no significant for the data set as a whole (Tajima D = -1.027, p-value = 0.2; Fu´s test = -1.906, p-value = 0.2). However, both indexes were both negative and significant for most of the individual islands. Signy, Deception and Doumer Islands were the exception, with Signy island showing positive but not significant values. The Bayesian Skyline plot analysis, where Signy, Penguin, King George, Greenwich, Deception Islands were included (start-like topology in the haplotype network) supports the hypothesis of a recent population expansion of K. subquadrata (Fig. 4). Based on this analysis, the time of the most recent common ancestor (trmca) was 5500 years ago, while the onset of population expansion occurred approximately 65,000 years ago.
The pre-evaluation step of the ABC procedure performed in DIYABC allowed to reduce the number of scenrios to test and improve DIYABC ability to reveal the true demographic model (Fig. S2-S4). The final analysis (Fig. 5) suggested both scenarios were realistic. The scenario 1 and 2 supports an initial divergence between Doumer and Signy islands and posterior admixture in which the population of Livingston island originates. The scenario 1 reveals an admixture event between Signy and Greenwich Islands giving birth to admixed populations www.nature.com/scientificreports/  www.nature.com/scientificreports/ in SSIs (King George, Penguin and Deception Islands). The scenario 2 supports a recent divergence in SSIs (Fig. 5). When scenarios were compared, posterior probabilities for the scenario 1 based on direct estimations reached the higest value of 0.8220 (0.4867,1.0000) while based en logistic approach, the scenario 2 reached the value of 0.9212 (0.9067, 0.9358). Type I error for both scenarios were low (scenario 1: Type I error = 0.038 for direct estimation; 0.009 for logistic approach; scenario 2: Type I error = 0.0039 for direct estimation; 0.010 for logistic approach).

Discussion
Our analysis supports the genetic differentiation among sampled populations of K. subquadrata, in accordance with the prediction that in the WAP, this brooder species have limited dispersal potential and would therefore display a highly structured population 2,49,50 . Reproductive life history might have a crucial influence in the genetic structure patterns and consequently in the degree of connectivity among populations in Antarctic organisms. For example, a study of two common Antarctic benthic surface-grazing gastropods with contrasting development strategies revealed a high spatial structure among populations in the brooder Margarella antarctica while no spatial population structure was found in the broadcaster Nacella concinna (planktonic larva) throughout the geographical range studied 51 . Similar results were reported in other works on this Antarctic limpet; these  www.nature.com/scientificreports/ studies, whose sampling sites covered a large geographical range, suggest the existence of a single panmictic unit 17,52 . More recently, for sea star brooders, researchers have concluded that broadcasters are less spatially structured than brooders 53 . The local oceanographic dinamyc may have also played a relevant role in the spatial interaction among populations in the interglacial periods in the SO. The spatial genetic structure pattern of K. subquadrata may have been enhanced by the Bransfield Strait current system. The surface circulation in Bransfield Strait is influenced by the Antarctic Coastal Current (ACC) and the Antarctic Slope Front 54 . Seaward of the South Shetland Islands (SSIs), the southwestward slope currents persist over the long term, they are originated by the Antarctic Slope Current that flows around the Antarctic Peninsula 55 . The water masses of the ACC transport nutrients west of the SSIs and enter the area through the Bransfield Strait flowing strongly to the east-northeast along the southern margin of the SSIs 56 . The Weddell Sea water enters through the Bransfield Strait mainly from the north-east and flows in a south-west direction along the Antarctic Peninsula. A similar process has been proposed to explain a genetic barrier identified in the Antarctic annelids of the family Phyllodocidae family in the strait north of Deception Island 57 , which coincides with the oceanic front that is generated by the intrusion of seawater masses from the Weddel Sea in the Bransfield Strait 55 . Therefore, the structure in K. subquadrata populations could be explained by the system of currents in the Bransfield Strait that forms a cyclonic circulation 58,59 . In addition, Livingston Island showed the presence of unique haplotypes that are separated by more than 10 mutational steps from the ancestral polymorphic haplotype but lie geographically lying within the SSIs. This undoubtedly requires a more detailed explanation and further investigation. According to a modified map from Barlett (2018), Livingston Island has low bathymetry (100-200 m) 60 . This pattern could result in water retention processes that increase the retention of marine invertebrate in the island, thus acting as a sink that results in a strong spatial genetic differentiation of its K. subquadrata population. However, the number of individuals analyzed in this study is  www.nature.com/scientificreports/ small and an increased sample size would provide a more complete picture of the spatial process in this locality and could reveal the actual phylogeographic pattern. On the other hand, glacial activity in Antarctica has been suggested as a driver of phylogeographic patterns for species by causing isolation and the presence of in situ glacial refuges 2,12-14, 61 . It now appears that the continental shelf was not ice-covered equally across the Antarctic coastline, which allowed some ice-free refuges for fauna during the glacial maxima 62 . According to the literature, the last glacial maximum occurred 20,000 to 17,000 years ago and West Antarctica and the adjacent islands were apparently covered in ice 62,63 . The deglaciation of Signy Island and the SSIs began 14,000-11,000 years ago and spread during the Holocene until 8,000-6,000 years ago [64][65][66][67] . Strong signals of recent demographic expansion as well as founder effect associated with recolonization can be inferred from our results in the haplotype network and DIYABC analysis. Therefore, a potential historical demographic hypothesis could be that strong genetic subdivision is a consequence of a pattern of multiple glacial refugia with a Pleistocene post-glacial expansion 19,21 .
In particular, the most widely distributed and frequents haplotypes were shared by five islands: Signy, Penguin, King George, Greenwich and Deception. The Signy and Doumer Islands are the most geographically distance localities within the study's sampling range. The DIYABC reveals that Signy and Doumer are the ancestral populations and probably after the original split, Signy island constituted a source population (refuge) for the population of SSIs. Interestingly, the haplotype network showed that Signy also has a unique haplotype that is different from the ancestral polymorphic haplotype and more closely related to the unique Doumer haplotypes. Some individuals of K. subquadrata that were collected as part of our samples lived as epibionts on the macroalgae. www.nature.com/scientificreports/ Macroalgal rafting has been suggested to explain the low genetic differentiation of marine communities across the Subantarctic region 68 , and recent investigations have revealed the arrival of invasive species that reached the Antarctic continent alive on kelp rafts 69,70 . Thus, a potential explanation for this unusual spatial pattern of the genetic distribution in K. subquadrata populations could be that they have used macroalgae as a transport vector. This hypothesis implies some degree of successful colonization and posterior isolation of this unique haplotypes.
In the case of Doumer Island, the southernmost locality, we identified exclusively private haplotypes, they displayed low genetic diversity and the highest F ST values, suggesting a process known as leptokurtic dispersion [71][72][73] .
Since that time it has probably served as its own refuge and undergone divergence in isolation 19 . Similar patterns of genetic structure have been reported in the marine bivalves Arctica islandica 74 in the Northern Hemisphere. Non-pelagic development has been described in the Antarctic fauna, where brood protection appears dominant 75 . It is worth noting that cryptic species are commonly found in brooding species in the Antarctic Peninsula in, for example, echinoderms such as ophiuroids, echinoids and brooding pygnogonids 44,[76][77][78] . The study of the bivalve Lissarca notorcadensis was one of the first to show indications of cryptic speciation in Antarcic bivalves 27 . New evidence of two cryptic lineages in the Antarctic Peninsula was recently reported in the bivalve Aequiyoldia eightsii, with an estimation of genetic distance of 5.79% from CoxI 35 . Considering K. subquadrata's intrinsically low dispersal and the evidence of the probable occurrence of in situ refuges revealed in this study, a cryptic speciation process could potentially explain the genetic diversity in the population of Doumer Island. The samples obtained from Doumer Island were collected in a protected bay surrounded by glaciers, semi-isolated from the marine currents since there is only one entrance to the bay; these geographical features indicate that gene flow out of bay could be difficult. Here, based on CoxI we report a genetic distance of 2% among the central haplotype (H1) and the two unique haplotypes from Doumer Island; this value is within the limit established in the existing literature 79 for inferring a speciation process. Intraspecific divergence with CoxI is rarely greater than 2% in the phylogeographic analyses; even so it can serve as an effective tool in recognizing putative new lineages 61,73,80,81 . Moreover, there are reports of distinct glacial refugia in the SSIs and Antarctic Peninsula harboring cryptic species that have diverged recently in micro-allopatry [82][83][84] .
In conclusion, we report new evidence of a strong spatial genetic structure in a brooder microbivalve species unique to the WAP. We also suggest possible presence of in situ glacial refuges and we infer a cryptic speciation in progress. Mitochondrial DNA has been widely used in numerous studies of phylogeography 5,16,17,21,27,35,42 ; its greatest utility is in non-model species where there are no previous genetic data, since its use enables comparisons with previous studies. However, recent studies have reported significant dissimilarities between genetic diversity and structure encountered in marine species depending on the genetic markers used 85 . Differences in patterns of genetic structure have been linked to the fact that organelle DNA is more sensitive to introgression and/or rapid sweeps (due to selection or strong genetic drift) than is nuclear DNA 86 . Therefore, it is imperative to develop new nuclear markers (e.g., microsatellites or SNPs) in Antarctic marine invertebrates. Further studies could increase the number of populations sampled in the Southern Ocean to test for the existence of local adaptation using recently developed genomics and transcriptomics tools. DNA extraction, amplification, sequencing and alignment. The total DNA of each individual was isolated with Quick-DNA plus (Zymo Research) commercial kit following the procedures described by the manufacturer. Given the absence of genetic information for the genus Kidderia, the strategy to generate molecular markers was to perform an Myseq Illumina sequencing using the genomic DNA and an in silico enrichment to obtain fragments of the mitochondrial genome. The total length of the partial mitochondrial genome recovered was 1733 bp. Using this genome data, specific primers for the Cytochrome oxidase I (CoxI) gene two pair primers were designed. Two pairs of primers were used: the first corresponds to kg1F (5′-TTG GGC TGG GTT AAT AGG TACA-3′) and kg7R (5′-GAA AAC CAG CAA ACA TAG CA-3′) flanking a fragment with a total length of 861 bp; the second pair corresponds to kg7F: (5′-TGC TAT GTT TGC TGG TTT TC -3′) and kg8R (5′-CCC AAA AAG ACA TTT GAC CC -3′), which were used to recover a fragment of 279 bp. The CoxI gene amplified a total length of 1140 bp, coding 380 amino acids.

Methods
All PCRs were performed in a final volume of 25 μL containing 1X PCR buffer, 3.5 mM MgCl 2 , 0.2 mM each dNTP, 0.25 μM each primer, 0.2X BSA, 1 μL DNA concentrate (10-50 ng), 0.6 units GoTaq DNA polymerase (Promega) and H 2 O to reach the final volume. The CoxI gene was amplified using the following thermocycling profile: an initial denaturation step (

Spatial pattern of the genetic diversity.
To estimate the levels of genetic polymorphism in populations of K. subquadrata we used standard diversity indexes: number of haplotypes (k), number of segregating sites (S), haplotype diversity (H) and the average number of paired differences between sequences (π) for each region using DnaSP v5 88 . Pairwise genetic differentiation (pairwise F ST ) was estimated to determine the genetic structure among populations using the ARLEQUIN program v.3.51 89 . To test the hypothesis of isolation by distance (IBD) we used the relationship between genetic (FST⁄1-FST) 90 and the geographic distances among Signy, King George, Penguin, Deception, Greenwich, Livingston and Doumer Islands to perform a Mantel test 91 executed in the R environment with the Vegan (2.5-2 version) package 92 .
To test for spatial genetic differentiation, we calculated a global FST by performing a global AMOVA analysis with the ARLEQUIN 89 . Next, we developed a set of analyses using the geographic coordinates of each sampled island. Here, we analyzed the population structure without an a priori cluster hypothesis using a spatial analysis of molecular variance (SAMOVA) 93 . In this approach, sample sites are clustered based on a simulated procedure that aims to maximize the proportion of total genetic variance caused by differences among populations groups.
To infer the spatial pattern of genetic diversity in K. subquadrata we estimated the number and composition of panmictic groups and the spatial boundaries among them using a clustering method based on the Bayesian model computed with the GENELAND package, version 4.0.0 94 in the R environment version 2.4.1 92 . This software implements a Markov Chain Monte Carlo (MCMC) procedure to determine the best clustering of samples with regard to genetic and geographic information. Geographic information is considered as the Bayesian prior level, so clusters corresponding to spatially structured groups are considered to be more probable than clusters that are randomly distributed in space. Five × 10 6 MCMC iterations were sampled every 1000 steps with a 500step burn-in period, and a maximum number of clusters K = 8 was run to estimate the model parameters and posterior probabilities of group membership.
Historical demography. To examine the historical population demography, the Tajima D test 96 and the Fu Fs statistic 97 were performed. The Tajima D test is based on the fact that in a neutral model, estimates of the number of sites segregating and the average number of nucleotide differences are correlated. Fu's test is based on the model of infinite sites without recombination; it gives the probability of observing a random sample with the number of alleles equal to or less than the value observed, given the observed level of diversity and the assumption that all alleles are selectively neutral. Both statistics were calculated in the ARLEQUIN. In addition, we estimated past population dynamics over time using the Bayesian skyline plot method implemented in BEAST 2 98 . To develop this analysis, we selected the five sample sites (Signy, Penguin, King George, Greenwich, Deception Islands) that shared the most frequently observed haplotype; as described earlier, these sample sites showed a star-like topology in the haplotype network. We conducted two independent Bayesian MCMC runs using the generalized time-reversible (GTR) model with a gamma distribution (G), previously estimated with JModelTest2 99 and mutation rates calibrated for CoxI sequences of bivalves 100,101 1.0% Myr -1 . Substitution rates were modified to a tenfold evolutionary rate (10% per million), considering the correction for time dependence of molecular rates at the population level 102,103 . The two independent MCMC calculations were run for 1.5 × 10 7 generations (sampled every 1000 iterations), and the first 10% of the trees were discarded as burn-in. The convergence of runs was confirmed with Tracer v1.6 104 , ensuring a minimum of 600 effective samples for each statistic (ESSs). The median and corresponding credibility intervals of the Bayesian skyline plot were depicted with Tracer v1. 6.
In order to better understand the population history of the species we used the program DIYABC v2.1 105 to tested for different population demographic scenarios. This software evaluates population histories using Approximate Bayesian Computation (ABC) with genetic data, by testing scenarios that are built through a combination of population divergence, admixture and population size changes. Following the recommendations of Cabrera and Palsbøll 106 to improve DIYABC's ability to reveal the true demographic model, we focused on simple contrasting models and reduced the number of candidate scenarios after a set of hierarchical analysis beginig with 5 preliminary scenarios (Fig. S1). Subsequently, and following the results of the previous genetic analysis, two models were evaluated at the last level (Fig. 5). The mutation rate was set with a mínimum of 1 × 10 -9 and máximum 1 × 10 -4 and the mutation model used was kimura 2 parameters. For the historical models, priors were set by default and in accordance with the recommendations of the authors of the software, we performed 2,000,000 simulations.