The evolutionary history of Stomatopoda (Crustacea: Malacostraca) inferred from molecular data

Taxon sampling and DNA sequencing

Thirty-eight stomatopod species, representing six superfamilies, 12 families, and 28 genera, were included in this study (Table S1). Sequences were either generated de novo or obtained from GenBank. Sequences obtained from GenBank were collected from single specimens in most cases; otherwise they were from multiple specimens collected from the same location. Tissue samples were collected from specimens provided by the Australian Museum (AM) and Muséum National d’Histoire Naturelle Paris (MNHN). Outgroup taxa were selected on the basis of previous analyses of stomatopods, sequence availability, and uncertainties over the phylogenetic position of Hoplocarida within Malacostraca (Richter & Scholtz, 2001; Miller & Austin, 2006; Ahyong & Jarman, 2009; Liu & Cui, 2010). Representatives of five other malacostracan lineages were included, representing Decapoda (Homarus americanus and Homarus gammarus), Anaspidacea (Anaspides tasmaniae), Euphausiacea (Meganyctiphanes norvegica), Peracarida (Neomysis americana and Neomysis integer), and Phyllocarida (Paranebalia longipes). GenBank accession numbers for all sequences are given in Table S1.

DNA was extracted using a modified version of the Chelex rapid-boiling procedure (Walsh, Metzger & Higuchi, 1991; Ahyong & Jarman, 2009). Three molecular markers were selected for amplification, based on the study by Ahyong & Jarman (2009). Regions of two mitochondrial genes (12S and 16S) and one nuclear gene (D1 region of 28S) were amplified using three sets of primers (Table S2). PCR cycle conditions differed for each primer set. The 12S cycle parameters followed those of Mokady et al. (1994) and Mokady & Brickner (2001). The 16S cycle parameters followed Ahyong & Jarman (2009), whereas the 28S cycle parameters followed Schnabel, Ahyong & Maas (2011).

Sanger sequencing was performed by Macrogen (Seoul, South Korea). Contigs were aligned using the default assembly parameters in Sequencher v5.0.1 (Gene Codes Corporation). These data were combined with published sequences for 16S, 18S, cytochrome c oxidase subunit 1 (CO1), and 28S (D2–D7b and D9–D10 regions). Sequences were aligned separately for each gene using MUSCLE v3.8.31 (Edgar, 2004), with poorly aligned sites removed following visual inspection. Additionally, owing to large inconsistencies with the sequences from other taxa, we excluded small sections of the 28S sequences from Squilla rugosa (sites 3,737–3,895) and Kempella mikado (sites 3,856–3,911). Sequence alignments are available online (https://github.com/caravanderwal/Stomatopoda).

To check for saturation of nucleotide substitutions, we analysed each gene separately using Xia’s test in DAMBE v6 (Xia, 2013). We also performed separate tests of saturation for each of the three codon positions in CO1. We found evidence of saturation at the third codon positions in CO1, but not in any other subsets of the data (Table S3). Therefore, we excluded the third codon sites of CO1 from our phylogenetic analyses.

Phylogenetic analysis

We performed phylogenetic analyses of our data set using both maximum likelihood and Bayesian inference. The best-fitting data-partitioning scheme and substitution models were selected using PartitionFinder v1.1.1 (Lanfear et al., 2012). The optimal partitioning scheme split the data into four subsets: (i) 12S and 16S; (ii) 18S; (iii) first and second codon positions of CO1; and (iv) 28S. For these four subsets, the best-fitting substitution models were GTR + G, GTR + I + G, GTR + I + G, and GTR + I + G, respectively.

The maximum-likelihood phylogenetic analysis was performed in RAxML v8.2.4 (Stamatakis, 2014). Given that this software has the constraint of using the same type of substitution model across all data subsets, we assigned a GTR + G model to each of the four subsets. Bootstrap support values were estimated using 1,000 pseudoreplicates. Two replicates of the analysis were performed to check for local optima, each using 10 random starts.

The phylogeny and divergence times were co-estimated using the Bayesian phylogenetic software BEAST v.1.8.4 (Drummond et al., 2012). To check for rate variation across branches, we ran analyses using a strict clock and an uncorrelated lognormal relaxed molecular clock (Drummond et al., 2006). To check the sensitivity of the results to the choice of tree prior, we conducted analyses using a Yule speciation model and using a birth-death model for the tree prior. All model comparisons were done using marginal likelihoods, calculated using the stepping-stone estimator (Xie et al., 2011). To examine the joint prior distribution of divergence times, we performed a BEAST analysis in which we drew samples from the prior (with the tree topology fixed to match the maximum-clade-credibility tree inferred using the sequence data).

Posterior distributions of parameters were estimated using Markov chain Monte Carlo (MCMC) sampling over a total of 10⁸ steps, with samples drawn every 10⁴ steps. The initial 10% of samples were discarded as burn-in. Convergence was checked by running the analysis in duplicate and by visualizing the results in the program Tracer v1.6 (Rambaut & Drummond, 2014), which showed that the effective sample sizes of all parameters were above 200. We used TreeAnnotator v1.8.4, part of the BEAST package, to identify the maximum-clade-credibility tree and to scale the node ages to the median posterior estimates.

Fossil calibrations

Fossil evidence was used to constrain the ages of seven nodes in the tree (Table 1). We implemented each of these in the form of a uniform prior, which allows the node an equal probability of taking any age within the specified bounds (Ho & Phillips, 2009). Additionally, we performed an analysis using exponential priors. The first calibration point provided a maximum age for the separation of Stomatopoda from the other malacostracans. Phyllocarids are the oldest malacostracans in the fossil record, appearing in the Cambrian at least 500 Ma (Sepkoski, 1998; Collette & Hagadorn, 2010). Therefore, the maximum age constraint for the age of the root was set to the beginning of the Cambrian period (541 Ma).

The second calibration point provided a minimum age for crown stomatopods. We note that a Devonian age for Palaeostomatopodea has been proposed based on interpretation of Eopteridae as probable palaeostomatopodeans (Schram et al., 2013). Given that Eopteridae is based on fragmentary fossils of uncertain identity, we take a more conservative approach and base our calibration on unambiguous hoplocarid fossils. This age was determined by the appearance of Palaeostomatopodea and the proto-stomatopods Archaeostomatopodea in the Carboniferous approximately 313 Ma, particularly the well-described Daidal acanthocercus fossil (Schram, 2007). These are regarded as representatives of the stem lineage leading to crown-group Stomatopoda.

Five calibration points were implemented as minimum bounds on the ages of unipeltatan stomatopod superfamilies, families, and genera. We selected reliable fossils of Gonodactyloidea, Lysiosquilloidea, and Squilloidea. The oldest putative gonodactyloid fossil, Paleosquilla brevicoxa from the Cretaceous, could not be used in our analysis. The incompleteness of Paleosquilla brevicoxa together with the non-monophyly of Gonodactyloidea, prevented a confident assignment of the fossil taxon to any particular gonodactyloid clade. However, three fossil calibrations were implemented for Gonodactyloidea. These provided minimum age constraints on extant genera within the families Gonodactylidae, Hemisquillidae, and Pseudosquillidae (Schram et al., 2013). One calibration for Lysiosquilloidea was selected on the basis of the sister relationship between Lysiosquilla and Lysiosquillina (Ahyong & Harling, 2000). The oldest squilloid fossil, Ursquilla yehoachi, which belongs to either the Squilla or Oratosquilla groups (Ahyong, 2005; Haug et al., 2013), was used to specify a minimum age constraint for the divergence between Squilloidea and Parasquilloidea. Finally, a minimum age constraint was placed on the root of the tree, representing the divergence between Paranebalia (Leptostraca) and all other taxa, based on the oldest known leptostracan fossil, Rhabdouraea bentzi, from the Upper Permian (Schram & Malzahn, 1984).

Topology tests and ancestral state reconstruction

We used likelihood and Bayesian methods to test the monophyly of three groups: Gonodactyloidea, smashers, and smashers and spearers. To do this, we repeated the phylogenetic analyses in RAxML and BEAST as described above, but with the monophyly of each of the three groups being enforced in turn. The likelihood scores of the constrained and unconstrained tree topologies were compared using the approximately unbiased (AU) test in CONSEL (Shimodaira & Hasegawa, 2001). Bayesian tests of topology were done by calculating Bayes factors, based on marginal likelihoods calculated using the stepping-stone estimator in BEAST. Bayes factors were interpreted according to the recommendations of Kass & Raftery (1995).

Ancestral states were reconstructed using maximum likelihood to study the evolution of the advanced stomatopod eye. The number of midband rows and the shape of midband ommatidia were reconstructed at internal nodes in the tree using the program Mesquite v3.03 (Maddison & Maddison, 2015). Character states for modern taxa were recognized as categorical data and sourced from the literature (Ahyong, 1997; Ahyong, 2001; Porter et al., 2010). To trace the evolution of these visual characters, we used the maximum-clade-credibility tree from our Bayesian phylogenetic analysis, which was pruned using TreeGraph 2 (Stöver & Müller, 2010) to include only the species with known states for each of the characters.

Stomatopod phylogeny

Bayesian and maximum-likelihood analyses of the three mitochondrial and two nuclear markers yielded almost identical estimates of the stomatopod phylogeny, although with differing levels of node support (Fig. 1). Hemisquilla was consistently placed as the sister group to all other extant stomatopods. The analyses yielded moderate to low support for nodes along the backbone of the tree. Each superfamily was inferred to be monophyletic with high posterior probability, with the exception of a polyphyletic Gonodactyloidea. Smashing stomatopods were placed in a relatively nested position within the tree. Although smashing had a monophyletic origin, a clade containing spearers, Pseudosquillidae and Bathysquillidae, was nested among the smashers (albeit with weak support).

Five of the six stomatopod superfamilies included in our study were found to be monophyletic, the exception being the polyphyletic Gonodactyloidea (Fig. 1). Topology tests were used to evaluate support for the monophyly of three groups: Gonodactyloidea (with and without Hemisquilla), smashers, and smashers and spearers. The AU test failed to reject the hypothesis of monophyly in all four cases (P > 0.05; Table 2). In contrast, the Bayesian analyses consistently favoured the unconstrained tree topology over the constrained topologies, except with respect to monophyletic smashers and spearers, with log Bayes factors ranging from −2015.6 to 4346.7 (Table 2). We note that this result has been obtained even though Bayes factors tend to be biased towards the hypotheses involving constrained tree topologies (Bergsten, Nilsson & Ronquist, 2013).

Estimates of divergence times

Comparison of marginal likelihoods revealed strong support for an uncorrelated lognormal relaxed clock (lnL = − 28, 550.8) over a strict clock (lnL = − 28, 668.3), indicating the presence of rate variation across branches. In addition, we found support for a Yule tree model (lnL = − 28, 550.8) over a birth-death model (lnL = − 28, 556.0). The choice of tree prior did not have a measurable impact on the resulting estimates of divergence times. Furthermore, the marginal prior distribution of node times differed from the posterior estimates, indicating the presence of an informative evolutionary signal in the sequence data (Fig. S1).

Stomatopoda diverged from its closest crustacean relatives 340 Ma (95% CRI [401–313 Ma]; Fig. 1). The crown-group unipeltatan stomatopods arose 193 Ma (95% CRI [267–131 Ma]), when the lineage leading to Hemisquilla split from remaining stomatopods. Squilloidea was found to be a relatively young superfamily, with the first divergence within the group occurring 70 Ma (95% CRI [113–40 Ma]). Stomatopods with specialized spearing claws arose about 155 Ma (95% CRI [212–111 Ma]), with specialized smashers appearing 126 Ma (95% CRI [174–87 Ma]).

Ancestral state reconstructions

Character state reconstructions indicate that the ancestral adult stomatopod eye contained six midband rows and hexagonal ommatidia (Fig. 2). Species in two superfamilies (Gonodactyloidea and Lysiosquilloidea) in this study possess six midbands, whereas the remaining superfamilies have lost complexity in this character (Fig. 2). The parasquilloid, Pseudosquillopsis, has lost three midbands, while Faughnia has lost four, leaving two midbands. Similarly, species in Squilloidea and some in Eurysquilloidea have lost four midbands, whereas bathysquilloids have lost all six midbands. Gonodactyloidea is the only superfamily with rectangular midband ommatidia; all other superfamilies have hexagonal midband ommatidia. The ancestral state reconstructions suggest that rectangular midband ommatidia evolved independently in Hemisquillidae and in the remaining Gonodactyloidea.

Phylogenetic relationships

This study presents the most comprehensive molecular phylogenetic analysis of Stomatopoda so far, based on representatives from six superfamilies. Our estimates of stomatopod relationships are congruent with the findings of Porter et al. (2010), who inferred a weakly supported sister relationship between Gonodactyloidea and Squilloidea using four markers (16S, 18S, CO1, and 28S). In contrast, Ahyong & Jarman (2009) found Lysiosquilloidea and Squilloidea to be sister clades, based on an analysis of three markers (12S, 16S, and 28S). Our analysis is the first to include representatives from the superfamilies Eurysquilloidea and Bathysquilloidea. The eurysquilloid representative is the sister lineage to a clade comprising Squilloidea and Parasquilloidea, which is consistent with its superfamily status. Additionally, the close relationship between Eurysquilloidea, Parasquilloidea, and Squilloidea is consistent with evidence from their morphology (Ahyong & Harling, 2000; Ahyong, 2001; Ahyong, 2005). However, our analysis placed the single bathysquilloid as the sister taxon to Pseudosquillidae, a position that is not supported by morphological evidence.

The placement of Hemisquilla as the sister lineage to all other stomatopods corroborates the results of previous molecular studies by Ahyong & Jarman (2009) and Porter et al. (2010). Collectively, these results support the hypothesis that hemisquillids are not a true gonodactyloid lineage. This suggests that Hemisquillidae should be moved to a separate superfamily (Porter et al., 2010), although such an action would be premature until rare morphologically intermediate taxa (Alainosquillidae) can be analysed. The phylogenetic placement of Hemisquilla is consistent with its seemingly plesiomorphic morphology, such as its subcylindrical body, and raptorial appendage that is neither highly optimized for powerful smashing nor efficient spearing, closely resembling that of the stem-lineage pseudosculdids (Ahyong & Harling, 2000; Ahyong, Garassino & Gironi, 2007; Ahyong & Jarman, 2009). Additionally, the preference of Hemisquilla for soft substrates is more similar to that of squilloids and other spearers than to gonodactyloids, which usually occupy hard substrate cavities (Schram et al., 2013).

The composition of the remaining Gonodactyloidea (apart from the anomalous inclusion of Bathysquilla, which represents a distinct superfamily) is congruent with previous molecular and morphological evidence, corroborating the removal of Eurysquillidae and Parasquillidae to separate superfamilies by Ahyong & Harling (2000). Additionally, the inferred placement of Pseudosquillidae within Gonodactyloidea in our study accords with morphological evidence (Ahyong & Harling, 2000) and differs from the near-basal position among stomatopods (albeit weakly supported) inferred by Ahyong & Jarman (2009) and Porter et al. (2010). However, the Bayesian topology test suggests a non-monophyletic Gonodactyloidea (excluding Hemisquilla), indicating that further exclusions might be necessary. The inferred placement of Odontodactylidae within Gonodactyloidea is consistent with those of both Porter et al. (2010) and Barber & Erdmann (2000), who estimated Odontodactylidae to be the sister group to all other smashing gonodactyloids. Such a result is consistent with morphological evidence and a monophyletic smashing clade (Ahyong & Harling, 2000).

Evolutionary timescale of stomatopods

Our molecular-clock analysis has provided an estimate of when the major stomatopod superfamilies diverged, shedding light on the geological events that might have been associated with the diversification of the order. The results of our analyses support fossil evidence for a Carboniferous divergence between Stomatopoda and other malacostracans (Ahyong, 1997; Hof, 1998; Ahyong & Jarman, 2009).

The estimates of divergence times are consistent with the notion that the break-up of the supercontinent Pangaea was associated with a key period in the evolution of Stomatopoda, which has been suggested previously for Hemisquilla (Schram, 1977; Schram et al., 2013). Diversification of extant stomatopod lineages can be traced to the beginnings of the break-up of Pangaea (i.e., ∼175 Ma; Fig. 1), with all of the divergences among superfamilies occurring from 193 to 95 Ma. Diversification might have coincided with the creation of new coastal habitats, which became available as landmasses were separated to form the Tethys Sea. As mentioned above, stomatopods predominantly occur in tropical coastal waters; therefore, the appearance of new habitats might have allowed the group to expand substantially.

The diversification of stomatopods with the opening of new habitats is consistent with findings from other crustaceans, especially decapods (George, 2006; Palero et al., 2009). George (2006) proposed that spiny lobsters (Palinuridae) radiated to new habitats as continental plates shifted, fragmenting the Tethys Sea and dividing Southern Ocean habitats. The age estimates from our analysis, along with the largely cosmopolitan distribution of stomatopod families and superfamilies, indicate that these lineages probably originated in the Tethys Sea when sea levels were high during the Cretaceous (Veron, 1995; George, 2006). The major superfamilies diverged before the closing of the seaway during the Paleogene/Neogene and the subsequent separation of the Atlantic and Pacific Oceans (Ekman, 1953; Reaka, Rodgers & Kudla, 2008; Schram et al., 2013).

Both Bathysquilloidea and Squilloidea have previously been suggested as “basal” lineages in the Unipeltata. These hypotheses were based on superficial morphological similarities with the extinct sculdids in the segmentation of the uropodal exopod in the bathysquilloid Indosquilla, and strong dorsal carination of squillids, said to somewhat resemble that of the Jurassic sculdids (Kemp, 1913; Manning, Kropp & Dominguez, 1990). Our analysis supports a “basal” position for neither Bathysquilloidea nor Squilloidea, instead suggesting that Hemisquillidae is likely to be the sister lineage to all other stomatopods, which is supported by morphological and other molecular evidence (Ahyong & Jarman, 2009; Schram et al., 2013). However, because only a single hemisquillid was included in the analysis, the results require further corroboration using additional exemplars.

Evolution of raptorial appendages and vision

Two hypotheses regarding the origin of the raptorial appendage have been proposed. The first suggests that specialized smashing evolved from within the spearers after a long history of specialized spearing (Caldwell, 1991; Ahyong, 1997). The second proposes that specialized smashing and specialized spearing evolved early in the history of Unipeltata, with more-or-less monophyletic spearing and smashing clades, perhaps diverging from a Hemisquilla-like ancestor (Ahyong & Harling, 2000; Ahyong & Jarman, 2009). Our results reprise aspects of both hypotheses, being similar to the first in finding smashers deeply nested among the spearers. However, our analysis infers a short branch on the ancestral node leading to the smashers. This indicates that smashing evolved relatively quickly and might also account for the weak molecular support for the grouping of odontodactylids with other smashers, at least among the markers used in this study. Porter et al. (2010) also recovered a similarly short branch leading to the smashers. Moreover, a link has been suggested between the evolution of the smashing appendage and the diversification of coral reefs in the Cretaceous; this could explain the rapid diversification of Gonodactyloidea during this time (Ahyong & Harling, 2000). Owing to the weak support for a monophyletic Gonodactyloidea, further molecular research into the boundaries of the superfamily is clearly required.

Ancestral state reconstructions of the complex stomatopod visual systems suggest that the eyes of stem-lineage unipeltatans had six midbands as also found by Porter et al. (2010), but with hexagonal rather than rectangular ommatidia. Groups, such as squilloids, parasquilloids, bathysquilloids and some eurysquilloids, with fewer or no midbands, are deeply nested within the crown-group. This suggests that the reduction in morphological complexity occurred independently in different lineages. This reduction in complexity might be linked to the visual environment, which is strongly influenced by depth and turbidity (Schram et al., 2013).

Further support for the influence of environment can be seen in the eyes of bathysquilloids, parasquilloids, and squilloids, which have varying degrees of ommatidial reduction. These groups inhabit deep or turbid waters and reduced eyes and loss of visual complexity are typical adaptations to these environments (MacDonald, 1975; Gaten, 1998; Lins et al., 2012). Additionally, squilloids and possibly parasquilloids are thought to have monochromatic vision, having lost some or all bands for colour (parasquilloids variously have two or three ommatidial midbands; Cronin, 1985; Ahyong, 2005); further investigation is needed to confirm this. Monochromatic vision is advantageous in the turbid, dark waters that parasquilloids and squilloids usually inhabit, where increased sensitivity to movement of prey and predators is more important than colour resolution. Bathysquilloids, living at great depth, have not only lost all ommatidial midbands, but the ommatidia on the cornea surface are themselves highly reduced.