Comparative Mitogenomic Analyses and New Insights into the Phylogeny of Thamnocephalidae (Branchiopoda: Anostraca)

Thamnocephalidae, a family of Anostraca which is widely distributed on all continents of the world except Antarctica, currently consists of six genera and approximately 63 recognized species. The relationships among genera in Thamnocephalidae and the monophyly of Thamnocephalidae, determined using morphological characteristics or gene markers, remain controversial. In order to address the relationships within Thamnocephalidae, we sequenced Branchinella kugenumaensis mitogenomes and conducted a comparative analysis to reveal the divergence across mitogenomes of B. kugenumaensis. Using newly obtained mitogenomes together with available Anostracan genomic sequences, we present the most complete phylogenomic understanding of Anostraca to date. We observed high divergence across mitogenomes of B. kugenumaensis. Meanwhile, phylogenetic analyses based on both amino acids and nucleotides of the protein-coding genes (PCG) provide significant support for a non-monophyletic Thamnocephalidae within Anostraca, with Asian Branchinella more closely related to Streptocephalidae than Australian Branchinella. The phylogenetic relationships within Anostraca were recovered as follows: Branchinectidae + Chirocephalidae as the basal group of Anostraca and halophilic Artemiidae as a sister to the clade Thamnocephalidae + Streptocephalidae. Both Bayesian inference (BI)- and maximum likelihood (ML)-based analyses produced identical topologies.

Thamnocephalidae currently comprises 63 valid species widespread in Europe, Asia, North and South America, southern Africa, and Australia [1]. Thamnocephalids typically have

Phylogenetic Analysis
In total, mitochondrial 16S rRNA genes of 41 specimens and mitochondrial cytochrome C oxidase subunit I (cox1) genes of 30 specimens were retrieved from GenBank, representing 38 species, 12 genera, and 6 families of Anostraca (Table S2). Three species belonging to the Branchinectidae clade and five species belonging to the Chirocephalidae clade were used as outgroups. We estimated the taxonomic status within genus Branchinella by reconstructing a phylogenetic tree using 16S rRNA datasets. To better resolve the internal phylogeny of Thamnocephalidae, we then reconstructed the phylogenetic tree using 16S + cox1 datasets. In order to confirm the relationships between Thamnocephalidae and the related families, we conducted phylogenetic analyses with 27 complete mitogenomes of Branchiopoda, representing 16 species, 7 genera, and 5 families of Anostraca (Table S3). Three Spinicaudatan species and six Notostracan species were used as outgroups. The amino acid sequences of 13 PCGs were aligned using MUSCLE implemented in the MEGA X software [33]. The corresponding nucleotide sequences for each PCG were aligned using the aligned amino acid sequences implemented in DAMBE 6 [31]. Ribosomal RNAs were aligned by published rRNAs secondary structures for Triops granarius (Lucas, 1864) [29] in order to improve both the alignment and the tree reconstruction processes of rRNA data sets. We estimated saturation for four subset types of the concatenated data set: the first codon positions (nt1), second codon positions (nt2), the first codon positions + second codon positions (nt12), and the third codon positions (nt3) of PCGs. The best-fit nucleotide substitution models were determined by jModelTest version 0.1.1 software [34]. ProtTest 3 [35] was used for the amino acid dataset. The best selected partition schemes and models are listed in Table  S4. Phylogenetic reconstructions were conducted using ML and BI methods. ML was performed with RAxML 7.0.3 [36], with the most appropriate substitution model for each of the separate partition (1000 bootstrap). BI was performed with MrBayes 3.2 [37], with two simultaneous runs (4 chains) for 10 million generations, sampled every 100 generations, of which the first 25% were discarded as burn-in. Convergence was assessed using the Tracer v1.5 software [38].

General Features of B. kugenumaensis Mitogenomes
The complete mitochondrial genome for B. kugenumaensis from Suqian City, Jiangsu Province of China, was sequenced, assembled, and deposited in GenBank (accession number: OP133270). The total length of the mitogenome was 14,126 base pairs (bp), consistent with the previously sequenced mitogenome of B. kugenumaensis from Naruto City, Tokushima Prefecture, Japan (MW136376) [28]: 14,123 bp; (Table 1, Figure 1). The complete mitogenome of B. kugenumaensis from Jinhua City, Zhejiang Province, China (MN660045) [39], had longer total length (15,127 bp) due to an expanded control region (1182 bp). Similarly, the gene arrangement in the mitogenomes of B. kugenumaensis from the 3 localities consisted of 13 PCGs, 22 tRNAs, two mitochondrial rRNAs (rrnS and rrnL), a putative control region, and a large number of intergenic sequences. Twenty genes were encoded by the majority strand (J-strand) and seventeen by the minority strand (N-strand). The canonical "ATN" start codon was the most commonly used start codon in the core PCGs of B. kugenumaensis, with the exception of the cox1 and NADH dehydrogenase, subunit 1 (nad1) genes in B. kugenumaensis of Naruto and Suqian (MW136376 and OP133270), which used TTG as start codons. In cox1, nad1, NADH dehydrogenase, subunit 5 (nad5), and ATP synthase subunits 6 (atp6) genes of B. kugenumaensis from Jinhua (MN660045), TTG and GTG were used as start codons (Table 1). TAA/TAG was most commonly used as a stop codon in core PCGs of B. kugenumaensis mitogenomes, except for cox1, cytochrome C oxidase subunits III (cox3), and NADH dehydrogenase; subunit 4 (nad4) in B. kugenumaensis from Naruto and Suqian, and cox1; and cytochrome C oxidase subunits II (cox2), cox3, nad4 and nad5 in B. kugenumaensis from Jinhua, which ended with truncated termination codons (T), see Table 1. synthase subunits 6 (atp6) genes of B. kugenumaensis from Jinhua (MN660045), TTG and GTG were used as start codons (Table 1). TAA/TAG was most commonly used as a stop codon in core PCGs of B. kugenumaensis mitogenomes, except for cox1, cytochrome C oxidase subunits III (cox3), and NADH dehydrogenase; subunit 4 (nad4) in B. kugenumaensis from Naruto and Suqian, and cox1; and cytochrome C oxidase subunits II (cox2), cox3, nad4 and nad5 in B. kugenumaensis from Jinhua, which ended with truncated termination codons (T), see Table 1. The A + T content and skewness levels for the major strands of B. kugenumaensis and Streptocephalus cafer (Lovén, 1847; sensu [41]) are summarized in Table S5. All were ATrich (67.78~68.17%) throughout the entire genomes of both species. Little heterogeneity was observed for the AT%, AT-and GC-skews in Naruto and Suqian samples. The A + T content for PCGs was the lowest in cox1 (62.88% for Jinhua sample and 63.46% for Suqian sample) and highest in NADH dehydrogenase, subunit 6 (nad6) gene (76.29% for Suqian sample and 73.6% for Jinhua sample). Among the three codon positions within the 13 PCGs, the highest AT% was found for the third codon position, which was in accordance with some pancrustacean groups (see, e.g., [42]). Positive AT skew values were observed The A + T content and skewness levels for the major strands of B. kugenumaensis and Streptocephalus cafer (Lovén, 1847; sensu [41]) are summarized in Table S5. All were AT-rich (67.78~68.17%) throughout the entire genomes of both species. Little heterogeneity was observed for the AT%, AT-and GC-skews in Naruto and Suqian samples. The A + T content for PCGs was the lowest in cox1 (62.88% for Jinhua sample and 63.46% for Suqian sample) and highest in NADH dehydrogenase, subunit 6 (nad6) gene (76.29% for Suqian sample and 73.6% for Jinhua sample). Among the three codon positions within the 13 PCGs, the highest AT% was found for the third codon position, which was in accordance with some pancrustacean groups (see, e.g., [42]). Positive AT skew values were observed only for the entire mitogenomes and tRNA, confirming the existence of more adenine than thymine in this organism. For GC-skew, the first codon position of PCGs, rRNA, tRNAs, cox3, and nad1 were markedly positive, whereas the other part of mitogenomes showed negative values.
The patterns of amino acid composition and synonymous codon usage were similar in the three mitogenomes of B. kugenumaensis ( Figure S1a

Sequence Divergence within B. kugenumaensis Mitogenomes
The two mitogenomes of B. kugenumaensis from Naruto and Suqian (MW136376 and OP133270) were similar with 1% nucleotide dissimilarity (141 SNPs) across the complete mitogenome alignment, with 106 SNPs distributed across the 13 PCGs, of which NADH dehydrogenase, subunit 4L (nad4L) gene (264 bp) comprised the highest proportion (1.89%) relative to the size of the gene, whereas NADH dehydrogenase, subunit 3 (nad3) gene (345 bp) had the lowest proportion (0.29%). Notably, a high level of variation was detected in B. kugenumaensis from Jinhua. It was unexpected that a total of 2322 SNPs distinguished B. kugenumaensis from Jinhua and Suqian in the PCGs. There were 255 SNPs (11%) detected in the cox1 gene, which has been proposed as the universal barcode locus for Anostraca. High-sequence variation was also observed in the tRNA genes, excluding tRNA-Met, tRNA-Glu, and tRNA-Ser1, which were perfectly conserved ( Figure S2). A total of 115 SNPs occurred in the 22 mtDNA-encoded tRNA genes between B. kugenumaensis from Jinhua and Suqian (MN660045 and OP133270) and 45 SNPs in the stems. Instead, only three tRNA SNPs were identified in B. kugenumaensis from Suqian and Naruto (OP133270 and MW136376).
The overall level of mitogenome diversity, represented by the mean pairwise differences per site (π), between the B. kugenumaensis mitogenome sequences OP133270 and MN660045 was 22.0% for the 13 PCGs concatenated, which is much higher than that of observed between the B. kugenumaensis mitogenome sequences OP133270 and MW136376 (π = 0.99% for 13 PCGs). The high value was on the same order of magnitude as for species of Streptocephalus (π = 21.5%). For each of the 13 PCGs, the diversity of nucleotides ranged from 0.029% (nad3) to 1.92% (nad4L) between OP133270 and MW136376. Nevertheless, the average values of nucleotide diversity for B. kugenumaensis from Jinhua (MW136376) was much higher, ranging from 16.4% (cox1) to 29.64% (nad6), comparable to interspecific variation for species of Streptocephalus (27.46% ≥ π ≥ 16.%; Figure 2a).
The proportion of amino acid change for each of the PCGs ranged from 0% (cox3 and atp6) to 2.88% (nad4L) in the two B. kugenumaensis mitogenome sequences MW136376 and OP133270 (Figure 2b; Table S6). Meanwhile, more divergence was observed in B. kugenumaensis mitogenome sequence MN660045: 10.1-23.64% for complex I genes and ATP synthase (complex V) genes and 1.75-3.41% for complex IV genes (Table S6), which was also comparable to interspecific comparison for species of Streptocephalus (Figure 2b). ences per site (π), between the B. kugenumaensis mitogenome sequences OP133270 and MN660045 was 22.0% for the 13 PCGs concatenated, which is much higher than that of observed between the B. kugenumaensis mitogenome sequences OP133270 and MW136376 (π =0.99% for 13 PCGs). The high value was on the same order of magnitude as for species of Streptocephalus (π = 21.5%). For each of the 13 PCGs, the diversity of nucleotides ranged from 0.029% (nad3) to 1.92% (nad4L) between OP133270 and MW136376. Nevertheless, the average values of nucleotide diversity for B. kugenumaensis from Jinhua (MW136376) was much higher, ranging from 16.4% (cox1) to 29.64% (nad6), comparable to interspecific variation for species of Streptocephalus (27.46% ≥ π ≥ 16.%; Figure 2a). Comparing the rate of non-synonymous (Ka) and synonymous nucleotide substitutions (Ks) is important in understanding the dynamics of molecular sequence evolution [43][44][45][46]. Here, among the 13 genes of the B. kugenumaensis mitogenome sequences OP133270 and MW136376, ATPase synthase 8 (atp8) and nad4L genes had a Ka/Ks value greater than 1, suggesting that they were under positive selection during evolution and remained fixed in the population. Interspecific comparison for species of Streptocephalus presented an average Ka/Ks ratio ranging from 0.0063 in cox2 to 0.26 in atp8, and B. kugenumaensis from Jinhua (MN660045) had an average Ka/Ks ratio ranging from 0.0067 in cox2 to 0.25 in atp8, signifying lesser amino acid changes and indicating purifying selection. Phylogenetic analyses using 16S rRNA sequence indicated that Branchinella was initially divided into three major clades: Australian Branchinella, B. maduraiensis, and other Asian Branchinella (Figure 3). In the clade of Australian Branchinella, our phylogenetic reconstruction was largely congruent with the results of previous studies [47,48]. However, our phylogenetic analysis placed B. campbelli at a basal position within Australian Branchinella, which was the major difference from those results. Interestingly, B. maduraiensis was closely related to Streptocephalus (Bayesian posterior probability, BPP = 1.0), suggesting that Branchinella is not a monophyletic group. B. kugenumaensis from Suqian was sister to B. kugenumaensis from Naruto, whereas B. kugenumaensis from Jinhua formed a separate branch from all the other Asian Branchinella and was basal within this clade (BPP = 0.90). Both the ML and BI trees were highly similar in topology (Figure 3). The sequence identity for each PCG in the pairwise mitochondrial genomes of B. kugenumaensis also indicated that sequences from Naruto had around 16.46-29.78% divergence from Jinhua sequences and 0.3-1.89% from the Suqian sequences (Table S5). The results of the phylogenetic and genetic analyses support the idea that B. kugenumaensis from Jinhua is a different species. This study is the first to include three species of subgenus Branchinellites. Whether Branchinellites is monophyletic remains an open question, and more molecular loci and taxa are needed in future surveys to improve the resolution within Branchinella.

Paraphyly of Thamnocephalidae
The monophyly of Thamnocephalidae has been questioned using different molecular datasets [23,24]. Increased sampling across thamnocephalid mtDNAs has provided valuable information for understanding the phylogeny of Thamnocephalidae. The reconstruction of the phylogenetic relationships within Thamnocephalidae based on 16S + cox1 is presented in Figure 4. Contrary to traditional morphological taxonomy, the BI tree indicated that Thamnocephalus was sister to the clade Branchinella + Streptocephalus (BPP = 0.90), while Phallocryptus constituted the sister group of all other clades in Thamnocephalidae (BPP = 1.0). The non-monophyly of Thamnocephalidae was revealed with good support in this study, which was congruent with previous studies [23,24]. Furthermore, we identified that Asian Branchinella clustered to Streptocephalus (Streptocephalidae) as a sister group (BPP = 0.98), rather than Australian Branchinella. Although we clarified some relationships in Branchinella, it must be taken with caution, as only species of the genera Thamnocephalus, Phallocryptus, and Branchinella have been investigated so far. Future research is still needed with increased species sampling from the genera of Carinophallus, Dendro-

Paraphyly of Thamnocephalidae
The monophyly of Thamnocephalidae has been questioned using different molecular datasets [23,24]. Increased sampling across thamnocephalid mtDNAs has provided valuable information for understanding the phylogeny of Thamnocephalidae. The reconstruction of the phylogenetic relationships within Thamnocephalidae based on 16S + cox1 is presented in Figure 4. Contrary to traditional morphological taxonomy, the BI tree indicated that Thamnocephalus was sister to the clade Branchinella + Streptocephalus (BPP = 0.90), while Phallocryptus constituted the sister group of all other clades in Thamnocephalidae (BPP = 1.0). The non-monophyly of Thamnocephalidae was revealed with good support in this study, which was congruent with previous studies [23,24]. Furthermore, we identified that Asian Branchinella clustered to Streptocephalus (Streptocephalidae) as a sister group (BPP = 0.98), rather than Australian Branchinella. Although we clarified some relationships in Branchinella, it must be taken with caution, as only species of the genera Thamnocephalus, Phallocryptus, and Branchinella have been investigated so far. Future research is still needed with increased species sampling from the genera of Carinophallus, Dendrocephalus, and Spiralifrons in order to explore the detailed relationships among genera of Thamnocephalidae and Streptocephalidae.

Phylogenetic Relationships among Anostracan Families
With respect to 11 additional branchiopod mitogenomes (Tables S3 and S4), both BI and ML analyses based on nucleotide sequence data from 9 PCGs of the 16 Anostracan mitogenomes produced identical topologies with similar branch lengths, and most clades had high ML bootstrap support or BI posterior probabilities (Figure 5a). Both the BI and ML trees strongly supported the monophyly of the order Anostraca (Bootstrap support value, BP = 100, BPP = 1.0). Meanwhile, the phylogenetic topology of monophyletic Streptocephalus nested within a paraphyletic Thamnocephalidae was confirmed (BP = 100, BPP = 1.0).
The monophyly of Anostraca and the major clade assignments of Anostraca (halophilic Artemiidae + Paratemia separate from the other Anostracan families) have been broadly accepted according to molecular-based phylogenies [14,17,20,24]. However, this traditional classification of Anostraca has been challenged by molecular analysis carried out on large-genome-scale molecular sequences, supporting Branchinectidae as the basal clade of Anostraca [23]. In this study, Anostraca was clearly separated into two major groups-1) Branchinectidae, together with Chirocephalidae, as the basal group of Anostraca; and 2) halophilic Artemiidae as a sister to the clade Thamnocephalidae + Streptocephalidae-with high support values (BP = 84, BPP = 1.0). Such a relationship has also been supported previously by a morphological study [9] and molecular analysis based on large-genome-scale molecular sequencing [23]. In addition, when they were applied to amino acid sequences of nine PCGs, the BI and ML methods supported almost identical phylogenetic relationships, except for the relationships within Artemiidae (Figure 5b).
The present phylogenetic analyses shed new light on the evolution of morphological traits. Our phylogenetic analysis confirmed the paraphyletic status of Thamnocephalidae, suggested that the diagnostic characteristics seem to be plesiomorphic (e.g., eversible portion becoming explanate distally; second maxillae with single apical setae). One or more longitudinal rows of spines on the genitalia and females with an elongated brood pouch, extending to the base of the fourth to seventh abdominal segments, may be ancestral characteristics of Thamnocephalidae+ Streptocephalidae. The present phylogenetic analyses

Phylogenetic Relationships among Anostracan Families
With respect to 11 additional branchiopod mitogenomes (Tables S3 and S4), both BI and ML analyses based on nucleotide sequence data from 9 PCGs of the 16 Anostracan mitogenomes produced identical topologies with similar branch lengths, and most clades had high ML bootstrap support or BI posterior probabilities (Figure 5a). Both the BI and ML trees strongly supported the monophyly of the order Anostraca (Bootstrap support value, BP = 100, BPP = 1.0). Meanwhile, the phylogenetic topology of monophyletic Streptocephalus nested within a paraphyletic Thamnocephalidae was confirmed (BP = 100, BPP = 1.0).
The monophyly of Anostraca and the major clade assignments of Anostraca (halophilic Artemiidae + Paratemia separate from the other Anostracan families) have been broadly accepted according to molecular-based phylogenies [14,17,20,24]. However, this traditional classification of Anostraca has been challenged by molecular analysis carried out on large-genome-scale molecular sequences, supporting Branchinectidae as the basal clade of Anostraca [23]. In this study, Anostraca was clearly separated into two major groups-(1) Branchinectidae, together with Chirocephalidae, as the basal group of Anostraca; and (2) halophilic Artemiidae as a sister to the clade Thamnocephalidae + Streptocephalidae-with high support values (BP = 84, BPP = 1.0). Such a relationship has also been supported previously by a morphological study [9] and molecular analysis based on large-genome-scale molecular sequencing [23]. In addition, when they were applied to amino acid sequences of nine PCGs, the BI and ML methods supported almost identical phylogenetic relationships, except for the relationships within Artemiidae (Figure 5b).
The present phylogenetic analyses shed new light on the evolution of morphological traits. Our phylogenetic analysis confirmed the paraphyletic status of Thamnocephalidae, suggested that the diagnostic characteristics seem to be plesiomorphic (e.g., eversible portion becoming explanate distally; second maxillae with single apical setae). One or more longitudinal rows of spines on the genitalia and females with an elongated brood pouch, extending to the base of the fourth to seventh abdominal segments, may be ancestral characteristics of Thamnocephalidae+ Streptocephalidae. The present phylogenetic analyses support that Branchinectidae + Chirocephalidae were the first to diverge from the Anostraca, which implies that each thoracopod bearing a single pre-epipodite and distinct seminal vesicle were already present in the ground pattern of Anostraca. Furthermore, the origin of the halophilic Anostracans is more complex than expected and needs further investigation.
Genes 2022, 13, x FOR PEER REVIEW 10 of 14 support that Branchinectidae + Chirocephalidae were the first to diverge from the Anostraca, which implies that each thoracopod bearing a single pre-epipodite and distinct seminal vesicle were already present in the ground pattern of Anostraca. Furthermore, the origin of the halophilic Anostracans is more complex than expected and needs further investigation.  Supplementary Materials: The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/genes13101765/s1, Figure S1: Amino acid composition (a) and relative synonymous codon usage (b) in the mitogenomes of B. kugenumaensis; Figure S2: Polymorphisms in inferred secondary structures of 22 transfer RNAs found in the mitogenomes of B. kugenumaensis; Table S1: List of primer combinations used to amplify the mitochondrial genome of B. kugenumaensis; Table S2: COI & 16S rRNA sequences of Anostraca retrieved from GenBank used in this study; Table S3: Details of species and mitogenomes of Branchiopoda used in this study; Table S4: Partition schemes and best-fitting models for phylogenetic analyses; Table S5: A + T content (%) and skewness levels calculated for major strand of B. kugenumaensis and S. cafer; Table S6: Nucleotide identity (%) and number of codon substitutions for each PCG of three B. kugenumaensis mtgenomes. The references [49][50][51][52][53][54][55][56][57][58][59][60][61][62][63][64] are cited in the Supplementary Materials.