Poles Apart: The “Bipolar” Pteropod Species Limacina helicina Is Genetically Distinct Between the Arctic and Antarctic Oceans

The shelled pteropod (sea butterfly) Limacina helicina is currently recognised as a species complex comprising two sub-species and at least five “forma”. However, at the species level it is considered to be bipolar, occurring in both the Arctic and Antarctic oceans. Due to its aragonite shell and polar distribution L. helicina is particularly vulnerable to ocean acidification. As a key indicator of the acidification process, and a major component of polar ecosystems, L. helicina has become a focus for acidification research. New observations that taxonomic groups may respond quite differently to acidification prompted us to reassess the taxonomic status of this important species. We found a 33.56% (±0.09) difference in cytochrome c oxidase subunit I (COI) gene sequences between L. helicina collected from the Arctic and Antarctic oceans. This degree of separation is sufficient for ordinal level taxonomic separation in other organisms and provides strong evidence for the Arctic and Antarctic populations of L. helicina differing at least at the species level. Recent research has highlighted substantial physiological differences between the poles for another supposedly bipolar pteropod species, Clione limacina. Given the large genetic divergence between Arctic and Antarctic L. helicina populations shown here, similarly large physiological differences may exist between the poles for the L. helicina species group. Therefore, in addition to indicating that L. helicina is in fact not bipolar, our study demonstrates the need for acidification research to take into account the possibility that the L. helicina species group may not respond in the same way to ocean acidification in Arctic and Antarctic ecosystems.


Introduction
Over the past 200 years the world's oceans have absorbed approximately one third of the total carbon dioxide (CO 2 ) released into the atmosphere by human activities [1]. This CO 2 uptake is causing profound changes to seawater chemistry, including a reduction in pH (i.e., ocean acidification) and a reduction in the saturation state of calcium carbonate (CaCO 3 ) [2]. The latter is critical to the formation of CaCO 3 skeletal structures by a wide range of marine organisms, including molluscs, corals, echinoderms and crustaceans, as their calcification rates are directly related to the CaCO 3 saturation state of seawater [3]. Decreasing CaCO 3 saturation levels are of particular concern for organisms that build their skeletons out of aragonite, a metastable form of CaCO 3 that is ,50% more soluble than calcite, and for organisms in the polar regions where CaCO 3 undersaturation, and hence skeletal dissolution, is expected to occur first [4]. Recent projections are that localised aragonite undersaturation of Arctic surface waters may occur within a decade [5], while the surface waters of the Southern Ocean (Antarctic) may begin to become aragonite undersaturated by as early as 2030 [6].
Aragonite-shelled (thecosome) pteropods, pelagic swimming sea snails sometimes referred to as sea butterflies, occur in all oceans but are particularly abundant in the polar regions [7,8]. Here they are principally represented by what is considered to be a bipolar species, Limacina helicina (Phipps 1774) (Figure 1a). Because of its aragonite shell and polar distribution, L. helicina may be one of the first organisms affected by ocean acidification, and it is therefore a key indicator species of this process [3]. L. helicina is a major component of the polar zooplankton. It can comprise .50% of total zooplankton abundance (number of individuals per unit volume) and it plays a significant ecological role as a phytoplankton grazer and prey species for zooplankton and fish, while also contributing substantially to carbonate and organic carbon flux [8]. As one of the organisms most vulnerable to ocean acidification, and a key component of polar ecosystems, L. helicina has become a focal point for research on acidification impacts [3,9].
Currently, northern and southern hemisphere L. helicina are listed as the sub-species L. helicina helicina and L. helicina antarctica respectively. In addition, the taxonomic category ''forma'' has been applied to designate at least three morphotypes of L. h. helicina (acuta, helicina and pacifica) and two morphotypes of L. h. antarctica (antarctica and rangi). These forms typically have different geographical ranges but it remains unclear as to whether ''forma'' represent morphological responses to different environmental conditions or are indeed taxonomically distinct, and if the latter, their level of taxonomic separation [10]. Recent findings show that the response of organisms' calcification rates to acidification can vary markedly between taxonomic groups [11]. It is hypothesised that this varied response is due to physiological differences, occurring even at the species level. In the absence of a detailed understanding of, and ability to measure, the physiological processes controlling calcification rates, correct taxonomic data are critical for quantifying acidification impacts. Antarctic (L. helicina antarctica) respectively, but 33.5660.09% between poles. Support is indicated as posterior probabilities above nodes (* indicate 1.0 support) and bootstraps from a maximum likelihood analysis below (* indicate 100% support). The scale bar represents substitutions per site. GQ861824 and GQ861825 from the Amundsen Sea; GQ861831, GQ861832 and GU732830 from the vicinity of South Georgia; GQ861826/27/28/30 from the Beaufort Sea; AY22739 and AY227378 from [12]. doi:10.1371/journal.pone.0009835.g001 Given the importance of L. helicina as an indicator of ocean acidification there is an urgent need for research that will resolve the taxonomic status of the L. helicina group. Molecular techniques represent an appealing route to take as they avoid the potential confusion resulting from environmentally induced morphological plasticity. The cytochrome c oxidase subunit I (COI) gene has been demonstrated to be well suited to gastropod phylogenetics [12]. Based on a single specimen of each, Remigio and Hebert [12] provided initial evidence for the genetic separation of L. h. helicina and L. h. antarctica. Here, we build upon their study and use COI sequences from multiple specimens of the Arctic L. h. helicina forma helicina and the Antarctic L. h. antarctica to quantify genetic distance within and between these regions with the specific aim to assess the bipolar status of the L. helicina species group.

Results and Discussion
We found a 33.56% (60.09) difference in COI sequences between the Arctic L. h. helicina forma helicina and the Antarctic L. h. antarctica (Figure 1b). This degree of separation is sufficient for ordinal level taxonomic separation in other organisms [13] and convincingly demonstrates that L. helicina is not bipolar, but that the Arctic and Antarctic populations differ at least at the species level. Our results support Remigio and Hebert [12] in identifying L. helicina as a rate-accelerated lineage within pteropods (Figure 1b). A conservative divergence time estimate of 31 Ma (95% HPD interval 12-53 Ma) for Arctic and Antarctic L. helicina, indicates that they have undergone rapid independent evolution since the establishment of cold water provinces in the early Oligocene.
Our results show the need for a revision of the taxonomic status of the L. helicina species group. The high degree of separation at what is considered the sub-species level, suggests that COI sequences analysis may also provide an effective means to clarify the relationships between the ''forma'' of both L. h. helicina and L. h. antarctica. Our study only included one form of L. h. helicina (forma helicina). In the case of the Antarctic sub-species, forma was not determined for any of the specimens analysed. Based on known biogeographic distributions the Amundsen Sea specimens were most likely forma antarctica [10], while analysis of the South Georgia net samples indicated that only forma antarctica were present. Although it is possible that the South Georgia specimens sequenced were forma rangi, the high COI sequence similarity between Antarctic samples demonstrated that specimens were closely aligned. It remains to be determined whether this similarity was form specific, or whether forma are indeed morphotypes and not genetically distinct.
As highlighted in the introduction, due to unique physiologies, the response of organisms to ocean acidification may vary even at the species level. A recent study comparing the locomotor abilities of another supposedly bipolar pteropod species, Clione limacina, identified significant differences in the aerobic capacity of Arctic and Antarctic forms, associated with neuromuscular and mitochondrial composition [14]. Given the substantial genetic divergence between Arctic and Antarctic L. helicina populations observed in our study, similarly large physiological differences may exist between the poles for the L. helicina species group. Therefore, in addition to the taxonomic implications, our study demonstrates the need for acidification research to take into account the possibility that the L. helicina species group may not respond in the same way to ocean acidification in the Arctic and Antarctic. Physiological differences between taxa coupled with differences in the processes and rates of acidification at the poles [4,5,6], brings to light the possibility that differences in acidification impacts in the Arctic and Antarctic may extend beyond species to the ecosystem level.

Materials and Methods
Specimens of the Antarctic Limacina helicina antarctica were obtained from the Amundsen Sea and the vicinity of South Georgia Island (Figure 1b). The forma of these specimens was not determined. Specimens of the Arctic Limacina helicina helicina were identified as forma helicina and were obtained from the Beaufort Sea. Full locations and station data are available on Barcode of Life Data systems (BOLD)/GenBank. Extraction, amplification and sequencing followed standard DNA barcoding protocols [15,16]. DNA was also extracted using the high salt method [17]. PCR amplifications were performed using the standard Folmer [18] primers and sequencing was carried out by Macrogen (Korea). New sequences have been deposited in BOLD/GenBank (Accession numbers GQ861824-861832, GU7328230). The length of L. helicina sequences varied from 528 bp to 618 bp. This variation reflects difficulty in amplifying the fragments. Alternative COI primers, a combination of standard primers [18] and minibarcode primers yielding two overlapping fragments [19], had to be used in addition to recover these shorter sequences. The published sequences of Remigio and Hebert [12] for single specimens of L. h. helicina and L. h. antarctica were included in subsequent calculations of genetic distance.
The K2P model [20] of sequence evolution was used to calculate the genetic distance for L. h. helicina and L. h. antarctica both within and between regions, (i.e., Arctic and Antarctic) using PAUP 4.0b10 [21]. The genetic distance between COI sequences of five individuals collected from the Arctic was 0.1560.06%, whilst the genetic distance between COI sequences of six individuals collected from the Antarctic was 0.6060.07%. Genetic distance between individuals collected from the two regions was 33.5660.09%.
Bayesian analyses were conducted using BEAST v1.4.8 [22], using a SRD06 nucleotide model [23]. Analyses were run with both strict clock and uncorrelated log-normal relaxed clock [24] models, with the mean substitution rate fixed to 1.0. A Yule prior on branching rates was employed [24]. Gymnosomata and Thecosomata were assumed to be reciprocally monophyletic [25]. Two independent MCMC analyses were run for each parameter set. Acceptable mixing and an appropriate 'burnin' was determined using Tracer v1.4.1 [26]. Each analysis was conducted for 20 million generations sampling every 1000 generations. The Bayes factor [27] was used to compare strict and relaxed clock models as implemented in Tracer v1.4.1. The uncorrelated lognormal relaxed clock model was preferred with a Bayes Factor (natural log) of 20.960.2.
Phylogenetic maximum likelihood analyses were performed with RAxML v.7.0.4 [28]. All searches were completed with the GTRMIX option and bootstraps were calculated with 1000 replicates. To obtain a minimum divergence time estimate of L. helicina from Arctic and Antarctic regions we also analysed the data within BEAST v1.4.8 (using a SRD06 nucleotide model and an uncorrelated log-normal relaxed clock model) using a fixed calibration date of 58.7 Ma on the divergence of Limacina (Limacinidae) and Hyalocylis (Creseinae) [29]. L. mercinensis is the oldest known limacinoid fossil from the Thanetian (58.760.2-55.860.2 Ma) [30]. The oldest known Creseinae fossils are from the Ypresian (55.860.2-48.660.2 Ma) [29]. Therefore the Limacinoidea and Cresinae lineages must have diverged prior to the Thanetian. The age of Thecosomata was constrained to be less than 65 Ma as the group is understood to have evolved in the Cenozoic [31].
In recent classification [32] the family Cavoliniidae contains the subfamilies Cavoliinae, Clioinae, Cuvierininae and Creseinae. In contrast, in our topology the Cavoliniidae is paraphyletic, with a sister taxon relationship between Hyalocylis (Creseinae) and Limacina. This relationship was further supported by the possession of a shared indel by both taxa not present in any of the other species sequenced.