The question of how symbionts acquired hosts and diversified phylogenetically during their evolutionary history is a focus of attention in many symbiotic taxa. Marine polyclad flatworms are usually free-living, but some are symbiotic, using animals as hosts. However, the history of their acquisition of symbiotic systems is not well understood. Therefore, we focused on mollusc symbiotic flatworms in the suborder Acotylea and investigated the host specificity and phylogenetic history of symbiotic acquisition. Field surveys revealed that symbiotic flatworms utilized certain molluscs as hosts. In particular, Stylochoplana pusilla and Stylochoplana parasitica utilized different molluscan species as hosts sympatrically. Also, the phylogenetic analysis and the ancestral state reconstruction indicate that the mollusc symbiotic flatworms formed a monophyletic group and that their common ancestor shifted from free-living to mollusc symbiosis. These results suggest that each of the flatworms did not independently acquire a symbiotic system with molluscan hosts during its phylogenetic history, but that their common ancestor acquired a mollusc symbiotic system, which then underwent acquisition of host specificity and speciation. This study emphasises that multiple host use can be a driving force for niche advancement and speciation in the symbionts.

To investigate the phylogenetic relationships of the symbiotic flatworms, the whole mitochondrial genome and the nuclear ribosomal sequences of S. pusilla, S. parasitica and Stylochoplana sp. were determined. DNA extraction was conducted using the DNeasy Blood and Tissue Kit (Qiagen, Hilden, Germany). Partial sequences for CO1 and 28S were PCR-amplified following the PCR conditions described by the previous study ¹, and were obtained through sequencing on ABI 3130xl automated DNA sequencer at Atmosphere and Ocean Research Institute (AORI), The University of Tokyo (Kashiwa, Chiba). The extracted DNA was quantified using a Nanodrop One (Thermo Scientific, Wilmington, DEL) and used for the library preparation for next-generation sequencing analysis. The DNA shotgun sequencing was conducted using a Miseq System (Illumina, San Diego, CA) with paired-end 300 bp reads at the National Institute for Environmental Studies (Tsukuba, Japan) for S. parasitica, and using DNBSEQ-T7 system (MGI Tech, Shenzhen, China) with pair-end 150 bp reads at Bioengineering Lab. Co., Ltd. (Sagamihara, Japan) for S. pusilla and Stylochoplana sp., respectively. The raw reads (Table S2) for each species were assembled into mitogenome and nuclear ribosomal DNA sequences using the NovoPlasty ver. 4.3.4 ². Partial CO1 and 28S sequences obtained from the same specimens were used as seed sequences for the NovoPlasty assembly of the mitogenome and nuclear ribosomal DNA sequences for three Stylochoplana species. The annotation of the mitochondrial genomes was initially performed using the MITOS 2 web server ³, followed by manual verification of the annotation in Mesquite 3.11 ⁴. Due to the difficulty in automatic annotating of the ATP8 gene of Platyhelminthes species with MITOS ⁵, ATP8 was re-annotated through manual inspection methods (e.g. ^5–7). The putative open reading frames (ORFs) of ATP8 were expected to be located in gap regions without annotations. To identify the ATP8 gene, ORFfinder (https://www.ncbi.nlm.nih.gov/orffinder/) was used to analyze long gap regions. The gap regions lacking annotations, which contained putative ORFs, were re-annotated using the MITOS 2 web server, resulting in the prediction of the ATP8 gene. The predicted ATP8 sequences were validated through BLAST searches ⁸ and manually compared among the three Stylochoplana species analysed. Circular mitogenome maps were constructed using Proksee ⁹. The 28S and 18S sequences of nuclear ribosomal DNA were also annotated using BLAST. We constructed phylogenetic trees using 18S, 28S, 16S and CO1, which are commonly used in previous studies, and the full length of the mitochondrial genome, respectively. The mitochondrial and nuclear genome sequences newly determined in this study have been deposited in the DDBJ (DNA Data Bank of Japan) (Table S2). The sequences of 18S, 28S, 16S, and CO1 and the mitogenome sequences used for phylogenetic analysis have been summarized (Table S3 and Table S4). All sequences containing gaps were trimmed using trimAL after alignment ¹⁰, with the -gt option set to 0.8. The final trimmed alignments resulted in 2,657 bp for the 18S, 28S, 16S and CO1 dataset and 11,710 bp for the mitochondrial genome. Phylogenetic trees were constructed using Bayesian inference (BI) in MrBayes and maximum likelihood (ML) methods in IQ-TREE. ^11,12. Evolutionary models were selected using Kakusan and ModelFinder respectively, and the phylogenetic trees were finally constructed using the best model based on BIC (Table S5 and Table S6) ^13,14. Two Markov chains were run for Bayesian estimation. The BI tree was constructed using 18S, 28S, 16S, and CO1 with the following settings: ngen = 15,000,000, samplefreq = 1000, burn-in = 375,000. The construction of a BI tree of the whole mitochondrial genome was constructed with the following settings: ngen = 60,000,000, samplefreq = 1,000, burn-in = 15,000,000. Convergence of the model was verified using Tracer v1.7.1 ¹⁵. The model was considered to have converged if the ESS (MCMC effective sample size) of all statistics from the MCMC simulation was greater than 100. For ML tree, ultrafast bootstrapping was performed with 1,000 replicates ¹⁶.

              Ancestral trait reconstruction was performed using BayesTraitsV3 to estimate whether the common ancestor was a free-living or a symbiont ¹⁷. Branch lengths were scaled to an average of 0.1 with reference to the BayesTraitsV3 manual. The simulation was conducted using reversible-jump MCMC with the following settings: ngen = 15,000,000, samplefreq = 1000, burn-in = 5,000,000. The exponential reversible jump prior distribution (RevJump exp) was set to 10. Marginal likelihood estimation was performed using the stepping-stone sampler, with the number of stones set at 100 and the number of iterations at 1000. Ancestral trait reconstruction was based on the BI tree, considering nodes with high support (Bayesian posterior probability > 0.9 or ultrafast bootstrapping > 90%). Acotylean flatworms are typically free-living that attach to rocks and algae, but they also interact with other organisms in various ecological contexts, such as predator-prey relationships ^18–20. Therefore, based on clearer information, we defined only those species for which a symbiotic relationship (symbiosis or parasitism) was mentioned in previous studies as symbiotic species, and incorporated their host organisms as traits in the model. Specifically, I. zebra, which utilises hermit crabs ²¹, and E. lysiosquillae, which utilises the mantis shrimp Lysiosquilla maculata ²², were defined as symbionts, in addition to the mollusc-utilising flatworms, S. pusilla, S. parasitic, and Stylochoplana sp.. Their respective host organisms (molluscs, hermit crabs and the mantis shrimp) were considered as traits (Table S1). Species that were shown in previous studies to be free-living, or to have been attached to rocky or boulder sites, or that were not mentioned as being symbiotic, were classified as free-living.

1.  Oya, Y. & Kajihara, H. Molecular phylogenetic analysis of Acotylea (Platyhelminthes: Polycladida). Zoolog. Sci. 37, 271–279 (2020).
2.  Dierckxsens, N., Mardulyn, P. & Smits, G. NOVOPlasty: de novo assembly of organelle genomes from whole genome data. Nucleic acids research 45, e18–e18 (2017).
3.  Bernt, M. et al. MITOS: improved de novo metazoan mitochondrial genome annotation. Mol. Phylogenet. Evol. 69, 313–319 (2013).
4.  Maddison, W. P. & Maddison, D. R. Mesquite: a modular system for evolutionary analysis. Version 3.11. http://www.mesquiteproject.org. (2016).
5.  Egger, B., Bachmann, L. & Fromm, B. Atp8 is in the ground pattern of flatworm mitochondrial genomes. BMC Genomics 18, 414 (2017).
6.  Monnens, M. et al. The first mitochondrial genomes of endosymbiotic rhabdocoels illustrate evolutionary relaxation of atp8 and genome plasticity in flatworms. Int. J. Biol. Macromol. 162, 454–469 (2020).
7. Shimada, D., Hiruta, S. F., Takahoshi, K. & Kajihara, H. Does atp8 exist in the mitochondrial genome of Proseriata (Metazoa: Platyhelminthes)? Anim. Gene 30, 200161 (2023).
8. Altschul, S. F., Gish, W., Miller, W., Myers, E. W. & Lipman, D. J. Basic local alignment search tool. J. Mol. Biol. 215, 403–410 (1990).
9. Grant, J. R. et al. Proksee: in-depth characterization and visualization of bacterial genomes. Nucleic Acids Res. 51, W484–W492 (2023).
10. Capella-Gutiérrez, S., Silla-Martínez, J. M. & Gabaldón, T. trimAl: a tool for automated alignment trimming in large-scale phylogenetic analyses. Bioinformatics 25, 1972–1973 (2009).
11. Ronquist, F. et al. MrBayes 3.2: efficient Bayesian phylogenetic inference and model choice across a large model space. Syst. Biol. 61, 539–542 (2012).
12. 12.  Minh, B. Q. et al. IQ-TREE 2: new models and efficient methods for phylogenetic inference in the genomic era. Mol. Biol. Evol. 37, 1530–1534 (2020).
13. 13.  Tanabe, A. S. Kakusan4 and Aminosan: two programs for comparing nonpartitioned, proportional and separate models for combined molecular phylogenetic analyses of multilocus sequence data. Mol. Ecol. Resour. 11, 914–921 (2011).
14. 14.  Kalyaanamoorthy, S., Minh, B. Q., Wong, T. K. F., von Haeseler, A. & Jermiin, L. S. ModelFinder: fast model selection for accurate phylogenetic estimates. Nat. Methods 14, 587–589 (2017).
15. 15.  Rambaut, A. Tracer v1. http://beast.community/tracer (2007).
16. 16.  Minh, B. Q., Nguyen, M. A. T. & von Haeseler, A. Ultrafast approximation for phylogenetic bootstrap. Mol. Biol. Evol. 30, 1188–1195 (2013).
17. 17.  Meade, A. & Pagel, M. BayesTraits V3. https://www.evolution.reading.ac.uk (2017).
18. 18.  Gammoudi, M., Ben Ahmed, R., Bouriga, N., Ben-Attia, M. & Harrath, A. H. Predation by the polyclad flatworm Imogine mediterranea on the cultivated mussel Mytilus galloprovincialis in Bizerta Lagoon (northern Tunisia). Aquac. Res. 48, 1608–1617 (2017).
19. 19.  Landers, W. S. & Rhodes, E. W. Some factors influencing predation by the flatworm,Stylochus ellipticus (Girard), on oysters. Chesapeake Science 11, 55–60 (1970).
20. 20.  Rawlinson, K. A. Biodiversity of coastal polyclad flatworm assemblages in the wider Caribbean. Mar. Biol. 153, 769–778 (2008).
21. 21.  Lytwyn, M. W. & McDermott, J. J. Incidence, reproduction and feeding of Stylochus zebra, a polyclad turbellarian symbiont of hermit crabs. Mar. Biol. 38, 365–372 (1976).
22. 22.  Oya, Y., Nakajima, H. & Kajihara, H. A new symbiotic relationship between a polyclad flatworm and a mantis shrimp: description of a new species of Emprosthopharynx (Polycladida: Acotylea) associated with Lysiosquilla maculata (Crustacea: Stomatopoda). Mar. Biodivers. 52, (2022).