Phylogenetic relationship within Cumacea (Crustacea: Peracarida) and genetic variability of two Antarctic species of the family Leuconidae

Summary: Phylogenetic hypotheses for the peracarid order Cumacea are scarce and have not provided a solution to the full extent. In the present study, a fragment of the mitochondrial 16S rDNA was used to erect a phylogenetic hypothesis for three cumacean families, Diastylidae, Bodotriidae and Leuconidae, along with intra-family relationships of the latter. The Cumacea resolved monophyletic with tanaids and isopods as outgroup taxa. The Diastylidae were the only family with good support for monophyly. The genus Leucon resolved paraphyletic, whereas the subgenus Crymoleucon was monophyletic. Furthermore, the genetic structure was analysed for two leuconid species, Leucon antarcticus Zimmer, 1907 and L. intermedius Mühlenhardt-Siegel, 1996, from the Weddell Sea and the Ross Sea. The two species showed different patterns of in- traspecific genetic variability. In contrast to L. intermedius , a bimodal distribution of pairwise genetic distances was observed for L. antarcticus , which is correlated with geographical and depth distributions between the Ross Sea and the Weddell Sea. Although a clear evaluation of cryptic speciation in these species requires additional work on more specimens from more geographic regions and broader depth ranges, differences shown in the sequences of 16S rDNA can only be explained by genetic separation of populations between the Weddell Sea and the Ross Sea for an extended period of time.


INTRODUCTION
Cumaceans are a group of peracarid crustaceans predominantly inhabiting marine soft bottom habitats. They can occur in large numbers (e.g. San Vicente at al. 1997, Rehm et al. 2007 and are an essential component of the benthic fauna, thus being an important food source for demersal fish and other macrofauna (e.g. Cartes 1993, Schlacher andWoolbridge 1996). The first report of an Antarctic cumacean was published by Sars (1873). Additional descriptions of five Antarctic cumaceans followed during the next decade (Sars 1887). Today, about 100 cumacean species from all extant families of Cumacea (Bodotriidae, Ceratocumatidae, Diastylidae, Gynodiastylidae, Lampropidae, Leuconidae, Nannastacidae and Pseudocumatidae) are described for the Antarctic and sub-Antarctic (Błażewicz and Heard 1999, Mühlenhardt-Siegel 1999, Pabis and Błażewicz-Paszkowycz 2011. However, knowledge about Antarctic cumaceans is still incomplete and restricted to species inventory, diversity and biogeography. Hypotheses on the evolutionary history of cumacean families have been proposed by Zimmer (1941) and Lomakina (1968). Both regard the Lampropidae and Diastylidae as basal taxa, but their interpretations differ regarding the more derived families. Nevertheless, both authors are of the opinion that the pleotelsonbearing families are most derived. Testing phylogenetic hypotheses has been difficult for cumaceans because characters used for the taxonomy of this peracarid taxon are inconsistent within and often among families. Haye et al. (2004) discuss the monophyly of the pleotelsonbearing Bodotriidae, Leuconidae and Nannastacidae, as indicated by the phylogenetic analysis of amino acid sequences from the mitochondrial cytochrome oxidase I gene and morphological characters. With respect to the "pleotelson clade", their findings are in accordance with Zimmer and Lomakina, but monophyly (with reasonable support) was confirmed only for the families Gynodiastylidae and Lampropidae by molecular data. The present study aimed to investigate the phylogenetic relationship of three cumacean families and within the family Leuconidae using a fragment of the mitochondrial LSU gene (16S rDNA).
Furthermore, genetic variation in Antarctic species of the genus Leucon is studied to reveal possible patterns of cryptic speciation, which have been demonstrated for Antarctic isopod (Held 2003, Held and Wägele 2005, Raupach and Wägele 2006, amphipod (Baird et al. 2011), mollusc (Allcock et al. 1997, Linse et al. 2007, and crinoid species (Wilson et al. 2007). The discoveries of cryptic speciation indicated that diversity of Antarctic zoobenthic organisms in terms of species richness is much higher than previously believed. Therefore, circum-Antarctic distribution, which was postulated for many taxa, is not valid for a variety of taxa because many nominal species with circumpolar distributions may in reality consist of a series of cryptic species with more local distributions (Arango et al. 2011, Dornburg et al. 2016, Beermann et al. 2018. Patterns of cryptic speciation observed in shallow water species inhabiting the Antarctic continental shelf are assumed to be caused by geographic isolation and mainly glaciation processes on a Milankovitch timescale, which might have led to isolated shelters on the Antarctic shelf (Clarke and Crame 1989, Thatje et al. 2005. Species with pelagic larvae or dispersal life stages are commonly assumed to overcome the barriers separating 'biogeographic islands' on the Antarctic shelf, and thus ensuring gene flow between isolated populations (Thatje 2012), although unexpected evidence for strong population connectivity has also been found in benthic species without obvious dispersal capabilites (Leese et al 2010, Fraser et al. 2013. As cumaceans belong to the brooding crustacean supraorder Peracarida the hypotheses presented above is tested for evidence in this taxon.

Source and preservation of material
Antarctic Cumacea were collected during the 19 th Italian Antarctic expedition with RV Italica along the coast of Victoria Land in the Ross Sea (Rehm et al. 2007). Further material was obtained from the BEN-DEX (ANT XXI-2) expedition and ANDEEP cruises I and II to the Scotia-Arc region, the Antarctic Peninsula and the Weddell Sea carried out with RV Polarstern in the years 2002 and 2004. The species Diastylis rathkei was sampled in the Kiel Bay in the Baltic Sea (Table 1). The material was sorted by hand from trawled gear (Rauschert dredge and epibenthos sledge) using Antarctic Peninsula 3685 59°39.9'S 57°53.9'E 1 HQ450554.1 a dissecting microscope. Samples were preserved in pre-chilled 80% (0°C, -80°C, resp.) ethanol. Samples were obtained from depths of between 15 and 3685 m. Samples were stored at -30°C for at least 4 months and were kept at 5°C until further processing. During the cruise with RV Italica, samples were stored at -80°C during the first four days. Extraction and sequencing of cumacean material collected during the BENDEX expedition was less successful than treating the material from the campaign with RV Italica. In contrast to sample processing during the BENDEX expedition, during which samples were fixed with 0°C cold ethanol, samples were fixed at -80°C onboard RV Italica. Deep temperatures at the beginning of the fixation might be the reason for better results during molecular work, so we suggest cooling newly collected material at -80°C during the first weeks of fixation.

Molecular work
DNA was extracted from individual legs, from the pleon without telson and uropods, or from entire smaller specimens. The following modifications were applied to the protocol of the QIAamp DNA Mini Kit, which was used for DNA extraction: the spin column loaded with elution buffer was incubated for 5 min at 70°C before elution of the DNA and the volume of the elution buffer was decreased from 200 to 50 µl in order to increase the concentration of eluted DNA.

Primer choice and creation
For DNA amplification the broadly applicable primers 16Sar 5'-CGCCTGTTTATCAAAAACAT-3' and 16Sbr 5'-CCGGTCTGAACTCAGATCACGT-3' (Palumbi et al. 1991) were used. Despite the general application of these primers on arthropod taxa, amplification of cumacean DNA was weak. Therefore, cumacean-specific primers were designed on the basis of the sequences obtained in our pilot study and from GenBank. The program 'Fast PCR' (Kalendar 2003) was used to construct primers. Primers ALh (5'-GTAC-TAAGGTAGCATA-3') and CLr (5'-ACGCTGT-TAYCCCTAAAGTAATT-3') were designed for the cumacean family Leuconidae in highly conserved regions of the 16S gene and used during this study. The amplification protocol for ALh and CLr was 2 min at 94°C for initial denaturing, 38 cycles of 20 s at 94°C, 10 s at 46°C, and 1 min at 65°C, followed by 8 min for final extension.

DNA sequencing
PCR products were purified with the QIAquick PCR-purification kit of Qiagen, Hilden, Germany. To achieve higher concentrations of purified DNA, only 30 µl of elution buffer were used. DNA purity and amount of DNA were controlled on an ethidium bromide-stained 1.5% agarose gel. Cycle sequencing was performed according to the manufacturer's instructions of the BigDye Terminator v3.1 kit of Applied Biosystems (ABI) using a ABI 3130 sequencer (96°C 1 min initial denaturing, 30 cycles of 10 s 96°C, 50 s 50°C, 4 min 60°C). In general, 1 to 3 µl of purified DNA was used for cycle sequencing on an Eppendorf Master Cycler (4 µl was used for samples with very low DNA concentrations). Surplus dye was removed with the DyeEx 2.0 spin kit (Qiagen) and 10 µl samples were denatured for 3 min at 95°C with 10 µl ABI HighDiye formamide (Applied Biosystems). The samples were kept on ice prior to sequencing.

Sequence alignment and phylogenetic analysis
Raw electropherograms from the sequencer were assembled using the programs Pregap4 and Gap4 of the Staden package (Staden et al. 1989). For the first alignment of the contig sequences and additional sequences of the mitochondrial 16S rDNA gene, downloaded from GenBank (Table 2), the "ClustalW Multiple alignment" option of the program BioEdit (Hall 1999) was used. The alignments were improved further manually by identifying secondary structure elements of the homologous molecules in Drosophila melanogaster (mitochondrial ribosomal LSU, Accession No. X53506, Gutell et al. 1993). Loop regions were locally re-aligned using a hidden Markov model implemented in the program ProAlign version 0.5 (Löytynoja and Milinkovitch 2003). Default parameters were used for alignment sampling with 1000 replicates, if not stated otherwise. Estimated nucleotide frequencies were A=0.366, C=0.149, G=0.171, T=0.3131. The analysis included sites which could only be aligned in the ingroup or within the family Leuconidae. Corresponding sites of the outgroup or cumaceans other than Leuconidae, respectively, were substituted with gaps. Sites that were still ambiguously aligned at this stage were excluded from analysis.
Phylogenetic analyses were performed using maximum parsimony, maximum likelihood and Bayesian approaches. Bayesian analyses were performed with MrBayes 3.1.2 (Huelsenbeck and Ronquist 2001) on preset parameters, whereas for maximum likelihood and maximum parsimony analyses the programs RAxML 7.0.4 (Stamatakis 2006) and PAUP* 4.0b10 (Swofford 2003), respectively, were used. We used the general time-reversible model with invariable sites and gamma distribution (GTR+I+Γ), whose parameters were estimated using the program ModelTest version 3.7 (Posada and Crandall 1998) and the Akaike information criterion. The ratio of invariable sites was 0.1935, the gamma distribution shape parameter was 0.7813 and the base frequencies were A=0.3629, C=0.1332, G=0.1742 and T=0.3297. Rates for the six substitution types estimated from the dataset were AC=3.2491, AG=13.2695, AT=5.3873, CG=2.2114, CT=21.8755, and GT=1.0000).
The settings for maximum likelihood and maximum parsimony were a heuristic search with random sequence addition (10 replicates) and tree bisection reconnection. The robustness of the tree topologies was assessed with bootstrapping with 1000 and 10000 replicates for maximum likelihood and maximum parsimony, respectively.

Outgroup selection
Malacostracan phylogeny is far from being solved at the moment, as is the phylogeny of the Peracarida (for review see: Davis 2001, Richter andScholz 2001). Both Malacostraca (Wills 1998) and Peracarida (Schram andHof 1998, Watling et al. 1991) are even suggested to be polyphyletic. Siewing (1963) suggested Tanaidacea and Isopoda as a sister group to Cumacea, whereas Watling (1999) and Schram (1986) placed the isopods more basally. After Schram and Hof (1998), the tanaids are a more basal group, but the tree shown is not fully resolved, with a possible close relationship between Tanaidacea and Cumacea, while there is no closer relationship between Isopoda and Cumacea. Resent phylogenetic analysis by Richter and Scholz (2001) indicate a possible sister relationship of Tanaidacea and Isopoda to Cumacea. Following these authors, the exact position of the Cumacea and related taxa in the phylogenetic tree is still not resolved. Still, several morphological characters, such as the carapace as a respiratory structure (Watling 1999), the dorsally folded embryo, in addition to the manca stage and similar formation of the midgut (Hessler 1983), support a closer linkage of Cumacea, Isopoda and Tanaidacea. Following these common characters, Isopoda and Tanaidacea were selected as an outgroup to Cumacea.

RESULTS
The fragment amplified with the primers ALh/CLr varied between 255 and 256 bp in length, while those amplified with the primers 16Sa/16Sb ranged from 470 to 472 bp. The alignment is based on sequences obtained with the primers 16Sa/16Sb; sequences of Leucon antarcticus and Leucon rossi were solely obtained using the primers ALh/CLr. The total length of the alignment was 523bp. After the exclusion of ambiguously aligned positions, 382 remained, of which 106 were constant and 54 were parsimony-uninformative.
Maximum parsimony resulted in a tree with most taxa included in only one polytomy. Transition/transversion ratios from 0 to 10 (values presented in Fig.  1 were calculated with a ratio of 3) were tested, all yielding similar trees with differences only in the boot- Fig. 1. -Bayesian analysis consensus tree (50% majority rule) based on 16S rDNA. The GTR+I+Γ model was used according to the Akaike information criterion test. First numbers represent the portions of sampled trees, in which the corresponding node was found (posterior probability). Second numbers and third numbers represent bootstrap values of maximum likelihood (1000 iterations) and maximum parsimony (10000 iterations) analysis, respectively (Outgroup taxa see Table 2). Values below 0.5 or 50 are not shown. Species are grouped by colour shading.
strapping support. The cumacean family Diastylidae was the only well-supported monophylum (bootstrap support 83%) and Cumacea were not monophyletic because Apseudes latreiellei (Tanaidacea) was placed on the same branch. Trees obtained with the Bayesian and maximum likelihood methods were similar, but in general nodes had a weaker maximum likelihood support (Fig. 1). The Bayesian analysis indicated that the Cumacea are monophyletic, supported by a Bayesian score (BS) of one. Furthermore, the Diastylidae are well supported (0.99 BS), while Bodotriidae were not resolved as a monophylum. The Leuconidae are weakly supported as monophyletic family, but with a BS of only 0.67. At this node the tree is trichotomous with Leucon assimilis, Eudorella pusilla and the remaining Leuconidae. The latter are fairly well supported (0.97 BS). The subgenus Crymoleucon is monophyletic and also well supported (0.88 BS), but fails to resolve in the maximum likelihood tree. Species pairs, which exhibit high BS, are Leucon antarcticus/L. rossi and Diastylis sculpta/D. rathkei.
The sequence belonging to species of the genus Leucon are split into three groups (Fig. 2) when compared with pairwise p-distances. The first group comprises within-species comparison with p-distances from 0 to 0.05, whereas the second group gives the minimum distance (0.20-0.21) of interspecific variation of the two closely related species Leucon antarcticus and L. rossi (compare Fig. 1). Interspecific distances of the remaining species are confined to the third group (p-distance 0.30-0.36). Intraspecific variation in the 16S rDNA of the two species L. antarcticus and L. intermedius follows different patterns. Intraspecific pdistances of L. intermedius (Fig. 3A) range from 0 to 0.033, while sequence similarity of L. antarcticus (Fig.  3B) shows higher variation (0-0.052) and a bimodal distribution with no intermediate sequence, correlating to geographical distance and depth distribution. Pairs of sequences with p-values from 0 to 0.014 were obtained from specimens collected either in the Ross Sea (depth ranging from 316 to 358 m) or in the Weddell Sea (900 m), whereas p-distances from 0.038 to 0.052 were observed between these groups.

Phylogenetic analysis of 16S rDNA
The model test recommended a complex model (GTR+I+G) for the present dataset. This is reflected by the fact that maximum likelihood and Bayesian methods lead to more resolved tree topologies than maximum parsimony. Maximum parsimony describes observed changes of characters and the method does not consider complex evolutionary assumptions, which are contained in the GTR model. According to the rescaled consistency index (0.0980) calculated with the program PAUP*, a certain homoplasy is indicated for the data set. Consequently, the result of maximum parsimony is regarded as less informative and will not be discussed further.
Tree topologies observed from Bayesian and likelihood analyses both show that cumaceans including the  families Diastylidae, Bodotriidae and Leuconidae are monophyletic with regard to the outgroup. Additionally, the Diastylidae appear monophyletic. In the phylogenetic analysis of molecular data from the cytochrome oxidase I (COI) by Haye et al. (2004), the Bayesian and maximum likelihood methods, in contrast to maximum parsimony, could not confirm monophyly for the Cumacea. The authors assumed that this is due to the low taxon number of Pseudocumatidae represented in their study, which do not group with other cumaceans. COI data suggest that the Diastylidae may be paraphyletic. As the number of diastylid taxa was less than half in the present study, we cannot rule out that 16S data might also prove paraphyly for a greater number of Diastylidae. Nevertheless, Haye et al. (2004) point out that constraining the Diastylidae to be monophyletic results in a tree that is not significantly longer than the Bayesian tree.
Bodotriidae (two branches) and Leuconidae (weakly supported) are part of the same polytomous node during the present study. Therefore, the result of the COI data could not be confirmed, showing that Bodotriidae were paraphyletic with the other pleotelsonbearing families, Leuconidae and Nannastacidae, nested within. A "pleotelson clade" has very low support in both studies. On the other hand, this clade is confirmed by morphological data with the three families monophyletic each and the Nannastacidae as a possible intermediate taxon between the more basal Leuconidae and derived Bodotriidae (Haye et al. 2004).
The genus Atlantocuma was originally placed in the family Bodotriidae (Băcescu and Muradian 1974). Jones (1984) mentioned the nannastacid-like character of the species, but preferred to leave it as an aberrant form within the Bodotriidae, while Haye (2002) used the taxon as an outgroup in the phylogenetic analysis of the Bodotriidae because it grouped as a sister taxon to the nannastacid genera Cumellopsis and Scherocumella. The recent morphological analysis of Bodotriidae (Haye 2007) does not include Atlantocuma in the Bodotriidae. In the present study Atlantocuma is a sister taxon to Cyclaspis, so a close relationship of Atlantocuma to the Bodotriidae is highlighted. Nevertheless, the placement of Atlantocuma cannot be solved finally since no sequences of 16S rDNA for the family Nannastacidae were available.
Monophyly of the Leuconidae is only weakly supported by the data presented here, but within the family a monophyletic group comprises the monophyletic subgenus Crymoleucon and an undescribed species of the subgenus Leucon (pers. comm. Mühlenhardt-Siegel). The tree topology suggests good evidence that the subgenus Leucon is paraphyletic because L. assimilis also belongs to the subgenus Leucon. The species L. antarcticus and L. rossi, which represent a monophyletic group, are also closely related morphologically. Besides decreasing size of the dorsomedial teeth to the posterior end of the carapace, the species can be distinguished by the shape of the pseudorostrum, which is blunt in L. rossi and tipped and slightly upturned in L. antarcticus, as well as by a spine present on the first article of the exopod of the first pereopod (Rehm and Heard 2008).
Phylogenetic information provided during this study is reliable partially within the Leuconidae, in delimiting Cumacea from the outgroup, and in the monophyly of the Diastylidae with respect to the other ingroup taxa. Some authors have suggested that the Diastylidae are the most derived cumacean family (Băcescu and Petrescu 1999), while Lomakina (1968) and Zimmer (1941) placed this family next to the Lampropidae at the basis of the Cumacea. A derived position of the Diastylidae was not confirmed by the results of the present study. Moreover, phylogenetic analyses of morphological characters and the cytochrome oxidase I gene presented by Haye et al. (2004) both indicate a more basal position of the Diastylidae. Therefore, the assumption of Zimmer and Lomakina considering the general position of the Diastylidae has to be regarded as confirmed.
For a well-founded analysis of cumacean families, more taxa of all families need to be analysed. Since cumaceans represent a relatively old group, more genes including more slowly evolving ones than the mitochondrial 16S rRNA gene are needed to fully resolve the phylogeny of higher cumacean taxa. The more slowly evolving 18S gene is a possible candidate for further investigations; in addition more genes should be included to enhance the resolution of cumacean phylogeny (Hillis et al. 1996) .

Variation in 16S rDNA of Antarctic Leuconidae
It should be highlighted that no previous molecular information on the Cumacea from Antarctica is available. The mitochondrial 16S sequences of Leucon antarcticus show a pronounced barcoding gap, i.e. a bimodal distribution pattern of pairwise genetic distances with no intermediate values, indicating the possible existence of cryptic species (Fig. 3B; see Held 2003). The sequences fall into two groups: one was obtained from the Weddell Sea at a depth of 900 m, whereas the other was obtained from the Ross Sea at about 350 m water depth. From the Weddell Sea only two sequences were available for genetic analysis; therefore, it is possible that intermediate sequences exist. Even if the results represent true haplotype distribution, intermediate sequences could exist in geographically intermediate populations of L. antarctica. Nevertheless, 14 sequences of the 'Ross Sea haplotype', sampled at two stations with a distance of 340 km, vary only in one position in the alignment, whereas nine positions are different to the "Weddell Sea haplotype".
The second criterion for distinguishing cryptic species is the differentiation level of the gene, which should be in the range of clearly separated but closely related species. The differentiation between L. antarcticus and L. rossi (Fig. 2), which are closely related species (see phylogenetic analysis), is less than that between L. antarcticus and other leuconid species, but still five times higher than that within the two observed haplotypes of L. antarcticus. The study of 16S rDNA of brachyuran crabs from Jamaica has shown that cryptic speciation may take place at lower levels than revealed for L. antarcticus (Schubart and Koller 2005). On the other hand, p-values observed for the differentiation of cryptic Antarctic isopod species (Held 2003, Held andWägele 2005) is at the upper range or even higher than in L. antarcticus. A further indication for cryptic speciation might be the different distance pattern observed in L. intermedius, with the upper limit of p-values at 0.033 and intermediate values (Fig. 3B). The third criterion mentioned by Held is not applicable, as it demands constantly high level of differentiation in sympatry.
Morphological records of L. antarctica are ambiguous. The species was first described by Zimmer (1907) from the East Antarctic (compare also Zimmer 1913) and by Calman (L. australis) from the Ross Sea in the same year. Ledoyer described the species for a third time from the Weddell Sea. Zimmer presented a more detailed description, whereas the descriptions of Calman and Ledoyer are vague in several aspects. Zimmer mentioned five lateral spines on the carapace, while no spine is mentioned in Calman's description. For L. antarctica (sensu Ledoyer 1993, Calman 1907 no spine is mentioned either, but in the drawing one spine is depicted. All descriptions cover only a part of the appendages. Moreover, due to low quality of the drawings and insufficient descriptions given in the text, it is not possible to judge possible geographical differences reflected in morphology. Specimens from the Weddell Sea used for the present study bear a similar spine pattern on the carapace as specimens from the Ross Sea. Both populations show some variation, which does not allow a differentiation according to the lateral spines of the carapace.
In conclusion, morphological descriptions of L. antarctica are indistinct, the number of samples of the 16S rDNA gene and the geographical distribution of sample sites are not sufficient to allow a final evaluation of genetic variability and cryptic speciation. Still, differences exist in the sequences of 16S rDNA, which can only be explained by genetic separation of populations from the Weddell Sea and the Ross Sea for an extended period of time. Further studies with more sequences and extended geographical range of samples will provide a more detailed image of the genetic diversity of this species and finally bring the stage of speciation to light.