Polyclad phylogeny persists to be problematic

Two conflicting morphological approaches to polyclad systematics highlight the relevance of molecular data for resolving the interrelationships of Polycladida. In the present study, phylogenetic trees were reconstructed based on a short alignment of the 28S rDNA marker gene with 118 polyclad terminals (24 new) including 100 different polyclad species from 44 genera and 22 families, as well as on a combined dataset using 18S and 28S rDNA genes with 27 polyclad terminals (19 new) covering 26 different polyclad species. In both approaches, Theamatidae and Cestoplanidae were included, two families that have previously been shown to switch from Acotylea to Cotylea. Three different alignment methods were used, both with and without alignment curation by Gblocks, and all alignments were subjected to Bayesian inference and maximum likelihood tree calculations. Over all trees of the combined dataset, an extended majority-rule consensus tree had weak support for Theamatidae and Cestoplanidae as acotyleans, and also the cotylean genera Boninia, Chromyella and Pericelis appeared as acotyleans. With the most inclusive short 28S dataset, on the other hand, there is good support for the aforementioned taxa as cotyleans. Especially with the short 28S matrix, taxon sampling, outgroup selection, alignment method and curation, as well as model choice were all decisive for tree topology. Well-supported parts of the phylogeny over all trees include Pseudocerotoidea, Prosthiostomoidea, Stylochoidea, Leptoplanoidea and Cryptoceloidea, the latter three with new definitions. Unstable positions in the tree were found not only for Theamatidae, Cestoplanidae, Boninia, Chromyella and Pericelis, but also for Anonymus, Chromoplana and Cycloporus.


Introduction
Due to their colourful appearance, polyclad flatworms are among the most conspicuous members of the phylum Platyhelminthes, yet these animals are relatively poorly studied . Usually, polyclads occur in diverse marine habitats, such as under coastal stones, on reefs and in interstitial spaces (Hyman 1951;Prudhoe 1985;Curini-Galletti et al. 2008). About 800 to 1000 species of polyclads are currently recognised (Rawlinson 2008;Martín-Durán and Egger 2012).
The phylogenetic position of Polycladida within Platyhelminthes used to be very controversial . Only recently, Polycladida have been consistently recovered as sister group to Prorhynchida (a group harbouring only freshwater dwellers), forming the Amplimatricata, which is the sister group of all other Trepaxonemata (Egger et al. 2015;Laumer et al. 2015;Laumer and Giribet 2017). Lang (1884) was the first to distinguish between two groups of 'marine planarians', the Tricladida and the Polycladida. He further grouped the Polycladida into forms with a ventral sucker behind the genital openings (Cotylea), and those without (Acotylea). This classification system persists after some modifications (e.g. Laidlaw 1903;Bock 1913;Hyman 1953;Marcus and Marcus 1966) until today, and in the 1980s, Faubel (1983Faubel ( , 1984 and Prudhoe (1985) separately published monographs attempting to further clarify the Electronic supplementary material The online version of this article (https://doi.org/10.1007/s13127-019-00415-1) contains supplementary material, which is available to authorized users. interrelationships of polyclads on morphological grounds, using genital organs, especially the organisation of the prostatic vesicle (Faubel 1983(Faubel , 1984, the position of eyes and tentacles (Prudhoe 1985), or the pharynx organisation (Faubel 1983(Faubel , 1984Prudhoe 1985) as main systematic characters-however, the resulting classifications were largely incongruent. Interestingly, Faubel (1984) considered both Cotylea and Acotylea as not being monophyletic, but retained the names for taxonomic consistency. Prudhoe (1985) was also aware of problems with the classification and he cited several cases, where some families, such as Enantiidae and Boniniidae, have features fitting both to Cotylea and Acotylea.
For more than 30 years, these two conflicting systems have been in use by polycladologists (a term coined by J. Bahia, personal communication), stressing the need of a unifying system, based on morphology, on molecules, or both. The first molecular phylogenetic reconstruction of polyclad interrelationships was using a partial sequence of about 350 nucleotides of the marker molecule 28S (large nuclear ribosomal subunit) and was focussed on the family Pseudocerotidae, w i t h P e r i c e l i s a s t h e c o t y l e a n s i s t e r g r o u p o f Pseudocerotidae (Litvaitis and Newman 2001). Another molecular phylogenetic analysis of Polycladida based on partial 28S sequences (about 900 nt long) included just eight cotylean and six acotyleans-Cotylea was not recovered as monophyletic, since the cotylean species Pericelis cata appeared outside the other Cotylea as sister group of Acotylea, while Cestoplana rubrocincta emerged as an acotylean as in Faubel's and Prudhoe's systems (Rawlinson et al. 2011). With a very similar dataset, Rawlinson and Stella (2012) recovered both, Pericelis and Cestoplana, as basally branching cotyleans, thereby stressing the problematic position of these taxa. In a flatworm-wide phylogenetic study based on four genes, the acotylean Theama was grouped with the remaining Cotylea, not with the Acotylea (Laumer and Giribet 2014), which was corrobarated in a transcriptomic study in the following year (Laumer et al. 2015).
In 2017, three large molecular phylogenies of polyclads were published, two with different stretches of the 28S marker gene Tsunashima et al. 2017), and one with mitochondrial genes (Aguado et al. 2017). Of these studies, only Bahia et al. dealt with the aforementioned problematic taxa, namely Pericelis, Cestoplana and Theama-all of them showing up as cotyleans in their tree . However, this study only used a single alignment method and a single model with relatively low bootstrap support levels, so the reliability of the provided reconstruction remained unclear. During the review phase of this manuscript, another publication using the 28S marker gene was published (Litvaitis et al. 2019).
In the present study, we also have used partial 28S rDNA sequences, as well as a combined dataset of longer 18S and 28S sequences of a wide systematic range of polyclads. Most importantly, we have applied three different, widely used alignment algorithms and two different statistical approaches for tree reconstruction to test the stability and reliability of molecular phylogenies using one or two genes, and also, when possible, to infer relationships between groups based on a bigger data set.

Material and methods
Animal collection, identification of species and transcriptome data An overview of newly generated and published sequences is provided in Table 1. For most collected material, tissue was stored in 99% ethanol, and histological sections were made as described by Aguado et al. (2017) and Dittmann et al. (2019). Several published polyclad transcriptomes (Egger et al. 2015, Laumer et al. 2015 were searched for 18S and 28S sequences (see Table 1) using BLAST (Altschul et al. 1990).

DNA extraction, PCR amplification and sequencing
For all specimens, DNA was extracted from a small piece of ethanol-preserved marginal tissue. DNA extraction was performed following phenol-chloroform protocols (Sambrook et al. 1989;Chen et al. 2010). Concentration and possible contamination of extracted DNA were checked using NanoDrop (NanoDrop Fluorospectrometer Thermo Fisher Scientific, USA). PCRs were performed in a total volume of 25 μl or 50 μl. 18S rDNA was amplified either in two overlapping fragments using the published primer combinations 4fb + 1806R (ca. 1200 bp) and 5fk + S30 (ca. 900 bp) (Larsson and Jondelius 2008) or in one approach using 18S-1F + 18S9R (ca. 1800 bp) (Álvarez-Presas et al. 2008). 28S rDNA was amplified with the primers 28_LSU5_fw + L1642R (ca.1450 bp) or 28S_1F + 28S_6R (ca. 1600 bp) (Larsson and Jondelius 2008). PCR was performed using a 'Touch Down' protocol using the following protocol: 5 min of initial denaturation at 94°C; 30 s of denaturation at 94°C, annealing at 68-45°C for 30 s, extension at 72°C for 2 min; 12 cycles; 30 s of denaturation at 94°C, annealing at 45°C for 30 s, extension at 72°C for 2 min; 23 cycles; final extension at 72°C for 10 min, hold at 4°C. Successful products were purified using ExoSAP-IT (Affymetrix, USA), following manufacturer's protocol, or with the Wizard® SV gel and PCR clean-up system (Promega, USA) according to the manufacturer's quick protocol. PCR products were sequenced by CBMSO (Spain) or by Microsynth Austria GmbH, respectively. Sequences were assembled and edited by hand or using the software CLC Main Workbench 7 (Qiagen, Germany). Table 1 List of all species used in this study, including authorities, sample locations and accession or SRA numbers. In trees using only a single sequence of the same species, the first listed sequence was included, with the number omitted

Datasets for phylogenetic analyses
We made eight different single gene sequence collections of 'short' 28S sequences (see Table 1 for accession numbers of all newly generated and used published data). In general, we only used one sequence per species from the same authors. The first sequence collection used 108 polyclad terminals (including the first, gappy version of sequences published by Bahia et al. 2017 on NCBI, which was corrected and reuploaded by Bahia et al. in 2019 with a non-gappy version), 20 of which were generated by us, and Macrostomum lignano as an outgroup ('28Sshort1'), while all subsequent 'short' 28S sequence collections worked with the updated second sequence versions of Bahia et al. (2017): '28Sshort2' added Cycloporus japonicus, two Pericelis and four pseudocerotoid sequences, while '28Sshort3' only included all (updated) sequences of '28Sshort1'.
Sequences for each gene were separately aligned using three methods: MUSCLE v3.8.31 (Edgar 2004), MAFFT Q-INS-i and MAFFT E-INS-i v7.310 (Katoh and Standley 2013). They were manually trimmed, and in the case of the combined dataset, concatenated. For several alignments, we also used Gblocks with the least stringent settings (Castresana 2000). Conversion of fasta alignments to Nexus and Phylip formats was done using ALTER (Glez-Peña et al. 2010).
For ML trees, between 500 and 10,000 tree searches were performed, and between 500 and 1000 separate bootstrap replicates. At least 5-10 million generations were calculated for BI trees, or more until convergence (average standard deviation of split frequences < 0.01) was reached. For extended majority-rule consensus trees, we used RAxML with the concatenated trees of BI and ML analyses of the 28Sshort6 dataset (see Table 2). Phylogenetic trees were visualised in Figtree (http://tree.bio.ed.ac.uk/software/figtree/) and adapted in Inkscape (https://inkscape.org/) and Adobe Illustrator CS4.
The sequences generated during and/or analysed during the current study are available in the GenBank repository, under the accession numbers listed in Table 1. The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request.

Results
Effects of model choice, alignment, outgroup selection and taxon sampling on tree topology We recovered varied results with our combined 18S28Slong matrix (see Table 2, Suppl. Figs. S1-12): without using Gblocks for alignment curation, three of the six phylogenetic reconstructions supported Cestoplanoidea, Chromoplanoidea, Periceloidea, Anonymidae and Chromoplana as cotyleans (Suppl. Figs. S5, S10, S11), while three trees rendered most of these taxa as acotyleans or polytomic (Suppl. Figs. S4, S6, S12). We have visualised these changes caused by model choice in Fig. 1. Using Gblocks, only the two trees based on a Q-INS-i alignment recovered Cestoplanoidea, Chromoplanoidea, Periceloidea, Anonymidae and Chromoplana as cotyleans (Suppl. Figs. S2,8). In both E-INS-i trees and the ML MUSCLE tree, Anonymidae and Chromoplana are sister group of all other Polycladida, while in the BI MUSCLE tree, they are polytomic with Cotylea and Acotylea (see Table 2).
We continued our analyses with the first 28S-only dataset (28Sshort1) with many more taxa than available in the 18S28Slong dataset, including the first version of sequences published by Bahia et al. (2017). With this dataset, we calculated BI and ML trees based on three different alignments, and consiste ntly (100%) recovered Cestoplanoidea, Chromoplanoidea, Periceloidea, Anonymidae and Chromoplana as Cotylea. The corresponding MRE tree exhibited an identical topology as the BI MUSCLE tree shown in Fig. 2a. After obtaining the new sequence versions of Bahia et al. (2017) in January 2019, we recalculated all trees with the new sequences (and adding additional sequences, see 28Sshort2) and consistently (100%) recovered the aforementioned groups as Acotylea, regardless of outgroup selection or alignment curation (Fig. 3).
We tested different parameters, always using a short 28S dataset with BI MUSCLE with and without Gblocks for tree reconstruction. Using Gblocks, outgroup selection markedly c h a n g e d o t h e r p a r t s o f t h e t o p o l o g y, s u c h a s Prosthiostomoidea alternating between Acotylea and Cotylea ( Fig. 3). With only Xenoprorhynchus as outgroup, Chromoplana is the sister group of all other Polycladida (Fig. 3b), while with only Macrostomum as outgroup, Cycloporus variegatus takes the place of sister group of all other Polycladida (Fig. 3c). Using the same ingroup and outgroup taxa as in Fig. 3c, but without Gblocks, we recovered a topology with many basal polytomies (Fig. 3d). With both non-polyclad outgroups, a basal polytomy between Cycloporus variegatus, Euryleptidae and all other Polycladida was recovered (Fig. 3a). Prosthiostomidae are basally branching Acotylea with Macrostomum + Xenoprorhynchus, and only Macrostomum as outgroups. Xenoprorhynchus alone as outgroup provides a basal polytomy of Anonymus, Cotylea and Acotylea, except Chromoplana (Fig. 3b).
Consequently, we tested if the newly added sequences were responsible for the change in tree topology, especially of Cestoplanoidea, Chromoplanoidea, Periceloidea, Anonymidae and Chromoplana. We therefore removed all additional sequences compared to our first dataset leading to the 28Sshort3 alignment, and with the same alignment and 7. Chromoplana and Anonymus recover as clade 1 x 8. Clade 1 appears as sister group to a clade including Prosthiostomoidea and Pseudocerotoidea 16. Prosthiostomoidea appears as monophyletic 17. Prosthiostomoidea is sister group to Pseudocerotoidea      Now we also removed Chromoplana from the dataset (28Sshort4) and once more had C. variegatus as sister group to all other Polycladida. Also, Cestoplanoidea, Chromoplanoidea, Periceloidea and Anonymidae emerged as Acotylea (Fig. 4a). With the additional removal of C. variegatus from the sequence collection (28Sshort5), we recovered Cestoplanoidea, Chromoplanoidea, Periceloidea and Anonymidae as Cotylea once more-but only after alignment curation with Gblocks (Fig. 4b).
In a last change, we returned to the full dataset with updated sequences, but also used all available sequence variations for Cycloporus variegatus and Cycloporus gabriellae, instead of only using one sequence per species from the same authors (28Sshort6, Suppl. Figs. S13-24). We now recovered Cycloporus again within Cotylea, and present the detailed results using this dataset in the following section.
Comparative tree topology using 28Sshort6 and 18S28Slong matrices All 12 trees using the 28Sshort6 matrix (Suppl. Figs. S13-24), and most of the 12 trees using the 18S28Slong matrix (Suppl. Figs. S1-12) are different from each other. We have analysed the tree topologies to identify stable and unstable taxa ( Table 2). This table also gives an overview of which tree supports which topology. Additionally, we computed extended majority-rule consensus (MRE) trees from all 12 trees of the 18S28Slong matrix (Fig. 5), and all 12 trees of the 28Sshort6 matrix (Fig. 6). We also calculated separate 28Sshort6 and 18S28Slong matrix-based MREs for all alignments treated with or without Gblocks, respectively (Suppl. Figs. S25-28). If not otherwise stated, the MRE tree always refers to the MRE calculated from all twelve trees of each matrix. 'Trees' refers to both BI and ML trees, unless it is preceded by 'MRE'. Instead of citing all trees supporting a particular placement of a taxonomic group, we provide this information in Table 2 for better accessibility and overview.
In the following, we focus our comparisons on already defined families and superfamilies, mainly of the systems established by Faubel (1983Faubel ( , 1984 and Bahia et al. (2017).
The majority (63%) of our trees, and the 28Sshort6 MRE tree support Cotylea and Acotylea sensu Bahia et al. (2017) and in the following we use these terms according to their definition: in brief, Theama and Cestoplana are cotyleans instead of acotyleans.
Cestoplanoidea  and thereby its only family, Cestoplanidae, appear monophyletic in all our trees, even if its position within the trees differs widely. The majority (63%) of our trees, and the 28Sshort6 MRE tree, support Cestoplanoidea within (and as sister group to all other) Cotylea, but in the remaining trees, it is sister group to Acotylea (33%) or, in one case, polytomic. Also 63% of our trees, and the 28Sshort6 MRE tree, support the phylogenetic position of Chromoplanoidea (Bahia et al. 2017, including Theama, Chromyella and Boninia) within Cotylea. Only 29% of our trees (all of them 28Sshort6 trees), as well as the 28Sshort6 MRE tree, place Chromoplanoidea as sister group to all other Cotylea, except Cestoplanoidea. In 88% of our trees, and in both MRE trees, Theamatidae is sister group to a clade consisting of Boninia and Chromyella.
Periceloidea  and thereby its only family, Pericelidae, is also monophyletic in all of our phylogenetic reconstructions and both MRE trees. They are most often either sister group to all Cotylea except Cestoplanoidea and Chromoplanoidea (25% and the 28Sshort6 MRE tree), or sister group to Chromoplanoidea within Acotylea (25% of all trees, but 100% of the 18S28Slong Gblocks trees, and the 18S28Slong MRE tree). However, in 21% of the trees, Periceloidea is placed as sister group to all Cotylea except Cestoplanoidea, or, also in 21% of the trees, Periceloidea is sister group to all Acotylea and Cestoplanoidea.
All but one of our trees, and both MRE trees recover Chromoplana and Anonymus as clade 1 and this clade mostly (58% and both MRE trees) appears as sister group to a clade including Prosthiostomoidea (with the single family Prosthiostomidae) and Pseudocerotoidea (consisting of Pseudocerotidae, Euryleptidae and clade 2, see paragraph below).
Pseudocerotoidea and Pseudocerotidae sensu Bahia et al. 2017 are monophyletic in all but two trees, and in both MRE trees. Within Pseudocerotidae, all of our 28Sshort6 trees show that neither Pseudoceros, nor Pseudobiceros, nor Thysanozoon are monophyletic. The traditional family Euryleptidae sensu Faubel 1984 does not appear monophyletic in any of our trees, including the MRE trees. It is split into two clades (21 trees) or three clades (3 trees). In this work, we termed one of these clades 'clade 2' (while retaining the name Euryleptidae for the larger clade). The larger clade includes Cycloporus japonicus, Cycloporus variegatus, Prostheceraeus and Maritigrella in the 28Sshort6 trees, while Cycloporus is lacking in the 18S28Slong trees. Clade 2 consists of Euryleptodes galikias, Cycloporus gabriellae and Stylostomum ellipse in the 28Sshort6 trees, and only Stylostomum ellipse in the 18S28Slong trees. In the three trees, where the Euryleptidae sensu Faubel 1984 are split into three clades, even the clade still called Euryleptidae is recovered as paraphyletic. The genera Cycloporus and Prostheceraeus are never recovered as monophyletic in any of the 28Sshort6 trees, and also Maritigrella is only recovered as monophyletic in one third of the 28Sshort6 trees.
Prosthiostomoidea ) appears monophyletic in all trees, and in all but two trees as sister group to Pseudocerotoidea. Within Prosthiostomoidea, Enchiridium appears as sister group to a clade consisting of Prosthiostomum, Lurymare and Amakusaplana in all trees. Prosthiostomum appears polyphyletic in 83% of our 28Sshort trees and also the 28Sshort MRE tree, as Amakusaplana and Lurymare cluster within.
Leptoplanoidea sensu Faubel 1983 (in whose definition Hoploplana and Theama are included) is not supported in any of our trees. We have termed Leptoplanoidea sensu Faubel 1983, but without Hoploplana and Theama, clade 3 (supported by all trees), which we further subdivided into two clades (clades 5 and 6).
In all 18S28Slong trees, clade 5 is synonymous with Leptoplanoidea sensu Bahia et al. (2017), as the genus Pseudostylochus is not available in these datasets. In all 28Sshort6 trees, Pseudostylochus is part of clade 5, but Pseudostylochus is not included in the superfamily's definition given by Bahia et al. (2017). All of our 28Sshort6 trees show that the genus Leptoplana is monophyletic, and that Notoplanidae as a whole and also its genera Notoplana and Notocomplana are not monophyletic.
In 58% of the 28Sshort6 trees and also the corresponding MRE tree, clade 6 is monophyletic, appears as sister group to clade 5 and includes Discocelis, Adenoplana, Ilyella, Phaenocelis and Amemiyaia. In our 18S28Slong trees and the corresponding MRE tree, clade 6 is represented only by Discocelis and always recovered as the sister group of clade 5.
Clade 4 can be subdivided into two clades, clades 7 + 8 and their sister group Callioplana. This topology is supported by two of twelve 28Sshort6 trees, as well as the respective MRE tree, while in ten 28Sshort6 trees, either a polytomy exists          In half of the 28Sshort6 trees, and also in the corresponding MRE tree, clade 7 is formed by Hoploplana clustering as sister group to Idioplana, but in one third of 28Sshort6 trees, Hoploplana is sister group to Planoceridae in clade 8. In our 18S28Slong trees, clade 7 is only represented by Hoploplana californica, in nine of twelve trees (and also in the 18S28Sshort MRE tree) as sister group to Planocera pellucida. In two cases, Planocera pellucida is clustering within clade 8, while in one case, it is unresolved as a polytomy.
Clade 8 resembles Stylochoidea sensu Faubel (1983) and in the 28Sshort6 dataset, includes Leptostylochus, Stylochus, Imogine, Paraplanocera and Planocera, but excludes Faubel's Callioplana, Idioplana and Pseudostylochus. This clade is monophyletic in a minority (two of twelve) 28Sshort6 trees and in the 28Sshort6 MRE tree. Clade 8 is polytomic in seven 28Sshort6 reconstructions and paraphyletic (including Idioplana and/or Hoploplana) in the remaining three 28Sshort6 trees. In the 18S28Slong dataset, clade 8 is represented by Stylochus, Paraplanocera, Imogine and Planocera and thus conforming to Faubel's (1983) definition. Clade 8 is supported by only two of twelve of the 18S28Slong trees as well, but not in the 18S28Slong MRE tree, where Hoploplana is sister group to Planocera.
In 75% of our 28Sshort6 trees and also the corresponding MRE tree, Planocera is not monophyletic as Paraplanocera      Faubel 1983 and1984 are written in blue and red fonts, respectively. Species recovered as Acotylea or Cotylea in our trees are displayed with blue and red background, respectively. Branches and nodes are given the same colour as their respective taxon sp. clusters within. Similarly, Paraplanocera is not monophyletic in any of our 28Sshort6 trees, and Planoceridae sensu Faubel 1983 is not recovered as monophyletic in any tree. In the majority (75%) of all trees, as well as in both MRE trees, the genus Stylochus is not monophyletic and in all 28Sshort6 trees, Imogine is not monophyletic.

Discussion
Taxon sampling and outgroup selection, as well as the choice of marker genes, the alignment method and the analysing statistical models affect the resulting phylogenetic reconstructions significantly (see e.g. Lockyer et al. 2003;Puslednik and Serb 2008;Aguado and Bleidorn 2010;Laumer and Giribet 2017). For polyclad interrelationships using mainly a rather short stretch of the 28S rDNA marker gene, but also a longer sequence comprised of both partial 18S and 28S rDNA, we show that the change of any of these parameters can vastly change the resulting tree topology (Figs. 1, 2, 3 and 4). A strong hypothesis about valid polyclad interrelationships is thus challenging, and we have therefore used majority-rule consensus trees to help us decide between different topologies (Figs. 5 and 6) and also manually analysed the support of different hypotheses (Table 2). To our knowledge, this is the first time that these difficulties and inconsistencies are discussed or even mentioned in regard to polyclad interrelationships.

Alignment is important
MUSCLE (Edgar 2004) was the alignment method of choice in both recently published polyclad phylogenies based on partial 28S sequences Tsunashima et al. 2017), and was also used in one of the best-scoring trees in both datasets shown here (Table 2, Suppl. Figs. S10, 13). We have also used two different variants of MAFFT (Katoh and Standley 2013); previously, MAFFT E-INS-i was selected for the polyclad phylogeny based on mitochondrial sequences (Aguado et al. 2017) and for an all-flatworm phylogeny working with the nearly complete nuclear ribosomal marker genes, 18S and 28S (Laumer and Giribet 2017). The other best-scoring 28Sshort6 tree according to our scoring in Table 2 is MAFFT E-INS-i aligned (Suppl. Figs. S15), and another MAFFT E-INS-i tree (Suppl. Fig. S18) is also closest to the topology shown in the MRE 28Sshort6 tree (Fig. 6). MAFFT Q-INS-i is by far the most computationally demanding alignment method, and was also employed quite extensively for resolving flatworm interrelationships on the level of orders based on partial 18S and 28S, e.g. macrostomorphs (Janssen et al. 2015), rhabdocoels (van Steenkiste et al. 2013;Tessens et al. 2014) and proseriates (Casu et al. 2014;Scarpa et al. 2015Scarpa et al. , 2016Scarpa et al. , 2017. The two best-scoring 18S28Slong trees are both based on a MAFFT Q-INS-i alignment (    The tree topologies resulting from these widely used alignment methods are not consistent (Fig. 1, Suppl. Figs. S1-24), corroborating the findings of Laumer and Giribet (2017), in which they reanalysed their differently aligned dataset from their earlier publication (Laumer and Giribet 2014) and also recovered trees with several major differences. In their re-analysis, they used MAFFT E-INS-i instead of RNAsalsa (Stocsits et al. 2009), and then recovered a tree very similar to two independently made transcriptomic analyses of flatworm interrelationships (Egger et al. 2015, Laumer et al. 2015, suggesting that MAFFT E-INSi provided a more robust alignment than RNAsalsa.
From this work, we cannot give an unambigious recommendation for the most suitable alignment method, but recommend to use at least two different methods to check for consistency.

Model choice is important
In the work presented here, we consistently recovered inconsistent BI and ML topologies using the same datasets and alignments (Table 2, Fig. 1). In the most recently published polyclad phylogeny, both BI and ML trees gave congruent results (Litvaitis et al. 2019). In other recent polyclad phylogenies based on partial 28S, only either BI (Rawlinson and Stella 2012) or ML  were used, so no comparisons between different models can be made. In two polyclad phylogenies, both BI and ML analyses were run, and the trees show the same   Faubel 1983 and1984 are written in blue and red fonts, respectively. Species recovered as Acotylea or Cotylea in our trees are displayed with blue and red background, respectively. Branches and nodes are given the same colour as their respective taxon topology in Rawlinson et al. (2011) and are 'highly congruent' in a mitochondrial gene analysis with several switches within families, but not of the overall topology (Aguado et al. 2017). Both models were used to resolve interrelationships within other flatworm orders, and reported with very similar or identical results using combined 18S and 28S datasets (Casu et al. 2014;Tessens et al. 2014;Janssen et al. 2015;Scarpa et al. 2015Scarpa et al. , 2016Scarpa et al. , 2017, in one case also including mitochondrial markers (Janssen et al. 2015). While these studies usually use a matrix of more than 3000 nt, our own large matrix with more than 3000 nt positions gives less congruent results among different models and alignments than our short matrix (ca. 800 nt) (Table 2, Figs. 5 and 6), indicating that taxon sampling may be even more important than matrix length.
Again, we recommend to use both models (BI and ML) to check for consistency between the models. In our case, results were not consistent, indicating that taxon sampling and matrix length were not sufficient yet.

Outgroup selection is important
We have tested the influence of outgroup choice on tree topology with Macrostomum lignano, a basally branching rhabditophoran, and Xenoprorhynchus sp., a basally branching prorhynchid-Prorhynchida being sister group of Polycladida (Egger et al. 2015, Laumer et al. 2015, using the same alignment (MUSCLE) and alignment curation (Gblocks), as well as the same model (BI) and the same dataset (28Sshort2). We found markedly different tree topolo g i e s b e t w e e n u s i n g b o t h M a c ro s t o m u m a n d Xenoprorhynchus, only Xenoprorhynchus or only Macrostomum as ougroups (Fig. 3a-c). Especially the sister group relationships of either Chromoplana sp. or Cycloporus variegatus with all other polyclads (Fig. 3b, c) were the reason to also test the influence of taxon sampling on the polyclad tree topology (Fig. 4).
An almost identical dataset, aligned with the same algorithm and tree reconstruction done with the same model and by the same leading author yielded two different topologies: in the first account, both Cestoplana and Pericelis are basally branching Acotylea (Rawlinson et al. 2011), while these two taxa switch to basally branching Cotylea in the second account (Rawlinson and Stella 2012). The only two differences in the reconstructions are one instead of two outgroups and a third sequence of Amakusaplana acroporae in the second paper (Rawlinson and Stella 2012), indicating that a higher number of outgroups gives more reliable results in their case. In our own datasets, we found no clear preference for outgroup selection (Fig. 3), making us default on a single, basally branching outgroup (Macrostomum lignano) for our main datasets (28Sshort6 and 18S28Slong).

Taxon sampling is important
Not only the long-branching Chromoplana (therefore excluded from the analysis in Bahia et al. 2017), but also Cycloporus variegatus was prone to upend the tree topology in the 28S trees (Figs. 3b, c and 4). Interestingly, both the complete removal of Chromoplana and all Cycloporus sequences, and the addition of more variants of Cycloporus species yielded similar tree topologies (Figs. 4b and 6). We have not tested removing taxa from the 18S28Slong dataset, but at least in theory, it should be more robust to taxon sampling artefacts than the much shorter 28Sshort dataset. In general, and as stated above, taxon sampling seems to be more important for resolving a stable polyclad phylogeny than matrix length at this point.

Correct determination is important
The correct identification of species is far-reaching for the interpretation of phylogenetic trees. During our analysis, we realised several inconsistencies in species determination of so far published sequences. In several of our 28Sshort6 trees, as well as in the corresponding MRE (Fig. 6), a sequence tagged as Paraplanocera sp. (KY263699.2) on GenBank clusters within Planocera. Therefore, Planocera does not appear monophyletic (Table 2, Fig. 6). However, according to Bahia et al. (2017), this sequence and the associated accession number belongs to Planocera sp.; hence, Planocera would be monophyletic also in our trees. We found several similar problems with sequences listed as 'Leptoplana sp. or Notoplana sp.' in Table 1 of Bahia et al. (2017). In their table, these sequences have the accession numbers KY263695, KY263650, KY262696, KY263698 and KY263651. KY262696 is apparently a typo and should read KY263696, which together with KY263695 and KY263698 is tagged as 'Leptoplana tremellaris' on GenBank, while KY263650 and KY263651 are labelled as 'Notoplana sp.' on GenBank. In their tree, Bahia et al. (2017) also show an unlisted Notocomplana sp., but it is not clear to which accession number this species refers to. As usual, we only took one sequence of the same species from the same authors, and we have used KY263695 (Leptoplana tremellaris) and KY263650 (Notoplana sp.) in our phylogenetic reconstructions (Table 1). Interestingly, in our 28Sshort6 MRE tree (Fig. 6), this Notoplana sp. by Bahia et al. (2017) does not cluster with any other Notoplana, Notocomplana (Notoplanidae) or Leptoplana (Leptoplanidae) species, but with Melloplana (Pleioplanidae) and Echinoplana (Gnesiocerotidae).
Also Pseudoceros is not monophyletic in our analyses, as two species, Pseudoceros harrisi and Pseudoceros cf. maximus are clustering outside the other 13 included Pseudoceros species (Fig. 6). Pseudoceros harrisi is consistently recovered as sister group to all other Pseudocerotidae in our trees and also by Bahia et al. (2017) and Tsunashima et al. (2017). In its species description, which is based on a single damaged specimen, it is stated that 'This species does not resemble any other species of Pseudoceros. However, P. harrisi may be confused with members of Cycloporus [...]' (Bolaños et al. 2007). Hence, the phylogenetic position of Pseudoceros harrisi might be the result of a misdetermination of its genus in the original description. The Pseudoceros cf maximus sequence (KY263708) we used was published by Bahia et al. (2017) and it appears with high support within Pseudobiceros in our reconstructions (Fig. 6). We noticed that the species name Pseudoceros cf maximus does not appear in Bahia et al.'s tree. On the other hand, they show two branches labelled 'Pseudobiceros spp.' in their tree, but only list a single Pseudobiceros sp. sequence in their Table 1. Taking into account our own results, we believe it is possible that the sequence published as Pseudoceros cf maximus on GenBank is one of the 'Pseudobiceros sp.' in their tree.
Several sequences have undergone name changes after redetermination efforts by the authors, or have dubious affiliations. For example, Cestoplana rubrocincta from Australia (C. rubrocincta 2 in our tree, HQ659009.1) is labelled as C. australis in the tree provided by Rawlinson et al. (2011), but called C. rubrocincta in their table, and also on GenBank. Other sequence names were updated without changing their accession number versions. We originally downloaded the following sequences published in Tsunashima et al. (2017) from GenBank in June 2017, but

Sequence problems
When we started with this study in 2017, we noticed gaps in all newly generated sequences uploaded to GenBank by Bahia et al. (2017). The first set of 28Sshort trees we made was based on a dataset including these sequences. We later realised that the gaps in the sequences were caused by alignment curation using Gblocks (J. Bahia, pers. comm.), and all other trees (using the 28Sshort2-6 sequence collections) were based on the updated sequences (version 2 on GenBank). We provided reconstructions based on both, Gblocks curated and original alignments, and often recovered different topologies if all other parameters stayed the same (Table 2, Fig. 3). According to a recent publication, phylogeny may be even better without using Gblocks or similar alignment curation programs (Tan et al. 2015). In our own study, however, we find that the bestscoring trees were made with datasets using Gblocks for alignment curation ( Table 2).
Some of the sequences published by Tsunashima et al. (2017) appear to be quite different to all other polyclad sequences published, especially in the 5′ region: among these are the above-mentioned Cycloporus japonicus (LC100092), Thysanozoon brocchii (LC100093) and Thysanozoon japonicum (LC100094). We initially removed all of these sequences from further analyses, but later added Cycloporus japonicus (28Sshort2 and 28Sshort6) despite the divergent sequence. Also Chromoplana sp. from Laumer and Giribet (2014) was an unusual sequence and was therefore removed from the tree of Bahia et al. (2017), but is included in most of our reconstructions (except 28Sshort4-5).
Although termed as 'clones' on GenBank, there is a considerable difference between the four published Cycloporus variegatus sequences by Bahia et al. (2017); we believe these sequences are not derived from clones, but from different specimens of the same species.
Polyclad phylogenies based on partial 28S rDNA published by different authors used different primers, making the integration of all sequences a challenge, as the overlapping regions get smaller. Especially Tsunashima et al. (2017) used a region of the 28S gene more towards the 3′ end than all other studies, but we have still included most of their sequences, because they provide important taxa not covered by our own or other previously published sequences. For future studies, we recommend amplifying 28S starting with expansion segment D1 and stretching as long as possible, to maximise compatibility with published sequences.

Classification on suborder and superfamily level
On suborder level, our 28Sshort6 trees are mostly compatible with the molecular phylogenetic hypothesis of Bahia et al. (2017), supporting their redefinition of Cotylea and Acotylea (see Table 2 and Fig. 6). There, two traditional actoylean families, Cestoplanidae and Theamatidae, switched from Acotylea to Cotylea.
The majority of the 18S28Slong trees, on the other hand, support Cestoplanidae and Theamatidae as acotyleans. Also, the traditionally cotylean genera Pericelis, Boninia and Chromyella are recovered as acotyleans (Table 2, Fig. 5). In this scenario, a sucker would be a character at the base of Polycladida and would have been lost at least five times: in the traditional Acotylea, in some Cestoplanidae, in the anonymid Simpliciplana marginata (Kaburaki 1923), in Theamatidae, in Amakusaplana (Rawlinson et al. 2011), and possibly in Chromyella (Fig. 5 and Faubel 1983, 1984, Prudhoe 1985. In the 28Sshort6 scenario, a sucker would be present at the base of Cotylea and would have been lost one time less, i.e. not in Acotylea (Fig. 6). According to Bahia et al. (2017), a 'true sucker' may have gradually evolved and may be an apomorphy of Prosthiostomoidea and Pseudocerotoidea. A true sucker is muscular and characterised by a modified epithelium with a thin basement membrane, while the adhesive disc or pad found in Boninia and Cestoplana is just a shallow depression of the epithelium not differentiated from the parenchyma (Prudhoe 1985;Rawlinson and Litvaitis 2008). Both true sucker and adhesive disc/pad are always located posterior of the genital openings. Several Pericelis species (excluded from having a true sucker in Bahia et al. 2017, but listed as having a true sucker in Rawlinson and Litvaitis 2008) are described with a 'distinct sucker' (Dittmann et al. 2019), so we suggest that the true sucker behind the genital openings already is an apomorphy for the unnamed group including Periceloidea, Anonymus, Chromoplana, Prosthiostomoidea and Pseudocerotoidea (Fig. 6). The acotylean genus Leptoplana has a sucker (a socalled genital pit) between the genital openings (Prudhoe 1985); therefore, it is excluded from the definition of a cotylean sucker.
Based on this scenario of sucker evolution in polyclads, it is more parsimonious to support the 28Sshort6 tree topology, although the 18S28Slong alignment with ca. 3000 nt is almost four times as long as the 28Sshort6 alignment with ca. 900 nt. Also, the support values of the trees rejecting Cotylea and Acotylea sensu Bahia et al. (2017) are consistently lower than those supporting them (Suppl. Figs. S1-24). In five of the twelve 18S28Slong trees, Cotylea and Acotylea sensu Bahia et al. (2017) are actually supported, and also in the 18S28Slong MRE tree without Gblocks (Suppl. Fig. S26). Only the 18S28Slong dataset using Gblocks skews the picture towards a weakly supported topology making Cestoplanidae, Theamatidae, Pericelis, Boninia and Chromyella acotyleans (Suppl. Fig. S25), also in the combined 18S28Slong MRE tree (Fig. 5).
In all but two 28Sshort6 trees, Cotylea and Acotylea sensu Bahia et al. (2017) are well supported (Fig. 6,(18)(19)(20)(21)(22)24). On the other hand, we have shown that this topology is very much dependant on taxon sampling, outgroup selection, alignment method and curation, and model choice (Figs. 1, 2, 3 and 4). Possibly, the most important parameter is taxon sampling, and this would explain why a much larger alignment (18S28Slong) with 27 polyclad terminals and 26 different polyclad species gives less consistent results than the shorter matrix (28Short6) with 118 polyclad terminals and 100 different polyclad species. Bahia et al. (2017) show 136 polyclad terminals, but only 55 different polyclad species, and Tsunashima et al. (2017) use 53 polyclad terminals and 50 polyclad species in their phylogenetic trees. While we have not tested their original datasets with different parameters here, their results suggest that neither the number of taxa, nor sequences are decisive for tree topology, but that some sequences are prone to change tree topology, among them Chromoplana, Cycloporus variegatus and Cycloporus japonicus (Figs. 3 and 4). As long as single taxa included or excluded can drastically change tree topology even in the overall more consistent 28S-only trees, polyclad phylogeny remains only preliminarily resolved, calling for larger datasets like in transcriptomic phylogenies.
However, apart from the position of Cestoplanidae, Theamatidae, Pericelis, Boninia, Chromyella, Anonymidae and Chromoplana in the tree, we find that most polyclad taxa are included in very well-supported clades.
Our data support the following new superfamilies sensu Bahia et al. (2017): Pseudocerotoidea sensu Bahia et al. (2017); this superfamily includes Pseudocerotidae and two clades of Euryleptidae in their reconstruction. In this work, we termed one of these clades 'clade 2' as all relevant trees show this non-monophyly (Table 2). This division can also be observed in the study of Bahia et al. (2017), where Cycloporus gabriellae represents our clade 2 of Euryleptidae, while Cycloporus variegatus and Cycloporus japonicus are part of the remaining Euryleptidae. Also in a cladistic analysis, Euryleptidae was not recovered as monophyletic (Rawlinson and Litvaitis 2008). As already suggested before, the genus Cycloporus needs to be revised, but no obvious characters to distinguish between described Cycloporus species could be determined so far . Our data show that the separation of the Cycloporus species not only results from potential inconsistencies within the genus Cycloporus, as also Stylostomum and Euryleptodes appear within clade 2. Therefore, we propose the revision of the whole family of Euryleptidae. As Eurylepta has been shown to cluster as sister group of other Euryleptidae in a phylogeny based on mitochondrial genes (Aguado et al. 2017), the family name Euryleptidae should be retained for the group containing Maritigrella, Prostheceraeus, Cycloporus variegatus and Cycloporus japonicus (Fig. 7). Cycloporus japonicus has been shown to group with Maritigrella in Tsunashima et al. (2017) as well. We propose the new family name Stylostomidae fam. nov. for clade 2, including at least Stylostomum, Euryleptodes and Cycloporus gabriellae. In the recently published work by Litvaitis et al. (2019), both Euryleptidae and Cycloporus appear as monophyletic, but neither Stylostomum, nor Euryleptodes are included in their study. As in our study, Litvaitis et al. (2019) have recovered both, Prostheceraeus and Maritigrella, as non-monophyletic and consequently they have synonymised Maritigrella as junior synonym with Prostheceraeus.
Pseudoceros, Pseudobiceros and Thysanozoon are not recovered as monophyletic in our study, agreeing with Bahia et al. (2017) and Tsunashima et al. (2017), stressing the need of further revision of the family Pseudocerotidae (Litvaitis and Newman 2001;Rawlinson and Litvaitis 2008).
Prosthiostomoidea was erected by Bahia et al. (2017) and only contains a single family, Prosthiostomidae. All our data support the monophyly of this family/superfamily and most data ( Table 2) also their sister group relationship to Pseudocerotoidea as described by Bahia et al. (2017). Similar to our results, also in the study of Tsunashima et al. (2017), Prosthiostomum is not monophyletic, as Amakusaplana (and in our case, also Lurymare) clusters within. In Aguado et al. (2017), two different species of Lurymare do not form an adelphotaxon. The fact that Amakusaplana clusters within Prosthiostomum is not very surprising, as Faubel (1984) remarks that the genus Amakusaplana has to be eliminated, as it is too similar to Prosthiostomum. The genus Amakusaplana is distinguished from Prosthiostomum mainly by body shape and the arrangement of eyes (Kato 1938;Faubel 1984), and also by the absence of the ventral sucker (Kato 1938;Rawlinson et al. 2011). Only in two of twelve 28Sshort6 reconstructions is Prosthiostomum monophyletic (Suppl. Figs. S17, 23), and Litvaitis et al. (2019) synonymise Amakusaplana with Prosthiostomum. Our data support this decision. The position of Lurymare within Prosthiostomum was already assumed by Poulter (1975). He proposed a subdivision of the genus Prosthiostomum into the subgenera P. (Lurymare) and P. (Prosthiostomum), distinguishable by the constitution of the prostatic vesicle (Poulter 1975). Faubel (1984) remarks that this definition also includes Enchiridium and elevates both subgenera back as genera. At least Enchiridium may be monophyletic, as suggested by Bahia et al. (2017), Litvaitis et al. (2019) and our own trees. Together, the molecular phylogenies do not support any of the previously proposed genera (Kato 1938;Poulter 1975;Faubel 1984) except Enchiridium, i.e. the revision of the genera Prosthiostomum and Lurymare is required.
Our clade 1, consisting of Anonymus and Chromoplana, is extremely well supported and always monophyletic, except in one case, where it appears polytomic (Suppl. Fig. S21). We propose a new superfamily Anonymoidea superfam. nov.   Fig. 6, but with newly defined and named groups indicated. Acotylea and Cotylea sensu Faubel 1983 and 1984 are written in blue and red fonts, respectively. Species recovered as Acotylea or Cotylea in our trees are displayed with blue and red background, respectively. Branches and nodes are given the same colour as their respective taxon (Fig. 7), including the families Anonymidae (Lang 1884) and Chromoplanidae (Bock 1922).
Cestoplanoidea was defined by Poche (1926), emended by Prudhoe (1985) and supported by Bahia et al. (2017) and Litvaitis et al. (2019); a majority of the 28Sshort6, but only a minority of the 18S28Slong analyses support its sister group relationship with all other Cotylea, as suggested by Rawlinson and Stella (2012) and Bahia et al. (2017), even if the family was originally assigned to Acotylea (Lang 1884;Faubel 1983;Prudhoe 1985) and appears as an acotylean in a majority of the 18S28Slong analyses, as well as in Rawlinson et al. (2011). In his monograph, Faubel remarks that organisation features of Cestoplana, like the forward direction of the male complex, the multiplication of the female apparatus in Cestoplana polypora, or the presence of an adhesive organ in some (but not all) Cestoplana species, could imply that Cestoplanidae have possibly arisen from a cotylean ancestor (Laidlaw 1903;Bock 1922;Faubel 1983;Bahia et al. 2017).
In all of our reconstructions, Cestoplanoidea are monophyletic (Table 2), although the only representing genus is Cestoplana.
Periceloidea was also erected by Bahia et al. (2017) and also contains a single, monotypic family, Pericelidae. Our data support the monophyly of this group. Additionally, our 28Sshort6 MRE tree (Fig. 6) supports its sister group relationship with all remaining Cotylea except Cestoplanidae, as already assumed by Bahia et al. (2017) and Rawlinson and Stella (2012). In Tsunashima et al. (2017), Pericelis is also recovered as a cotylean, but as sister group to Boninia + Chromyella (Theamatidae), although Cestoplana is absent in Tsunashima et al.'s reconstruction. In Rawlinson et al. (2011) and our 18S28Slong MRE tree (Fig. 5), Periceloidea are grouping with Acotylea, however. Litvaitis et al. (2019) include Diposthus in their phylogenetic reconstruction, which emerges as sister group of Pericelis, and they argue for abolishing both Pericelidae and Periceloidea in favour of the family Diposthidae.
Our data do not support the following superfamilies sensu Bahia et al. (2017): The position of Chromoplanoidea within Cotylea is supported by most of our analyses (Table 2), although in the 18S28Slong MRE tree, Chromoplanoidea is recovered as acotylean (Fig. 5). The superfamily always is monophyletic, but the interrelationships between the three included chromoplanoid genera are differently resolved. In Bahia et al. (2017), Theama + Chromyella form a sister group to Boninia, while in almost all of our trees, including the MRE trees, Chromyella + Boninia are sister group to Theama. Curiously, in the only trees of our dataset supporting Theama + Chromyella (Suppl. Figs. S13, 15, 21), we used the same alignment method (MUSCLE), the same reconstruction method (RAxML), a partial 28S matrix and Gblocks, just like Bahia et al. (2017). In Giribet (2014, 2017), the remaining possibility is realised, i.e. Theama + Boninia are sister group to Chromyella.