Molecular evidence on evolutionary switching from particle-feeding to sophisticated carnivory in the calanoid copepod family Heterorhabdidae: drastic and rapid changes in functions of homologues

ABSTRACT The oceanic calanoid copepod family Heterorhabdidae is unique in that it comprises both particle-feeding and carnivorous genera with some intermediate taxa. Both morphological and molecular (nuclear 18S and 28S rRNA) phylogenetic analyses of the family suggest that an evolutionary switch in feeding strategy, from transitions of typical particle-feeding through some intermediate modes to sophisticated carnivory, might have occurred with the development of a specially designed ‘poison’ injection system. In view of the small amount of genetic differentiation among genera in this family, the switching of feeding modes and re-colonization from deep to shallow waters might have occurred over a short geological period since the Middle to Late Miocene.

The aim of the present study was to test the phylogenetic hypotheses of Ohtsuka et al. (1997) regarding the evolutionary switch from particle-feeding to carnivory in the family Heterorhabdidae using molecular-based phylogenetic analyses employing the nuclear 18S rRNA (18S) and 28S rRNA (28S) genes. These genes have proven to be informative in constructing phylogenetic trees of the higher-level taxa of copepods (cf. Huys et al. 2006Huys et al. , 2007Huys et al. , 2012Blanco-Bercial et al. 2011a, 2011b. Furthermore, we estimated the initial divergence timing from particle-feeding to carnivory based on analysis of mitochondrial cytochrome c oxidase subunit I (COI) sequences, which are often used as a molecular clock (Thum and Harrison 2009;Milligan et al. 2011;Blanco-Bercial et al. 2011b). This gene is frequently used to detect differences among copepod populations or speciation (cf. Bucklin et al. 1999;Makino and Tanabe 2009). Therefore, we used these three genes to infer evolutionary trends in feeding and colonization of the family Heterorhabidae.

Sampling
Specimens of heterorhabdid species were collected off the Nansei Islands, southwestern Japan, by oblique towing of an Ocean Research Institute (ORI) net (diameter 1.6 m; mesh size 0.33 mm) between 0 and 1056 m depth during cruises of the training and research vessel (TRV) Toyoshio-maru, Hiroshima University. The specimens were fixed and preserved with 99.5% ethanol immediately after collection (Table 2).   (Folmer et al. 1994) were used. Polymerase chain reaction was carried out with a KOD FX system (Toyobo, Osaka, Japan) according to the manufacturer's standard protocol. Cycling conditions followed a three-step cycle, as follows: initial denaturation for 2 min at 94°C, followed by 30 cycles of 10 s each at 98°C, 30 s at 50°C, and 1 min at 68°C; and a final extension for 5 min at 68°C. The polymerase chain reaction products were purified by the silica method (Boom et al. 1990). Fragments were sent to FASMAC (Kanagawa, Japan) for sequencing using a standard Applied Biosystems protocol (Applied Biosystems, Foster City, CA). The sequences were aligned with MAFFT v. 7.127b (Katoh et al. 2002(Katoh et al. , 2005 using the E-INS-i algorithm. The aligned 18S and 28S sequences were concatenated. Ambiguous regions of the aligned sequences were detected using Gblocks 0.91b (Castresana 2000) and were removed manually.
The phylogenetic trees based on the 18S and 28S sequences were reconstructed using two independent strategies. A maximum likelihood-based reconstruction was performed using raxml GUI ver 1.3.1 (Stamatakis 2006) under the GTRGAMMA option. Bootstrap values (Felsenstein 1985) were calculated with 1000 replications. In addition, a Bayesian inference analysis was carried out using MrBayes v3. 2. 2 (Ronquist and Huelsenbeck 2003). For the Bayesian inference method, the best-fit evolutionary model for the present data was estimated using jModelTest 2 (Darriba et al. 2012). According to the Bayesian information criterion, K80 (Kimura 1980) and K80 + G (Kimura 1980) were chosen for 18S and 28S, respectively. Four Markov chains were run for 1,000,0000 generations, and were sampled every 100 generations. The first 25,000 samples from each run were discarded as burn-ins. Pairwise genetic distances between species, derived from the concatenated sequences of the 18S and 28S genes, were estimated using MEGA v. 5.2 (Tamura et al. 2011). The genetic distances of the concatenated 18S and 28S sequences of the family Metridinidae G.O. Sars, 1902 in the superfamily Arietelloidea G.O. Sars, 1902 (including Heterorhabdidae) were also calculated for comparison.
The divergence time between Disseta scopularis (Brady 1883) and Mesorhabdus brevicaudatus (Wolfenden 1905) was estimated using BEAST v1.8.2 (Drummond et al. 2012) from the mitochondrial COI sequences of eight species (Table 3). Although there is no molecular clock determined for copepods, there are several reports of mutation rates for crustaceans (Thum and Harrison 2009;Milligan et al. 2011;Blanco-Bercial et al. 2011b). In accordance with a previous study (Knowlton and Weigt 1998) we applied the mutation rate of mitochondrial COI for alpheid shrimp (1.4% per million years) for the divergence time estimation. Analyses were run for this species pair using an uncorrelated log-normal molecular clock (Drummond et al. 2006) under the Hasegawa, Kishino, and Yano model, which was the best-fit model chosen by jModelTest 2 (Darriba et al. 2012), and the unweighted pair-group method using an arithmetic averages (UPGMA) starting tree. Tracer v1.6 (Rambaut et al. 2014) was used for determination of the divergence times between taxon pairs and their 95% confidence interval values.
Previously registered sequences of the 18S, 28S and COI genes of Heterorhabdidae and Metridinidae were downloaded from GenBank and are listed in Table 3. The 18S and 28S sequences of Metridia effusa (Accession Nos HM997054 and HM997024) and Pleuromamma abdominalis (Accession Nos AB625963 and AB753524) were used as outgroups, because the sequence lengths are sufficient for use in our phylogenetic analysis and they belong to the same superfamily as Heterorhabdidae.

Results
The sequenced lengths of the 18S and 28S genes ranged from 908 to 1743 base pairs (bp) and from 606 to 786 bp, respectively. The concatenated data set of 18S and 28S sequences was filtered to a final length of 1416 bp (906 bp for 18S and 510 bp for 28S). A total of 148 nucleotide sites were variable, 74 of which were parsimony-informative.
For the 18S and 28S sequences, the pairwise genetic distances within the Heterorhabdidae ranged from 0.07% to 2.27% (Table 4). In particular, the genetic distance values were low (0.14-0.92%) within the four genera Hemirhabdus, Neorhabdus, Paraheterorhabdus and Heterorhabdus. Pairwise genetic distance values for members of the Metridinidae, belonging to the superfamily Arietelloidea, as well as among members of the Heterorhabdidae ranged from 0.21% to 2.21% (Table 4). The intergeneric distances within Heterorhabdidae were nearly the same as the interspecies distance values observed in the Metridinidae.
The maximum likelihood and Bayesian inference trees showed basically the same topology ( Figure 2). The genera Disseta, Mesorhabdus and Heterostylites formed the first three offshoots with high and moderate bootstrap support, respectively. The remaining four genera split into two clades: Hemirhabdus + Neorhabdus and Heterorhabdus + Paraheterorhabdus. The sequenced lengths of the COI genes were 562-655 bp, and the multiple-aligned 557-bp sequences, not including indels, were used for the divergence time estimation. BEAST analyses showed that the divergence between Disseta scopularis and Mesorhabdus brevicaudatus probably occurred 11.8 million years ago (median value; 95% confidence interval 5.5-18.6 million years ago), corresponding to the middle to late Miocene.

Discussion
Genetic analyses based on molecular markers are highly informative for inferring phylogenetic relationships among higher copepod taxa (e.g. Huys et al. 2006Huys et al. , 2007Huys et al. , 2012Blanco-Bercial 2011a, 2011b. These molecules have also shed light on the phylogeny and evolution of the family Heterorhabdidae in the present study to test the phylogenetic hypotheses proposed by Ohtsuka et al. (1997). The phylogeny of the family Heterorhabdidae has so far been analysed using only morphological characters (Ohtsuka et al. 1997;Park 2001). The present molecular-based phylogeny corresponds with morphology-based studies to suggest that the earliest and most recently derived groups are the primitive particle-feeder Disseta and the highly modified carnivores Heterorhabdus + Paraheterorhabdus, respectively, with the genus Mesorhabdus in an intermediate position. However, the two previous morphology-based studies had difficulty resolving the position of Heterostylites in the phylogeny (Ohtsuka et al. 1997;Park 2001) (Figure 1). The cladistic analysis of Ohtsuka et al. (1997) produced two phylogenetic hypotheses: (1) Heterostylites is positioned between Mesorhabdus and the more derived carnivorous genera ( Figure 1A), or (2) Heterostylites forms a clade with Hemirhabdus + Neorhabdus ( Figure 1B). Park (2001) also produced two potential cladograms on more intuitive grounds: the first showing topology (2) of Ohtsuka et al. (1997) (Figure 1B), and the other with a different topology in which Heterostylites forms a clade with Paraheterorhabdus + Heterorhabdus ( Figure 1C). However, Park's (2001) analysis was not based on the principle of parsimony. Moreover, in the present study, Heterostylites did not emerge as a sister to the Paraheterorhabdus + Heterorhabdus clade in the maximum likelihood tree (Figure 2). Therefore, a monophyly of (Heterostylites, Paraheterorhabdus, Heterorhabdus) does not have strong support. Our molecular phylogenetic analyses, based on ML methods, showed the same topology as tree (1) of Ohtsuka et al. (1997; Figure 1A); and the divergence between Heterostylites and the other four genera was robustly supported with high bootstrap values (98% in the ML tree; Figure 2). During the evolutionary switch in feeding strategy of the Heterorhabdidae, the following morphological transformations are inferred to have occurred (Nishida and Ohtsuka 1996;Ohtsuka et al. 1997): (1) modifications of segmental and elemental homologues, including tubularization of the mandibular ventral-most tooth, reductions and chitinization of elements on the maxilla, and elongation of the maxillary praecoxa; (2) transformations of a labral gland system to a poison-injection system; (3) body miniaturization. Considering the molecular-based phylogenetic trees of the present study (Figure 2), these three transformations appear to have occurred independently. The maxillary segments and elements seem to have gradually transformed from those in the primitive particle-feeder Disseta, through the intermediate forms and primitive carnivores such as Neorhabdus and Hemirhabdus, and finally to the most developed carnivore Heterorhabdus. It seems that the innovation of a poison-injection system in the carnivorous genera (i.e. conversion of the labral gland system and mandibular ventralmost teeth to inject poison into prey), occurred simultaneously with the divergence of these genera (the node marked with a star in Figure 2). Finally, body size seems to have decreased since the divergence of Heterorhabdus. This evolutionary event is hypothesized to have been caused by its re-colonization into the highly productive epipelagic zone, where certain anti-predation strategies such as body transparency and miniaturization are required for survival (Hamner 1995). Indeed, habitat shifts occurring in concert with an increase in productivity have been reported for the fossil-rich ostracods (Bradford-Grieve 2002), and Roe (1972) clearly showed that the body sizes of calanoid copepods are relatively smaller in the epipelagic zone than in deeper waters. These observations support our hypothesis that the originally deep-sea genus Heterorhabdus has successfully re-invaded the epipelagic zone through anti-predatory miniaturization.
Unfortunately, it is difficult to infer the geological timings of such evolutionary events because of the paucity of fossil copepods. Bradford-Grieve (2002) first considered the origin and colonization of calanoid copepods over geological time scales, and proposed that their origin might date back to the Palaeozoic, which is supported by a recent discovery of the Cambrian fossils of copepods (Harvey and Pedder 2013). According to Bradford-Grieve (2002), the primitive superfamily Arietelloidea, containing the family Heterorhabdidae, might have first colonized the deep sea in the Palaeozoic. However, the effect of the end-Permian anoxia in the oceans, corresponding to the most drastic mass extinction in the history of the Earth (Wignall and Twitchett 1996;Bottjer et al. 2008;Payne and Clapham 2012), was not considered in the above paper.
Hence, the evolutionary habitat shift in copepods is presumably related to spatial and temporal changes in productivity (Bradford-Grieve 2002). The present results can also provide insight into the relationships between such evolutionary and geological events based on the molecular phylogeny of the genera of the Heterorhabdidae. The fact that the 18S and 28S pairwise genetic distances between the heterorhabdid genera were similar to those between the metridinid species suggests that the timing of generic divergences in the former group most likely corresponds to the timing of speciation in the latter. Although the molecular clock is a powerful tool for the reconstruction of evolutionary history, it is difficult to apply this method to the analysis of poorly fossilized organisms. We estimated that divergence between Mesorhabdus brevicaudatus and Disseta scopularis occurred in the middle to late Miocene. Colonization of new habitats by planktonic copepods could have become possible through a combination of physiochemical and biological innovations, as Bradford-Grieve (2002) suggested. In the Cenozoic era, deep waters were well-ventilated, and so productivity became enormously high, and the abundance of zooplankters such as copepods and euphausiids has been suggested to have increased in concert with the phytoplankton blooms (Salamy and Zachos 1999). We suggest that colonization of the Heterorhabdidae in the deep sea might have also occurred during this era. This hypothesis is supported by considering the scenario of the colonization of oceanic waters by the terrestrial ancestors of the Cetacea (Fordyce 1980;Fordyce and Barnes 1994;Uhen 2010). Accelerated cetacean evolution is attributed to high primary productivity following deep-ocean ventilation since the Eocene-Oligocene boundary (Kennett and Stott 1990;Salamy and Zachos 1999;Zachos et al. 2001). Moreover, large body-sized baleen whales first appeared in the early Oligocene and then became the most taxonomically diverse group by the late Miocene (Uhen 2010), corresponding to an increase in planktivorous fish (Eiting and Smith 2007). Recently derived baleen whales and fish feed on large amounts of zooplankton; therefore, predation pressure on zooplankters must have been high in the late Miocene. Therefore, it can be speculated that the colonization of planktonic copepods into the deep sea might have been driven by a combination of ventilation, a high to moderate food supply in deep waters, and high predation pressure in shallow waters (Bradford-Grieve 2002). The ancestor of the Heterorhabdidae would have colonized the deep sea so as to escape from predators such as baleen whales and fish, and then developed carnivory with a specially designed 'poison' injection system because of the low productivity in the deep sea. This evolutionary switch of feeding modes might have occurred in the Miocene epoch. According to Fleminger (1986), the speciation rates of calanoid copepods may fall between 10,000 and 1,000,000 years. If this figure is applicable to the evolution of the Heterorhabdidae, re-colonization from the deep to shallow waters might have occurred soon after the Miocene.