Three nuclear protein-coding genes corroborate a recent phylogenomic model of the Branchiopoda (Crustacea) and provide estimates of the divergence times of the major branchiopodan taxa

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 clarify-ing 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 conﬁrming 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 ﬁrst divergence of extant branchiopods occurred about 534 Ma during the early Cambrian period and that diversiﬁcation within the extant branchiopod lineages started in or after the late Permian.

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 cur-rently 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 Lozano-Fernandez et al., 2019). For a review of this topic, see Giribet and Edgecombe (2012).
As for molecular evidence, mitochondrial 12S rRNA data supported Notostraca as the sister group to Laevicaudata (Braband et al., 2002) (Fig. 1D) (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 ( 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 morphologybased 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).

MATERIALS AND METHODS
Taxonomic sampling 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.
RNA isolation, reverse transcription, PCR and sequencing 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 analyses
Nucleotide 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 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.
Divergence time estimation An 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 Stecher et al., 2020) based on the 22-sample dataset using the RelTime method (Tamura et al., 2012) 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).

RESULTS AND DISCUSSION
Sequence data and alignment 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 Branchiopoda
Using 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, 2007Olesen, , 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 Branchiopoda
Both 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, 2007Olesen, , 2009; see also Olesen and Richter, 2013), and furthermore agree with the results of a recent phylogenomic study (Schwentner et al., 2018).

Phylogeny of Cladocera As shown in
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,  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 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. 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.
Conclusion The aa sequences of three protein-coding genes (RPB1, RPB2 and DPD1) proved to be highly reliable in inferring the phylogenetic relationships of higherlevel 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.