Population substructure and signals of divergent adaptive selection despite admixture in the sponge Dendrilla antarctica from shallow waters surrounding the Antarctic Peninsula

Antarctic shallow‐water invertebrates are exceptional candidates to study population genetics and evolution, because of their peculiar evolutionary history and adaptation to extreme habitats that expand and retreat with the ice sheets. Among them, sponges are one of the major components, yet population connectivity of none of their many Antarctic species has been studied. To investigate gene flow, local adaptation and resilience to near‐future changes caused by global warming, we sequenced 62 individuals of the sponge Dendrilla antarctica along the Western Antarctic Peninsula (WAP) and the South Shetlands (spanning ~900 km). We obtained information from 577 double digest restriction site‐associated DNA sequencing (ddRADseq)‐derived single nucleotide polymorphism (SNP), using RADseq techniques for the first time with shallow‐water sponges. In contrast to other studies in sponges, our 389 neutral SNPs data set showed high levels of gene flow, with a subtle substructure driven by the circulation system of the studied area. However, the 140 outlier SNPs under positive selection showed signals of population differentiation, separating the central–southern WAP from the Bransfield Strait area, indicating a divergent selection process in the study area despite panmixia. Fourteen of these outliers were annotated, being mostly involved in immune and stress responses. We suggest that the main selective pressure on D. antarctica might be the difference in the planktonic communities present in the central–southern WAP compared to the Bransfield Strait area, ultimately depending on sea‐ice control of phytoplankton blooms. Our study unveils an unexpectedly long‐distance larval dispersal exceptional in Porifera, broadening the use of genome‐wide markers within nonmodel Antarctic organisms.


| 3153
LEIVA Et AL. playing a key role by providing shelter and food for many other marine invertebrates (e.g., Moles et al., 2017). Its distribution spans along the Antarctic Peninsula and its associated islands, to the South Orkney Archipelago as the northernmost point of its range (data from World Porifera Database: www.marin espec ies.org/porif era/ porif era.php?p=taxde tails &id=164875). D. antarctica is a brooding sponge, with yolky lecithotrophic larvae that are released during the Antarctic summer . In the present study, we aim to assess the genetic diversity, demographic history, and genetic connectivity of D. antarctica at a regional scale in the WAP and South Shetland Islands using double digest (dd)RADseq-derived SNPs. We also evaluate the suitability of the full mitochondrial genome in D. antarctica to assess genetic diversity and connectivity.
Finally, we test for genetic signatures of divergent selection using SNPs identified in an F ST outlier test, and measure the expression levels of the genes identified under selection in three transcriptome samples spanning the whole latitudinal range of our sampling area.

| Sample collection, preservation and DNA extraction
For the population genomics study with ddRADseq, we collected ~1 cm 3 of tissue from 67 specimens of Dendrilla antarctica during the 2015-2016 austral summer in seven locations across the WAP and the South Shetland Islands (Figure 1; Table 1). Sampling was performed by SCUBA diving at 5-25 m depth. Sponge fragments were preserved in 96% ethanol, with the ethanol being replaced three times, and stored at −20°C until further processing. We extracted DNA from all samples using the DNeasy Blood & Tissue kit (Qiagen) following the manufacturer's protocol, with minor modifications in the cell lysis time (which was conducted with an overnight incubation) and the final DNA elution step (performed twice using 50 µl of elution buffer each time). DNA quantity was assessed with a Q ubit dsDNA HS assay (Life Technologies).
For mitogenome reconstruction, we collected a fragment (~1 cm 3 ) of a specimen in 96% ethanol from Deception Island to perform draft-level genomic sequencing, with genomic DNA (gDNA) extracted as described above. Furthermore, we subsampled tissue fragments (~1 cm 3 ) of three individuals for additional mitogenome reconstruction and transcriptomic analysis, from three different sampling stations (O'Higgins Bay, n = 1; Deception Island, n = 1; and Adelaide Island, n = 1), covering the whole latitudinal range of the sponge in our study. We preserved the subsampled tissue fragments in RNAlater (Life Technologies) immediately after collection, stored them for 24 hr at 4°C, replaced the RNAlater once, and then stored samples at −80°C until further processing.

| Transcriptomic and genomic library preparation and sequencing
For transcriptomics, total RNA was extracted using a standard trizol-based method using TRI Reagent (Life Sciences) following the manufacturer's instructions. Subsequent mRNA purification was

| Transcriptome and mitochondrial assembly, and mitogenome screening
A total of 123,782,845 paired reads (Den_ROT_3 = 29,840,417 reads, Den_OH_2 = 39,434,819 reads, and Den_DEC_19 = 54,507,609 reads) were obtained in our transcriptomic run. A total of 9,644,983 paired raw reads were obtained in our gDNA run, 8,484,436 paired reads after trimming and filtering. Both transcriptomic and gDNA reads were cleaned using trimmomatic 0.33 (Bolger, Lohse, & Usadel, 2014) with the following settings: ILLUMINACLIP:../Adapters. gDNA reads from the sample from Deception Island were assembled using velvet 1.2.10 (Zerbino & Birney, 2008) at k-mer sizes of 71 and 91, which were the best k-mers after optimization trials. A local blast database was made from these gDNA assemblies using the makeblastdb command (Altschul et al., 1997 above. In addition to the complete mitochondrial genome recovered from the gDNA sample, three more mitochondrial genomes were obtained from the transcriptomic reads or the assembled transcriptomes using the pipeline Trimitomics (Plese et al., 2018). Subsequently, all four mitochondrial genomes (one coming from the gDNA sample and three from the transcriptomes) were aligned in geneious 8.1.8 (Kearse et al., 2012) using the Q-INS-I algorithm of mafft version 7 (Katoh & Standley, 2013), which is used as the default algorithm for rRNA alignments because it considers secondary structure information, as a form of base-pairing probability (Katoh & Toh, 2008). The software dnasp 5.10.01 (Librado & Rozas, 2009) was used to calculate the number of segregating sites (S), haplotype number (H), haplotype diversity (Hd) and nucleotide diversity (π).

| ddRADseq library preparation and sequencing
Library preparation was conducted following Peterson, Weber, Weber, Kay, Fisher, and Hoekstra (2012) with the following modifications (as in Combosch et al., 2017). gDNA (100-1,000 ng) was digested using the high-fidelity restriction enzymes Sbf1 and EcoR1 (New England Biolabs

| ddRADseq locus assembly and outlier detection
Quality filtering and locus assembly was conducted with the stacks pipeline, version 1.44 (Catchen, Hohenlohe, Hohenlohe, Bassham, Amores, & Cresko, 2013). RAD-tags (DNA fragments with the two appropriate restriction enzyme cut sites that were selected, amplified and sequenced) were processed using process_radtags, where raw reads were quality-trimmed to remove low-quality reads, reads with uncalled bases, and reads without a complete barcode or restriction cut site. The process_radtags rescue feature (-r) was used to recover minimally diverged barcodes and RAD-tags (--barcode_dist 3; --adapter_mm 2). The process_radtags trimming feature (-t) was used to trim remaining reads to 120 bp, in order to increase confidence in SNP calling. After performing these filtering steps in we finally retained loci present in at least 40% of the individuals, and filtered out individuals with less than 30% of the final loci, resulting in a data set containing 577 SNPs and 62 individuals.
In order to differentiate neutral SNPs from putative SNPs under positive selection, the database containing 577 SNPs was analysed using default parameters in lositan (Antao, Lopes, Lopes, Lopes, Beja-Pereira, & Luikart, 2008). We used lositan because it implements the FDIST2 approach of Beaumont and Nichols (1996), which provides a robust method when populations deviate from the island model of migration (Tigano, Shultz, Shultz, Edwards, Robertson, & Friesen, 2017). Also, it incorporates heterozygosity and simulates a distribution for neutrally distributed markers (De Mita et al., 2013;Narum & Hess, 2011). We considered that these features were more appropriate for our model species studied herein, D. antarctica, than the characteristics of other F ST -outlier methods. Moreover, we also run bayescan (Foll & Gaggiotti, 2008) with default parameters, which only detected one locus under selection already detected by lositan.
We discuss lositan results given the reasons stated above.
Population structure was assessed using the function snapclust in the R package adegenet (Beugin, Gayet, Gayet, Pontier, Devillard, & Jombart, 2018), the discriminant analysis of principal components (DAPC) as implemented in the adegenet R package (Jombart, Devillard, & Balloux, 2010), and structure version 2.3 (Pritchard, Stephens, & Donnelly, 2000). structure and snapclust may produce similar individual membership probability plots, but they have totally different approaches to the genetic clustering problem: while structure uses a Bayesian approach with Markov chain Monte Carlo (MCMC) method to estimate allele frequencies in each clus- Both analyses were run for 200,000 MCMC iterations using the admixture model, with a burn-in of 100,000 iterations, setting the putative K (number of clusters) from 1 to 10 with 20 replicates for each run. structure harvester (Earl & vonHoldt, 2012) and clumpp version 1.1.2 (Jakobsson & Rosenberg, 2007) were used to determine the most likely number of clusters and to average each individual's membership coefficient across the K value replicates, respectively.
Pairwise F ST values were estimated to measure the differentiation between pairs of sampling stations using the pairwise.fst function in the hierfstat R package. Their significance was tested with 1,000 permutations using the ade4 R package (Dray & Dufour, 2007

| Annotation, expression values and structure of putative loci under selection
To improve the annotation step, all RAD-loci containing an outlier SNP under positive selection were mapped back to our de novo reference transcriptome using CLC Genomics Workbench 5.1 local blast (Altschul et al., 1997) to obtain the contig for each RAD-locus. All contigs with uniquely mapped loci were then subjected to a blastx and blastn search against nr (default parameters, Altschul et al., 1997), and this annotation was retained for the mapped RAD-locus.
A functional annotation analysis was then performed using david 6.8 (Huang Sherman, & Lempicki, 2008a, 2008b. Expression values for putative contigs under selection in each of the three RNA-seq data sets were determined by mapping reads from individual samples to the de novo reference transcriptome using the rsem package (Li & Dewey, 2011;--aln_method bowtie2) in Trinity (Grabherr et al., 2011 was also used on the loci under positive selection and, subsequently, two AMOVAs were also performed in the poppr R package to test whether the DAPC grouping or the barrier's genetic break explained a significant part of the genetic variation.

| Mitogenome diversity
The total length of the mitochondrial genome obtained from our gDNA reads from Deception Island was 19,498 bp, of which 10,949 bp comprised protein-encoding sequences. From our de novo transcriptomes, we were able to recover 10,859 bp for the sample from Adelaide Island, 10,682 bp for the sample from

| Population genetics analyses using neutral SNPs
Population genetics statistics are shown in Table 2  The AMOVA results for the neutral SNP data set are shown in Table 4a. Both the DAPC and the barrier groupings appeared to be nonsignificant portions of the genetic variance (p = 0.468 and 0.248, respectively), the two of them representing less than 1% of the total variation (Table 4a).  (Table 5).

| Putative loci under selection
One of them corresponded to an uncharacterized protein, and another one matched a bacterial aminotransferase (

| Mitogenome diversity
Our study unveiled an extremely low mitochondrial diversity in pair mechanisms (Huang, Meier, et al., 2008).

| Population genomic analyses using neutral SNPs
In contrast with the low variability of the mitochondrial genome in  (Table 2), indicating a deviation in the haplotype frequencies from the neutrality model (Tajima, 1989). These results support the existence of a recent and rapid demographic expansion of D. antarctica in the WAP and the South Shetlands, which could have started after the last glacial-interglacial alternation (~20,000-10,000 years ago) when the last Antarctic shelf recolonization took place (see Allcock & Strugnell, 2012). This hypothesis has been suggested for other shallow-water Antarctic invertebrates (Díaz, Féral, Féral, David, Saucède, & Poulin, 2011;González-Wevar et al., 2011;Leiva, Riesgo, Riesgo, Avila, Rouse, & Taboada, 2018;Thornhill, Mahon, Norenburg, & Halanych, 2008), which could have migrated northwards to sub-Antarctic islands during glacial periods and recolonized the Antarctic shelf during interglacial periods.
This expansion-contraction model has already been tested for the Antarctic limpet Nacella concinna (Strebel, 1908), demonstrating its glacial survival in the sub-Antarctic South Georgia Island, followed F I G U R E 4 Functional annotation analysis showing the 14 annotated loci under selection, their cellular function, and relative expression levels presented in heatmaps (actual values can be seen in Table 5). Red indicates higher expression, white moderate expression and blue low expression, based on the values reported in Table 5  Taking expected heterozygosity (H E ) as a measure of genetic diversity, as originally defined (Nei, 1973), we found significantly lower genetic diversity values for D. antarctica (  (Brown et al., 2017).
In other examples using SNPs in different animal phyla, H E ranged from 0.298 to 0.312 in the salmon louse Lepeophtheirus salmonis (Jacobs et al., 2018), from 0.211 to 0.214 in the Galapagos shark population decimations, as has already been reported for other shallow-water Antarctic fauna (see Allcock & Strugnell, 2012 (Table 3), suggesting high connectivity and dispersal capability of D. antarctica throughout the sampling area, which covered most of the species distribution. We propose that this could be due to the relatively long planktonic life of D. antarctica larvae as a result of the great amount of proteinaceous yolk that they contain  in comparison with sponge larvae from congeneric species from lower latitudes (e.g., Ereskovsky & Tokina, 2004). Furthermore, the strong oceanic currents in the study area (Moffat, Beardsley, Owens, & Van Lipzig, 2008;Zhou, Niiler, & Hu, 2002; see Figure 1b) may increase the dispersal ability of D. antarctica larvae. Remarkably, our results differ from most previous population genetic studies on sponges, which generally report highly structured and differentiated populations, even at local and regional scales (e.g., Brown et al., 2017;DeBiasse et al., 2010;Pérez-Portela, Noyer, & Becerro, 2015;Riesgo et al., 2016). Even compared to some oviparous sponges such as Cliona delitrix which appears to disperse along the ~315 km of the Florida reef track (Chaves-Fonnegra, Feldheim, Secord, & Lopez, 2015), our results suggested an unprecedent ~900 km contemporary migration. However, although this long-distance connectivity is unusual in sponges, it is common in other Antarctic marine invertebrates. Examples of high gene flow are shown in many Antarctic species, such as the brittle stars A. agassizii (Galaska et al., 2017a) and Ophionotus victoriae (Galaska et al., 2017b), the Antarctic limpet N. concinna (González-Wevar et al., 2013), the nemertean Parborlasia corrugatus (Thornhill et al., 2008) and the annelid Pterocirrus giribeti (Leiva et al., 2018).
However, both DAPC and barrier groupings appeared as a nonsignificant part of the total genetic variation in the AMOVAs (Table 4a).
This may be due to the high admixture and migration detected in the neutral data set (Figures 2a,b and 3), and hence we suggest that both groupings should be understood as permeable barriers.

| Signals of divergent adaptive selection
In our SNP data set we identified 140 outlier SNPs as candidates for positive selection. Based on this data set, we recovered a high with a similar pattern are the Atlantic cod Gadus morhua (Barth et al., 2017) and the Atlantic herring Clupea harengus (Limborg et al., 2012).
These results are particularly relevant for the Southern Ocean, as they challenge the classic consideration of Antarctic organisms as stable and homogeneous along their distributions.
The structure and snapclust results from the 140 SNPs under positive selection also showed the effects of the admixture and high migration discussed above in the neutral data set (Figure 5a,b).
However, we observed two unique genetic clusters at the central and southern WAP (Paradise Bay and Adelaide Island), one of them appearing at both stations (purple individuals in Figure 5b) and the other one exclusively present at Adelaide Island (blue individuals in Figure 5b). In agreement, the DAPC analysis clustered together all the stations from the Bransfield Strait area, separating Adelaide Island (as the most divergent sampling station) and Paradise Bay ( Figure 5c). The significant 6.72% of the variance explained by this DAPC grouping (Table 4b) (Asai & Koonce, 2001;Carter, 2013;McKenzie et al., 2015).
As with most other sponges, D. antarctica is a filter-feeding sponge that relies on flagellar beating to modulate the inflow current for particle feeding, and therefore we suggest that the selection signatures in the previously mentioned genes might be related to divergent filtering abilities between the Bransfield Strait area and the central and  (Ciobanasu, Faivre, & Le Clainche, 2013), which is essential for immune responses (Wickramarachchi, Theofilopoulos, & Kono, 2010). In addition, heavy chain dyneins such as DNAH3 have also been reported to aid in the formation of stress granules (SGs) (Kwon, Zhang, & Matthias, 2007;Loschi, Leishman, Berardone, & Boccaccio, 2009). SGs are cytosolic aggregations consisting of RNAs and RNA-binding proteins which appear in response to different stressors, with important function in preserving mRNA and regulating its translation during stress responses . In addition, SGs also prevent apoptosis (e.g., Buchan & Parker, 2009), contain antioxidant machinery (Takahashi et al., 2013), and are involved in cellular recovery after stress exposure . Gleason and Burton (2016) (Figure 4). Apoptosis is a conserved mechanism that occurs during antibacterial responses in sponges (Wiens et al., 2006). Da Fonseca, Kosiol, Kosiol, Vinař, Siepel, and Nielsen (2010) reviewed previous studies on selective pressure in apoptosis-related genes, concluding that positive selection in apoptotic genes is caused by their immune function. Indeed, some of the other genes under positive selection in our data set were related to ubiquitination ( Figure 4), a function also related to the immune system, as it regulates the pattern-recognition receptor signalling that mediates immune responses (Hu & Sun, 2016). Moreover, ubiquitination has been related to local adaptation in corals, where it responds to different environmental factors that cause stress (Bay & Palumbi, 2014;Jin et al., 2016;van Oppen et al., 2018). This is due to its role in removing macromolecular debris such as reactive oxygen species (ROS) generated by cellular stress (Kültz, 2003). Interestingly, The signatures of selection in stress and immune responses that we detected are mostly related to the molecular toolkit that sponges, which are generally filter-feeders, use to discriminate between, and react to, food, pathogens and symbionts in the seawater that they filter and that runs through their bodies (Pita, Hoeppner, Ribes, & Hentschel, 2018). Hence, different microbiome complements in the seawater in different areas would elicit divergent adaptive strategies in sponges in the particular genes that we detected here as under positive selection. Interestingly, differences in sea-ice duration in the Antarctic Peninsula's shallow waters usually translate into highly divergent seawater microbiota, both in composition and in abundance (Ducklow et al., 2013;Vernet et al., 2008). While total sea-ice duration in the vicinity of Adelaide Island is around 250 days a year, with a summer sea-ice retreat, total sea-ice duration is below 150 days in the other sampling stations, generally with spring retreats . This difference in sea-ice duration is key to maintaining vastly different planktonic communities between the southern WAP and the Bransfield Strait area, because the presence and magnitude of phytoplankton blooms in the Southern Ocean are regulated by the timing of sea-ice retreat (Ducklow et al., 2013;Luria, Ducklow, & Amaral-Zettler, 2014;Vernet et al., 2008). Generally, the later the sea-ice retreats, the higher the phytoplankton productivity, as a consequence of sea-ice inhibition of the formation of a spring deep mixed layer, which in turn inhibits phytoplankton (Ducklow et al., 2013 & Estrada, 2001) and by the bacterial dependence on dissolved organic matter, which in turn depends on phytoplankton (Ducklow et al., 2013). Apart from the effects on the planktonic communities, a later sea-ice retreat produces fresher and colder summer surface waters in the southern WAP, due to more recent or ongoing seasonal ice melting (Ducklow et al., 2013).
Thus, due to the ice-plankton interaction outlined above, phytoplankton and bacterial communities, as well as summer surface water temperature, differ widely between the southern WAP, and the Bransfield Strait area. Because sponges are able to feed on both diatoms and bacteria, the different composition of these communities across our study area could potentially drive local adaptation of D. antarctica populations, not only because of their relevant role in food availability, but also as potential agents of diseases and other stresses. Further studies will be directed to test whether this local adaptation hypothesis we suggest for D. antarctica is a general pattern also present in other benthic filter-feeding invertebrates sampled in the same studied area, or even whether microplankton composition generally drives the adaptation of sponges.  (Baldwin & Smith, 2003) and chemicals from local geothermal activity (Deheyn, Gendreau, Baldwin, & Latz, 2005;Elderfield, 1972). Moreover, the fumarolic emissions and geothermal springs spotting the sedimentary seafloor confer upon Port Foster unusually high bottom-water temperatures of 2-3°C (Ortiz et al., 1992). These features undoubtedly affect and stress benthic filter feeders such as D. antarctica, and may have contributed to the upregulation of genes related to different stresses.
Proper differential gene expression analyses will be conducted to test whether particular physicochemical water features at Deception Island are determinant in shaping gene expression in a wide array of shallow-water invertebrates, thus testing their adaptation potential at the transcriptome level.

| CON CLUS IONS
Overall, the current gene flow scenario for Dendrilla antarctica is characterized by high migration and low population differentiation, with a subtle population substructure driven by the oceanic features of the region. Remarkably, despite this background of population admixture, we identified divergent selective pressures along the studied region that could be explained by the sea-icebenthos coupling via planktonic communities. Local adaptation was long assumed to be erased when high population connectivity was present in marine organisms. However, recent investigations indicate that even though few larvae might suffice to maintain genetic homogeneity between populations, that is hardly possible for loci under selection (Sanford & Kelly, 2011). The implications of our results are therefore vast. Our relatively slight patterns of local adaptation are indicative of the potential for plastic physiological responses to environmental shifts. In addition, and in contrast to previous studies of shallow-water sponges, we report a well-connected network of populations across ~1,000 km. Our study therefore corroborates that populations that appear homogeneous for neutral loci may also exhibit local adaptation. In this sense, our study suggests a finely tuned physiological response to current conditions but high resilience to future changes for D. antarctica in the Antarctic Peninsula. However, due to larval reliance on oceanic currents to maintain high dispersal abilities, this exceptional gene flow might be threatened by changes that increasing sea temperature could create in Southern Ocean oceanographic circulation patterns, which are not completely understood yet (Meijers, 2014). Moreover, a general reduction of planktonic larval duration is expected for all larvae in the near future, because their metabolic, developmental and growth rates are determined by water temperature (O'Connor et al., 2007). Thus, a shorter larval stage could imply a reduction of the dispersal capabilities of D. antarctica, with implications for its gene flow and resilience, due to a putatively higher proportion of larvae dying before reaching a suitable settlement site, as has been proposed for fish larvae (Kendall, Poti, Poti, Wynne, Kinlan, & Bauer, 2013;O'Connor et al., 2007). Therefore, our results can be used as a baseline for future assessments of the effects of a changing Southern Ocean on the population connectivity and resilience of D. antarctica.