http://www.nps.gov/cave/planyourvisit/maps.htm [Accessed on

Background The phylogenetic position of pycnogonids is a long-standing and controversial issue in arthropod phylogeny. This controversy has recently been rekindled by differences in the conclusions based on neuroanatomical data concerning the chelifore and the patterns of Hox expression. The mitochondrial genome of a sea spider, Nymphon gracile (Pycnogonida, Nymphonidae), was recently reported in an attempt to address this issue. However, N. gracile appears to be a long-branch taxon on the phylogenetic tree and exhibits a number of peculiar features, such as 10 tRNA translocations and even an inversion of several protein-coding genes. Sequences of other pycnogonid mitochondrial genomes are needed if the position of pycnogonids is to be elucidated on this basis. Results The complete mitochondrial genome (15,474 bp) of a sea spider (Achelia bituberculata) belonging to the family Ammotheidae, which combines a number of anatomical features considered plesiomorphic with respect to other pycnogonids, was sequenced and characterized. The genome organization shows the features typical of most metazoan animal genomes (37 tightly-packed genes). The overall gene arrangement is completely identical to the arthropod ground pattern, with one exception: the position of the trnQ gene between the rrnS gene and the control region. Maximum likelihood and Bayesian inference trees inferred from the amino acid sequences of mitochondrial protein-coding genes consistently indicate that the pycnogonids (A. bituberculata and N. gracile) may be closely related to the clade of Acari and Araneae. Conclusion The complete mitochondrial genome sequence of A. bituberculata (Family Ammotheidae) and the previously-reported partial sequence of Endeis spinosa show the gene arrangement patterns typical of arthropods (Limulus-like), but they differ markedly from that of N. gracile. Phylogenetic analyses based on mitochondrial protein-coding genes showed that Pycnogonida may be authentic arachnids (= aquatic arachnids) within Chelicerata sensu lato, as indicated by the name 'sea spider,' and suggest that the Cormogonida theory – that the pycnogonids are a sister group of all other arthropods – should be rejected. However, in view of the relatively weak node confidence, strand-biased nucleotide composition and long-branch attraction artifact, further more intensive studies seem necessary to resolve the exact position of the pycnogonids.


Background
Pycnogonids, or sea spiders, comprising about 1163 extant species [1], are a special group of exclusively marine arthropods that are distributed from the intertidal zone to the abyssal depths in all the seas around the world. The pycogonids are traditionally classified as a sister group or even an ingroup taxon of euchelicerates, which comprises the arachnids, xiphosurans and extinct eurypterids [2,3]. However, the phylogenetic position of pycnogonids has been vigorously debated during the last two centuries, as summarized in a recent review of pycnogonid affinities by Dunlop and Arango [4]. According to this review, there have been a number of controversial hypotheses related to this issue. The pycnogonids have been regarded as either degraded crustaceans [5,6], a transitional form between crustaceans and arachnids [7][8][9], an isolated group unrelated to other arthropods [10][11][12], a sister group of aquatic arachnids or mites within Chelicerata [13][14][15], a sister group of euchelicerates [3,4,16,17], or a sister group of all other extant euarthropods (encompassing chelicerates, myriapods, crustaceans and hexapods) [18][19][20]. Because of rapid methodological improvements in phylogenetic systematics and related fields such as developmental biology, molecular biology and computational biology (bioinformatics), recent phylogenies have continuously attempted to resolve this issue and have eliminated many of the possible hypotheses, leaving only the two most plausible to be examined: namely that pycnogonids are either a sister group of euchelicerates or a sister group of euarthropods.
On the other hand, the hypothesis that the pycnogonids are a sister group of euarthropods (all extant arthropods except for pycnogonids) was initially formally proposed by Zrzavy et al. on the basis of all the available morphological and gene sequence data [18]. They named the clade of euchelicerates and mandibulates (excluding Pycnogonida) "Cormogonida;" they share the putative autapomorphic characteristic of a gonopore on the trunk. Pycnogonids have multiple gonopores on the bases of their legs. Although some authors have recovered this relationship [14,19,20,33,34], the presence of gonopores on the trunk in Cormogonida has been doubted as a plesiomorphic character [4].
Interestingly, a detailed neuroanatomical study of the Anoplodactylus sp. protonymphon lava [20] suggested that the chelifore of euchelicerates is not homologous to the chelicerae of pycnogonids, but rather that it is homologous to the protocerebral appendages ("great appendages") proposed for ancestral arthropods [35,36]. These authors suggested that the pycnogonids are a sister group of all the euarthropods. However, more recent studies of Hox expression patterns in Endeis spinosa [29,30] refuted the result of the neuroanatomical study [20]. The Hox expression data suggested that all the anterior-most appendages of all extant arthropods (chelifore of pycnogonids, chelicerae of euchelicerates and the first antennae of mandibulates) are homologous and deutocerebral, indicating that protocerebral appendages have been lost in all extant arthropods.
In recent years, comparison of complete mitochondrial genomes has become a very powerful tool for reconstructing arthropod phylogeny [15,[37][38][39][40][41]. Mitochondrial genomes contain a variety of useful phylogenetic information, such as gene orders and orientations, alternative start codons of protein-coding genes, transfer RNA and ribosomal RNA secondary structures, genetic code variations, and features of the control region for genome replication and transcription [42]. Typical metazoan mitochondrial genomes are circular, 14-16 kb in size, and encode 13 proteins, large and small subunit ribosomal RNAs (rrnL and rrnS) and 22 tRNAs [38,43]. The 13 polypeptides are involved in ATP synthesis coupled to electron transfer during O 2 consumption [ATP synthetase subunits (atp6 and atp8), cytochrome C oxidase subunits (cox1-cox3), apocytochrome b (cob) and NADH dehydrogenase subunits (nad1-nad6 and nad4L)].
The N. gracile mitochondrial genome has a peculiar gene order with extensive inversion of protein-coding genes and translocations of 10 tRNA genes, which are not typically found in arthropods [44]. These peculiar characteristics are probably restricted to the N. gracile lineage, since the partial E. spinosa mitochondrial genome possesses the arthropod ground pattern of mitochondrial gene arrangements. Furthermore, the N. gracile branch on the phylogenetic trees presented by Podsiadlowski and Braband [44] was very long, as was also the case for Acari. This may have given rise to an artifactual relationship between Acari and Pycnogonida, which suggests that N. gracile may not be an appropriate representative of the pycnogonids. Thus, additional pycnogonid mitochondrial genomes must be sequenced if such information is to be used to address the problem of the phylogenetic position of pycnogonids [45,46].
In this study, we present a complete new pycnogonid mitochondrial genome from a sea spider, Achelia bituberculata, belonging to the family Ammotheidae, which combines a number of anatomical features considered plesiomorphic with respect to other pycnogonids. The genome is characterized and compared to those of other arthropods, including two pycnogonids, N. gracile and E. spinosa. We attempt to use the data obtained in this study to address the long-standing and hotly-debated issue of the phylogenetic position of pycnogonids.

Genome organization
The mitochondrial genome of A. bituberculata is 15,474 bp in length, which is similar to those of arthropod mitochondrial genomes (Table 1). It includes 13 protein-coding genes (cox1-cox3, nad1-nad6, nad4L, atp6, atp8 and cob), 22 tRNA genes (trnA, trnR etc.), two ribosomal RNA genes (rrnL and rrnS) and one large non-coding region [= a plausible control region (CR); 977 bp], as typically found in metazoa ( Fig. 1 and Table 2). The components are compactly juxtaposed. There are even some overlapping nucleotides between trnC and trnW (1 bp), nad3 and trnA (3 bp), and trnE and trnF (1 bp) ( Table 2). The order and orientation of the gene arrangement pattern is identical to that of a horseshoe crab, Limulus polyphemus, except for the position of the trnQ gene (Fig. 2).

Comparison of gene arrangements
As shown in Fig. 2, the mitochondrial gene arrangement of A. bituberculata was compared with that of L. polyphemus, which is considered a putative ground pattern for arthropods [38], and with those of two sea spiders, N. gracile and E. spinosa. The identical gene arrangement pat-tern shown in A. bituberculata and L. polyphemus except for the trnQ gene is an interesting feature because the other sea spider, N. gracile, has an extensive inversion (from trnI to cox2) and ten tRNA gene translocations (Fig. 2). The partial genome of E. spinosa has the same gene arrangement pattern as those of L. polyphemus and A. bituberculata.

Base composition and codon usage pattern
As shown in Table 1, the overall A+T content in the mitochondrial genome of A. bituberculata is 77.0% (+ strand: A = 38%; C = 12%; G = 11%; T = 39%). This A+T content is much higher than those in L. polyphemus (67.6%) and the scorpions (64.5 and 67.0%), but is similar to those in Acari and Araneae (69.8-80.0%), including N. gracile (76.6%). Metazoan mitochondrial genomes usually present a clear strand bias in nucleotide composition [15,44]. Table 3 shows the AT-and CG-skews of each of the 13 proteincoding and 2 ribosomal RNA genes and those of the whole genome (total) in A. bituberculata mitochondria.
The results show no marked strand-bias in nucleotide composition. The AT-skew for all genes is negative whether they are on the (+) or the (-) strand. The CG-skew is positive for four genes and negative for five on the (+) strand. The CG-skew for all six genes on the (-) strand is negative. This is totally different from the clear strand-bias observed in N. gracile.
Investigation of the anticodon and codon usage patterns in the protein-coding genes of the A. bituberculata mitochondrial genome indicates that the nucleotides in the 3rd codon position are primarily 'A' or 'T' rather than 'C' or 'G'. Codons composed of only 'A' and 'T' are predominantly used, which seems to reflect the high A+T bias composition of the entire genome. A strong preference for A+T-rich codons is found in A. bituberculata as in mitochondrial protein-coding genes in general (Table 4). For example, the most frequently-used codons are L(UUA) (367 times) and I(AUU) (365 times). Additional file 1 shows the usage of leucine and serine codons in chelicer-  (363) are used than in the other chelicerates. The A+T content of the third codon positions of the 13 protein-coding genes is 85.99%, which is lower than in N. gracile (90.50%) and higher than in L. polyphemus (74.70%).
The start and stop codons of the 13 protein-coding genes of the A. bituberculata mitochondrial genome were identified by comparison with those of other known arthropods; particular attention was given to the amino acid sequence alignment in order to identify the most likely codon in ambiguous situations. As shown in Table 2, the start codons of 12 mitochondrial of the protein-coding genes (excluding cox1) are 'ATN', which is typical of most metazoan mitochondria [43]: 'ATG' is the start codon for nad2, nad4L, nad6, cox2, cox3 and cob; 'ATT' is the start codon found in nad5 and atp6; and 'ATA' is the start codon for nad1, nad3, nad4 and atp8. A less typical 'GAT' start codon appears in the remaining cox1 gene. On the other hand, only four of the protein-coding genes are terminated with a complete stop codon: 'TAA' in nad4L and cox3 and 'TAG' in nad1 and nad3. The remaining eight genes are terminated using truncated stop codons that are presumably completed by post-transcriptional polyadenylation [51]: 'TA' for nad4 and nad6, and 'T' for nad2, nad5, cox1, cox2, cob, atp6 and atp8.

Transfer RNA genes
All 22 of the tRNAs typically found in metazoan mitochondrial genomes were identified in A. bituberculata. All these genes are located among protein-and/or rRNA-coding genes ( Fig. 1 and Table 2). Twenty-two putative secondary structures were predicted from the tRNA gene sequences (Fig. 4). Most of the tRNAs are capable of forming the typical clover-leaf structure, with the exception of tRNA Ala (which lacks both the DHU and TψC arms), tRNA Ser(AGN) , tRNA Tyr and tRNA Val (which lacks only the DHU arm). The lengths of the TψC and DHU arms in the 18 tRNAs with a stable clover-leaf shape range from 4 to A circular map of the complete mitochondrial genome of the sea spider Achelia bituberculata (Family Ammotheidae) Figure 1 A circular map of the complete mitochondrial genome of the sea spider Achelia bituberculata (Family Ammotheidae). Protein and rRNA genes are abbreviated as follows: atp6 and atp8 (genes for ATPase subunits 6 and 8), cox1-cox3 (genes for cytochrome C oxidase subunits I-III), cob (gene for apocytochrome b), nad1-nad6 and nad4L (genes for NADH dehydrogenase subunits 1-6 and 4L), and rrnS and rrnL (genes for 12S and 16S rRNAs). All the 22 tRNA genes are located among protein-and/or tRNA-coding genes. The tRNA genes are named using single-letter amino acid abbreviations, with the exception of those coding for leucine and serine, which are named L1 for the tRNA Leu(UUR) (anticodon TAA) gene, L2 for the tRNA Leu(CUN) (anticodon TAG) gene, S1 for the tRNA Ser(UCN) (anticodon TGA) gene and S2 for the tRNA Ser(AGN) (anticodon GCT) gene. The arrows indicate the orientations of the gene components. Translocation of trnQ between the 12S rRNA (rrnS) and con-trol region (CR) observed in a sea spider, Achelia bitubercu-lata, two terrestrial spiders, Ornithictinus huwena and Habronattus oregonensis, and a collembolan, Tetrodontophora bielanensis Comparison of the mitochondrial gene arrangements of three sea spiders, Achelia bituberculata, Endeis spinosa and Nymphon gracile, and a horseshoe crab, Limulus polyphemus, which is considered to be an arthropod ground pattern Figure 2 Comparison of the mitochondrial gene arrangements of three sea spiders, Achelia bituberculata, Endeis spinosa and Nymphon gracile, and a horseshoe crab, Limulus polyphemus, which is considered to be an arthropod ground pattern. The mitochondrial gene arrangements in A. bituberculata (complete) and E. spinosa (partial) are the same as the arthropod ground pattern shown in L. ployphemus, except for the unexpected placement of trnQ (Q) in A. bituberculata. In N. gracile, however, the gene arrangement differs markedly from the other two pycnogonids. The two trnL1 and trnL2 genes are marked by L1 and L2, respectively. L1 and L2 in A. bituberculata and N. gracile have the orientation and position generally shown in the chelicerates and myriapods among the arthropods. Refer to the Fig. 1 legend for abbreviations of components.

Control region
A large noncoding region [= a plausible control region (CR); 977 bp] is found between trnQ and trnI. It includes three tandem repeats, the first and second of which are completely identical in length (315 bp) and sequence (Fig. 5A). However, the third repeat (282 bp) is 33 bp shorter than the other two at the 3'-end, and one nucleotide is different in the 315-bp repeat (Fig. 5A). The repeats have unique stem and loop structures that may play key roles in controlling the replication and transcription of the mitochondrial genome (Fig. 5B). The A+T content of the CR is 79%, which is similar to those of most chelicerates (Table 1).
Putative secondary structures of the 22 tRNAs identified in the mitochondrial genome of the sea spider Achelia bituberculata Figure 4 Putative secondary structures of the 22 tRNAs identified in the mitochondrial genome of the sea spider Achelia bituberculata.

The inclusion of Pycnogonida within Arachnida
As shown in Fig. 6, phylogenetic analyses based on the amino acid sequences of 12 protein-coding genes for 30 arthropods and 4 outgroups (1 onychophoran, 2 annelids and 1 mollusk) as shown in Table 1  Two previous studies based on the mitochondrial genome sequences of E. spinosa [15] and N. gracile [44] indicated that the Pycnogonida-Acari affinity might have resulted from biased features such as the strand-biased nucleotide contents and long-branch attraction. Interestingly, the A. bituberculata mitochondrial genome shows no marked strand-bias nucleotide composition (Table 3). Since most of the other chelicerates and N. gracile have strand-biased nucleotide contents [44]as well as extremely long branches on the phylogenetic tree (Fig. 6), the possibility remains that the appearance of Pycnogonida within Arachnida may result from artifacts of analysis due to the strand-biased and long-branching patterns.
To test whether the inclusion of Pycnogonida within Arthropoda results from long-branch attraction artifacts, we removed 9 long-branch taxa -all 7 Acari members and 2 terrestrial spiders (H. The present findings conflict with the two major hypotheses that Pycnogonida is a sister group of the monophyletic Euchelicerata (= Arachnida + Xiphosura) or that it is a sister group of the monophyletic Cormogonida (= Euchelicerata + Myriapoda + Pancrustacea). Nevertheless, the inclusion of Pycnogonida within Arachnida is compatible with some previous analyses based on combinations of 253 morphological characteristics and molecular data (complete 18S rDNA and the D3 region of 28S rDNA) [14], and the mitochondrial protein-coding genes [15,44], which indicated that Pycnogonida appears within Arachnida. More recently, Mallatt and Giribet [31] on the basis of nearly-complete 18S and 28S rDNA sequences, and Giribet et al. [32] on the basis of parsimony analyses of combined data for 7 nuclear genes, 2 mitochondrial genes and 375 morphological characters, both showed that Pycnogonida is united with Chelicerata with very weak support and Cormogonida is strongly rejected. The present result is consistent with theirs in that Cormogonida is rejected and Pycnogonida has an affinity with Chelicerata. However, they did not suggest that Pycnogonida may possibly be included within Arachnida; many morphological taxonomists may not agree with this point of view. Thus, it is likely that the inclusion of Pycnogonida within Arachnida, deduced only from mitochondrial genome data [15,44], will remain an open question.
In addition, Scorpiones are often considered to be "one of the most primitive arachnids" because they possess a number of plesiomorphic character states, even though their phylogenetic position remains controversial. Recent publications based on mitochondrial protein-coding genes [56,57] show that Scorpiones is a sister group of all other arachnids (Acari and Araneae). The present result (Fig. 6) is consistent with previous reports [56,57] showing that Scorpiones is a sister group of the clade of Pycnogonida plus Araneae plus Acari.
Primary and secondary structures of the control region in the mitochondrial genome of the sea spider Achelia bituberculata Maximum likelihood tree inferred from amino acid residues deduced from 12 mitochondrial protein-coding genes of 30 arthro-pods and 4 non-arthropods

Conclusion
This study presents the complete sequence of the mitochondrial genome from a pycnogonid, A. bituberculata (Family Ammotheidae). This complete sequence and a previously-reported partial sequence of Endeis spinosa show the typical gene arrangement patterns of arthropods (Limulus-like), but they differ from that of N. gracile. Thus, it is most likely that the peculiar arrangements found in N. gracile are only present in the family Nymphonidae or even in the genus Nymphon. Phylogenetic trees based on amino acid sequences deduced from mitochondrial protein-coding genes showed that Pycnogonida may be authentic arachnids (= aquatic arachnids) within Chelicerata, as indicated by the name 'sea spider', and suggest that the Cormogonida theory, which considers the pycnogonids a sister of all other arthropods, should be rejected. However, since node confidence values supporting the Pycnogonida inclusion within Arachnida are slightly low and the possibility of long-branch attraction artifact still remains, further more intensive studies seem necessary to resolve the problem.

Specimen collection and DNA extraction
The pycnogonid A. bituberculata was collected from the intertidal zone of the East Sea (Pohang) in Korea. This specimen was immediately fixed in 70% ethanol and brought to our laboratory where it was preserved in absolute ethanol and stored at -20°C. The specimen was washed with distilled water, and total cellular DNA was extracted from a single individual using the DNeasy Tissue Kit (QIAGEN Co., Germany).

PCR amplification and cloning
The entire mitochondrial genome was amplified by two polymerase chain reactions (PCR) using the following previously-reported primers: The PCR reactants were loaded on a 0.7% agarose gel and stained with ethidium bromide to visualize the bands using an ultraviolet transilluminator. PCR products of approximately 5 kb (rrnL-cox1) and 11 kb (cox1-rrnL) covering the entire mitochondrial genome were observed. The two PCR products were purified using the PCR Purification Kit (QIAGEN Co., Germany) and were ligated with the pGEM T-easy vector (Promega Co., USA). The ligation mixtures were transformed into Escherichia coli host strain DH5-α. Correct recombinant clones were selected by the blue/white colony selection method using X-gal and IPTG. Plasmid DNA was extracted using the AtmanBio plasmid purification kit (Takara Co., Japan), which was used for all subsequent sequencing reactions.

Sequencing and sequence analysis
Both strands of the purified plasmid DNA were sequenced using a primer-walking strategy. The sequencing reactions were performed using the dideoxy-nucleotide termination method with the big-dye terminator system and an ABI3700 model automatic sequencer. The complete pycnogonid mitochondrial genome [GenBank: AY457170] was annotated using several stand-alone programs including Clustal X [63] and GeneJockey II v. 1.6 (Biosoft Inc., UK) and the BLAST web-based program. Thirteen proteincoding genes and two rRNA genes were identified by their sequence similarities to previously published genes. Potential secondary structures of 18 of the 22 tRNAs were identified using tRNAscan-SE v. 1.1 [64], and those of the other 4 by eye (trnA, trnS2, trnY and trnV).

Phylogenetic analysis
For the present phylogenetic analyses we employed 30 arthropod ingroup taxa and 4 outgroup taxa, as listed in Table 1. The amino acid sequences of the 13 protein-coding genes were used. The 13 multiple alignment subsets of these sequences were created using a Clustal X multiple alignment program [63] under the default option. Only well-aligned, conserved alignment sites were extracted from each alignment subset using the Gblock program [65] with the default option. The atp8 gene was excluded from the following analyses because it is very variable in sequence and too short in length, so no region was selected by the Gblock program. The extracted conserved blocks were subsequently concatenated into a unified, single large alignment set with the Gblock program. Phylogenetic analyses inferred from nucleotide sequences are not presented because the resultant trees are unresolved in deep branches in spite of removing 3 rd -codon positions. It is likely that amino acid sequences are considered more appropriate for resolving deep branching patterns in this study.
The refined alignment (2232 aa positions in length) was subjected to two different tree-making algorithms: the maximum likelihood (ML) and Bayesian inference (BI) methods. Rather than using hierarchical likelihood ratio tests to select the best fitting model for the evolution of sequences, and to calculate the related parameter values (I and Γ), ProtTest ver. 1.3 [66] was used under Akaike Information Criterion (AIC) because it has several important advantages [67]. Among the 36 models implemented in this program, the best-fitting model selected was mtREV [68] with among-site substitution-rate heterogeneity described by a gamma distribution (Γ = 0.86) and a fraction of sites constrained to be invariable (I = 0.19). The ML analysis was carried out using PHYML v2.4.4 [69]. The bootstrap proportions (BP ML ; 100 replicates) of the ML tree were obtained by the fast-ML method using PHYML. BI analysis was carried out using the MrBayes v3.0b4 program [70] with the following options: 1,000,000 generations, 8 chains (2 hot and 6 cold) and a burn-in step of the first 50,000. Node confidence values of the BI tree were presented with Bayesian posterior probabilities (BPP). We applied bootstrap re-sampling procedures to the Bayesian approach since it is known that posterior probabilities place excessive confidence on a given phylogenetic hypothesis. Bayesian bootstrap proportions (BP BI ) were calculated as follows: 100 bootstrapped datasets were generated with the SEQBOOT program under the PHYLIP package Version 3.6b [71], and then each dataset was analyzed by MrBayes v3.0b4 [70], with four independent Markov chains run for 500,000 Metropolis-coupled MCMC generations, with tree sampling every 100 generations and a burn-in of 10,000 trees. Finally, the 100 Bayesian majority rule consensus trees were used to construct the BBP consensus tree with the CONSENSUS program using the PHYLIP package Version 3.6b [71].