Contrasting molecular diversity and demography patterns in two intertidal amphipod crustaceans reflect Atlantification of High Arctic

The distribution of two common intertidal amphipod species Gammarus oceanicus and Gammarus setosus was studied along the coast of Svalbard Archipelago. Genetic analysis showed geographical homogeneity of G. oceanicus with only one molecular operational taxonomic unit (MOTU) and much higher diversification of G. setosus (5 MOTUs) in the studied area. Only two MOTUs of G. setosus are widespread along the whole studied Svalbard coastline, whereas the remaining three MOTUs are present mainly along the northern and eastern parts of archipelago’s largest island, Spitsbergen. Distribution analysis indicates that the demographic and spatial expansion of G. oceanicus in the northern Atlantic has started already during the Last Glacial Maximum (LGM, ca. 18 ka), while G. setosus seems to be a long-persistent inhabitant of the Arctic, possibly even through the LGM, with slower distribution dynamics. Combining the results of our molecular study with previous field observations and the knowledge upon the direction of ocean currents around the Svalbard Archipelago, it can be assumed that G. oceanicus is a typical boreal Atlantic species that is still continuing its postglacial expansion northwards. In recent decades it colonized High Arctic due to the climate warming and has partly displaced G. setosus, that used to be the only common gammarid of the Svalbard intertidal zone.


Introduction
The Arctic has been warming at a rate about three times faster than the global mean (IPCC 2014). This results in the melting of glaciers at large spatial scale as well as in the reduction of the duration, area and thickness of sea ice, including both fast ice on the coast and the drifting pack ice (ACIA 2005;Stroeve et al. 2012). Changes in the average air temperature are accompanied by the rise of ocean temperatures and increased northward transport of relatively warm Atlantic waters with the West Spitsbergen Current. These warm waters penetrate more northern parts of the European Arctic (Walczowski et al. 2012), intensifying the recently observed changes in the environment. All of the above described processes are sometimes called the "Atlantification" of the Arctic (Polyakov et al. 2017).
The increasing intrusion of Atlantic waters into the Arctic, consequently changes the distribution ranges of both the Atlantic and the Arctic species (Greene et al. 2008). The northward distribution shifts have been observed at all trophic levels, including bloom-forming coccolithophores (Neukermans et al. 2018), zooplankton (Beaugrand et al. 2002;Weydmann et al. 2014), and benthic mollusks (Berge et al. 2005). The success of arriving species depends on whether they are able to survive in the new environment. Nonetheless, the so far observed changes in the distribution shifts have moved up to several degrees north from their former occurrence ranges (Beaugrand et al. 2002;Neukermans et al. 2018). In consequence, this may lead to the reorganization of both marine and terrestrial ecosystems of the Arctic (Stempniewicz et al. 2007;Węsławski et al. 2009;Vihtakari et al. 2018).
Arrival of new species, reappearance of the former inhabitants from the last climate optimum (7000-4000 years ago) and range expansions of boreal species from the Norwegian Sea shelf have been observed also in Svalbard waters (Blacker 1957;Hopkins 2002;Berge et al. 2005;Fleischer et al. 2007;Beuchel and Gulliksen 2008;Węsławski et al. 2018). The warming of the Arctic and disappearance of land-fast ice, together with intensified erosion of the shores and coastal water turbidity, is responsible for the structural changes in the rocky-bottom and intertidal faunal communities. The most obvious are extreme fluctuations of salinity and temperature, as well as the devastating disturbance caused by scouring of land-fast ice.
Marine and brackishwater amphipods, particularly gammarids, due to their usually high abundance and biomass, combined with predominantly detritivorous feeding, are among keystone organisms maintaining the high rate of energy transfer in the intertidal ecosystems of the temperate and boreal climatic zones (Gerhardt et al. 2011). However, in contrast to their ecological importance, the genetic diversity, population structure and demography of these organisms have been still insufficiently studied. In comparison, high cryptic diversity, well pronounced spatial population genetic structuring and local endemism have been observed in several freshwater gammarids, particularly those with large distribution ranges (e.g. Mamos et al. 2016;Copilaş-Ciocianu and Petrusek 2017;Grabowski et al. 2017). Even the few published studies upon the marine taxa revealed that divergent cryptic lineages or even haplogroups may differ by geographical distributions, various habitats requirements, altitudinal and depth ranges, demographic histories or tolerance to environmental factors (Krebes et al. 2011;Hupało et al. 2018). In addition, a more general pattern of higher degree of genetic differentiation in non-invasive versus invasive gammarids was recently suggested (Baltazar-Soares et al. 2017). This underlines the importance of historical factors for the present composition and dynamics of benthic communities. Thus, the general scarcity of such studies upon the marine coastal gammarids is surprising, particularly in areas undergoing fast environmental changes, such as the High Arctic.
Both occur at localities covered with loose stones and rocks, reaching highest densities up to 2000 ind/m 2 , just below the mean low water mark (Węsławski 1994;Ikko and Lyubina 2010). They prefer moderately to low exposed sites, where large loose stones serve as a shelter against birds during low tide. According to Węsławski (1994), each of them seems to be very tolerant to changes in salinity (0 to 35 PSU) and temperature (− 1.5 to 8 °C on Spitsbergen). Both species are typical inhabitants of coastal marine and brackish waters in the North Atlantic with very similar reproductive potentials and wide, largely sympatric, distribution areas. Gammarus oceanicus ranges from Long Island on the eastern coast of North America, through the Arctic, to the southern North Sea and the Baltic Seas, while G. setosus ranges from the Bay of Fundy on the Atlantic coast of North America, to the Arctic Steele 1970, 1972;Lincoln 1979;Węsławski 1983Węsławski , 1991Gaston and Spicer 2001). Recent molecular studies have suggested presence of geographic structuring in both species, with diversified populations and even divergent lineages along the boreal North American coast versus the Arctic and boreal European coastlines (Costa et al. 2009;Krebes et al. 2011).
Although the habitat preferences of both littoral Gammarus species are very similar elsewhere (Tzvetkova 1975;Steele and Steele 1975;Węsławski 1994), the Arctic G. setosus is expected to perform better in lower temperature, while boreal G. oceanicus in warmer conditions (Steele and Steele 1970, 1972, 1974Whiteley et al. 2011). Large scale comparison of Gammarus species occurrence in Svalbard was done between 1980(Węsławski et al. 2018. It shows that G. oceanicus has apparently expanded its distribution range along the west and north coasts of Spitsbergen and has increased its dominance on the expense of G. setosus, which has generally remained where it was observed 20 years ago.
In the present paper, we aim to reveal whether the observed extension of distribution range and increasing domination of G. oceanicus in Svalbard versus the generally decreasing abundance of G. setosus is reflected in genetic diversity of the local populations of both species. We hypothesised that G. oceanicus, that is presumably a recent newcomer to the High Arctic, would be characterised by rather uniform population genetic structure, without a clear spatial pattern, yet with strong signal of recent demographic and spatial expansion. In contrast, the populations of G. setosus, being a more cold-water and presumably typical Arctic species, would show higher diversification at a molecular level and more pronounced spatial patterns of genetic diversity, but much less dynamic demography.

Study area
The Archipelago of Svalbard is situated in the European sector of the Arctic Ocean. In terms of water masses, this area is a transition zone between the Atlantic and Arctic domains. The west coast of its largest island, Spitsbergen, is under the influence of the West Spitsbergen Current (WSC), which transports relatively warm and saline Atlantic water northwards to the Arctic Ocean. The south-eastern part of the archipelago is influenced by the coastal South Cape Current (SCC) that carries cold, less saline Arctic-type water (Piechura et al. 2001;Cottier et al. 2005). Presence, amounts and salinity of water masses in Svalbard fjords result from the local processes and balance between the flows of main Atlantic and Arctic water masses. Therefore, the fjords and their fauna, are potentially sensitive indicators of environmental change (Svendsen et al. 2002;Cottier et al. 2005;Nilsen et al. 2008).
The Svalbard coast consists mainly of low gravel and stony beaches, with several tidal glacier cliffs, moraine lagoons, and frequent rocky peninsulas with half-submerged skerries (Węsławski et al. 1993). The intertidal zone of Svalbard was described by Węsławski et al. (1993), and its recent changes by Węsławski et al. (2010). The source material for the Gammarus occurrence observations in Svalbard is presented on the Svalbard Intertidal Project website http://www.iopan .gda. pl/proje cts/SIP/.

Sampling and morphological identification
The material was gathered in the intertidal zone around Svalbard Archipelago, in late July and early August of 2015 and 2016, when neither fast ice nor sea ice was present on the shore (Fig. 1, Table 1). The amphipods were collected at 15 sampling sites, under the stones, each time from the surface of c.a. 1 m 2 . The collection took place during low tide in lower intertidal, with the use of hand net to catch individuals under the loose, flat stones. The samples were sorted at a site and immediately fixed in 96% ethanol. In total 389 specimens of gammarids were collected. They were observed and identified to species level under the stereomicroscope NIKON SMZ800 with the use of available keys (Lincoln 1979;Stephensen 1940). The material has been kept in the permanent collection of the Institute of Oceanology, Polish Academy of Sciences.

Molecular analysis
Three pleopods were cut out with a sharp-edged tweezers from each individual. The total DNA extraction was performed with phenol-chloroform (Hillis et al. 1996) and/ or Chelex (Casquet et al. 2012) procedures. In case of phenol-chloroform method the muscle tissue was incubated overnight at 55 °C in a 1.5 ml tube containing 200 μl of Queen's lysis buffer (Seutin et al. 1991) and 5 μl of proteinase K (20 mg/ml). The collection of extracts is kept in the Department of Invertebrate Zoology and Hydrobiology, University of Lodz and in the Department of Ichthyology and Hydrobiology, St.-Petersburg University, Russia. The mitochondrial marker COI (Cytochrome c oxidase subunit I), ca. 670 bp in length, was chosen for PCR amplification according to the Hou et al. (2007) procedure with HCO2198/LCO1490 (Folmer et al. 1994) and HCO JJ/LCO JJ (Astrin and Stuben 2008) primer pairs. The PCR was performed in a total volume of 10 µl containing 5 µl of DreamTaq Master Mix (2×) Polymerase (Fermentas), 0.8 µl of each primer (concentration 5 lM), 2.4 µl ultrapure H 2 O and 1 µl of DNA template. A 2 µl aliquot of the PCR products was visualized in Midori Green-stained (Nippon Genetics) 1.0% agarose gels to check PCR product quality and length. The PCR products were purified with Exo-FastAP (Thermo Scientific) (Werle et al. 1994) and sequenced with BigDye terminator technology by Macrogen Inc., Korea.
We also used sequence data deposited already in BOLD Systems by other investigators for comparison with our results.
The haplotypes were identified with DnaSP 6 (Librado and Rozas 2009) while the relationships between them were illustrated with Minimum Spanning Network (MSN) in Pop-ART 1.7 (Bandelt et al. 1999).
In order to examine the level of cryptic diversity within the investigated species, the Molecular Operational Taxonomic Units (MOTU) were identified with the distancebased barcode-gap approach using both the Automatic Barcode Gap Discovery (ABGD) method (Puillandre et al. 2012) and the Barcode Index Numbers (BINs) ascribed to the sequences by BOLD (Ratnasingham and Hebert 2013).
To illustrate haplotype and/or MOTU distribution around Svalbard, a series of non-metric multidimensional scaling plots (nMDS) was performed in PRIMER 7.0 (Clarke and Warwick 1994). The analysis was conducted on the standardized proportion of haplotypes for G. oceanicus, while in the case of G. setosus separate plots were performed for all MOTUs and, separately, for the most diverse and haplotyperich MOTU A and MOTU B. To illustrate which haplotype and/or MOTU influenced the distribution of samples in a given plot, overlay vectors were shown based on multiple correlations of haplotype or MOTU proportions among sampling sites.

Phylogeny
The time calibrated within-species phylogeny of G. setosus was reconstructed, using all identified haplotypes, with Bayesian methodology in BEAST v. 1.8.2 software package . A haplotype of G. oceanicus from the present study was chosen as an outgroup. Three, 10 M iterations long runs of Monte Carlo Markov Chains (MCMC), sampled each 1000 iterations, were performed. The strict molecular clock was calibrated with the standard mitochondrial rate for arthropod COI equal to 0.0115 substitution per million years per base pair (Brower 1994). Such rate was accepted as a reasonable midpoint approximation for Gammarus taking into account the recently published cross-validations for other Gammarus species, based on the independent calibration points (Grabowski et al. 2017). The runs were examined in Tracer v1.6 (Rambaut et al. 2014), so they reached the values of effective sampling size (ESS) of at least 200 for each parameter. The Yule speciation process

Demography
The demography of both, G. oceanicus and G. setosus, populations in Spitsbergen was analyzed with mismatch distribution, under a sudden demographic and spatial expansion models, with 1000 bootstrap replicates, in Arlequin 3.5.1.3 software with Rcmd plug-in Lischer 2010, 2011). We verified the validity of both models with the sum of squared deviations (SSD) between observed and expected mismatches and by the Harpending's Raggedness statistics (HRI). We also tested for recent demographic expansion in both species with two mutation neutrality tests, Tajima's D (Tajima 1989) and Fu's Fs (1997), also implemented in Arlequin, with 1000 simulation as evaluation for test significance. We used the same software to test for the genetic dissimilarities between populations of both species, sampled at different localities around Spitsbergen with the F st estimator, as well as by the exact test for population differentiation based on haplotype frequencies (Raymond and Rousset 1995), using default program set-ups. The fluctuations of demography over time was illustrated with the extended Bayesian Skyline Plot (eBSP), reconstructed in BEAST v1.8.2 and based on three runs of MCMC 20 M iterations long, and sampled every 1000 iterations. Taking into account that even small differences in accepted mutation rate may produce large shifts in time estimations for very recent demographic events, we performed preliminary analyses, each with different mutation rate based on different calibration schemes for Gammarus: 0.0113, 0.0129 (Grabowski et al. 2017), 0.0165 (Copilaş-Ciocianu and Petrusek 2017) and the standard rate of 0.0115 (Brower 1994). The runs did not differ substantially and eventually the standard rate of 0.0115 was applied for analysis as a reasonable midpoint approximation. The runs were examined in Tracer v1.6 (Rambaut et al. 2014). The Hagesava-Kishino-Yano (HKY) model of evolution chosen in jModel Test 2 (Darriba et al. 2012) was set as prior. The eBSP plot was produced in MS Excel.

Results
In the material collected from a total of 15 sampling localities around Svalbard, we found 110 individuals of G. setosus at 14 sites and 279 individuals of G. oceanicus at 12 sites. Both species co-occurred at 11 sampling sites (Table 1). The morphological identification of all the collected individuals was confirmed with DNA barcodes.
The 389 COI sequences were ascribed to 20 haplotypes within G. setosus and 22 haplotypes within G. oceanicus. The reticulate topology of the Minimum Spanning Network illustrating phylogenetic relationships between haplotypes of G. setosus suggests rather long history of divergence within the species around Svalbard and presence of five distinct lineages (Fig. 2). In contrast, the populations of G. oceanicus appear genetically uniform and the star-like topology of the network implies a recent population expansion in the area (Fig. 2).
The delimitation of MOTUs, both with the ABGD and BIN algorithms, indicated the presence of cryptic diversity with five MOTUs and BINs within G. setosus (MOTU: BIN): A: ABY7004, B: ABZ4044, C: ACE3104, D: ACE3105, E: ACE4619. In contrast, all the individuals of G. oceanicus from Svalbard were grouped within the same MOTU/BIN (AAA1262) (Fig. 2). The topology of the COIbased Bayesian chronogram suggests that the G. setosus MOTUs started to diverge already in early Pleistocene, at ca. 1.5 Ma (Fig. 3).
Illustration of the distribution of G. setosus haplotypes and MOTUs around Svalbard implied at least some geographical structuring, despite presence and high share of the most common MOTU A in all the sites where the species was found (Figs. 2 and 4). Interestingly, MOTU A was the only one recorded in the southern Spitsbergen. Another common and highly frequent MOTU B was present elsewhere with the exception of North-East Spitsbergen. However, we cannot exclude that it may be an effect of undersampling. The remaining MOTUs were present in a few disjunct localities around Svalbard, with the exception of southern Spitsbergen. A closer look into the haplotype distribution within the most common MOTUs shows high frequency of the abundant haplotype 1 in samples from the West, North-West and South Spitsbergen and of haplotype 4 in the West Spitsbergen, both within MOTU A (Fig. 5). There was also a high frequency of haplotype 3, belonging to MOTU B, at stations located on Prins Karls Forland and West Spitsbergen (Fig. 5).
In G. oceanicus, the illustration of geographic distribution of haplotypes showed some, however little pronounced, spatial genetic population structure around Svalbard (Figs. 2,  4). Only the nMDS analysis (Fig. 5) allowed for distinction of two main sampling site groups differing by the haplotype domination structure: first with sites (from Isfjorden and East Hinlopen Strait) dominated by haplotype 1, and second with other sites (from Prins Karls Forland, South and North-West Spitsbergen) where the haplotype 2 was most frequent. The only exception was the sample from Kongsfjorden, where the haplotypes 2 and 14 co-dominated. The F ST values were generally low (usually much lower than 0.3) (Fig. 6). It is also important to stress that the significant F ST values measured for G. oceanicus overlapped (with one exception) with the significant differences provided by exact test for population differentiation (Table Suppl. 1, Fig. 6). On the other hand, regarding G. setosus, the significant F ST values only in about half of the cases pointed to the same pairs of sampling areas as did the results of the exact test for population differentiation (Table Suppl. 2, Fig. 6).
Mismatch analysis and the associated neutrality tests showed no traces of recent demographic expansion of G. setosus around Svalbard, however some degree of spatial expansion could not be excluded. The results of eBSP show that the population of G. setosus in Svalbard could contract in Late Pleistocene, ca. 30 ka. The population have possibly started to expand at ca. 20 ka, however the very high span of the HPD intervals confounds unequivocal interpretation (Fig. 7).
In case of G. oceanicus, both, the mismatch analysis and the neutrality tests clearly suggest that its population around Svalbard is in the stage of demographic and spatial expansion ( Table 2). The eBSP revealed that the rapid  (Fig. 7).

Discussion
Our study shows contrasting, to some extent, patterns in molecular diversity and demography of the two sibling gammarids, G. oceanicus and G. setosus, dominating intertidal communities in Svalbard. Most recently, Węsławski et al. (2018) observed that, during the last 20-30 years, the former species has expanded its distribution range by over 1300 km along the west and north coasts of Spitsbergen. They reported also a dominating position of G. oceanicus on the number of sites, where it was previously just an occasional species, while G. setosus has remained where it was observed already three decades ago. We observed absence of cryptic diversity within G. oceanicus from Svalbard, star-like topology of haplotype network and lack of spatially structured genetic differentiation, yet with relatively high number of weakly divergent haplotypes (Fig. 2). Such features are frequently observed in recent colonizers with high propagule pressure (e.g. Haydar et al. 2011;Rewicz et al. 2015). Also the set of demographic analyses we performed, points out to a recent demographic and spatial expansion of the species. However, our eBSP analysis suggests, that G. oceanicus has started to expand already in    (Fig. 5) and not only during the recent decades, as observed by Węsławski et al. (2018). Such results were observed for other North Atlantic intertidal macroinvertebrates with typically boreal distribution. They include e.g. the snail Littorina obtusata, crustaceans such as Semibalanus balanoides and Idotea balthica (Wares and Cunningham 2001) as well as the Arctic planktonic copepods Calanus glacialis (Weydmann et al. 2018). Such an early onset of expansion of boreal marine species may be explained by the relatively milder conditions along the northern European shores, if compared to Western Atlantic, that enabled the existence of vast refugial areas extending from as much north as Faroe Islands to Northern Africa (e.g. Ingólfsson 1992;Kettle et al. 2011;Waltari and Hickerson 2013). Thus, it seems that the mass colonisation of Svalbard by G. oceanicus in the last few decades can be considered as another step forward in the ongoing post-glacial expansion of this species. Apparently, this expansion has been greatly intensified by the recent climatic changes, but it has not left any specific signature in the genetic structure of the local populations.
A question arises, where to look for the source of this colonisation and for initial refugia from which the expansion northwards started. The Svalbard population represents the same MOTU as the populations inhabiting North Sea/ European Atlantic coast, Baltic Sea and the White Sea. On the other hand, the North American boreal coast can be  excluded as another, genetically divergent, MOTU inhabits this area. Then, the Canadian Arctic population contains just a small subset of the haplotypes found in European Arctic, suggesting recent westward colonization (BOLD Systems, Krebes et al. 2011, Radulovici 2012. Of the European seas, the Baltic emerged as a freshwater Baltic Ice Lake (BIL) at ca. 16 ka BP, after the retreat of the ice-sheet. Then it passed through brackish and purely freshwater phases, and only since ca. 8 ka regained its connection to the North Sea and present brackish water conditions (Björck 1995;Andrén et al. 2011). Thus it is highly doubtful as the glacial refugium for G. oceanicus, and the population present in this basin, is presumably also of a very recent origin. This leaves either the North Sea/European Atlantic Coast or the White Sea as the potential refugia. However, at this stage we cannot discriminate between the two for the absence of reliable reference data and insufficient sampling in these areas. In most already studied cases of species that expanded post-glacially to the Arctic, the North Sea/European Atlantic Coast is the most plausible refugium (Wares and Cunningham 2001;Maggs et al., 2008). Such direction of colonization would be also favored by the local system of sea currents (Beszczynska-Möller et al. 2011; Fig. 1). Interestingly, fossil data show that in late Pleistocene/early Holocene the ocean quahog (Arctica islandica), showing at present similar distribution as G. oceanicus, occurred in the Bay of Biscayne and Mediterranean as well as along the Siberian coast adjacent to the present White Sea (Dahlgren et al. 2000). Thus, however no fossil data is available in this case, we cannot exclude a similar situation and colonization of Arctic from two refugia in case of G. oceanicus. The climatic conditions in High Arctic during the LGM are still disputable. The reconstructions of the ice extent over Svalbard and the adjacent sea in late Weichselian range from almost no ice (Troitsky et al. 1979, Boulton et al. 1982 to thick ice-cover extending to the western shelf edge of the Barents Sea as well as north of Svalbard (Schytt et al. 1968;Mangerud et al. 1998;Landvik et al. 1998). More recent studies, both geological and phylogeographical, agree that at least limited areas in High Arctic remained unglaciated through LGM, providing microrefugia for Pleistocene terrestrial biota (Funder 1979;Brochmann et al. 2003;Landvik et al. 2003;Klütsch et al. 2017). Concerning marine shelf refugia, it is known that continental shelf areas reduced greatly all over the world, due to the general eustatic sea level regression during LGM (Clark and Mix 2002;Boehme et al. 2012). In the Arctic, vast coastal areas were presumably covered by ice sheets and glaciers (Clark and Mix 2002;Lambeck et al. 2002). Yet, several studies pointed out to presence of Arctic and periglacial refugia for various marine taxa, e.g. along the Norwegian coast, Island and Beringia (for summary and discussion see Hardy et al. 2011). The high cryptic diversity of G. setosus in the Svalbard archipelago, with five MOTUs of Pleistocene age (Fig. 3) distributed all over this relatively limited area (Fig. 4), is a very surprising discovery. However, it must be stressed that the illustrated distribution patterns may, to some extent, be an artefact resulting from the generally low frequency and abundance of G. setosus in samples dominated by G. oceanicus.
It suggest a rather long and complex evolutionary history of the species. It must be also taken into account, that G. setosus is a cold-water stenotherm, nowadays typical of Arctic and sub-Arctic waters (Steele and Steele 1970;Węsławski 1991). Thus, two main scenarios may be drawn to explain such diversification pattern: (1) long persistence and local diversification of G. setosus in the High Arctic, including Svalbard, during Pleistocene; (2) multiple late glacial and post-glacial colonization of High Arctic by the already diverged MOTUs of G. setosus from several, so far unknown, Arctic or periglacial refugia. Concerning the first scenario, so far, there is no evidence for the existence of marine shelf refugia in High Arctic, their presence is also disputable in the adjacent Arctic areas (e.g. Ingólfsson 2009), however the lack of such evidence cannot be treated as an evidence that they did not exist. In situ divergence in High Arctic archipelagos and only geographically limited subsequent recolonizations of formerly glaciated coastlines would be most-parsimonious explanation. Geological studies show that even if, ge, most of the Svalbard coastline and fjords were covered with glaciers, there were also prolonged episodes with much milder and iceless conditions, both on land and in the ocean (Troitsky et al. 1979;Landvik et al. 1998). Thus, presence of patchy distributed ice-free shelf and/or intertidal habitats on Svalbard and other High Arctic archipelagos is quite probable. The first scenario is also supported by our results of the demographic analyses that do not indicate explicit population expansion in Late Pleistocene (Fig. 5), as in the case of other boreal taxa. In such instance, G. setosus in Svalbard would be a Pleistocenic relict. The second scenario should be associated with much more prominent expansion signals, such as in cases of other recently expanded species (Berge et al. 2005;Ware et al. 2013;Berge et al. 2015). However, a few decades ago G. setosus was apparently much more abundant in Svalbard, until the very recent expansion of the boreal G. oceanicus (Węsławski et al. 2018). Thus, it cannot be excluded that the latter species outcompeted G. setosus what caused a significant reduction of its population size and masked the signal of former demographic/spatial expansion. At this moment there is no available data on molecular diversity and identity of populations in possible refugial areas, such as Iceland, Lofoten or Siberian coast, that could help to back the second scenario.
In case of both species, the F ST values between regions were usually lower than 0.3, which is much less than the average level for native gammarids and close to the range of fast spreading invasive gammarids (Baltazar-Soares et al. 2017). Such generally low F ST values suggest high connectivity between different collection sites (Excoffier and Lischer 2011). The discrepancies between the distribution of significant F ST values and significant sample differences measured with the use of the exact test for population differentiation (Fig. 6) could be explained by disproportions in the sample sizes deriving from different regions (Table 1).
Concluding, we have shown for the first time contrasting molecular diversity and demographic patterns in two crustaceans, with different thermal optima, dominant in the Arctic intertidal zone, that currently co-exist in the coastal habitats of the Svalbard Archipelago in High Arctic. The generally boreal G. oceanicus, that has been observed to colonise Svalbard in recent decades, is apparently still in the phase of the post-glacial mass population expansion from extra-Arctic refugia, that is likely to be enhanced by the Atlantification of High Arctic. On the other side, the typically cold-water G. setosus seems to have a much longer evolutionary history in the Arctic, where it apparently survived the LGM. Yet, its demography is much less dynamic than that of the former species. Very low number of haplotypes in some MOTUs may even reflect population decline associated with mass expansion of G. oceanicus. This points out to the necessity for further studies upon the still understudied molecular biogeography patterns in High Arctic biota, that may be heavily influenced by the recent climatic changes.
Acknowledgements Molecular work at the University of Lodz was financed by the NCN Project No. 2014/15/B/NZ8/00266 and by the university statutory funds. The Centre for Molecular and Cell Technologies of St. Petersburg State University is acknowledged for technical help. Sampling has been performed during Svalbard Intertidal Project expeditions, supported by Sysselmannen pa Svalbard Miljovernfond RIS 3423 and Polish Ministry of Science support fund SPUB-SIP 2017, https ://www.iopan .pl/proje cts/SIP/SIP_2015/index .html.

Compliance with ethical standards
Conflict of interest The authors declare that they have no conflict of interest.
Ethical approval All applicable international, national, and institutional guidelines for the use of animals were followed. The studied species neither are threatened by extinction. They are also not subjected to any other national or international legal restrictions. The samples were collected within the official sampling campaigns supported by Sysselmannen pa Svalbard Miljovernfond (Norway) and Ministry of Science and Higher Education (Poland), with no harm to the local environment and in compliance with all the legal regulations and approvals, both from the Norwegian and Polish authorities.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creat iveco mmons .org/licen ses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.