Global Kinetoplastea phylogeny inferred from a large-scale multigene alignment including parasitic species for better understanding transitions from a free-living to a parasitic lifestyle

Euki Yazaki, Sohta A. Ishikawa*, Keitaro Kume, Akira Kumagai, Takashi Kamaishi, Goro Tanifuji, Tetsuo Hashimoto and Yuji Inagaki Graduate School of Life and Environmental Sciences, University of Tsukuba, 1-1-1 Tennodai, Tsukuba, Ibaraki 305-8572, Japan Department of Biological Sciences, Graduate School of Sciences, University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan Graduate School of Systems and Information Engineering, University of Tsukuba, 1-1-1 Tennodai, Tsukuba, Ibaraki 305-8573, Japan Miyagi Prefecture Fisheries Technology Institute, 97-6 Sodenohama, Watanoha, Ishinomaki, Miyagi 986-2135, Japan National Research Institute of Aquaculture, Fisheries Research Agency, 422-1 Nakatsuhamaura, Minami-Ise, Mie 516-0913, Japan Department of Zoology, National Museum of Nature and Science, 4-1-1 Amakubo, Tsukuba, Ibaraki 305-0005, Japan Center for Computational Sciences, University of Tsukuba, 1-1-1 Tennodai, Tsukuba, Ibaraki 305-8577, Japan


INTRODUCTION
Trypanosomatid flagellates have been studied extensively, as some of them are causative agents of human African trypanosomiasis (sleeping sickness), Chagas disease, and leishmaniasis.Besides their clinical importance, these flagellates possess intriguing properties that are shared by only a few or no other eukaryotes.Trypanosomatids are known to possess a unique peroxisomederived organelle, the glycosome, that encloses glycolytic enzymes (Opperdoes and Borst, 1977;Gualdrón-López et al., 2012).Mitochondria of trypanosomatids contain a complex network of two types of circular DNA molecules, maxicircles and minicircles, and their mitochondrial mRNAs undergo intricate and distinctive editing prior to translation (Lukeš et al., 2002(Lukeš et al., , 2005)).The 5′ termini of mRNAs from trypanosomatid nuclear genomes also undergo post-transcriptional modification (Campbell et al., 2003;Michaeli, 2011).Although the properties described above are also observed in phylogenetic relatives of trypanosomatids (i.e., members of other orders in the class Kinetoplastea; see below), trypanosomatids, for which various experimental techniques in molecular and cell biology (e.g., genetic modification) are available, have been the center of research on Kinetoplastea.
Trypanosomatida, together with Eubodonida, Parabodoida, Neobodonida and Prokinetoplastida, comprise the class Kinetoplastea (Moreira et al., 2004;Simpson et al., 2006).All known members of Trypanosomatida and of Prokinetoplastida are parasites (Simpson et al., 2006;Lukeš et al., 2014).However, the remaining three orders are dominated by free-living members, and only a few or none of the members are known to be parasitic.Previously published phylogenies of small subunit ribosomal DNA (SSU rDNA) sequences have constantly and robustly united Neobodonida, Parabodonida, Eubodonida and Trypanosomatida, excluding Prokinetoplastida (Simpson et al., 2002;Moreira et al., 2004;von der Heyden et al., 2004).This tree topology suggested that Trypanosomatida and Prokinetoplastida acquired parasitic lifestyles separately (Moreira et al., 2004;Simpson et al., 2006;Lukeš et al., 2014).Owing to their importance in public health, the origin of parasitism in the extant trypanosomatids is one of the major questions in the evolution of Kinetoplastea.To address this question, the precise relationship among Neobodonida, Parabodonida, Eubodonida and Trypanosomatida has been explored mainly by analyzing SSU rRNA genes or genes encoding highly conserved proteins, but has not been resolved with high statistical support (Dolezel et al., 2000;Simpson et al., 2002;Moreira et al., 2004;Simpson et al., 2004;von der Heyden et al., 2004;Deschamps et al., 2011).A phylogenetic analysis of 64 genes encoding highly conserved proteins successfully designated Eubodonida as the closest relative of Trypanosomatida (Deschamps et al., 2011).
Compared to pathogenic trypanosomatids, other parasitic members in Kinetoplastea have received less research attention.Two types of parasites belong to Prokinetoplastida: fish parasites (e.g., the causative agent of ichthyobodosis, Ichthyobodo necator: Callahan et al., 2002) and intracellular parasites in the amoebozoan Paramoeba pemaquidensis (i.e., Perkinsela sp. or Ichthyobodo-related organism: Dyková et al., 2003;Caraguel et al., 2007;Dyková et al., 2008;Feehan et al., 2013;Lukeš et al., 2014).Parabodonida includes fish parasites that cause cryptobiosis in salmonid and cyprinid fishes (e.g., Cryptobia salmositica and Trypanoplasma borreli: Woo and Poynton, 1995), as well as the snail parasite C. helicis (Leidy, 1846).Among the known neobodonids, there is a single parasitic member, Azumiobodo hoyamushi, which infects ascidians and causes soft tunic syndrome (Hirose et al., 2012).As the transition from a free-living to a parasitic lifestyle occurred after the divergence of the extant parabodonids/neobodonids, the parasites in Parabodonida and Neobodonida are potentially useful to retrace the evolutionary path from a free-living to a parasitic lifestyle.
For a deeper understanding of the evolution of parasitism in Kinetoplastea, a well-resolved, taxon-rich phylogeny is indispensable.Deschamps et al. (2011) analyzed an alignment of 64 proteins and elucidated the relationship among Trypanosomatida, Eubodonida, Neobodonida and Parabodonida.However, their analyses contained two potential limitations.First, the alignment analyzed in Deschamps et al. (2011) contained no prokinetoplastid species.Second, each of Parabodonida and Neobodonida was represented by only a single free-living species but no parasitic member.In this study, we overcame these limitations by analyzing a new alignment of 43 proteins (43-gene alignment), which covered all five orders in Kinetoplastea, and Parabodonida and Neobodonida were represented by both free-living and parasitic members.Combining the global Kinetoplastea phylogeny, updated by analyzing the 43-gene alignment, with a taxon-rich SSU rDNA phylogeny, we discuss the transition from a free-living to a parasitic lifestyle in the evolution of Kinetoplastea.

MATERIALS AND METHODS
Cultures, RNA extraction and sequencing The laboratory culture of A. hoyamushi established by Hirose et al. (2012) was grown and maintained in sea water containing 2% heat-inactivated fetal bovine serum (Gibco, Thermo Fisher Scientific, Waltham, Massachusetts, USA) at 17 °C.Trypanoplasma borreli ATCC50836 was purchased from the American Type Culture Collection, and grown in live-infusion-tryptose medium (Fernandes and Castellani, 1966) at 17 °C.Total RNA was extracted from the harvested cells using Trizol (Thermo Fisher Scientific, Waltham, Massachusetts, USA), following the manufacturer's protocol.Construction of a cDNA library and subsequent sequencing by the Illumina HiSeq2500 system were performed at Hokkaido System Science (Sapporo, Hokkaido, Japan).We generated 401,725,240 and 433,374,224 paired-end, 100-base reads from the Azumiobodo and Trypanoplasma libraries, respectively (deposited in NCBI Sequence Read Archive under accession numbers SRX2210809 and SRX22109115, respectively).The two sets of initial reads were separately assembled into 28,134 contigs (Azumiobodo) and 18,591 contigs (Trypanoplasma) by Trinity (Grabherr et al., 2011;Haas et al., 2013).Transcriptomic data generated from Perkinsela by the Illumina HiSeq2000 system were retrieved from NCBI Sequence Read Archive (accession number ARX255943), and assembled into 18,600 contigs by Trinity.The contig data of Azumiobodo and Trypanoplasma are freely available at https://sites.google.com/site/eukiyazaki/home/data-archive/kinetoplastida.

Phylogenetic alignments
We prepared a multigene alignment including sequences from Azumiobodo, Trypanoplasma and Perkinsela by following Deschamps et al. (2011), which assessed the phylogenetic relationship among seven kinetoplastids based on 64 protein-coding genes.First, we prepared four sets of the 64 genesthose of Trypanosoma brucei, T. cruzi, Leishmania major and L. infantum-by surveying in NCBI (http://www.ncbi.nlm.nih.gov/),TriTrypDB (http://tritrypdb.org/tritrypdb), and Sanger Institute databases (http://www.sanger.ac.uk/resources/databases/). The putative Trypanosoma/ Leishmania proteins identified in the first survey were then used as queries for TBLASTX against the contig data of Azumiobodo, Trypanoplasma and Perkinsela.We recovered Azumiobodo, Trypanoplasma and Perkinsela transcripts that matched query sequences with E-values below 10 − 10 as candidates encoding the proteins of interest.The candidate transcripts were subjected to BLASTX against the NCBI nr protein database to confirm that these transcripts encode proteins that were homologous to the queries in the initial BLAST analysis (TBLASTX).Finally, the selected transcripts of Azumiobodo, Trypanoplasma and Perkinsela were conceptually translated into amino acid sequences by EMBOSS Transeq (http://www.ebi.ac.uk/Tools/st/emboss_transeq/).We repeated the same procedure above for three kinetoplastids (Rhynchomonas nasuta, Procryptobia sorokini and Bodo saltans: GenBank accession numbers HO651677-HO651928), the amoebozoan Dictyostelium discoideum (NCBI BioProject: PRJNA201), the heterolobosean Naegleria gruberi (NCBI BioProject: PRJNA43691), the diplonemid Diplonema papillatum (TBestDB organism ID: DP) and the euglenid Euglena gracilis (NCBI SRA: ERX324117).Dictyostelium and Naegleria were included in the alignments to identify sequences originated from the amoebozoan Paramoeba pemaquidensis (i.e., the host organism of Perkinsela), which potentially contaminate the Perkinsela transcriptomic data.If the Paramoeba sequences were misassigned as Perkinsela sequences, we anticipated that such sequences would probably group with the Naegleria and Dictyostelium sequences, with the specific affinity to the latter (i.e., amoebozoan) sequences, in the phylogenetic analyses described below.
Sixty-four single-gene alignments (each containing the 14 species described above) were automatically generated by MAFFT 7.205 (Katoh et al., 2002;Katoh and Standley, 2013).Ambiguously aligned positions were manually excluded prior to the phylogenetic analyses described below.The alignments were subjected separately to maximum-likelihood (ML) phylogenetic analyses with the LG model (Le and Gascuel, 2008), incorporating empirical amino acid frequencies and among-site rate variation approximated by a discrete gamma distribution with four categories (LG + Γ + F model).The tree search was started from ten maximum-parsimony (MP) trees, each of which was generated by random stepwise addition of sequences by RAXML 8.0.3 (Stamatakis, 2014).Although not shown here, the ML trees inferred from 21 of the 64 single-gene alignments failed to reconstruct the monophyly of kinetoplastids, and were omitted from the multigene alignment described below.The remaining 43 single-gene alignments listed in Supplementary Table S1 were then concatenated into a single alignment containing 6,842 amino acid positions (43-gene alignment).Prior to phylogenetic analyses, we excluded both Dictyostelium and Naegleria sequences, which were used as "probes" to detect Paramoeba sequences potentially contaminating the Perkinsela transcriptomic data.

Phylogenetic analyses
The 43-gene and SSU rDNA alignments were subjected to ML phylogenetic analyses by RAXML 8.0.3.The detailed settings of tree searching were the same as those used in the single-gene analyses described above, but we assigned the LG + Γ + F model for ML analyses of the 43-gene alignment, and the GTR + Γ model (Rodríguez et al., 1990) for those of the SSU rDNA alignment.One hundred bootstrap replicates were generated from each alignment, and subjected to tree searching as described above.The resultant bootstrap trees were used to calculate ML bootstrap support values (MLBPs).
The two alignments were also analyzed with Bayesian method.The 43-gene alignment was subjected to Phy-loBayes 3.3 (Lartillot and Philippe, 2004;Lartillot et al., 2009) with the CAT + Poisson model.Two Markov chain Monte Carlo (MCMC) runs were conducted for 10,000 cycles with burn-in of 2,500 (the maxdiff value were 0.00013).Subsequently, the consensus tree with branch lengths and Bayesian posterior probabilities (BPPs) was calculated from the remaining trees.The SSU rDNA alignment was subjected to MrBayes 3.2.3(Huelsenbeck and Ronquist, 2001;Huelsenbeck and Ronquist, 2004).We assigned the same substitution model as in the ML method described above.The MCMC run was performed with one cold and three heated chains with default chain temperatures.We ran 1,000,000 generations, and sampled log-likelihood scores and trees with branch lengths every 1,000 generations.The first 25% generations were discarded as burn-in.The consensus tree with branch lengths and BPPs were calculated from the remaining trees.

RESULTS AND DISCUSSION
Global Kinetoplastea phylogeny revisited We updated the large-scale multigene (phylogenomic) alignment used in Deschamps et al. (2011) by adding Trypanoplasma, Azumiobodo and Perkinsela, and reexamined the global phylogeny of Kinetoplastea.The ML tree and MLBPs inferred from a 43-gene alignment are presented in Fig. 1.As the topology inferred from Bayesian method was identical to that from the ML method, we only mapped BPPs on the ML tree shown in Fig. 1.Although the 43-gene phylogeny includes three parasitic kinetoplastids that were absent in the previous phylogenomic analysis, the overall tree topology (Fig. 1) was similar to that presented by Deschamps et al. (2011).Trypanosomatida was found to be tied with the other four taxons/clades in the following order: (i) Bodo saltans representing Eubodonida, (ii) the Parabodonida clade comprising Procryptobia sorokini and Trypanoplasma borreli, (iii) the Neobodonida clade comprising Rhynchomonas nasuta and Azumiobodo hoyamushi, and (iv) Perkinsela sp.representing Prokinetoplastida.All the internal nodes in the ingroup were supported by MLBPs of 100% and BPPs of 1.0, except the clade of Rhynchomonas and Azumiobodo, which received an MLBP of 97% and a BPP of 0.99 (Fig. 1).Prior to this study, the phylogenetic position of Prokinetoplastida has been assessed only by single-gene alignments including kinetoplastid species sampled from all five orders (Callahan et al., 2002;Moreira et al., 2004;Simpson et al., 2004;Breglia et al., 2007;Hirose et al., 2012) or by four-gene and 11-gene alignments with restricted taxon samplings (Tanifuji et al., 2011;Cenci et al., 2016).Thus, the current study provides the first phylogenomic support for the earliest-branching status of Prokinetoplastida in Kinetoplastea (Fig. 1).

Evolution of parasitism in Kinetoplastea
We analyzed an SSU rDNA alignment, in which taxon sampling was much richer than that of the 43-gene alignment, to illustrate the sporadic distribution of parasitic and freeliving species in Kinetoplastea (Fig. 2).Prokinetoplastida was excluded from the clade of Trypanosomatida, Eubodonida, Parabodonida and Neobodonida with an MPBP of 100% and a BPP of 1.0.The monophylies of Eubodonida, Parabodonida and Neobodonida were recovered with MLBPs of 63-91% and BPPs of 0.71-1.0,while that of Trypanosomatida was not positively supported.Although Paratrypanosoma confusum was excluded from the clade of other trypanosomatids in the SSU rDNA phylogeny (Fig. 2), we regard Paratrypanosoma as an ancestral branch in Trypanosomatida based on a series of multigene phylogenetic analyses presented in Flegontov et al. (2013).We assume that the SSU rDNA phylogeny presented here failed to place Paratrypanosoma in the correct position due to lack of phylogenetic signal.Thus, the combination of the 43-gene phylogeny, which resolved the backbone of the tree of Kinetoplastea (Fig. 1), and the taxon-rich SSU rDNA phylogeny (Fig. 2) enables us to depict how and when parasitic species emerged during the evolution of Kinetoplastea.We schematically illustrate the evolution of lifestyles in Kinetoplastea in Fig. 3 (see below for details).As the life cycles of most "freeliving" kinetoplastids are not well understood, we cannot exclude the possibility that some of them have parasitic stages.
The 43-gene phylogeny united the obligatory parasitic order Trypanosomatida with Eubodonida (Fig. 1), which comprises free-living members (see Fig. 2 for the details).Thus, as discussed in Simpson et al. (2006), Deschamps et al. (2011) and Flegontov et al. (2013), a parasitic lifestyle was most likely established after the separation of Trypanosomatida and Eubodonida, but before

Rhynchomonas nasuta
Diplonema papillatum Fig. 1.Global Kinetoplastea phylogeny inferred from an alignment comprising 43 genes encoding highly conserved proteins.Tree topology was inferred from the maximum-likelihood (ML) method.Nodes marked by dots were supported by ML bootstrap values of 100% and Bayesian posterior probabilities of 1.0.For each taxon, the percentage of missing data is presented in parentheses.Parasites are marked by diamonds.Phylogenomic analysis of Kinetoplastea the divergence of the extant trypanosomatids including Paratrypanosoma (Fig. 3).All trypanosomatids known so far are extracellular parasites, but an intracellular stage has been also reported for Leishmania spp.and Trypanosoma cruzi (Tyler and Engman, 2001;Handman and Bullen, 2002).As Leishmania spp.and T. cruzi are distantly related in the SSU rDNA phylogeny (Fig. 2), the ability to invade host cytoplasm was likely acquired by the two separate lineages after the divergence of trypanosmatids (Fig. 3).Parabodonida contains parasitic members, namely Trypanoplasma borreli and Cryptobia spp., as well as free-living members (Fig. 2).In the SSU rDNA phylogeny (Fig. 2), T. borreli, C. catostomi, C. bullocki and C. salmositica formed a robust clade (MLBP of 100% and BPP of 1.0).The four parabodonids are commonly found in the blood stream, although ectoparasitic forms have also been reported for C. bullocki and C. salmositica (Bower and Margolis, 1983;Woo and Wehnert, 1983).We can conclude that an extracellular parasitic lifestyle was established on the branch leading to the Trypanoplasma-Cryptobia clade, as proposed in Simpson et al. (2006).However, we currently have no evidence to determine whether the ectoparasitic form is the ances-tral trait of this clade.Cryptobia helicis (Leidy, 1846), which was found in the seminal receptacle of snails, most likely acquired an extracellular parasitic lifestyle independent from the Trypanoplasma-Cryptobia clade, as the SSU rDNA phylogeny united C. helicis with freeliving Parabodo nitrophilus and P. caudaus (MLBP of 96% and BPP of 1.0; Fig. 2).Altogether, we propose that parasitism emerged at least twice in Parabodonida (Fig. 3).The potential third parasitic lineage in Parabodonida is Jarrellia atramenti, found in the mucus of the respiratory tract of the pygmy sperm whale (Poynton et al., 2001).As only morphological and no molecular data are available for this species, it will be necessary to assess the phylogenetic relationship between Jarrellia and other parabodonids to better understand the evolution of parasitism in this order.
Among the diversity of neobodonids, Azumiobodo is the sole parasitic member known to date (Hirose et al., 2012).This flagellate was found in the tunics of ascidians with soft tunic syndrome by histopathology (Kumagai et al., 2010).Intriguingly, Azumiobodo appeared to bear a phylogenetic affinity to a commensal in the ascidian intestine, Cruzella marina (Frolov and Malysheva, 2002), in the SSU rDNA analysis (MLBP of 79% and BPP of 0.98; Fig. 2).The affinity between Azumiobodo and Cruzella prompts us to propose serial lifestyle changes in Neobodonida, as follows: the common ancestor of Azumiobodo and Cruzella established a symbiotic relationship with ascidians, and, after separation of the two species, the former became a pathogenic extracellular parasite causing soft tunic syndrome.Although the above scenario needs to be examined in future studies, we are confident that the transition from a free-living to a parasitic lifestyle occurred once after the divergence of the extant neobodonids (Fig. 3).
Independent from the lifestyle changes discussed above, the transition from a free-living to a parasitic lifestyle was proposed for Prokinetoplastida (Simpson et al., 2006).As all the members belonging to Prokinetoplastida known to date are parasitic, it is straightforward to assume that they emerged from a parasitic ancestor (Fig. 3).In the SSU rDNA phylogeny (Fig. 2), prokinetoplastids were split into two subclades; one is of parasites that infect the skin, fins and gills of fishes (e.g., Ichthyobodo necator) and the other is of intracellular parasites of amoebozoans (e.g., Perkinsela sp.).The conspicuous difference in parasitic mode between the two subclades in Prokinetoplastida-ectoparasitism and intracellular parasitism-demands a more complex scenario than that assuming a single transition from a free-living to a parasitic lifestyle prior to the divergence of prokinetoplastids (Fig. 3).Notably, the hosts of Perkinsela sp. and its relatives are amoebozoans, which parasitize marine animals (Munday et al., 2001;Young et al., 2007).Thus, it is attractive to hypothesize that the common ancestor of Perkinsela sp. and its relatives was an ectoparasite of a marine animal, then switched its host to an amoebozoan parasitizing marine animals, and finally invaded and settled in the cytoplasm of the amoebozoan host.The above scenario will be favored if future surveys find a novel ectoparasitic species that branches at the base of the clade of Perkinsela sp. and its relatives.It is also important to pursue the possibility that some free-living prokinetoplastids have been overlooked in natural environments, as environmental sequence data hinted that the full diversity of Prokinetoplastida has yet to be unveiled (Moreira et al., 2004;von der Heyden et al., 2004).We will certainly need to reevaluate how parasitism was established in Prokinetoplastida after its organismal diversity is sufficiently depicted in the future.

Future perspectives
The phylogenetic relationship among the five orders in Kinetoplastea was found to be largely unchanged before and after incorporating three parasitic kinetoplastids, Azumiobodo, Trypanoplasma and Perkinsela, into a phylogenomic alignment (Deschamps et al., 2011;this study).However, we still need to improve taxon sampling in phylogenomic alignments to re-examine the global Kinetoplastea phylogeny in the future.For instance, the monophyly of Neobodonida received only weak statistical support in the SSU rDNA analysis (Fig. 2).Thus, future phylogenomic analyses including additional neobodonids are required to strengthen their monophyly inferred from the SSU rDNA alignment.Another potential concern is the diversity of Eubodonida.It is not clear whether the eubodonids identified so far represent the true diversity of this order, as the SSU rDNA phylogeny (Fig. 2) implies that the diversity of Eubodonida is considerably lower than that of any other order in Kinetoplastea.Thus, future studies may identify novel kinetoplastid flagellates that bear phylogenetic affinities to the currently known eubodonids.If such kinetoplastid flagellates exist, it would be worthwhile to incorporate them into phylogenomic analyses to reflect the proper diversity of Eubodonida.
Because of the serious threats to public health they pose, members of Trypanosomatida have attracted more research attention than other kinetoplastids.From an evolutionary biological perspective, trypanosomatids are regarded as model organisms for studying how and when kinetoplastid flagellates acquired a parasitic lifestyle and pathogenicity (Deschamps et al., 2011;Lukeš et al., 2014).Nevertheless, we also anticipate that Parabodonida and Neobodonida will provide insights into the mechanism involved in the transformation of a free-living species into a parasite.In Parabodonida, the snail parasite Cryptobia helicis showed a specific affinity to the free-living members Parabodo caudatus and P. nitrophilus (Fig. 2).This change in lifestyle took place after the divergence of parabodonids (Fig. 3), and we may have a chance to pinpoint the set of genes that played a pivotal role in the change by comparing the genomic and transcriptomic data of the parasite and its closest freeliving relatives.Similarly, the comparison between two closely related neobodonids, Azumiobodo and Cruzella, a parasite and a commensal of ascidians, respectively, may provide insights into the origin of parasitism and pathogenicity.

Fig. 2 .
Fig.2.Global Kinetoplastea phylogeny inferred from an alignment of small subunit ribosomal DNA sequences.Tree topology was inferred from the maximum-likelihood (ML) method.Nodes marked by dots were supported by ML bootstrap values (MLBPs) of 100% and Bayesian posterior probabilities (BPPs) of 1.0.MLBPs and BPPs are presented only for the nodes that are critical for discussing the evolution of lifestyle in Kinetoplastea.Parasites are marked by diamonds.

Fig. 3 .
Fig. 3. Transition from a free-living to a parasitic lifestyle in Kinetoplastea.The branching order among Trypanosomatida, Eubodonida, Parabodonida, Neobodonida and Prokinetoplastida is based on Fig. 1.Blue branches/clades indicate free-living species, while red clades indicate parasitic species.Putative changes in lifestyle are marked by stars.Due to the lack of a precise phylogenetic position, the parasitic parabodonid Jarrellia atramenti is omitted from this figure.Orange triangles indicate the lineages that acquired the ability to invade host cytoplasm, namely Leishmania spp.and Trypanosoma cruzi in Trypanosomatida and the intracellular parasites of marine amoebozoans in Prokinetoplastida.
S. A. I. was supported by research fellowships from the Japan Society for the Promotion of Science (JSPS) for Young Scientists (Nos.2400007 and 2802725).This work was supported in part by grants from JSPS awarded to Y. I. (23117006 and 16H04826) and T. H. (23117005 and 15H05231).The phylogenetic analyses conducted in this work have been carried out under the "Interdisciplinary Computational Science Program" in the Center for Computational Sciences, University of Tsukuba.