Multi-Gene Phylogeny of the Ciliate Genus Trachelostyla (Ciliophora, Hypotrichia), With Integrative Description of Two Species, Trachelostyla multinucleata Spec. nov. and T. pediculiformis (Cohn, 1866)

Many hypotrich genera, including Trachelostyla, are taxonomically challenging and in a need of integrative revision. Using morphological data, molecular phylogenetic analyses, and internal transcribed spacer 2 (ITS2) secondary structures, we attempt to cast more light on species relationships within the genus Trachelostyla. The present multifaceted approach reveals that (1) a large-sized species with numerous macronuclear nodules, isolated from sandy littoral sediments in southern China, is new to science and is endowed here with a name, T. multinucleata spec. nov.; (2) two other Chinese populations previously identified as T. pediculiformis represent undescribed species; and (3) multigene phylogeny is more robust than single-gene trees, recovering the monophyly of the genus Trachelostyla with high bootstrap frequency. Additionally, ITS2 secondary structures and the presence of compensatory base changes in helices A and B indicate the presence of four distinct taxa within the molecularly studied members of the genus Trachelostyla. Molecular data are more suitable for delimitation of Trachelostyla species than morphological characters as interspecific pairwise genetic distances of small subunit (18S) rDNA, ITS1-5.8S-ITS2, and large subunit (28S) rDNA sequences do not overlap, whereas ranges of multiple morphometric features might transcend species boundaries.

Trachelostyla Borror, 1972, is one of many taxonomically "difficult" hypotrich genera. Distinguishing Trachelostyla species is difficult because the proposed taxa have a rather similar body size and shape as well as overlapping numbers of macronuclear nodules. This genus was introduced by Kahl (1932) and originally included two species, namely, T. pediculiformis (Cohn, 1866) and T. caudata Kahl, 1932. As none of them was fixed as type species, the erection of Trachelostyla was invalid at that time. Borror (1972) made the genus name available by fixing T. pediculiformis as the type species (Berger, 1999;Aescht, 2001). Due to the lack of type locality and information on some taxonomically important morphological features, Gong et al. (2006) neotypified T. pediculiformis and revised the diagnosis of Trachelostyla.
Currently, this genus is defined by a combination of the following features: a dorsoventrally flattened, nonspirally twisted, and elongated body with a peristomial region conspicuously narrowed; 11 cirri in the peristomial region (i.e., three frontal cirri, four frontoventral cirri, one buccal cirrus, and three postoral ventral cirri); one left and one right marginal cirral row not confluent posteriorly; and three caudal cirri. Gong et al. (2006) also established a phylogenetically related genus, Spirotrachelostyla, which differs from Trachelostyla by the spirally twisted body and about 13 cirri scattered in the anterior peristomial region. Three species (Stichotricha simplex Kahl, 1932;Trachelostyla spiralis Dragesco and Dragesco-Kernéis, 1986; and Trachelostyla tani Hu and Song, 2002) were originally included in Spirotrachelostyla.
Trachelostyla was reported by comparatively many researchers (for reviews, see Berger, 2008). However, the vast majority of reports did not provide detailed descriptions of the ciliary pattern based on silver-impregnated material. Some morphometric differences should not be, therefore, over-interpreted. In this light, Berger (2008) recognized only three valid species, namely, T. caudata (Kahl, 1932) Berger, 2008;T. pediculiformis (Cohn, 1866) Borror, 1972;and T. rostrata (Lepsi, 1962) Berger, 2008. Because of their overall morphological similarity as well as the lack of protargol-impregnated specimens and molecular data, some Trachelostyla populations were very likely misidentified, or their taxonomic status was at least questionable. For instance, as already mentioned by Berger (2008), a Chinese population of T. pediculiformis studied by Xu and Song (1999) conspicuously differs from the neotype population of Gong et al. (2006) in body size and the location of the rearmost postoral ventral cirrus, which is distinctly set off from the other two postoral cirri. Kahl's (1928) population of T. pediculiformis was preliminarily classified as incertae sedis by Berger (2008) because it deviates from the typical T. pediculiformis by having two ventral cirral rows and two macronuclear nodules. On the other hand, T. pediculiformis does not exhibit any ventral cirral rows and possesses many scattered macronuclear nodules. Besides the interphase morphology, ontogenetic data might provide some clues for clearing taxonomic confusion as well. However, a detailed description of the whole ontogenetic process was reported only for a single T. pediculiformis population by Shao et al. (2007).
Previous phylogenetic studies suggest that 18S, ITS1-5.8S-ITS2, and 28S rDNA data potentially have the power to separate closely related hypotrich species (Huang et al., 2016;Li et al., 2017;Chen X. et al., 2021;Fan et al., 2021;Wang et al., 2021b). In the present study, we assessed the validity of phylogenetic relationships among Trachelostyla species, using a combination of morphological and multigene data as well as the secondary structure of the ITS2 molecule.

Sampling and Cultivation
Trachelostyla multinucleata spec. nov. was collected from the top 5 cm of sandy littoral sediments on the Xiaolajia Island in Daya Bay, near the city of Huizhou, southern China (22 • 36 40 N,114 • 38 02 E) on April 1, 2018, when the water temperature was 26 • C and the salinity was 33 . Trachelostyla pediculiformis population 1 (pop. 1) was isolated from sediments and seawater collected from Tangdao Bay, Qingdao, China (35 • 55 57 N, 120 • 11 43 E) on November 5, 2017, when the water temperature was 15 • C and the salinity was 32 . Both samples were divided into aliquots that were used to establish raw cultures in Petri dishes at room temperature (23 • C). Some rice grains were added to stimulate the growth of bacteria, which served as prey organisms for ciliates. Specimens from cultures were used for the subsequent molecular analyses.

Taxonomic Methods and Terminology
Trachelostyla multinucleata spec. nov. and T. pediculiformis pop. 1 were investigated using a combination of detailed in vivo observation and protargol impregnation as described by Zhang et al. (2020b). Morphometric data were obtained from living and protargol-impregnated specimens. Illustrations of living cells were based on free-hand sketches and photographs. To distinguish parental and daughter structures during the morphogenetic processes, new (daughter) structures are painted solid, whereas old (parental) ciliary structures are depicted by contour. General terminology and classification follow Berger (2008).
Thirteen taxonomically important features were counted and measured on 21 protargol-impregnated specimens of T. multinucleata spec. nov., 11 individuals of T. pediculiformis pop. 1, and 16 cells of T. pediculiformis population 3 (pop. 3 from Huang et al., 2016; Table 1 and Supplementary Table 1). Morphometric data were processed in Python ver. 3.6.6 with the libraries NumPy (Oliphant, 2015) and Pandas (McKinney, 2010) to calculate pairwise similarities between all specimens using Gower's coefficient. The pairwise Gower's similarity matrix served as an input for the metric multidimensional scaling (MDS), which was conducted with the help of the SMACOF algorithm and the scikit-learn ver. 1.0 package 1 (Pedregosa et al., 2011). MDS included 250 initializations each with 20,000 iterations. The MDS diagram was graphically prepared using the Matplotlib ver. 3.4.3 package (Hunter, 2007).

DNA Extraction, PCR Amplification, and Sequencing
Extraction of the genomic DNA and PCR amplification of 18S, ITS1-5.8S-ITS2, and 28S rDNA as well as sequencing of T. multinucleata spec. nov. followed our previous study (Zhang et al., 2020b). Conspecific sequences obtained from multiple samples were identical, and therefore, only those derived from the single-cell sample were included in the subsequent phylogenetic analyses.
As for T. pediculiformis pop. 1, the single-cell genome isolation and amplification were carried out with the Repli-g single-cell kit (Qiagen, Hilden, Germany). PCR products were purified with AMPure XP beads (Beckman Coulter, IN, United States) and quantified by Qubit 3.0 Fluorometer (Invitrogen, Waltham, MA, United States). The final library was sequenced on the Illumina NovaSeq 6000 platform (Illumina, San Diego, CA, United States) with a 2 × 150 bp paired-end run in the Novogene Company (Beijing, China). Finally, rDNA sequences were extracted from assembled contigs. The quality of the rDNA contig was measured by Bowtie2 version 2.4.4 (Langmead and Salzberg, 2012) and Samtools version 1.14 (Danecek et al., 2021).

Molecular Phylogenetic Methods
The newly obtained sequences were blasted against the nucleotide NCBI database 2 . The BLASTn algorithm revealed that they belong to the family Trachelostylidae (subclass Hypotrichia). All available trachelostylid sequences were carefully checked; some sequences designated as Trachelostyla sp. but without associated publications were excluded from the phylogenetic analyses and genetic distance analyses. Taxon sampling in the single-gene data set (18S rDNA) and multigene data set (18S + ITS1-5.8S-ITS2 + 28S rDNA) mostly followed Huang et al. (2016). Additionally, some other sequences published after the study of Huang et al. (2016) and related to Trachelostyla according to the BLASTn search were included in phylogenetic analyses as well. Four species (Apodiophrys ovalis Jiang and Song, 2010;Diophrys scutum (Dujardin, 1841) Kahl, 1932Paradiophrys zhangi Jiang et al., 2011;and Uronychia multicirrus Song, 1997) were used as outgroup taxa, which follows Huang et al. (2016). GenBank accession numbers are shown in the respective figure and Supplementary Table 2. The sequences of T. pediculiformis populations 2 and 3 are from Huang et al. (2016). However, no 18S rDNA sequences are available from population 2, and thus it could not be included in the phylogenetic analyses.
Sequences were aligned online on the MAFFT ver. 7 server 3 (Katoh et al., 2019), using the iterative refinement G-INS-i method, the gap opening penalty at 1.53, and the 200PAM/κ = 2 scoring matrix for nucleotide sequences. No masking strategy was used. The 5 and 3 ends of the resulting alignments were trimmed manually in the program BioEdit ver. 7.0 (Hall, 1999). SeaView ver. 4 was used to prepare the concatenated data set (Galtier et al., 1996;Gouy et al., 2010). The single-gene alignment contained 1,798 nucleotide positions, and the multigene data set comprises 3,678 positions. The reference alignments are provided in the Supplementary Material. The number of unmatched nucleotides and the pairwise p-distances of Trachelostyla and Spirotrachelostyla species were calculated with the help of the program BioEdit ver. 7.0 (Hall, 1999), using the sequence difference count matrix and sequence identity matrix options.
Maximum likelihood (ML) analyses were computed with the program IQ-TREE ver. 1.6.10 (Nguyen et al., 2015) on the IQ-TREE server 4 (Trifinopoulos et al., 2016). Each molecular marker (i.e., 18S rDNA, ITS1-5.8S-ITS2 region, 28S rDNA) was assigned the best evolutionary substitution model as chosen by the in-built program under the Bayesian information criterion (Supplementary Table 3). Nodal support was assessed with 1,000 ultrafast bootstrap pseudo-replicates. All other parameters were left default. Bayesian inference was carried out in MrBayes ver. 3.2.7 (Ronquist et al., 2012) on the CIPRES portal ver. 3.3 5 (Miller et al., 2010). Prior parameters of evolutionary models as estimated with IQ-TREE were implemented in Bayesian analyses with the "prset" command. Four Markov chain Monte Carlo simulations were run for 5,000,000 generations with a sampling frequency of 100 and a relative burn-in fraction of 25% (first 12,500 trees). Convergence of the MCMC analyses was confirmed in that the average standard deviation of split frequencies was well below 0.01, the potential scale reduction factor approached 1, effective sample sizes were >200, and no obvious trends were in the plots of generations vs. log probability. ML and BI trees were computed as unrooted and were rooted using the outgroup taxa in FigTree ver. 1.2.3 6 .

Prediction of Internal Transcribed Spacer 2 Secondary Structure
Boundaries of ITS2 were determined according to Weisse et al. (2008) as well as Obert and Vd'ačný (2020). More specifically, we aligned trachelostylid sequences against Meseres corlissi (EU399522-29) and Plagiotoma lumbrici (MN176618-23), and then searched for the 5.8S-28S rRNA proximal stem of the ITS2 molecule. Predictions of the putative secondary structure of the ITS2 molecules included the formation of the imperfect 5.8S-28S rRNA helix, homology modeling, and free-energy minimization approach as implemented on the Mfold server ver. 3.0 7 (Zuker, 2003). Thermodynamically optimal secondary ITS2 structures were manually processed in VARNA ver. 3.93 (Darty et al., 2009). The number of nucleotides in bulges and loops was counted and evaluated for each structural domain of the ITS2 molecules. Compensatory base changes (CBCs) were determined with the CBCAnalyzer option (Wolf et al., 2005) implemented in 4SALE ver. 1.7.1 (Seibel et al., 2006), and hemi-CBCs were searched for manually as recommended by Shazib et al. (2016). The tertiary structure of the ITS2 molecule was modeled using the online program RNAComposer ver. 1.0 8 (Popenda et al., 2012).  Small andLynn, 1985 Genus Trachelostyla Borror, 1972 Trachelostyla multinucleata Spec. nov.

Diagnosis
Size in vivo 200-310 × 35-60 µm. Body elongated, distinctly constricted in anterior fifth and slightly tail-like narrowed in posterior fifth. Macronuclear apparatus composed of about 60-90 nodules. Contractile vacuole located in posterior body fifth. Three frontal, one buccal, and four frontoventral cirri; three postoral ventral, two pretransverse, and five transverse cirri; three inconspicuous caudal cirri. One marginal cirral row on each side, composed of about 25-35 left and 30-45 right marginal cirri. Adoral zone slightly bipartite, composed of usually four apical ordinarily spaced membranelles and about 100 lapel narrowly spaced membranelles. Seven dorsal kineties with extremely long, spine-like dorsal cilia. Marine sandy sediment habitat.

Type Locality
Sandy littoral sediments from the Xiaolajia Island in the Daya Bay, near the city of Huizhou, southern China (22 • 36 40 N, 114 • 38 02 E).

Type Material
A protargol slide (no. ZTY2018040106_1) with the holotype specimen marked with an ink circle and nine paratype slides (no. ZTY2018040106_2-10) as well as the DNA sample (ID Collection Code: C287) of a voucher specimen have been deposited in the Laboratory of Protozoology, Ocean University of China (OUC).

Etymology
The species-group name multinucleata is a composite of the stem of the Latin quantifier mult·us (many), the thematic vowel ·i-, and the Latin adjective nȗcleāt·us, -a, -um ([m; f; n], containing a nucleus), referring to the many macronuclear nodules whose number significantly exceeds that in other congeners.
Size of specimens from fresh raw cultures 205-310 × 35-60 µm, usually about 260 × 45 µm, length:width ratio ca. 5-6:1; body rather flexible but not contractile. Body shape elongated with narrowly rounded anterior end and broadly rounded posterior end; outline highly characteristic because tripartite: anterior body fifth distinctly constricted to a snoutlike structure, trunk region cylindrical with more or less parallel margins, and posterior body fifth slightly tail-like narrowed ( Figures 1A,B, 2I); posterior body constriction indistinct or missing in cultivated specimens, causing the body to appear bipartite ( Figures 1B, 2J-L). Ventral side flat or slightly concave, dorsal side often arched in mid-body; dorsoventrally flatted 2:1 ( Figures 1C, 2M); As many as 58-89 macronuclear nodules scattered throughout cytoplasm in trunk and tail-like region; individual nodules ovoid or ellipsoidal, 3.8-7.9 µm long after protargol impregnation, containing small chromatin bodies; two to five roughly globular micronuclei, 2.6-4.2 µm in largest diameter, scattered among macronuclei, sometimes difficult to recognize from smaller macronuclear nodules and exact number, thus, difficult to determine (Figures 1H, 2A,B,E and Table 1). Contractile vacuole located in posterior fifth of body length, approximately 18 µm across during diastole, usually very difficult to observe because cytoplasm studded with innumerable dark cytoplasmic inclusions and, thus, seen in about two out of the 13 specimens examined (Figures 1B, 2J-L). Cytoplasm colorless, packed with many granules 3-5 µm across, inclusions, and food vacuoles, rendering the cell with a dark and opaque appearance at low magnifications (Figures 2I-O). Cortex highly fragile, often does not withstand coverslip pressure, causing the cell to burst; no cortical granules recognizable. Crawls rapidly on debris particles, sometimes swims by rotation about main body axis.    Dorsal kineties, number 7 7 7 7.0 7.0 0 0 10 -6 6 6.0 6.0 0 0 6 -6 6 6.0 6.0 0 0 16 All data based on protargol-impregnated specimens. Measurements in µm. H, holotype specimen; Min, minimum; Max, maximum; Mean, arithmetic mean; SD, standard deviation; CV, coefficient of variation in %; n, number of cells investigated.  Frontal-ventral-transverse cirral pattern and number of cirri constant, viz., cirri arranged in a 3:1:4:3:2:5 pattern (Figures 1D,G, 2A-C,F and Table 1). Cilia of frontal, buccal, frontoventral, postoral ventral, and pretransverse cirri 16-19 µm long; cilia of marginal cirri about 16 µm long; and cilia of transverse cirri 25-30 µm long. Three frontal cirri same in size as remaining cirri, obliquely arranged behind apical adoral membranelles. One buccal cirrus situated at level of mid-portion of adoral zone. Four frontoventral cirri located posterior to rearmost frontal cirrus and anterior to level of buccal cirrus, form a roughly rectangular pattern. Three postoral ventral cirri displaced distinctly anteriorly compared with most oxytrichids, i.e., located right anterior to paroral membrane, two anterior cirri conspicuously separated from third posterior cirrus (Figures 1D,G, 2C). Five transverse cirri distinctly enlarged and arranged in a J-shaped pattern caudally; two pretransverse cirri situated anterior to leftmost and rightmost transverse cirrus (Figures 1G, 2F). Right marginal cirral row commences slightly above anteriormost frontoventral cirrus and terminates slightly behind level of posteriormost transverse cirrus; two anteriormost cirri remarkably more widely spaced than remaining cirri (Figures 1D,G, 2A-C). Left marginal cirral row begins slightly above buccal vertex and terminates ahead of level of posteriormost transverse cirrus, marginal cirral rows, thus, not confluent posteriorly (Figures 1G, 2A,B); one extra cirrus (likely the first left marginal cirrus) situated at right of beginning of the left marginal cirral row and above the proximal end of the adoral zone of the membranelles, sometimes difficult to recognize in deeply impregnated specimens because they are easily hidden by deeply impregnated membranelles (Figures 1D,F,G, 2B,D, red arrowheads).
Adoral zone of membranelles occupies about 40% of body length in vivo and approximately 48% in protargol preparations, divided into an apical and a lapel region. Apical region composed of four ordinarily spaced and radially extending membranelles, cilia 25-28 µm long (Figures 1A,D Table 1). Lapel region, composed of 86-113 narrowly spaced membranelles, cilia up to 15 µm long, forms a Gonostomumlike pattern, i.e., runs along left body margin to level of midportion of paroral membrane, where it bends rather abruptly rightward, plunging into a comparatively long buccal cavity; a zigzag structure very likely composed of membranellar fibers recognizable right of lapel membranelles in some protargolimpregnated specimens (Figures 1D,G, 2A-D). Undulating membranes arranged almost in parallel, hence resembling a Gonostomum-like pattern. Paroral membrane commences left of postoral ventral cirri and runs to buccal vertex, distinctly longer than endoral membrane, about 35 µm long. Endoral membrane starts about in mid-portion of paroral membrane and terminates more or less at same level as paroral membrane, about 16 µm long (Figures 1D,F,G, 2D).

Morphogenesis
Only one early and one late divider were found in protargol preparations (Figures 2H, 3A-C). Morphogenesis commences apokinetally with the formation of two fields of densely arranged basal bodies, viz., the proter's oral primordium situated near the buccal vertex and the opisthe's oral primordium located below the mid-body (Figures 2H, 3A). All parental cirri appear intact at this stage and are, thus, apparently not involved in the formation of new oral primordia. In the late division stage, parental structures (i.e., adoral membranelles, undulating membranes, frontal-ventral-transverse cirri, marginal cirri, and caudal cirri) are almost resorbed and new structures are nearly completely developed both in the proter and the opisthe (Figures 3B,C).

Physiological Reorganization
Physiological regeneration resembles divisional morphogenesis, but there is only one oral primordium and a single set of six FVT cirral anlagen. Two early stage (Figures 4A,B, 5A-C) and three middle-stage (Figures 4C-F, 5D-H) reorganizers as well as one late-stage reorganizer (Figures 4G, 5I,J) were found in protargol preparations.
The oral primordium is formed de novo as a longitudinal field of densely arranged basal bodies posterior to the buccal vertex (Figures 4A, 5A, blue triangles). Parental structures are not involved in the formation of new FVT cirral anlagen. These emerge as short streaks on the right side of the anterior half of the oral primordium (Figures 4B, 5B blue arrowheads). At this stage, also new adoral membranelles begin to organize at the anterior end of the oral primordium in a posteriad direction (Figures 4B, 5B). The new FVT cirral anlagen grow posteriorly along the differentiating adoral zone of membranelles to become longitudinal streaks. In middlestage reorganizers, the six FVT streaks split new cirri in the following pattern 1, 3, 3, 3, 4, 4 ( Figures 4C-G, 5D,F,H,J). The formation of 18 FVT cirri in reorganizers, thus, corresponds to the divisional morphogenesis shown in Figure 3B. The new undulating membranes are developed from the FVT cirral anlage I (Figures 4G, 5J, red asterisks).
Both new marginal cirral anlagen are generated within the mid-portion of the parental rows (Figures 4B-F, 5B,C,E,F,H, red arrowheads and red double-arrowheads). Subsequently, the new marginal cirral anlagen differentiate into marginal cirri in a posteriad direction to completely replace the parental ones (Figures 4G, 5J). The reorganization of the dorsal ciliature occurs in a unique way. Two parallel dorsal kinety anlagen are generated intrakinetally within the mid-portion of the parental dorsal kinety 7 (Figures 4B, 5G). They grow in a posteriad direction to become the new dorsal kineties 6 and 7 (Figures 4C-G, 5I). On the other hand, only a single anlage develops in the middle of the parental dorsal kinety 1 (Figures 4B-E, 5B,E,F,H). Later on, this anlage generates five dorsal kineties in the late reorganization stage (Figures 4G, 5J). Finally, a single caudal cirrus is differentiated at the posterior end of dorsal kinety 5 (Figures 4G, 5J, red arrows) and two caudal cirri are generated at the posterior end of dorsal kineties 6 and 7 (Figures 4G,  5I, red arrows).

Morphological Redescription of the Qingdao Population of Trachelostyla pediculiformis (Cohn, 1866)
Body size 105-150 × 20-40 µm in vivo, length:width ratio 4.0-5.5:1; body rather flexible but not contractile. Body elongated and bipartite, i.e., anterior body fifth distinctly snout-like constricted and trunk broadly cylindrical; anterior body end narrowly rounded while rear body end broadly rounded (Figures 6A,B,  7A,B,D-F). Ventral side flat, dorsal side slightly arched in mid-body; dorsoventrally flatted approximately 2:1 ( Figure 7C). Seven to 16 ovoid or ellipsoidal macronuclear nodules, 3.5-12.0 × 3.0-8.7 µm in size after protargol impregnation, arranged in a ring-like pattern in trunk; two more or less globular micronuclei, 2.8-3.5 µm in diameter, usually one micronucleus located near anteriormost macronuclear nodules and one micronucleus near posteriormost nodules or in cell center (Figures 6B, 7H-J,N). Contractile vacuole located in posterior body fifth, approximately 15 µm across during diastole, difficult to recognize due to presence of numerous endoplasmic granules ( Figure 7D). Cytoplasm colorless and packed with innumerable irregularly shaped granules 0.5-3.0 µm across, providing the cell with a grayish and opaque appearance at low magnification (Figures 7A-G). Cortex highly fragile; no cortical granules recognizable. Crawls rapidly on debris particles.

17-19 dikinetids
Adoral zone of membranelles occupies about 40% of body length in vivo and about 50% in protargol preparations; composed of 44 membranelles on average; divided into an apical and a lapel region (Figures 6A, 7A, Table 1). Apical region composed of three to five, 18-20 µm long, strong, and radially extending membranelles. Lapel membranelles narrowly spaced and extend along left margin of anterior body portion to plunge into a comparatively long buccal cavity, forming a Gonostomum-like pattern; length of membranellar cilia up to 10 µm (Figures 6A, 7A-F,H-J and Table 1). Undulating membranes arranged almost in parallel; paroral membrane begins at level of postoral anteriormost ventral cirri, extends to buccal vertex, and about 14 µm long in protargol preparations; endoral membrane commences about in mid-portion of paroral membrane, terminates at buccal vertex, and only about 6 µm long after protargol impregnation (Figures 6A, 7N).

Morphometric Analyses
Morphological variation and boundaries of Trachelostyla species were analyzed using the multidimensional statistical approach and Gower's similarity coefficient. MDS generated two clusters that were well separated along the second ordination axis. One cluster contained only specimens of T. multinucleata spec. nov., and the other cluster comprised individuals of T. pediculiformis pop. 1 and 3. Although there was a tendency to separate these two populations of T. pediculiformis, morphometric data were not sufficient to unambiguously delimit their boundaries ( Figure 6C). This result, however, conflicts with the present phylogenetic, distance, and CBC analyses (see below) as both populations are genetically well separated and obviously represent distinct biological species.

Phylogenetic Analyses
The rDNA sequences of T. multinucleata spec. nov. were obtained by the bidirectional Sanger sequencing, whereas the rDNA contig of T. pediculiformis pop. 1 (7,495 bp) was obtained from the assembled contigs of high-throughput sequencing. There were 176,857,976 raw reads in total, and 313,070 of them could be mapped to our rDNA contig with 95.3% similarity. The reads depth for this contig is 3,607, further proving its high quality and reliability. GenBank accession numbers and length of the nuclear 18S rDNA, ITS1-5.8S-ITS2, and 28S rDNA sequences of T. multinucleata spec. nov. and T. pediculiformis pop. 1 are summarized in Table 2.
The phylogenetic positions of T. multinucleata spec. nov. and T. pediculiformis pop. 1 were determined with ML and BI. Topologies of ML and BI trees were nearly congruent, and hence, only the ML tree is shown with nodal supports from both methods (Figures 8, 9). Although the taxon sampling slightly differed between the single-and multigene data sets, Trachelostyla and Spirotrachelostyla consistently clustered with the urostylid Caudikeronopsis marina with full statistical support. The monophyletic origin of trachelostylids was strongly supported only in the multigene ML tree (97%), whereas it was very weakly supported in the single-gene ML tree (56%). Bayesian analyses provided only poor support (0.90) for monophyly of trachelostylids in the case of multigene data, whereas their common origin was not recognized in the case of single-gene data. The genus Trachelostyla was depicted as monophyletic in the multigene trees and paraphyletic in the single-gene trees. The reader is, therefore, referred to the fully resolved and strongly statistically supported multigene trees in which Trachelostyla multinucleata spec. nov. groups with T. pediculiformis pop. 3, and T. pediculiformis pop. 1 clusters with the neotype population of T. pediculiformis (Figure 9).
The number of unmatched nucleotide positions and the pairwise p-distances among members of the family Trachelostylidae are summarized in Table 3. The 18S rDNA sequence divergences among T. multinucleata, T. pediculiformis neotype pop., and T. pediculiformis pop. 1 and 3 ranged from 2.1% to 2.5%, corresponding to as many as 34-41 unmatched nucleotide positions. The observed genetic divergence in 18S rDNA sequences between T. pediculiformis pop. 1 and the neotype T. pediculiformis was 0.2% as there were only two unmatched nucleotide positions. The p-distances between Trachelostyla and Spirotrachelostyla species varied from 2.7% to 3.2%, i.e., there were 43-52 unmatched nucleotide positions ( Table 4). There was no difference in the ITS1-5.8S-ITS2 region between T. pediculiformis pop. 1 and the neotype T. pediculiformis. On the other hand, p-distances among T. multinucleata, T. pediculiformis neotype pop., T. pediculiformis pop. 1-3 varied from 4.0% to 14.9%, with 18-68 unmatched nucleotide positions. Distance between Trachelostyla and Spirotrachelostyla species ranged from 8.1% to 15.1% as there were 37-69 unmatched nucleotide positions ( Table 4).

Putative Internal Transcribed Spacer 2 Secondary Structure and Compensatory Base Change Analyses
The ITS2 secondary structures of five Trachelostyla taxa and Spirotrachelostyla tani were proposed using Mfold and are shown in Figures 10, 11A,B. The tertiary structure model of the T. multinucleata ITS2 molecule was predicted based on the secondary structure model and is shown in Figure 11C. The length of ITS2 molecules ranges from 188 nt in T. pediculiformis pop. 3 to 194 nt in the T. pediculiformis neotype, T. pediculiformis pop. 1, and T. pediculiformis pop. 2. The trachelostylid ITS2 molecules consist of a central loop radiating two helices (A and B) of unequal length. The central loop is composed of 54 nucleotides in all taxa. Helix A is distinctly shorter than helix B and invariably consists of 22 nucleotides, forming eight pairs and a terminal hexaloop. Helix B displays a much more complex structure as it is composed of a furcation loop and four subhelices. Subhelix B-1 contains seven nucleotide pairs in all taxa except for T. pediculiformis pop. 3, which possesses six pairs instead due to one pyrimidine-pyrimidine mismatch. Subhelix B-2 is composed of one unpaired nucleotide (G) and 16 or 20 paired nucleotides. It is separated from subhelix B-1 by an adenineadenine mismatch in all taxa except for T. pediculiformis pop. 3, which has an adenine-cytosine mismatch, and S. tani, which displays one unpaired uracil at the boundary of subhelices B-1 and B-2. Subhelix B-2 is separated from subhelices B-3 and B-4 by a furcation loop. Subhelix B-3 is the most variable among Trachelostyla and Spirotrachelostyla species in both length and nucleotide composition. This subhelix contains four nucleotide pairs in T. pediculiformis pop. 3, five pairs in T. multinucleata, six pairs in S. tani, and seven pairs in T. pediculiformis pop. 2 and T. pediculiformis (both the neotype and pop. 1). Subhelix B-4 is the longest with 16 nucleotide pairs, two bulges, and a terminal tetraloop. The numbers of unmatched nucleotides and the counts of compensatory base changes between Trachelostyla and Spirotrachelostyla species are summarized in Supplementary  Table 4. There was one CBC detected in helix B between T. pediculiformis (neotype pop. and pop. 1) and T. pediculiformis pop. 3 as well as between T. pediculiformis (neotype pop. and pop. 1) and S. tani (Supplementary Table 4). Hemi-CBCs were found in both helices and between all species (Figures 11A,B). The high number of unmatched nucleotides as well as the presence of CBCs and hemi-CBCs document that T. pediculiformis (neotype pop. and pop. 1), T. pediculiformis pop. 2, T. pediculiformis pop. 3, and T. multinucleata represent a distinct biological species each. As expected, no CBCs were found between the neotype population and population 1 of T. pediculiformis supporting their conspecificity.
Our new species undoubtedly belongs to the genus Trachelostyla as documented by the nonspirally twisted and elongated body with a conspicuously narrowed peristomial region, the 11:2:5:3 cirral pattern (i.e., 11 cirri in the frontal region, two pretransverse cirri, five transverse cirri, and three caudal cirri), and the Gonostomum-like oral apparatus. Berger (2008) based the taxonomy of the genus Trachelostyla on the morphology of the posterior body end and the structure of the nuclear apparatus. As mentioned, he recognized three valid species: T. caudata, T. pediculiformis, and T. rostrata.
Although Trachelostyla caudata has been reported several times, all descriptions are based only on live observations, and hence, no type material is available. As its specific epithet indicates, T. caudata has a tail-like posterior body end, which is also a property of T. multinucleata. Nevertheless, T. multinucleata can be distinguished from all populations of T. caudata by the much larger body (205-310 µm vs. 120-240 µm according to Kahl, 1932Kahl, , 1933 µm according to Maeda and Carey, 1984;Carey, 1992, and Al-Rasheid, 1999, and 80-120 µm according to Al-Rasheid, 2001) and the distinctly higher number of macronuclear nodules (58-89 vs. around 10 according to Kahl, 1932, 12 according to Maeda andCarey, 1984, more than 20 according to Al-Rasheid, 1999; Table 3).
Trachelostyla multinucleata can be separated from T. pediculiformis and T. rostrata by the much higher number of macronuclear nodules (58-89 nodules in T. multinucleata vs. 9-17 nodules in the T. pediculiformis neotype and two FIGURE 8 | Maximum likelihood tree inferred from SSU rDNA sequences, showing the phylogenetic positions of Trachelostyla multinucleata spec. nov. and Trachelostyla pediculiformis pop. 1 (indicated by bold face and red color). Numbers near branches denote bootstrap values for maximum likelihood (ML) and posterior probabilities for Bayesian inference (BI). Black asterisks indicate the disagreement between ML and BI trees. GenBank accession numbers are provided after species names. Trachelostyla pediculiformis pop. 3 marked by red asterisk was very likely misidentified by Huang et al. (2016). Scale bar corresponds to three substitutions per 100 nucleotide positions. nodules in T. rostrata). According to Borror (1963), the number of macronuclear nodules in T. pediculiformis is unusually variable ranging from 16 to 64. This wide range needs to be, however, taken with caution because Borror (1963) very likely mixed multiple species that are not conspecific with T. pediculiformis and that even originated from different habitats including sand, shell fragments, exoskeletons of isopods, dead oysters, and detritus. Further comparisons with Borror (1963) populations are, unfortunately, not possible as their morphological descriptions are very insufficient.
Finally, T. multinucleata has a tail-like narrowed posterior body end (fresh material is needed), and T. pediculiformis and T. rostrata display a broadly rounded rear body end. The morphometrical distinctness of T. multinucleata from T. pediculiformis is also shown in the MDS diagram ( Figure 6C). Their separation is strongly corroborated also by genetic data in that the pairwise p-distances between T. multinucleata and the T. pediculiformis neotype are 2.5% in 18S rDNA, 7.5% in ITS region, and 6.3% in 28S rDNA sequences ( Table 4).  Supplementary Table S2. Trachelostyla pediculiformis pop. 3 marked by red asterisk was a misidentified material in Huang et al. (2016). Scale bar corresponds to three substitutions per 100 nucleotide positions.

Species Identity of Previous
Trachelostyla pediculiformis Population 2 and 3 Huang et al. (2016) assigned their populations 2 and 3 to T. pediculiformis due to their high morphological similarity with the T. pediculiformis neotype population. Indeed, the present multivariate analyses show that population 3 cannot be unambiguously differentiated from T. pediculiformis even when a combination of 13 morphometric characters is considered ( Figure 6C and Supplementary Table 1). Unfortunately, detailed morphometric data are not available for population 2 and, therefore, its similarity with T. pediculiformis could not be statistically assessed. Despite that fact, populations 2 and 3 genetically differ so conspicuously from T. pediculiformis neotype pop. and T. multinucleata spec. nov. that they cannot be conspecific (Figures 8, 9 and Table 4).

Morphological Versus Molecular Taxonomy of Trachelostyla
Multiple Trachelostyla populations were very likely misidentified in the past (for a review, see Berger, 2008). Even rather recently, two genetically distinct species, Trachelostyla pediculiformis pop. 2 and T. pediculiformis pop. 3, were misassigned to T. pediculiformis based on their high morphological similarity to the neotype population of T. pediculiformis Huang et al., 2016). Indeed, the present multidimensional statistical approach could not unambiguously differentiate T. pediculiformis pop. 1 from T. pediculiformis pop. 3 even when a combination of 13 morphometric characters was employed. On the other hand, T. multinucleata spec. nov. could be morphometrically clearly separated from T. pediculiformis populations (Figure 6C). Despite this fact, multiple morphometric features transcend boundaries among Trachelostyla species (Tables 1, 3).
As molecular taxonomic methods are commonly available nowadays and molecular approaches are less complicated than morphological multivariate analyses, we propose a shift to molecular taxonomy and DNA barcoding to achieve more objective identification and delimitation of Trachelostyla species. Moreover, molecular data appear to be more suitable for delimitation of Trachelostyla species than morphological characters as interspecific pairwise genetic distances of 18S, ITS1-5.8S-ITS2, and 28S rDNA sequences do not overlap. On the other hand, multiple morphometric features might cross species boundaries (Tables 1, 3). Considering the p-distances between the T. pediculiformis neotype pop. and T. pediculiformis pop. 1, we suggest a 0.2% difference in 18S and a 0.6% divergence in 28S rDNA sequences as thresholds for discrimination of Trachelostyla species.
Also, CBC analysis of the ITS2 secondary structures might be of great help to more unambiguously delimit Trachelostyla species. The presence of CBCs strongly correlates with the existence of distinct biological species and, therefore, CBCs have been used as a reliable marker for species delimitation in various eukaryotes (e.g., Coleman, 2005Coleman, , 2009Wolf et al., 2005Wolf et al., , 2013Müller et al., 2007;Ruhl et al., 2010;Sorhannus et al., 2010). CBC-based species delimitation worked also well in multiple ciliate genera, for instance, in Anoplophrya, Metaradiophrya, Paramecium, Pseudokeronopsis, and Vorticella (Coleman, 2005;Sun et al., 2010Sun et al., , 2013Zhan et al., 2019;Obert et al., 2021). On the other hand, no CBCs were detected among distinct Spirostomum (Shazib et al., 2016(Shazib et al., , 2019, Trichodina (Rataj and Vd'ačný, 2021), and clevelandellid (Pecina and Vd'ačný, 2021) species, which very likely corresponds to the shortness and high GC content of their ITS2 molecules. As concerns Trachelostyla, CBC and hemi-CBCs were observed among all Trachelostyla populations studied except the conspecific T. pediculiformis neotype pop. and pop. 1 (Figures 11A,B and Supplementary  Table 4). This indicates that CBC analysis might be a powerful tool for delimitation of Trachelostyla species.
To summarize, the nuclear rDNA cistron (including 18S, ITS-5.8S-ITS2, and 28S rDNA sequences), the ITS2 secondary structure, and CBC analysis might be successfully utilized for identification and discrimination of Trachelostyla species as well as for uncovering new Trachelostyla species. Moreover, the combination of the multifaceted molecular approach might significantly reduce the risk of misidentification and subsequent errors in distributional and ecological analyses. On the other hand, some morphometric features might transcend species boundaries and even multivariate analyses might fail to separate genetically distinct taxa.

Phylogeny of the Family Trachelostylidae
According to Berger (2008), the family Trachelostylidae comprises only two marine genera, viz., Trachelostyla and Spirotrachelostyla. Monophyly of trachelostylids was corroborated not only by morphological data but also by multigene phylogenetic analyses (Figure 9), albeit not by single-gene analyses (Figure 8). Various studies have already shown that concatenation of 18S, ITS1-5.8S-ITS2, and 28S rDNA sequences might lead to better resolved phylogenetic trees in the subclass Hypotrichia than 18S rDNA sequences per se (e.g., Huang et al., 2016;Obert and Vd'ačný, 2020;Zhang et al., 2020b;Wang et al., 2021b). The same applies also to the monophyletic origin of the genus Trachelostyla, which was only indicated in 18S rDNA phylogenies (Figure 8) although well recognizable in the multigene trees (Figure 9). Thus, when multiple molecular markers were employed, Spirotrachelostyla was robustly placed outside Trachelostyla. This is also consistent with morphological classifications as Spirotrachelostyla possesses a spirally twisted and spindle-shaped body, and Trachelostyla displays a nontwisted and bi-or tripartite body.
Caudikeronopsis marina was revealed to be a sister-group taxon of the family Trachelostylidae with full statistical support in both the single-and multigene phylogenies (Figures 8, 9). This result is surprising from a morphological viewpoint because Caudikeronopsis exhibits an urostylid ventral cirral pattern, and trachelostylids have an oxytrichid ventral pattern composed of more or less 18 FVT cirri.
Nevertheless, the mixing of taxa with urostylid and oxytrichid cirral patterns in molecular phylogenies is well known and was reconciled by the CEUU (Convergent Evolution of Urostylids and Uroleptids) hypothesis proposed by Foissner et al. (2004). Later on, Foissner andStoeck (2006, 2008) and Kim et al. (2014) reported further taxa corroborating the CEUU hypothesis. The urostylid ventral cirral pattern very likely evolved or was lost multiple times independently during hypotrich evolution. The close kinship of Caudikeronopsis and trachelostylids is supported by the presence of caudal cirri, which are often lacking in the "core" urostylids. More importantly, Caudikeronopsis and trachelostylids share multiple unique deletions in the ITS2 molecule, which caused their helix B to be composed of only one furcation loop and four subhelices (Figures 10, 11A,B,D). By contrast, helix B contains two furcation loops and five subhelices in all other hypotrichs ( Figure 11D) and oligotrichs studied so far (e.g., Coleman, 2005;Weisse et al., 2008;Li et al., 2013;Zhan et al., 2019;Obert and Vd'ačný, 2020). The single furcation loop of Caudikeronopsis and trachelostylids corresponds to the second furcation loop of oligotrichs and other hypotrichs ( Figure 11D). The remnant of the first furcation loop is composed of an adenine-adenine mismatch in all trachelostylids except for T. pediculiformis pop. 3, which exhibits there an adenine-cytosine mismatch, and S. tani, which possesses only one unpaired uracil at the place of the former furcation loop (Figures 10, 11B). However, the vestige of the first loop is slightly larger and more complex in Caudikeronopsis ( Figure 11B). The first furcation loop was apparently lost due to the complete deletion of subhelix B-5 already in the last common ancestor of Caudikeronopsis and trachelostylids. This loss also caused subhelices B-1 and B-2 to appear confluent in Caudikeronopsis and trachelostylids, whereas both subhelices are well recognizable in oligotrichs and all other hypotrichs (e.g., Coleman, 2005;Weisse et al., 2008;Li et al., 2013;Zhan et al., 2019;Obert and Vd'ačný, 2020).

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/ Supplementary Material.

AUTHOR CONTRIBUTIONS
YW designed and supervised the research. TYZ and TTZ performed laboratory work. TYZ and PV analyzed data and wrote the manuscript. CS, WS, and SA-F revised the manuscript.