2021 Volume 96 Issue 1 Pages 13-24
The class Branchiopoda (Crustacea) shows great diversity in morphology and lifestyle among its constituent higher-level taxa: Anostraca, Notostraca, Laevicaudata, Spinicaudata, Cyclestherida and Cladocera. The phylogenetic relationships among these taxa have long been controversial. We sequenced three orthologous nuclear genes that encode the catalytic subunit of DNA polymerase delta and the largest and second-largest subunits of RNA polymerase II in the expectation that the amino acid sequences encoded by these genes might be effective in clarifying branchiopod phylogeny and estimating the times of divergence of the major branchiopodan taxa. The results of phylogenetic analyses based on these amino acid sequences support the monophyly of Branchiopoda and provide strong molecular evidence in support of the following phylogenetic relationships: (Anostraca, (Notostraca, (Laevicaudata, (Spinicaudata, (Cyclestherida, Cladocera))))). Within Cladocera, comparison of the nucleotide sequences of these same genes shows Ctenopoda to be the sister group of Haplopoda + Anomopoda. Three statistical tests based on the present amino acid sequence data—the approximately unbiased test, Kishino–Hasegawa test and weighted Shimodaira–Hasegawa test—tend to refute most of the previous molecular phylogenetic studies on Branchiopoda, which have placed Notostraca differently than here; however, our results corroborate those of one recent phylogenomic study, thus confirming the effectiveness of these three genes to investigate relationships among branchiopod higher taxa. Divergence time estimates calibrated on the basis of fossil evidence suggest that the first divergence of extant branchiopods occurred about 534 Ma during the early Cambrian period and that diversification within the extant branchiopod lineages started in or after the late Permian.
The class Branchiopoda (Crustacea) is a relatively small, primarily freshwater taxon comprising about 1,120 described species (Adamowicz and Purvis, 2005; Brendonck et al., 2008; Forró et al., 2008; Ahyong et al., 2011). Branchiopods display a wide diversity of morphology and lifestyle (Olesen, 2009) and comprise six main groups: Anostraca (fairy shrimps); Notostraca (tadpole shrimps); Laevicaudata, Spinicaudata and Cyclestherida (clam shrimps); and Cladocera (water fleas). Laevicaudata, Spinicaudata and Cyclestherida were previously classified together as ‘Conchostraca’, but that grouping is now recognized as paraphyletic (Olesen and Richter, 2013). All non-cladocerans (fairy shrimps, tadpole shrimps and clam shrimps) live in inland freshwaters or salt lakes and are generally known as “large branchiopods”. Although they share characteristics such as many serially similar phyllopodous trunk limbs (Olesen, 2009), morphological and molecular analyses suggest that “large branchiopods” are paraphyletic (see Olesen and Richter, 2013). About half of the known species of Branchiopoda are cladocerans, some of which are marine. The Cladocera are further divided into four subgroups, Anomopoda, Ctenopoda, Haplopoda and Onychopoda (Calman, 1909; Fryer, 1987; Walossek, 1993; Olesen et al., 1997; Olesen, 1998, 2004; Negrea et al., 1999).
The systematics of Branchiopoda has recently received much attention, not only because of interest in the intrinsic diversity within the group, but also on account of branchiopods’ possible importance for understanding the origin and evolution of insects (Hexapoda). Molecular systematics has greatly changed our traditional understanding of arthropod phylogeny by revealing that Crustacea and Hexapoda form a common clade (Pancrustacea), which is now widely accepted among zoologists; for a review, see Giribet and Edgecombe (2019). Among crustaceans, either Remipedia or Remipedia + Cephalocarida is currently considered by many to be the clade most closely related to Hexapoda (Regier et al., 2010; von Reumont et al., 2012; Oakley et al., 2013; Schwentner et al., 2018; Lozano-Fernandez et al., 2019; Noah et al., 2020), but Branchiopoda has also been considered to be the sister group of Hexapoda (Glenner et al., 2006; Lozano-Fernandez et al., 2019). For a review of this topic, see Giribet and Edgecombe (2012).
In recent years, many investigations of branchiopod phylogeny and evolution have been conducted, based on analyses of both morphological and molecular data (Hanner and Fugate, 1997; Olesen, 1998, 2000; Negrea et al., 1999; Taylor et al., 1999; Spears and Abele, 2000; Braband et al., 2002; deWaard et al., 2006; Stenderup et al., 2006; Richter et al., 2007; Olesen, 2009; Schwentner et al., 2018; Luchetti et al., 2019). Relationships among the higher-level taxa, especially the phylogenetic positions of Notostraca and Cyclestherida and the phylogeny within Cladocera, have long been disputed (Stenderup et al., 2006; Richter et al., 2007; Olesen, 2009) (see Fig. 1). Nonetheless, morphological analyses have succeeded in providing a model of branchiopod sister-group relations and phylogeny, namely (Anostraca, (Notostraca, (Laevicaudata, (Spinicaudata, (Cyclestherida, Cladocera))))) (Richter et al., 2007; Olesen, 2009) (Fig. 1A).
Phylogenetic hypotheses of Branchiopoda based on morphological and molecular analyses. A, mainly based on morphological studies (Richter et al., 2007; Olesen, 2009; see also Olesen and Richter, 2013) and supported by a recent phylogenomic study (Schwentner et al., 2018); B, C and D, based on various molecular markers with different analytical methods as described in the Introduction (for B: deWaard et al., 2006; Stenderup et al., 2006; Richter et al., 2007; for C: Stenderup et al., 2006; for D: Braband et al., 2002).
As for molecular evidence, mitochondrial 12S rRNA data supported Notostraca as the sister group to Laevicaudata (Braband et al., 2002) (Fig. 1D), whereas maximum likelihood (ML) analysis and Bayesian inference based on mitochondrial 16S rRNA and nuclear 28S rRNA data supported Notostraca as the sister group to Cyclestherida + Cladocera (Stenderup et al., 2006) (Fig. 1C) and maximum parsimony analysis supported Notostraca as the sister group to Spinicaudata + Cyclestherida + Cladocera (Stenderup et al., 2006) (Fig. 1B). This last phylogeny (Fig. 1B) was also supported by an analysis based on six molecular loci (COI, 16S rRNA, 18S rRNA, EF-1 alpha, 12S rRNA and 28S rRNA) (deWaard et al., 2006), but when the same authors restricted their analysis to only three of these loci (COI, 16S rRNA and 18S rRNA), four of the large branchiopod taxa, Anostraca, Notostraca, Laevicaudata and Spinicaudata, were clustered into one clade. In the analyses of Richter et al. (2007), who used sequence data from six loci and 65 morphological characters, the morphological data supported a sister-group relationship between Notostraca and Diplostraca (= Laevicaudata + Spinicaudata + Cyclestherida + Cladocera) (Fig. 1A), while the molecular data strongly supported Notostraca as the sister group to Spinicaudata + Cyclestherida + Cladocera (Fig. 1B). Most recently, mitochondrial genomic sequence data provided a topology in which Notostraca, Spinicaudata and Cladocera formed one group, but with weak bootstrap support and no data from Laevicaudata and Cyclestherida (Luchetti et al., 2019). Most analyses have supported a sister-group relationship between Cyclestherida and Cladocera, but some studies have suggested a position of the former within the latter (deWaard et al., 2006; Richter et al., 2007).
Phylogenomic studies have increased in recent years because it has become easier to obtain genome-wide sequence data. In one such study, Schwentner et al. (2018) provided a strongly supported topology that, for branchiopods, basically coincides with the morphology-based phylogenies of Richter et al. (2007) and Olesen (2009). However, phylogenomic datasets produced by different workers sometimes give contradictory results despite showing strong support values. For example, one phylogenomic dataset shows the Cephalocarida (a class of Crustacea) as sister to Remipedia (Regier et al., 2010; Noah et al., 2020), while another shows it as sister to Branchiopoda + Remipedia + Hexapoda (Schwentner et al., 2018). These differences are probably due to the inclusion of inappropriate genes in the datasets. It is important to select molecular markers that are known to be appropriate for phylogenetic analyses. Previous studies (Ishiwata et al., 2011; Sasaki et al., 2013; Miyazawa et al., 2014) have shown that the amino acid (aa) sequences encoded by three nuclear genes, RNA polymerase II largest subunit (RPB1), RNA polymerase II second-largest subunit (RPB2) and DNA polymerase delta catalytic subunit (DPD1), are effective for resolving the relationships among the higher-level taxa of Hexapoda and Myriapoda. Here, we used the sequences of these three genes (RPB1, RPB2, DPD1) obtained from 15 species of branchiopod representing all six higher taxa to address the questions raised above. Phylogenetic analyses were conducted based on the inferred aa sequences of the encoded proteins, and previously proposed phylogenetic hypotheses were evaluated statistically. In addition, the divergence times of branchiopod lineages were estimated based on fossil calibration points including a recently discovered early spinicaudatan fossil (Gueriau et al., 2016).
The nucleotide sequences reported here have been deposited in the DDBJ/EMBL/GenBank nucleotide sequence databases with accession numbers LC532181–LC532216 as shown in Supplementary Table S1.
Supplementary Table S1 lists all the taxa used in this study along with their collection localities and the respective sequence accession numbers for all three nuclear protein-coding genes (DPD1, RPB1, RPB2). Fifteen branchiopod species (Fig. 2) representing all major groups were sequenced for these genes, with the exception of Onychopoda (one of four infraorders within the suborder Cladocera), living material of which was not available. Twenty-two other arthropod species, including two chelicerates, one myriapod, 16 hexapods and three non-branchiopod crustaceans, were used as outgroups in the phylogenetic analysis (Supplementary Table S1). We would have liked to include representatives of the purportedly basal crustacean classes Cephalocarida and Remipedia in both sets of outgroups, but no living material of these rare animals was available and none of their relevant nuclear gene sequences has yet been archived.
Photographs of branchiopod species included in this study.
Live specimens were used for RNA isolation. Total RNA was extracted from a part of the specimen (legs and tails) for anostracans and notostracans and from the whole body for other branchiopod specimens using ISOGEN (Nippon Gene) following the manufacturer’s instructions. Total extracted RNA was dissolved in 22 μl DNase/RNase-Free Distilled Water (Invitrogen), and 1 μg was reverse-transcribed to cDNA using the SMART RACE cDNA Amplification Kit (Clontech) and SuperScript III Reverse Transcriptase (Invitrogen). The resulting cDNA was then used as a template to amplify each of the three genes by PCR with LA-Taq (Takara Bio) using sense and antisense degenerate primers (Supplementary Fig. S1 and Supplementary Table S2). The experimental procedures for PCR amplification and sequencing were according to Ishiwata et al. (2011) and Sasaki et al. (2013).
For Cyclestheria hislopi, total RNA (822 ng) was used to construct an RNA library using the TruSeq RNA Library Preparation Kit v2 (Illumina). The library was then used for next-generation sequencing (NGS) on MiSeq (Illumina). Sequence assembly was performed with CLC Genomics Workbench (Filgen). Most of the nucleotide sequences of the three genes used in this study were obtained from the NGS sequence data, and the missing parts were amplified by PCR as described above.
Sequence alignment and phylogenetic analysesNucleotide sequences were corrected and edited using 4Peaks (Netherlands Cancer Institute) and connected by ATSQ ver. 5.1 (GENETYX Corporation). The consensus sequences of the three genes were translated to the corresponding encoded aa sequences with GENETYX-MAC ver. 16.0.1 (GENETYX Corporation), aligned by MAFFT L-INS-i (Katoh et al., 2005) and then inspected by eye. Unambiguously aligned aa sites were selected using Gblocks ver. 0.91b (Castresana, 2000; Talavera and Castresana, 2007).
Selection of an appropriate outgroup often improves the results of phylogenetic analysis. Outgroup taxa that represent long branches may cause misplacement of long-branched ingroup taxa (Schneider and Cannarozzi, 2009), but the use of certain outgroups that are much more closely related to the ingroups may reduce the effects of such artifacts. With this in mind, we prepared two sequence datasets for phylogenetic analysis. One consisted of all 37 samples (37-sample dataset) used in this study and was intended to help confirm or refute the monophyly of Branchiopoda. The 22 outgroup taxa included representatives of the Chelicerata, Myriapoda, Hexapoda and two non-branchiopod crustacean classes; the hexapods included ten species of wingless insect (two proturans, three diplurans, two bristletails and three silverfishes) and six species of winged insect (Fig. 3, Supplementary Table S1). The other dataset comprised 22 samples (22-sample dataset) with only seven hexapod and non-branchiopod crustacean species used as outgroups. The 22-sample dataset was used to analyze the phylogenetic relationships among branchiopods, statistically test the previously proposed hypotheses of their relationships and estimate divergence times.
Phylogenetic tree of Branchiopoda based on aa sequences of DPD1, RPB1 and RPB2 genes with 22 outgroup species (37-sample dataset). Phylogenetic analyses were performed with RAxML and MrBayes. Bootstrap values and posterior probabilities are shown at each node.
Phylogenetic trees were inferred by the ML method and Bayesian analysis. The ML analysis was carried out on RAxML 7.2.8 (Stamatakis, 2006) with the best-fit model, namely the LG+G model that was determined for each dataset using ProtTest 3 under the AIC, BIC and AICc criteria (Darriba et al., 2011). To assess the node supports, 1,000 nonparametric bootstrap runs were performed for each ML analysis. Bayesian inference was performed with MCMC analysis in MrBayes v3.1.2 (Ronquist and Huelsenbeck, 2003) using the LG+G model. Each analysis was run for 2,000,000 generations, and the trees were sampled every 1,000 generations (burn-in = 500,000 generations). To test the convergence of chains, the log file of the MrBayes analyses was examined by calculating the effective sample sizes of all parameters using Tracer v1.5 (http://tree.bio.ed.ac.uk/software/tracer/). Bayesian posterior probabilities were obtained from the majority rule consensus tree sampled after the initial burn-in period.
In addition, phylogenetic analyses within Anostraca, Notostraca and Cladocera were conducted using the nucleotide sequences of the three genes. Two sequence datasets were constructed: one for Anostraca and Notostraca (11,225 nucleotide sites), and the other for Cladocera (11,202 nucleotide sites). For both datasets, sequence alignment was performed with MAFFT (Katoh et al., 2005) and all gaps were removed. Model testing was carried out using jModelTest 2.1.5 (Guindon and Gascuel, 2003; Darriba et al., 2012), and GTR+G+I was selected as the best model for the two datasets. The phylogenetic trees were inferred with RAxML 7.2.8 (Stamatakis, 2006) under this model, again with 1,000 bootstrap replicates.
Statistical testsPrevious molecular phylogenetic studies have suggested three possible placements of Notostraca (Braband et al., 2002; Stenderup et al., 2006; Richter et al., 2007). In this study, these hypotheses were statistically tested based on the 22-sample dataset using the CONSEL program (Shimodaira and Hasegawa, 2001) with the approximately unbiased test (AU), the Kishino–Hasegawa test (KH) and the weighted Shimodaira–Hasegawa test (wSH) (Kishino and Hasegawa, 1989; Shimodaira and Hasegawa, 2001; Shimodaira, 2002).
Divergence time estimationAn ephemeral-pool branchiopod community from the late Famennian (372.2–358.9 Ma) was recently discovered at the Strud locality in Belgium, and two fossil arthropods, Haltinnaias serrata Gueriau et al., 2016 and Gesvesia pernegrei Gueriau et al., 2016, were described as being the earliest known unequivocal representatives of the groups Anostraca and Spinicaudata, respectively (Gueriau et al., 2016). The minimum age of G. pernegrei (358.5 Ma) provides a calibration point for the time of divergence between Spinicaudata and Cladoceromorpha (= Cyclestherida + Cladocera). With this and four other fossil-based calibration points (Wolfe et al., 2016), a timetree was constructed on MEGA X (Kumar et al., 2018; Stecher et al., 2020) based on the 22-sample dataset using the RelTime method (Tamura et al., 2012, 2018) and the LG model (Le and Gascuel, 2008) with gamma distribution (+G) and evolutionarily invariable (+I) settings. The Tao et al. (2019) method was used to set minimum and maximum time boundaries on nodes for which calibration densities were provided. According to Wolfe et al. (2016), the minimum ages of the other four calibration points are 125.71 Ma for crown-group Anostraca, 162.5 Ma for crown-group Spinicaudata, 386.9 Ma for crown-group Diplostraca and 405 Ma for crown-group Hexapoda. For all calibration points, the maximum age was set to 521 Ma (Wolfe et al., 2016).
We determined the sequences of three nuclear protein-coding genes (DPD1, RPB1 and RPB2) from three species of Anostraca, three species of Notostraca, one species of Laevicaudata, three species of Spinicaudata, one species of Cyclestherida and four species of Cladocera. For DPD1, the complete encoded aa sequence was obtained from representatives of four higher taxa and ranged in length from 1,047 to 1,120 aa, although that obtained from the species of Notostraca and Laevicaudata lacked a short part of the N-terminal region (Supplementary Table S3). For RPB1, the C-terminal region (approx. 400 aa) had repeated sequences that were not appropriate for phylogenetic analysis, and this region was excluded from the analysis, while the N-terminal sequence was determined completely for ten species representing all six taxa and had a length ranging from 1,529 to 1,907 aa (Supplementary Table S3). For RPB2, the complete sequences obtained from nine species ranged in length from 1,174 to 1,179 aa. The remaining six species were each missing a short part of the N-terminal region and consequently had a sequence length ranging from 1,146 to 1,164 aa (Supplementary Table S3).
The 37-sample dataset including numerous outgroup taxa (22 species) yielded a total of 3,447 aa sites (DPD1, 874 aa; RPB1, 1,434 aa; RPB2, 1,139 aa) selected by Gblocks (Castresana, 2000; Talavera and Castresana, 2007) as suitable for phylogenetic analyses. The 22-sample dataset, including only seven outgroup species, had a total of 3,490 aa sites (DPD1, 886 aa; RPB1, 1,467 aa; RPB2, 1,137 aa). Neither of the two datasets contains missing aa sites. Sequence information including parsimony-informative and variable aa sites in the two alignment datasets is shown in Supplementary Table S4. The proportion of parsimony-informative sites was highest for the DPD1 gene, and the RPB2 gene was the most highly conserved. These results were consistent with those observed for these genes in insects (Ishiwata et al., 2011; Sasaki et al., 2013) and myriapods (Miyazawa et al., 2014).
Monophyly of BranchiopodaUsing the 37-sample dataset including 22 outgroup species, all analytical methods generated the same tree topology, with the exception of the relationships among the three species of Triops (Fig. 3). Monophyly of Branchiopoda was confirmed with robust support (bootstrap value of 100% and posterior probability of 1.00). Such monophyly had been indicated previously by morphological analyses including larval characters (e.g., Olesen, 2007, 2009), and by some molecular studies based on different molecular data and sampling sets (Regier et al., 2010; von Reumont et al., 2012; Schwentner et al., 2018).
When chelicerates and myriapods were included as outgroups, Hexapoda was most closely related to Branchiopoda with high posterior probability but low bootstrap support (Fig. 3). The other two crustaceans, belonging to Copepoda and Malacostraca, showed a close relationship to each other with high support (Fig. 3), in agreement with the results of von Reumont et al. (2012) and Schwentner et al. (2018). However, to understand the phylogeny and evolution of Crustacea, it is necessary to analyze more crustacean samples including representatives of the two classes Cephalocarida and Remipedia, which were not available in this study. Despite extensive genomic data-based studies, it remains unclear which taxa are most closely related to Branchiopoda; Schwentner et al. (2018) showed that Branchiopoda is the sister group of Remipedia + Hexapoda, while Noah et al. (2020) indicated a sister relationship between Branchiopoda and Multicrustacea (Copepoda + Thecostraca + Malacostraca). In light of the phylogenetic relationships suggested by these studies, the outgroups used in the present study (Hexapoda, Copepoda and Malacostraca) can be judged as appropriate, if not ideal, for phylogenetic analyses of Branchiopoda.
Phylogenetic relationships among higher-level taxa of BranchiopodaBoth the 37-sample and 22-sample datasets resulted in the same strongly supported inferences of phylogenetic relationships among the sequenced branchiopods considered in this study, except for minor differences among the compared species of Anostraca and Notostraca (Figs. 3 and 4). The choice of outgroups did not appreciably affect the outcome. From the topology obtained, Anostraca was identified as the sister group to all remaining branchiopods, with the latter constituting a monophyletic Phyllopoda. Within Phyllopoda, Notostraca branched off first, followed by Laevicaudata and Spinicaudata in that order, and Cyclestherida was the sister lineage to Cladocera (Figs. 3 and 4). These relationships were all supported by very high bootstrap values (≥ 89%) and posterior probabilities (1.00) (Fig. 4). The present results thus strongly support the monophyly of several higher-level taxa, namely Phyllopoda, Diplostraca, Onychocaudata, Cladoceromorpha and Cladocera (Fig. 4), which had mainly been proposed or corroborated on the basis of morphological studies (Richter et al., 2007; Olesen, 2007, 2009; see also Olesen and Richter, 2013), and furthermore agree with the results of a recent phylogenomic study (Schwentner et al., 2018).
Phylogenetic tree of Branchiopoda inferred from aa sequences of DPD1, RPB1 and RPB2 genes, with seven outgroup species including four hexapods and three non-branchiopods (22-sample dataset). Bootstrap values and posterior probabilities are shown at each node.
In most molecular analyses conducted up to now, except for Schwentner et al. (2018), the detailed results have differed from ours on account of the use of different datasets or analytical methods (Braband et al., 2002; deWaard et al., 2006; Stenderup et al., 2006; Richter et al., 2007). In general, Notostraca has fared the worst in terms of inconsistent phylogenetic position, appearing as either 1) the sister group of Laevicaudata (Braband et al., 2002) (Tree No. 2 in Fig. 5); 2) the sister group of Spinicaudata + Cyclestherida + Cladocera (deWaard et al., 2006; Stenderup et al., 2006; Richter et al., 2007) (Tree No. 3 in Fig. 5); or 3) the sister group of Cyclestherida + Cladocera (Stenderup et al., 2006) (Tree No. 4 in Fig. 5). It would simplify matters if we could reject one or more of these options with statistical confidence. Statistical analyses using the AU, KH and wSH tests as outlined above, based on the 22-sample dataset, resulted in P values for all tests that were less than 0.01 (Fig. 5B). All three tree topologies were thus rejected by these tests at the 1% significance level. The present results provide strong molecular evidence to support Notostraca as the sister group of all extant non-anostracan branchiopods, in agreement with the prevailing morphology-based model of branchiopod phylogeny (Richter et al., 2007; Olesen, 2009) and also with a recent phylogenomic study (Schwentner et al., 2018). This confirms the effectiveness of using the aa sequences encoded by these three genes—as opposed to more expensive, time-consuming and technically difficult phylogenomic analyses—to investigate relationships among crustaceans, as well as among the higher taxa of Hexapoda (Ishiwata et al., 2011; Sasaki et al., 2013) and Myriapoda (Miyazawa et al., 2014).
Statistical testing of hypotheses concerning the phylogenetic position of Notostraca. The AU test, KH test and wSH test were performed based on the same aa sequence dataset as in Fig. 4. A, tested tree topologies: the ML tree obtained in this study (Tree No. 1) and the three previously proposed candidate topologies (Tree Nos. 2–4). The candidate topologies were inferred by ML analysis with the same method and dataset used in Fig. 4, with Notostraca constrained into a single cluster with Laevicaudata (Tree No. 2; see Braband et al., 2002), with Onychocaudata (Tree No. 3; see deWaard et al., 2006; Stenderup et al., 2006; Richter et al., 2007) or with Cladoceromorpha (Tree No. 4; see Stenderup et al., 2006). The nodes within Notostraca, Onychocaudata (in Tree No. 3) and Cladoceromorpha (in Tree No. 4) were not fixed. B, results of statistical tests for each tree topology with tree and specimen numbers corresponding to those shown in A. ΔlnL, the difference of the log-likelihood between ML tree and test tree.
As shown in Figs. 3 and 4, the phylogenetic relationships within Anostraca, Notostraca and Cladocera as revealed by the encoded aa sequences of the three analyzed nuclear genes were not well resolved due to low node-supporting bootstrap values. This is probably due to the lower-level taxa within each group being too closely related, and thus having accumulated too few sequence changes in these three genes to resolve their phylogenetic relationships. Additional analyses for these taxa using nucleotide sequences that included synonymous substitutions, and not just the encoded aa sequences, showed well-resolved phylogenetic relationships with strong bootstrap support (Fig. 6). The branching patterns of the three anostracan species and the three species of Triops were (Artemia franciscana, (Branchinella kugenumaensis, Eubranchipus uchidai)) and (Triops cancriformis, (T. longicaudatus, T. granarius)) (Fig. 6A), respectively. With Cyclestherida as the outgroup, the branching pattern of the four cladoceran species (Fig. 6B) showed Diaphanosoma sp. (Ctenopoda) as the sister group of the remaining species (Haplopoda + Anomopoda), and the species of Anomopoda clustered together. Despite many studies based on both morphological and molecular data, the relationships among the four higher taxa of Cladocera, namely the Anomopoda, Ctenopoda, Haplopoda and Onychopoda, have remained unclear (Braband et al., 2002; deWaard et al., 2006; Richter et al., 2007; Olesen, 2009; Olesen and Richter, 2013). Similar to our results, Schwentner et al. (2018) showed that Ctenopoda is the sister group of the remaining groups of Cladocera, without, however, including any representative of Haplopoda in their study. The present results provide partial resolution, but only four cladoceran species were included in the analysis, with no representative of the Onychopoda, and the aa sequence analysis gave different results than the nucleotide sequence analysis, albeit with quite weak support for the former. Elucidation of the phylogeny of Cladocera still requires more study.
ML trees of Notostraca and Anostraca (A) and Cladocera (B) inferred from nucleotide sequences of DPD1, RPB1 and RPB2 genes. Phylogenetic analyses were performed with RAxML (GTR+G+I model). Bootstrap values are shown at each node.
In the timetree produced by using five branchiopod fossil ages as calibration constraints (see Materials and Methods), the divergence between Anostraca and Phyllopoda is estimated to have occurred 534 Ma, during the early Cambrian, and the divergence times between Notostraca and Diplostraca (= Laevicaudata + Spinicaudata + Cyclestherida + Cladocera) and between Laevicaudata and Onychocaudata (= Spinicaudata + Cyclestherida + Cladocera) are estimated to be 496 Ma (late Cambrian) and 419 Ma (beginning of Devonian), respectively (Fig. 7). The split between Cyclestherida and Cladocera is estimated to have occurred 305 Ma, during the late Carboniferous, and the Cladocera first diverged 261 Ma, during the late Permian (Fig. 7).
Timetree of branchiopod higher taxa cladogenesis inferred using the RelTime method and the LG+G+I model. Evolutionary analyses were conducted in MEGA X (Kumar et al., 2018; Stecher et al., 2020). Filled circles indicate the calibration points (see MATERIALS AND METHODS). The estimated log likelihood value is -34133.33. The divergence times between taxa are shown at each node, and the bars represent 95% confidence intervals. The chronology bar is according to the International Chronostratigraphic Chart v2020/01 of the International Commission on Stratigraphy.
Luchetti et al. (2019) estimated the divergence time of four higher branchiopod taxa using a phylogenetic tree based on mitochondrial genomic sequences. The divergence date between Spinicaudata and Cladocera and the internal divergence ages of various anostracan and cladoceran taxa were similar to or older than ours, but the times of divergence between Anostraca and Phyllopoda, and between Notostraca and Onychocaudata, were younger than our estimates. This may be due to the saturation of base substitutions accumulated in the mitochondrial DNA sequences of those lineages. If so, including additional fossil calibration points for higher taxa crown groups such as Onychocaudata (Gueriau et al., 2016) and Diplostraca (Wolfe et al., 2016) might have reduced the effect of such saturation. Another study (Sun et al., 2016) estimated the divergence times of the major lineages of Branchiopoda using various molecular clock methods, namely BEAST (Drummond and Rambaut, 2007), Multidivtime (Thorne et al., 1998), MCMCTree (Yang, 2007) and r8s (Sanderson, 2003). For the divergence times of Anostraca–Phyllopoda and Notostraca–Diplostraca, BEAST calculated maximum ages with averages of 606 Ma and 548 Ma, and r8s gave minimum ages with averages of 495 Ma and 465 Ma, respectively. Our estimate yielded divergence times within these bounds, 534 Ma and 496 Ma (Fig. 7), which are similar to Sun et al.’s (2016) Multidivtime and MCMCTree estimates. These divergence ages also roughly support the model of the origin and evolutionary process of Branchiopoda proposed by Negrea et al. (1999) on the basis of morphological data and the then-known fossil evidence.
Phylogenetic analysis of Anostraca (Weekers et al., 2002) has shown a division of this order into two distinct monophyletic groups. The first consists of the Artemiidae and Parartemiidae, representing saltwater taxa, whereas the second contains all freshwater taxa. Artemia franciscana, one of the three anostracan species used in the present study, belongs to the first group and our other two species belong to the second group. Our results suggest that the two clades of extant anostracans diverged 142 Ma, during the early Cretaceous (Fig. 7). This, together with the lack of fossil evidence of any additional stem- or crown-group anostracan diversity between the Middle Devonian and the present day (Negrea et al., 1999; Belk and Schram, 2001), suggests that the anostracan lineage was evolutionarily static up until this binary split. A similar case may be seen in Cyclestherida, which our data show splitting from the common ancestor of Cladocera 305 Ma, during the late Carboniferous (Fig. 7). Although the extant populations seem to represent the results of a split into three geographical lineages during the Middle to Late Cretaceous (Schwentner et al., 2013), all are still similar enough to be formally classified as a single species, Cyclestheria hislopi. The first divergence among extant spinicaudatans is estimated to have occurred 170 Ma, in the middle Jurassic, much later than that inferred in a previous study (Schwentner et al., 2020), either 294.6 Ma or 403 Ma depending on the choice of fossil-based minimum-age calibration point. In sum, the present results (Fig. 7) suggest that the major branchiopodan lineages were established at least before the break-up of Pangaea, and diversification within the extant branchiopod lineages started in or after the late Permian, perhaps in connection with historical vicariance events such as the break-up of Pangaea into Gondwana and Laurasia, until which time they had evolved minimally, and silently. This notion should be tested by further analyses using additional species from the lineages of the Notostraca and Laevicaudata.
ConclusionThe aa sequences of three protein-coding genes (RPB1, RPB2 and DPD1) proved to be highly reliable in inferring the phylogenetic relationships of higher-level taxa within Branchiopoda. Our results strongly support the monophyly of Branchiopoda and some of its constituent taxa, Phyllopoda, Diplostraca, Onychocaudata, Cladoceromorpha and Cladocera, and corroborate Schwentner et al.’s (2018) placement of the Notostraca as the sister to all extant non-anostracan branchiopods. Furthermore, the nucleotide sequence analyses showed that the suborder Ctenopoda is the sister group to Haplopoda + Anomopoda within Cladocera. The divergence time estimated on the basis of firm data including that provided by a recently discovered unequivocal spinicaudatan fossil suggests that the first divergence of crown-group branchiopods occurred during the early Cambrian period about 534 Ma, and that diversification within the extant branchiopod lineages started in or after the late Permian, before which they had evolved but little.
We are grateful to M. Goto, K. Igarashi, Y. Kusuoka, K. Masunaga, T. Ojika, Y. Okamoto, R. J. Smith and T. Yamamoto for their invaluable contribution in the supply of samples. We deeply thank P. Dabseepai and her student for sampling Cyclestheria hislopi in Khon Kaen, Thailand. Thanks are also due to A. Sasaki for sampling and sequencing of Cyclestheria hislopi and producing illustrations of the representative branchiopods. This work was supported in part by grants from JSPS KAKENHI (No. 17570191 and No. 16K07471). Some collections were made as part of Lake Biwa Museum Comprehensive Research Project S06-02, and M. J. G.’s later contributions were facilitated by support for National Taiwan Ocean University’s Center of Excellence for the Oceans under the Featured Areas Research Center Program within the R.O.C. (Taiwan) Ministry of Education’s Higher Education Sprout Project.
