Phylogeography of the brittle star Ophiura sarsii Lütken, 1855 (Echinodermata: Ophiuroidea) from the Barents Sea and East Atlantic

Ophiura sarsii is a common brittle star species across Arctic and subarctic regions of Atlantic and Pacific oceans. In the Barents Sea O. sarsii is among the dominant echinoderms. We studied genetic diversity of O. sarsii by sequencing the 548 bp fragment of mitochondrial COI gene. O.sarsii demonstrated high genetic diversity in the Barents Sea. Both major Atlantic mtDNA lineages were present in the Barents Sea and were evenly distributed between the northern waters around Svalbard archipelago and the southern part near Murmansk coast of Kola Peninsula. Both regions, as well as other parts of the O.sarsii range, were characterized by high haplotype diversity with a significant number of private haplotypes, being mostly satellites to the two dominant haplotypes, each belonging to a different mtDNA clade. Demographic analyses indicated that the demographic and spatial expansion of Ophiura sarsii in the Barents Sea most plausibly has started during the Bølling–Allerød interstadial, during the deglaciation of the western margin of the Barents Sea.


Introduction
Marine communities of the Northern Atlantic experienced remarkable changes during the Quaternary sea level oscillations, when ranges of species were contracting during the advances of glaciations and expanding during warmer periods (Wares & Cunningham, 2001 etc). Also, several taxa migrated from the Pacific by transarctic migrations (Laakonen et al., 2020 etc.). While historical ecology of nearhore and intertidal ecosystems are being extensively studied (Colson, & Huges, 2006;Rock et al., 2007;Maggs et al., 2008;Krebes et al., 2011;Genelt-Yanovskiy et al., 2019), our understanding of phylogeography and population structure in The Ophiuroidea are important components of marine benthic communities in the Northern Hemisphere, often dominating vast areas of marine shelves (Piepenburg & Schmid, 1996;Blachard et al., 2013;Ravelo et al., 2017). The brittle star Ophiura sarsii Lütken, 1855 is a circumpolar cold-water species widely distributed throughout the North Pacific, North Atlantic, and Arctic oceans where it can reach depths of ca. 3000 m (Djakonov 1954;Anisimova 2000). In the Barents Sea, O.sarsii represent up to 4% of total echinoderm biomass, typically occurring at muddy and muddy fine sand communities at a depth range between 3 and 500 meters (Anisimova 2000). O.sarsii is a predator-scavenger species, and an important food item for many fish species (Packer et al., 1994). It is also a typical marine broadcast spawner species, spawning in the Barents Sea in May-June (Kuznecov, 1963) and then having a feeding larva (ophiopluteus) with planktonic existence of up to 30 days (Dautov & Selina 2009).
Genetic structure of the O.sarsii were investigated in detail only in the Pacific region (Li et al., 2020). But data from the Atlantic are solitary, when this species were used in barcoding studies (Khodami et al., 2014;Laakmann et al., 2016;Layton et al., 2016). Here we presented the pilot study of the mtDNA variation in O. sarsii populations from the Barents Sea.

Samples
Sampling was conducted during the 17-th joint Barents Sea Ecosystem Survey (BESS) cruise of Institute of Marine Research (IMR, Norway) and Polar branch of the VNIRO (PINRO, Russia) (Michalsen et al. 2013) in September 2019 and 114 cruise RV "Vilnus" (Polar branch of the VNIRO) (Fig. 1)

DNA extraction, PCR, and sequencing
Total genomic DNA was extracted from the arm tissues using the QIAamp Fast DNA Tissue Kit (QiaGen). Newly designed forward primer OphS_COI_full_F (TTC TAC CAA CCA TAA GGA TAT AGG G) with HCO-2198 (TAA ACT TCA GGG TGA CCA AAA AAT CA) (0. Folmer et al. 1994) were used to amplify a 700 bp fragment of cytochrome oxidase (COI) gene. PCR amplifications were carried out with a BioRad T-100 cycler. Each PCR reaction mixture (volume 25 μl) contained 3 μl template DNA, 5 μl RCR master-mix ScreenMix-HS (Evrogen), 13 μl nuclease-free water, and 2 μl each primer (10 pmol). The amplification reaction was performed with an initial denaturation step of 2 min at 94°C, followed by 40 cycles of 94°C denaturation for 30 s, 52°C annealing for 60 s, 72°C extension at 60 s, and a final 7 min extension at 72°C during the last cycle. Amplified fragments of each individual were sequenced in both directions at Evrogen (Russia).
Following the pairwise alignments of 5'-3' and 3-'5' individual sequence reads, all polymorphisms were double-checked on the raw chromatograms. Sequences of this study are available from GenBank (accession numbers MW376210-MW376262).

Molecular analysis
Sequences were aligned by Muscle 3.8.425 (Edgar 2004) as implemented in Geneious Prime® 2020.1.2. The mtDNA diversity was firstly illustrated with Bayesian Inference phylogenetic tree reconstruction, for which the GTR+I+G model was selected by JModelTest 2.1.10 (Darriba et al., 2012). Bayesian tree reconstruction was performed in MrBayes 3.2.6 (Ronquist et al., 2012). Analysis started with random trees and performed with four independent Markov Chain Monte Carlo (MCMC) chains for 100000 generations with sampling every 100th generation, the standard deviations of split frequencies were below 0.01; potential scale reduction factors were equal to 1.0; stationarity was examined in Tracer v1.7 (Rambaut et al., 2018). A consensus tree was constructed based on the trees sampled after the 25% burn-in. Mean pairwise distances within and between the two species were calculated using MEGA X (Kumar et al., 2018). Sequences of Ophiura kinbergi (NCBI accession No MN183291) and Ophiura robusta (NCBI accession No MG935300).
Descriptive genetic statistics and haplotype analysis were used with R-package pegas (Paradis, 2010). Haplotype network was built using an infinite site model (uncorrected or Hamming distance) of DNA sequences and pairwise deletion of missing data as implemented in pegas R package (Paradis, 2010). The neutrality tests, Tajima' D (Tajima, 1989) and Fu' Fs (Fu, 1996), were conducted for each of the major haplotype groups present in the Barents Sea.
To test departures from a constant population size, three approaches were used. Firstly, the Tajima's D and Fu's Fs neutrality tests were calculated in DNAsp v6 (Rozas et al., 2017) for each of the major East Atlantic mtDNA lineages occuring in the Barents Sea. The demographic history was further assessed by constructing pairwise mismatch distributions in DNAsp. For the major mitochondrial lineages, the demographic parameter tau (τ) was calculated using DNAsp software to infer the time since regional population expansion. The time since the start of population expansion was calculated as t = τ/2μ, where t is a number of years since population expansion, τ represents a unit of mutational time and μ is the per locus per year mutation rate (Rogers & Harpending, 1992). Finally, the coalescent-based Bayesian skyline plots (BSP) using BEAUti v 2 and BEAST v 2.6.3 (Bouckaert et al., 2014). For BSP priors included the implementation of the GTR substitution model defined by jModelTest, a strict clock model, and the constant skyline model. For both mismatch and BSP analyses, a mutation rate of 2.48*10 -8 per lineage per year (2.48% per million years) was used for the COI, following other studies in Ophiuroidea (Naughton et al., 2014;Sands et al., 2015).

Results
The total aligned COI data set consisted of 548 nucleotide positions for 76 Ophiura sarsii sequences, including 53 newly obtained sequences from the Barents Sea and earlier published sequences deposited in GenBank (Khodami et al., 2014;Laakmann et al., 2016;Layton et al., 2016).
Bayesian tree reconstruction demonstrated that O.sarsii from Hudson Bay grouped with Pacific sequences from British Columbia (Fig.2  The haplotype network revealed the presence of three different haplotype clusters separated by more than five substitutions (Fig. 2). First cluster was mostly represented by a single star-shaped group of haplotypes from the Barents Sea (SVD, BAR), Iceland (ICE) and the North Sea, and was two mutation steps apart from the haplotypes from the Canadian Arctic (BAF, LAW). Second cluster also included haplotypes from the North Sea (NS), Barents Sea (SVD, BAR) and Iceland (ICE) and was more complex with several dominant haplotypes with adjacent satellite singletons connected by 1 substitution. Third cluster of haplotypes was fully attributed to the Hudson's Bay. First and second haplogroups were separated by 11 substitutions, while second and third clusters were separated by 33 mutation steps.
To test the presence of population structure and identify the homogenous groups for the subsequent demographic analyses, a series of AMOVA analyses was implemented. Firstly, AMOVA was implemented on a complete dataset, where two Barents Sea samples (BAR and SVD) were grouped together; four West Atlantic sequences (BAF and LAW) as well as Pacific and Hudson Bay (BC and HUD) were combined in corresponding samples; Iceland and the North Sea were left as individual samples. AMOVA showed that 57% of total variance was distributed among groups, 43% variance within populations with only 0.28% among populations within groups (Table 3). A separate AMOVA test (Table 3) did not show internal heterogeneity between brittle stars from SVD and BAR samples. As no clear population structure was revealed among East Atlantic samples, demographic analyses were implemented for each of the two major O.sarsii haplotype clades. Samples from the West Atlantic, Pacific and Hudson Bay were excluded from the subsequent analyses. Thus, clade 1 was represented by a star-shaped pattern formed by dominant Barents Sea haplotype, and satellite haplotypes from Barents Sea and Iceland; clade 2 was represented by two sub-dominant haplotypes shared between Barents Sea, North Sea and Iceland (see Fig.2 and right indent on Fig.4 for details).
The neutrality Tajima' D index was not significantly different from zero for both clade 1 (D = -1.596, p=0.56), and clade 2 (D = -0.950, p=0.59). Fu's Fs was significantly negative for both clade 1 (Fs = -6.26, p = 0.002) and clade 2 (Fs = -6.58, p = 0.001), indicating population expansion. The mismatch distribution for the clade 1 was unimodal and bimodal for the clade 2 showing peaks at 1 and 4 differences (Fig. 3). Star-shape of clade 1 with the unimodal mismatch distribution and the significantly negative values of Tajima' D and Fu's Fs tests indicated a recent expansion of this clade. Conversion of demographic parameter τ into years was performed using the 2,48% per MY divergence rate (Naughton et al., 2014) and sequence length of 547 bp (μ = 2.8*10 -8 *547 = 0.0000136). Considering τ = 1.039, the expansion time for clade (1) was calculated as 38295 years BP. Similarly, for the clade 2 demographic parameter τ was estimated as 2.350, and time since last population expansion was calculated as 86616 years BP. The clade 2 itself is complex and consists of several abundant haplotypes forming star-like patterns with one dominant haplotype and several satellites differing by one or two nucleotide substitutions. We thus additionally calculated time since past population expansion for the biggest star-like pattern within clade 2. Mismatch analysis (not shown) yielded a demographic parameter τ equal to 0.858, and time since last population expansion was calculated as 31624 years BP. Bayesian Skyline Plots further confirmed demographic expansion for both clades of O.sarsii, yielding in younger median estimates of time since last population expansion as ca. 18000 years BP for clade 1 and ca. 58000 years for clade 2.

Discussion
Our study revealed two major mitochondrial lineages in Atlantic Ophira sarsii, while no clear population structure was found across East Atlantic. O.sarsii sequences from Hudson Bay appeared to be closer to the Pacific than the Atlantic samples. Our findings corroborate the results of recent study of Pacific populations of O.sarsii indicating multiple mtDNA lineages across the species range with one from Yellow Sea identified as separate subspecies, O.sarsii vadicola previously described on the basis of morphological features (Li et al., 2020). Atlantic-Pacific lineage divergence is typical for most amphi-boreal and subarctic marine species (Laakkonen et al., 2020 and references therein). Bayesian tree reconstruction and haplotype network also indicated that West Atlantic O.sarsii sequences formed a compact cluster being sister to the East Atlantic sequences within the clade 1, yet only four sequences were available from the region. Several intertidal and subtidal amphi-Atlantic species demonstrates pronounced divergence between West and East Atlantic haplotypes (Wares & Cunningman 2001;Magnúsdóttir et al., 2019), but cases of shared haplotypes and lineages between two parts of the Atlantic are also common in e.g. sea star Asterias rubens (Wares & Cunningman 2001), quahog Arctica islandica (Dahlgren et al., 2000) and Modiolus modiolus (Halanych et al., 2013).  Naughton et al., 2014;Sands et al., 2015), and representing average divergence between sister species pairs of family Ophiocomidae occuring on two coasts of Isthmus of Panama (Naughton et al., 2014). While the method of calibrating the mutation rates by estimating the divergence rate between pairs of closelyrelated species on both shores of Isthmus of Panama, that emerged ca. 3.5 MY BP, is widely used for many taxa (Kolman & Bermingham, 1997;Donaldson & Wilson, 1999;Lessios, 2008), in this case it is related to species from different families within the Ophiuroidea. Recent review of demographic studies of marine species (Grant, 2015) indicates that choosing a mutation rate to place a historical reconstruction of a population remains the main source of errors in mismatch and BSP analyses.
During the Last Glacial Maximum, the Barents Sea was partly or nearly completely covered by ice (Ingólfsson & Landvik 2013). The analysis of benthic sediments and cores indicated that signatures of glaciation can be observed down to a depth of about 200 meters, and only central parts of the Barents Sea with depths more than 300 meters were ice-free and connected with the Norwegian Sea (Biryukov et al., 1988). The ice-sheet retreat most likely started during the Bølling-Allerød interstadial, i.e. ca. 15 kyr BP at western margin of the Barents Sea (Polyak et al., 1995), and this period thus most plausibly can indicate the beginning of eastward spread of subtidal marine species.
Until recently, the study of mtDNA diversity of O.sarsii, a widely distributed brittle star species in the North Atlantic, was limited to general barcoding projects being aimed on regional revisions of echinoderm faunas of Canada (Layton et al., 2016), Iceland (Khodami et al., 2014) and the North Sea (Laakmann et al., 2016). Our study represents an important contribution in the context of phylogeography of brittle stars in the region in line with recent studies of Atlantic-Mediterranean species (Pérez-Portela, Almada, Turon, 2013;Weber et al., 2014;Taboada & Pérez-Portela, 2016). As a result, we present the data on the northernmost European populations of O.sarsii from the Barents Sea, and since sampling from other parts of the Atlantic remains predominantly qualitative and not quantitative, further studies are needed to implement more broad phylogeographic analysis on this brittle star species.