Morphological versus molecular delimitation of ciliate species: a case study of the family Clevelandellidae (Protista, Ciliophora, Armophorea)

The endozoic ciliates of the family Clevelandellidae Kidder, 1938 typically inhabit the hindgut of wood-feeding panesthiine cockroaches. To assess the consistency of species delimitation in clevelandellids, we tested the utility of three sources of taxonomic data: morphometric measurements, cell geometrical information, and 18S rRNA gene sequences. The morphometric and geometrical data delimited the clevelandellid morphospecies consistently and unambiguously. However, only Paraclevelandia brevis Kidder, 1937 represented a homogenous taxon in both morphological and molecular analyses; the morphospecies Clevelandella constricta (Kidder, 1937) and C. hastula (Kidder, 1937) contained two or three distinct, more or less closely related genotypes each; and the genetic homogeneity of the morphospecies C. panesthiae (Kidder, 1937) and C. parapanesthiae (Kidder, 1937) was not corroborated by the 18S rRNA gene sequences at all. Moreover, the 18S rRNA gene phylogenies suggested the C. panesthiae-like morphotype to be the ancestral phenotype from which all other clevelandellid morphotypes arose. The only exception was the C. constricta-like morphotype, which very likely branched off before the diversification of the C. panesthiae-like progenitor. The present molecular analyses also suggested that a huge proportion of the clevelandellid diversity still waits to be discovered, since examination of only four panesthiine populations revealed 10 distinct clevelandellid genotypes/molecular species.


Introduction
At least 25 species concepts have been proposed so far. Interbreeding, genetic or phenotypic cohesion, evolutionary cohesion or evolutionary history are among the most often utilized indicators for delimitation of taxonomic units. Already the existence of multiple species concepts suggests that none of them could be universally used, and the problems of species discrimination might differ case by case (for a review, see Coyne & Orr 2004). There is no simple way out of the 'species problem' because species arise from other species and during the speciation process there are always some ambiguities (Hey 2001). Therefore, a combined approach involving the morphological, phylogenetic, and ecological methods is recently most preferred in protistology (Boenigk et al. 2012;Abraham et al. 2019). In the past, a morphological approach based on multiple morphological characteristics was massively utilized to define ciliate species (e.g., Foissner et al. 2002). With the advent of molecular methods and the availability of DNA sequencing, phylogenetic methods have begun to accompany the morphological approach in delimiting the species boundaries (e.g., Shazib et al. 2016Shazib et al. , 2019Zhao et al. 2016Zhao et al. , 2018Obert & Vďačný 2019). However, the morphological and genetic data may sometimes bring spurious results. Individuals morphologically identical at first glance do not necessarily have to form a monophyletic taxon in phylogenetic analyses, which could be caused by recent divergence or strong natural selection acting on morphological traits (Chenuil et al. 2019). On the other hand, individuals classifiable into distinct taxa based on morphology might be products of the same genotype. These differences in the morphology are considered as a response to various environmental conditions, the process also well known as phenotypic plasticity (Scheiner 1993;Via et al. 1995). In the present study, we attempt to address the consistency of morphological and molecular methods in species discrimination with the example of the ciliate family Clevelandellidae Kidder, 1938. We selected this endozoic family because (1) it contains only about a dozen morphospecies and multiple species co-occur in the same host; (2) all clevelandellid species exhibit a very distinct shape, which makes their identification comparatively simple (Kidder 1937(Kidder , 1938Mandal & nair 1974;Pecina & Vďačný 2020); and (3) some clevelandellid morphospecies may be a cryptic species complex, consisting of multiple genotypes (Lynn & Wright 2013). To address the consistency of morphological and molecular delimitation of clevelandellid species, we analyzed three sources of data: the morphometric measurements, the cell geometrical information, and the 18S rRna gene sequences. For the first time in ciliates, we thus utilized in addition to metric data and DNA sequences also the body shape information, which allowed us to avoid any subjective shape descriptions.

Material collection and processing
Four populations of wood-feeding cockroaches were acquired from commercial breeders in order to be inspected for the presence of clevelandellid ciliates in their hindguts. One population represented an undetermined panesthiine species and three populations belonged to Panesthia angustipennis (Illiger, 1801), specifically to the subspecies P. angustipennis angustipennis (Illiger, 1801), and P. angustipennis cognata Bey-Bienko, 1969. Panesthia angustipennis angustipennis originated from a tropical rain forest in Khao Yai, Nakhon Ratchasima Province, Thailand (referred to as 'Thai I population' henceforth); one population of P. angustipennis cognata came from a tropical rain forest on the Côn Sơn (Grande-Condore) island, Bŕ Rịa-Vũng Tŕu Province, Vietnam (referred to as 'Vietnamese population'); the other P. angustipennis cognata population was obtained from a not more closely specified locality in Cambodia (referred to as 'Cambodian population'); and Panesthia sp. originated from a tropical rain forest in Mae Taeng, Chiang Mai Province, Thailand (referred to as 'Thai II population'). Cockroaches from individual sampling sites were bred separately at room temperature in plastic boxes and were fed by decaying wood. Prior to dissection, cockroaches were euthanized in chloroform vapors. Consequently, their hindguts were extracted and their content was released in sterile Petri dishes containing Ringer's solution (0.63% NaCl solution) to obtain endozoic ciliates.
processed to generate the cell outlines. Altogether 94 outlines were statistically analyzed, using the Momocs package (http://CRAN.R-project.org/package=Momocs), as implemented in the statistical environment R ver. 3.5.1 (R Development Core Team 2020). Prior to the elliptic Fourier analysis, which serves to extract the geometrical information from the outlines, the number of harmonics needs to be determined. The Fourier analyses has its basis in the idea of Fourier series, that is, in decomposing a periodic function into a sum of more simple trigonometric functions such as sine and cosine. These simple functions have frequencies that are integer multiples, i.e., are harmonics, of one another. The lower harmonics provide approximations for the coarse-scale trends in the original periodic function, while the high-frequency harmonics fit its fine-scale variations (Bonhomme et al. 2014). The function 'calibrate_harmonicpower_efourier' was used to determine the number of harmonics that gather 99.9% of the total cumulative harmonic power. The elliptical Fourier analysis was then computed with the function 'efourier', whereby the number of harmonics was set to 17 and the Fourier coefficients were normalized.
The efourier method generates a 'Coe' class object that can be directly used for the principal component analysis (PCA) and hierarchical clustering. The PCA diagram was generated using the functions 'PCA' and 'plot.PCA'. The hierarchical clustering was performed with the 'CLUST' function and included a combination of two coefficients of distance (the Euclidean distance and the Manhattan city block distance) and five clustering algorithms (average linkage, complete linkage, median, centroid, and Ward's D2 method). Finally, multivariate analysis of variance (MANOVA) was performed on the PCA objects and served for testing of difference between shapes of individual clevelandellid morphospecies.

Molecular methods
For the purposes of molecular analyses, individual cells of each clevelandellid morphospecies were picked from the hindgut content of all four cockroach populations. Ciliates were then carefully washed in a series of drops of Ringer's solution. Thoroughly cleaned specimens were lysed separately in 180 µl CLD buffer (Promega, Fitchburg, WI, USA) and stored at 6°C until DNA extraction. Altogether 36 single cell samples were prepared. Genomic Dna was isolated from these samples using the ReliaPrep ™ Blood gDna Miniprep System (Promega, Fitchburg, WI, USa). The 18S rRna gene was amplified using the universal eukaryotic primers Euk a (5'-aaC CTG GTT GaT CCT GCC aGT-3') and Euk B (5'-TGa TCC TTC TGC aGG TTC aC-3') designed by Medlin et al. (1988). Individual polymerase chain reactions (PCR) included 5 µl of the extracted template DNA, 0.4 µl of the forward and reverse primers each (50 pmol/µl), and 10 µl of the GoTaq ® Long PCR Master Mix (Promega, Fitchburg, WI, USA). The final volume was adjusted to 20 µl by adding deionized distilled water. PCR conditions followed our previous protocol (Pecina & Vďačný 2020): initial hot start denaturation at 95°C for 15 min, 30 identical amplification cycles (denaturing at 95°C for 45 s, annealing at 55°C for 1 min, and extension at 72°C for 2.5 min), and final extension at 72°C for 10 min. The resultant PCR products were purified by calf intestinal alkaline phosphatase and exonuclease I (New England Biolabs ® Inc., Ipswich, MA, USA). The purified PCR products were subsequently sequenced on an aBI 3730 automatic sequencer at Macrogen Europe B.V. (Amsterdam, The Netherlands).

Phylogenetic methods
18S rRNA gene sequences acquired from six morphospecies (one belonging to the family Nyctotheridae amaro, 1972 and five to the family Clevelandellidae) were trimmed at the 5' and 3' ends in Chromas ver. 2.6.6 (Technelysium Pty Ltd., South Brisbane, Australia) and assembled into contigs with the aid of BioEdit ver. 7.2.5 (Hall 1999). Based on our previous studies Pecina & Vďačný 2020), all available sequences of the order Clevelandellida de Puytorac & Grain, 1976 and their closest relatives from the paraphyletic order Metopida Jankowski, 1980 were downloaded from GenBank (https://www.ncbi.nlm.nih.gov/genbank) and included in phylogenetic analyses. GenBank accession numbers of the analyzed sequences are provided in the phylogenetic trees. Two 18S rRNA gene alignments were prepared: one dataset contained taxa from the order Clevelandellida and some closely related members of the order Metopida, and the other dataset comprised only taxa from the family Clevelandellidae. Both alignments were constructed on the basis of the secondary structure of the 18S rRNA molecule on the online platform at https://www.arb-silva.de/aligner/ with the aid of the SILVA rRNA database (Quast et al. 2013;Beccati et al. 2017). Since this is a curated database, no masking strategy was used, and any editing was kept to a minimum as recommended by Pruesse et al. (2007).
Two maximum likelihood (IQTREE and PhyML) and two Bayesian (MrBayes and Phycas) methods were employed to reconstruct the phylogenetic relationships within the order Clevelandellida and to determine the phylogenetic positions of the clevelandellid morphotypes. The 50%-majority rule consensus IQtrees were constructed under the best selected evolutionary model using IQTREE (Nguyen et al. 2015) on XSEDE ver. 1.6.10 on the CIPRES portal ver. 3.3 (http://www.phylo.org/) (Miller et al. 2010). The reliability of the branching pattern was assessed with 1000 ultrafast bootstrap replicates (Hoang et al. 2018). The model selection was performed using the in-built ModelFinder program with the best model selected according to the BIC criterion. For the purpose of the three other phylogenetic methods, the best evolutionary substitution models were calculated and selected under the Akaike Information Criterion in jModelTest ver. 2.1.10 (Darriba et al. 2012) on the CIPRES portal. Consequently, all analyses were conducted under the best respective GTR + G + I evolutionary models with parameters as estimated in jModelTest. The PhyML trees were constructed in PHYML ver. 3.0 with the SPR swapping algorithm and 1000 non-parametric bootstrap replicates on the South of France bioinformatics platform (http://www.atgc-montpellier.fr/phyml/) (Guindon et al. 2010). Bayesian Markov Chain Monte Carlo analyses were carried out with the program MrBayes on XSEDE ver. 3.2.7 (Ronquist et al. 2012) on the CIPRES server. They included two runs each with four (one cold and three heated) simultaneous chains of five million generations. The sampling frequency was set to one hundred and the burn-in fraction was specified as 25% of the first sampled trees. Convergence in Bayesian analyses was confirmed in that the average standard deviation of split frequencies was well below 0.01, the potential scale reduction factor approached 1, and no obvious trends were observed in the plots of generations versus log probability. The other Bayesian method was conducted in Phycas ver. 2.2 (Lewis et al. 2015) with the following settings: (1) the GTR + Γ + I evolutionary model as estimated by jModelTest, (2) the inverse gamma hyperprior with mean 2.1 and variance 0.90909 assigned to the hyperparameter μ governing the mean of the branch length prior, (3) the polytomy prior with C = exp(1) favoring unresolved trees with polytomies over more-resolved trees when the difference is lower or equal to one likelihood unit on the log scale, (4) 100 000 burn-in cycles to autotune the updaters in MCMC simulations, (5) one million cycles in MCMC simulations, and (6) sampling trees and parameters every 100 cycles. Convergence and an adequate sample of the posterior were checked in Tracer ver. 1.7.1 (Rambaut et al. 2018).
All trees were computed as unrooted and were rooted a posteriori in FigTree ver. 1.2.3 (http://tree.bio.ed.ac.uk/software/figtree/). Specifically, the root was placed on a branch separating the order Clevelandellida from the selected members of the order Metopida in the case of the larger dataset. As concerns the smaller dataset, the root was positioned on a branch separating C. constricta from the remaining taxa of the family Clevelandellidae, following the topology of the trees derived from the larger dataset. As FigTree may display support values at incorrect nodes after re-rooting the tree, the branch support values were also checked in Mesquite ver. 2.73 (Maddison & Maddison 2007), a tree visualization program whose re-rooting behavior is correct according to the analyses of Czech et al. (2017).

Diversity, prevalence and distribution of clevelandellids
Overall, 24 specimens of Panesthia Serville, 1831 cockroaches were examined for the presence of clevelandellid ciliates: ten specimens from the Thai I population of P. angustipennis angustipennis, ten from the Vietnamese population of P. angustipennis cognata, three from the Cambodian population of P. angustipennis cognata and one specimen from an undetermined panesthiine species originated from the Thai II population. All cockroaches contained clevelandellid ciliates and there were more than 200 specimens per cockroach.
The morphospecies Clevelandella constricta, C. hastula and C. parapanesthiae were most common and were detected in all cockroach populations studied. The morphospecies C. panesthiae was detected in P. angustipennis angustipennis from the Thai I population and in P. angustipennis cognata from the Vietnamese population. The morphospecies C. lynni was noted only in the Thai I population of P. angustipennis angustipennis. Paraclevelandia brevis was recorded in all cockroach populations except for the Cambodian population of P. angustipennis cognata. Nyctotherus galerus Pecina & Vďačný, 2020 was found in all populations except for the Thai II population.
No clear geographic patterns were observed in the occurrence of the clevelandellid genotypes. Some genotypes were noted only at one sampling site, while others were present in various cockroach populations originating from geographically distant localities. Nevertheless, a distinctly increased sampling, covering the whole distribution area of the host cockroaches, is needed to more robustly address the occurrence of individual genotypes across the Asian-Australian region.
During the examination of the hindgut content of panesthiine cockroaches, resting cysts of clevelandellids were also detected. Unfortunately, their numbers were very low and hence insufficient for detailed morphological analyses and for the affiliation to individual genotypes. Vegetative cells of clevelandellids died soon after isolation from the host and wet chambers could therefore not be used for formation of resting cysts. The processing of morphologically studied cysts with molecular methods might be a promising tool for matching various cyst types with individual clevelandellid morphospecies/genotypes in the future.

Characterization of clevelandellid morphospecies
Based on the general morphology, morphometry and geometry, we recognized six clevelandellid morphospecies in the hindgut of the wood-feeding cockroaches P. angustipennis angustipennis from Thailand and P. angustipennis cognata from Vietnam: Clevelandella constricta, C. hastula, C. lynni, C. panesthiae, C. parapanesthiae, and Paraclevelandia brevis. If the same clevelandellid morphospecies was encountered in both cockroach subspecies, their descriptions are kept separate.

Description of Vietnamese population
Size in vivo about 100-155 × 25-50 μm, usually 120 × 35 μm, as calculated from some in vivo measurements and morphometric data; length:width ratio ranging from 3:1 to 4.2:1 after protargol impregnation (Table 1). Body spindle-shaped, more or less distinctly constricted in anterior third, usually widest at mid-portion, i.e., about at level of proximal end of adoral zone of membranelles; dorsoventrally flattened 1.4:1. anterior end bluntly pointed; posterior body portion differentiated into a short, inconspicuous peristomial projection, distinctly constricted at its base; left and right body margins slightly concave at level of macronucleus (Figs 1a-n, 3a, E-G). Macronucleus located in anterior second fifth of body length; ellipsoidal to almost spherical, with a length:width ratio of 1.1-1.8:1 in protargol preparations; 20-28 × 15-20 μm in size after protargol impregnation; filled with numerous globular structures (very likely nucleoli) 0.7-1.8 μm across after protargol impregnation, well recognizable in vivo and sometimes in lightly impregnated specimens. Karyophore attached to right and left cell's margins in middle third of body, usually at level of macronucleus, rarely in midbody. Micronucleus attached to anterior side of macronucleus, i.e., near the place where longitudinal cell axis crosses macronucleus or on its left side; globular and about 4.1-4.3 × 4.5-4.7 μm in size in vivo (Table 1; Figs 1a-L, 3E). Contractile vacuole near left body margin in posterior third of cell, i.e., close to canal leading to cytopyge (Fig. 1a, L). Cortex flexible, no cortical granules recognizable. Cytoplasm colorless; finely granulated; divided by karyophore into an anterior and a posterior part; cytoplasm anterior to macronucleus contains a frontal lamina transversely stretching slightly posterior to anterior body end and densely packed, oval, refractile bodies (probably paraglycogen platelets) observable in vivo and after protargol impregnation; compartment posterior to macronucleus contains some (symbiotic?) bacteria and/or archaea freely scattered throughout the main body portion and food  (Kidder, 1937). Vietnamese specimens isolated from Panesthia angustipennis cognata Bey-Bienko, 1969 from life (A) and after protargol impregnation (B-N).
A. Ventral view of a representative specimen, length 120 μm. B-K. Variability of body shape and size as well as of the nuclear (shaded grey) and oral (shaded yellow) apparatus. L. Semi-schematic diagram, showing the general body organization. Black double arrowhead marks densely packed, oval, refractile bodies (probably paraglycogen platelets). M-N. Ciliary pattern of ventral and dorsal sides. Arrow marks the right suture, black arrowheads indicate the position of the ciliary whorl (posterior suture). O. Prokaryotes freely scattered throughout the cytoplasm posterior to the macronucleus. P. Detail of oval, refractile bodies (probably paraglycogen platelets) anterior to the macronucleus. Scale bars = 50 μm. Somatic ciliature holotrichous; cilia about 4.5-6.5 μm long in vivo and narrowly arranged. approximately 80 ciliary rows narrowly spaced over entire body surface and about 10 ciliary rows encroaching onto peristomial projection. Almost all ciliary rows commence from a whorl (posterior suture) on left side near contractile vacuole to radiate over ventral and dorsal sides toward right body margin; some kineties shortened anteriorly or posteriorly (Figs 1M-n, 3F-G). Right suture extends from base of peristomial   (Kidder, 1937). Thai I specimens isolated from Panesthia angustipennis angustipennis (Illiger, 1801) after protargol impregnation. A-J. Variability of body shape and size as well as of the nuclear (shaded grey) and oral (shaded yellow) apparatus. Scale bar = 50 μm. (Kidder, 1937). Vietnamese (a, E-G) and Cambodian (D) specimens isolated from Panesthia angustipennis cognata Bey-Bienko, 1969, as well as Thai I specimens (B-C) isolated from Panesthia angustipennis angustipennis (Illiger, 1801) from life (a, D-G) and after protargol impregnation (B-C). A-C. Ventral view of specimens with well-preserved body shape. D-E. Ventral view, showing the general body organization. Arrows mark oval, refractile bodies anterior to the macronucleus, black arrowheads mark the proximal end of the adoral zone of membranelles, white arrowheads denote the karyophore attached to the right and left body margins and black double arrowhead marks the canal leading from the contractile vacuole to the cytopyge. F-G. Ciliary pattern of ventral and dorsal sides. Asterisks mark the position of the ciliary whorl (posterior suture), white double arrowhead denotes the right suture. Scale bars: a-C, E-G = 50 μm; D = 20 μm.
Peristomial projection occupies on average 16% of body length and measures 12-23 × 15-20 μm in protargol preparations. Peristomial opening situated on ventral side of peristomial projection, triangular and short, i.e., about 11% of body length and 10-15 × 13-18 μm in size after protargol impregnation (Figs 1a-M, 3a, E). Peristomial funnel approximately 50 μm long after protargol impregnation. Adoral zone extends obliquely from distal end of peristomial projection along left side of peristomial funnel to terminate about in mid-portion of cell; occupies 41% to 55% of body length; composed of on average 55 membranelles; cilia of distalmost membranelles about 9 μm long in vivo and projecting out of peristomial opening (Table 1; Figs 1A, L, 3E). Paroral membrane diplostichomonad, i.e., composed of two rows of basal bodies; extends in parallel with adoral zone but on opposite side of peristomial funnel; commences at level of proximal end of peristomial opening and terminates near cytostome at proximal end of peristomial funnel (Fig. 1L). Pharyngeal fibres spread from proximal end of adoral zone and paroral membrane, run transversely leftwards forming a conical funnel about 20 μm long in vivo.

Notes on Thai I population
The Thai population matches very well the Vietnamese population. However, the Thai population shows a slightly wider range in the body length (95-175 μm vs 100-155 μm) and is slightly more slender than the Vietnamese population (length:width ratio 3.4-5.0:1 vs 3.0-4.2:1). The variability of body size and shape of the Thai population is summarized in Table 1 and shown in Figs 2a-J, 3B-C.

Description of Vietnamese population
Size in vivo about 75-105 × 25-35 μm, usually 90 × 30 μm, as calculated from some in vivo measurements and morphometric data; length:width ratio ranging from 2.6:1 to 3.5:1 in protargol preparations (Table  2). Body spear-shaped, widest at mid-portion, i.e., about at level of contractile vacuole. Anterior end pointed; posterior body portion differentiated into a conspicuous, long peristomial projection; left side distinctly curved at level of proximal end of adoral zone of membranelles and hence forming a lobe above the base of peristomial projection (Figs 4A-N, 5A-H). Macronucleus located in anterior second fourth of body; ellipsoidal, with a length:width ratio of 1.3-1.7:1 in protargol preparations; 15-23 × 11-15 μm in size after protargol impregnation; filled with innumerable globular structures (presumably nucleoli) 0.7-1.8 μm in diameter after protargol impregnation, well observable in vivo and in some protargol preparations. Karyophore absent. Micronucleus invariably attached to anterior side of macronucleus; almost globular to broadly ellipsoidal with a length:width ratio of 1.0-1.7:1; about 4-6 × 3-4 μm in size after protargol impregnation (Table 2;  Peristomial projection conspicuous because it occupies on average 38% of body length and measures 23-34 × 9-15 μm after protargol impregnation. Peristomial opening situated on left ventral side of peristomial projection, roughly triangular and relatively large, i.e., about 17% of body length and 11-15 × 9-12 μm in size after protargol impregnation (Figs 4a-M, 5a-G). Peristomial funnel approximately 34 μm long in protargol preparations. adoral zone extends slightly obliquely from distal end of peristomial   projection across right side of peristomial funnel to terminate about at level of base of peristomial projection; occupies 36% to 50% of body length; composed of on average 32 membranelles; cilia of distalmost membranelles about 8 μm long in vivo and projecting out of peristomial funnel (Table 2; Figs 4a, L, 5B-G). Paroral membrane diplostichomonad, i.e., composed of two rows of basal bodies; runs in parallel with adoral zone on opposite side of peristomial funnel; commences about at level of proximal end of peristomial opening and terminates near cytostome at proximal end of peristomial funnel (Figs 4L, 5E). Pharyngeal fibres originate from proximal end of adoral zone and paroral membrane, run transversely leftwards forming a conical funnel about 12 μm long in vivo.

Description of Thai I population
A detailed morphological description, including morphometric characterization, is provided in Pecina & Vďačný (2020). The original description is to be published in Journal of Eukaryotic Microbiology.
To avoid nomenclatural problems, the name Clevelandella lynni is disclaimed here for nomenclatural purposes (Article 8.3 of the ICZN 1999).

Description of Thai I population
Size in vivo about 95-120 × 45-65 μm, usually 105 × 55 μm, as calculated from some in vivo measurements and morphometric data; length:width ratio rather stable, ranging from 1.7:1 to 2.2:1 after protargol impregnation (Table 3). Body shape obcordiform in ventral view, widest about at midportion, i.e., usually near level of proximal end of peristomial funnel. Anterior end pointed or bluntly pointed; posterior body portion differentiated into a short but conspicuous peristomial projection; left body margin markedly convex (Figs 6A-N, 8B-D). Macronucleus located in anterior half of body; tearshaped, with a length:width ratio of 1.8-2.3:1 in protargol preparations; anterior end broadly rounded, posterior end narrowly rounded and oriented toward right body margin; 34-44 × 17-23 μm in size after protargol impregnation; filled with numerous globular structures (very likely nucleoli) 0.9-1.9 μm across after protargol impregnation, well recognizable in vivo and after protargol impregnation. Karyophore attached to right and left cell's margins (Table 3; Figs 6A-L, 8B-D). Micronucleus not observed. Contractile vacuole at left side of body, slightly above base of peristomial projection (Figs 6A, L, 8B, D). Cortex flexible, no cortical granules recognizable. Cytoplasm colorless; finely granulated; divided by karyophore into an anterior and a posterior part; compartment anterior to macronucleus contains few oval, refractile bodies (very likely paraglycogen platelets) recognizable in vivo and after protargol impregnation; compartment posterior to macronucleus contains some free (symbiotic?) bacteria and/or archaea scattered throughout cytoplasm and few food vacuoles only about 1.8-2.9 μm in diameter with prey prokaryotes (Fig. 6A). Swims slowly; dies quickly on microscope slides, possibly due to presence of oxygen; body shape changed in dying and strongly squeezed cells, i.e., body becomes thicker and peristomial projection shortens.
Somatic ciliature holotrichous; cilia about 4.0-6.0 μm long in vivo and very narrowly arranged. Approximately 100 ciliary rows narrowly spaced over entire body surface and about 20 ciliary rows encroaching onto peristomial projection. Ciliary rows on peristomial projection appear as oblique lines in ventral view, while as very shallow arcs in dorsal view ( Fig. 6M-N). Almost all ciliary rows begin from a whorl (posterior suture) on left body side slightly above base of peristomial projection, i.e., near location of contractile vacuole ( Fig. 6M-N, arrowheads), radiate across ventral and dorsal sides toward right body margin; some kineties shortened anteriorly or posteriorly. Right suture spreads from base of Table 3. Morphometric data on Clevelandella panesthiae (Kidder, 1937) from Thailand and Vietnam.   peristomial projection to anterior body end; formed by obliquely abutting ventral and dorsal ciliary rows (Fig. 6M, arrow).
Peristomial projection short but conspicuous, occupying on average 21% of body length and measuring 16-22 × 12-19 μm in protargol preparations. Peristomial opening situated on the left ventral border of the peristomial projection, roughly triangular and long, i.e., occupying about 21% of body length and measuring 18-22 × 8-12 μm after protargol impregnation (Figs 6a-M, 8B-D). Peristomial funnel about 44 μm long after protargol impregnation. adoral zone extends slightly obliquely from distal end of peristomial projection across right side of peristomial funnel, terminating near posterior end of macronucleus; occupies from 41% to 55% of body length; composed of on average 45 membranelles; cilia of distalmost membranelles about 8 μm long in vivo and projecting out of peristomial funnel (Table 3; Fig. 6A, L). Paroral membrane diplostichomonad, i.e., composed of two rows of basal bodies; runs in parallel with adoral zone on opposite side of peristomial funnel; commences about at level of proximal end of peristomial opening and terminates near cytostome at proximal end of peristomial funnel (Fig. 6L).

Notes on Vietnamese population
The Vietnamese population displays a similar cell organization as the Thai population. However, the former population is more slender and has a longer peristomial projection. On the other hand, the macronucleus is markedly wider in the Vietnamese population than in the Thai population. The variability of body size and shape of the Vietnamese population is summarized in Table 3 and shown in Figs 7a-D, 8a, E-G. Fig. 7. Clevelandella panesthiae (Kidder, 1937). Vietnamese specimens isolated from Panesthia angustipennis cognata Bey-Bienko, 1969 after protargol impregnation. A-D. Variability of body shape and size as well as of the nuclear (shaded grey) and oral (shaded yellow) apparatus. Scale bar = 50 μm. Clevelandella parapanesthiae (Kidder, 1937) Figs 9-11

A-D.
Ventral views of specimens with well-preserved body shape. Black arrowheads mark the proximal end of the adoral zone of membranelles. E-G. A strongly squeezed specimen by pressure of the cover slip, causing the body to become markedly wider and the notch at the base of the peristomial projection to be lost. The general body organization is shown in (E), the ciliary pattern of ventral and dorsal sides is shown in (F) and (G). asterisks mark the position of the ciliary whorl (posterior suture), white arrowhead denotes the karyophore attaching to right body margin, white double arrowhead denotes the right suture. Scale bars = 30 μm.
Somatic ciliature holotrichous; cilia about 4.5-6.0 μm long in vivo and very narrowly arranged. about 90 ciliary rows narrowly spaced over entire body surface and about 25 ciliary rows running onto peristomial projection. Peristomial ciliary rows in a form of strongly oblique lines in ventral view, while in a form of shallow arcs in dorsal view (Fig. 9M-N). Almost all body ciliary rows begin from a whorl (posterior suture) on left body side slightly above base of peristomial projection, i.e., near location of contractile vacuole ( Fig. 9M-N, arrowheads) to radiate over ventral and dorsal sides toward right body margin; some kineties shortened anteriorly or posteriorly. Right suture extends from base of peristomial projection to anterior body end; formed by obliquely abutting ventral and dorsal ciliary rows (Fig. 9M, arrow).
Peristomial projection conspicuous because it occupies on average 24% of body length and measures 15-23 × 10-12 μm in protargol preparations. Peristomial opening situated on left side of peristomial projection, roughly triangular and slender, i.e., about 21% of body length and 13-20 × 7-9 μm in size after protargol impregnation (Figs 9a-M, 11C-D). Peristomial funnel about 37 μm long in protargol preparations. Adoral zone extends slightly obliquely from distal end of peristomial projection across right side of peristomial funnel to terminate about in half of body length; occupies 41% to 57% of body length; composed of on average 38 membranelles; cilia of distalmost membranelles about 8 μm long in vivo and projecting out of peristomial funnel (Table 4; Figs 9A, L, 11C-D). Paroral membrane diplostichomonad, i.e., composed of two rows of basal bodies; extends in parallel with adoral zone but on opposite side of peristomial funnel; commences about at level of proximal end of peristomial opening and terminates near cytostome at proximal end of peristomial funnel (Fig. 9L). Pharyngeal fibres originate from proximal end of adoral zone and paroral membrane, run transversely leftwards forming a conical funnel about 25-30 μm long in vivo.

Notes on Vietnamese population
The Vietnamese population exhibits similar morphological features as the Thai population. However, the Vietnamese population shows a slightly narrower range in the body length (70-90 μm vs 70-110 μm) and is slightly more slender (length:width ratio 1.8-2.1:1 vs 1.7-1.9:1) than the Thai population. The macronucleus is shorter in some Vietnamese specimens, but its variability completely falls within the range of the Thai population (length 26-35 μm vs 27-47 μm). The variability of body size and shape of the Vietnamese population is summarized in Table 4 and shown in Figs 10a-J, 11a-B, E-G.

Description of Vietnamese population
Size in vivo about 45-60 × 20-35 μm, usually 50 × 30 μm, as calculated from some in vivo measurements and morphometric data; length:width ratio ranging from 1.7:1 to 2.2:1 after protargol impregnation (Table 5). Body obcordiform, widest at beginning of posterior third, i.e., about at level of contractile vacuole; dorsoventrally flattened 1.4:1. anterior end usually bluntly pointed, rarely pointed; posterior body portion broadly rounded; left side slightly curved at level of peristomial opening and thus forming a more or less noticeable lobe; edge of dorsal side longer than that of ventral side and hence forming a distinct lip overhanging peristomial opening (Figs 12A-N, 14D-E). Macronucleus located in anterior half of body; ellipsoidal with a length:width ratio of 1.8-2.3:1 in protargol preparations; 18-23 × 9-12 μm in size after protargol impregnation; filled with innumerable globular structures (very likely nucleoli) 0.7-1.9 μm in diameter after protargol impregnation, well observable in vivo and after protargol impregnation. Karyophore completely surrounds macronucleus, attaches to anterior body pole forming a distinct funnel (Table 5;    Anterior body end to macronucleus, distance VN 6.6 6.5 flexible, no cortical granules recognizable. Cytoplasm colorless; finely granulated; contains some free (symbiotic?) bacteria and/or archaea and food vacuoles about 2.6-3.7 μm across with prey prokaryotes (Fig. 12A). Swims slowly; dies quickly on microscope slides, possibly due to presence of oxygen.
Somatic ciliature holotrichous; cilia about 4.0-5.0 μm long in vivo and very narrowly arranged. Approximately 50 ciliary rows narrowly spaced over entire body surface. Almost all ciliary rows begin from a whorl (posterior suture) on left body side, near location of contractile vacuole, to radiate across ventral and dorsal sides toward right body margin; some kineties shortened posteriorly (Fig. 12M-N, arrowheads). Right suture extends from level of peristomial opening to anterior body end; formed by obliquely abutting ventral and dorsal ciliary rows (Fig. 12M, arrow).
Peristomial opening situated on posterior body pole, occupies about 11% of body length and measures 4-6 × 10-17 μm after protargol impregnation (Figs 12a-M, 14D-E). Peristomial funnel about 17 μm long in protargol preparations. Adoral zone extends almost in parallel with main body axis along right side of peristomial funnel to terminate about at level of posterior end of macronucleus; occupies 34% to 41% of body length; composed of on average 18 membranelles; cilia of distalmost membranelles about 7 μm long in vivo and projecting out of peristomial funnel (Table 5; Figs 12a, L, 14D-E). Paroral membrane not observed.

Notes on Thai I population
The Thai population matches very well the Vietnamese population. However, the body size of the Thai population is smaller than in the Vietnamese population (30-45 μm vs 45-60 μm). On the other hand, Fig. 13. Paraclevelandia brevis Kidder, 1937. Thai I specimens isolated from Panesthia angustipennis angustipennis (Illiger, 1801) after protargol impregnation. A-J. Variability of body shape and size as well as of the nuclear (shaded grey) and oral (shaded yellow) apparatus. Scale bar = 20 μm. the peristomial funnel is longer in the Thai population than in the Vietnamese population (up to 55% vs 41% of body length). The variability of body size and shape of the Thai population is summarized in Table 5 and shown in Figs 13a-J, 14a-C, F-G.

Morphometric and shape analyses
To assess the morphometric variation, distinctness, and boundaries of the six clevelandellid morphospecies, we utilized a multivariate approach including cluster analyses (CA) and metric multidimensional scaling (MDS). Altogether 16 morphometric features were used to calculate pairwise similarities (measured by Gower's coefficient) of 20 Paraclevelandia and 74 Clevelandella specimens. The pairwise Gower's similarity matrix was then subjected to six hierarchical CA, employing the average linkage, the weighted average linkage, the complete linkage, the median, the centroid, and the Ward's clustering method. Since all clustering algorithms consistently depicted each morphospecies as a distinct group, we present here only one dendrogram that was produced by the weighted average linkage method (Fig. 15A). Likewise, MDS conducted on the Gower's similarity matrix generated six mutually well-isolated and homogenous groups each representing one morphospecies (Fig. 16A).
Shape analyses, including CA and PCA, brought very similar results as did morphometric analyses, i.e., each morphospecies formed a well-delimited and homogenous group. We present here only one illustrative dendrogram that was computed with the Ward's D2 clustering method in a combination with the Manhattan city block distance (Fig. 15B). This dendrogram was also completely consistent with the PCa diagram based on the Fourier coefficients (Fig. 16B). ManOVa performed on the PCa objects revealed statistically significant differences between shapes of individual clevelandellid morphospecies (Hotelling-Lawley trace = 48.25, approximate F 55, 382 = 67.03, P < 2.2 × 10 -16 ).

Phylogenetic analyses
Thirty-six new 18S rRNA gene sequences of endozoic ciliates belonging to the order Clevelandellida, isolated from four populations of Panesthia cockroaches, were obtained. Length, GC content, and GenBank accession numbers of the new sequences are provided in Table 6. Four phylogenetic methods (IQTree, PhyML, MrBayes, and Phycas) were used to reconstruct relationships within the order Clevelandellida and to test for monophyly of the clevelandellid morphotypes. With respect to individual datasets, all phylogenetic methods resulted in similar topologies except for some statistically poorly supported nodes that might be considered as soft polytomies. Therefore, only IQTrees are presented along with nodal supports from all statistical methods (Figs 17-18).
The larger dataset served, especially, to uncover the phylogenetic position of the family Clevelandellidae within the order Clevelandellida, to test for the monophyly of the family Clevelandellidae, and to reveal its fundamental bifurcation. Apart from the outgroup, the larger dataset included members of the family Sicuophoridae Amaro, 1972 (represented only by Sicuophora multigranularis Xiao et al., 2002), Nyctotheridae, and Clevelandellidae. In all analyses, S. multigranularis clustered with the Nyctotheridae + Clevelandellidae clade with full statistical support. The family nyctotheridae was depicted as paraphyletic and contained the monophyletic family Clevelandellidae (96% IQTrees, 78% PhyML, 1.00 MrBayes, 1.00 Phycas). The nyctotheridae + Clevelandellidae clade exhibited a clustering specific for higher taxa of their host organisms. One cluster received full statistical support in all analyses and comprised ciliates that had been isolated exclusively from the large intestine of amphibians (Nyctotherus cordiformis (Ehrenberg, 1831), Nyctotheroides deslierresae Affa'a, 1991, Nyctotheroides hubeiensis Li et al., 1998, Nyctotheroides parvus (Walker, 1909, Nyctotheroides pyriformis (Nie, 1932), and Nyctotheroides sp. AF147882). The other cluster obtained poor statistical support (57% IQTrees, 43% PhyML, 0.79 MrBayes) and contained ciliates that had been isolated exclusively from the hindgut of arthropods (all remaining members of the family Nyctotheridae and all members of the family    Pecina & Vďačný, 2020, which had been isolated from three Panesthia populations, were identical and grouped together with full statistical support. This species joined a cluster of all other ciliates isolated from cockroaches, forming a so-called 'cockroach clade' though with variable statistical support (85% IQTrees, 79% PhyML, 0.88 MrBayes, 1.00 Phycas). Nyctotherus velox Leidy, 1849, isolated from a myriapod, was placed in a sister position to the 'cockroach clade' in most analyses. However, this grouping received only negligible statistical support (57% IQTrees, 43% PhyML, 0.79 MrBayes), and was not recognized in Phycas analyses at all. Therefore, we consider the phylogenetic position of N. velox to be questionable (Fig. 17).
As concerns the family Clevelandellidae, its type genus Clevelandella was depicted paraphyletic because it contained members of the genus Paraclevelandia (Figs 17-18). All sequences of P. brevis were identical and therefore grouped together with full statistical support in all analyses. Paraclevelandia brevis grouped with the C. hastula clade with variable support in both datasets. All members of the C. constricta morphotype clustered together with full or strong statistical support in all trees inferred from both datasets. Clevelandella constricta was depicted as sister to all other clevelandellids with strong statistical support. By contrast, the Asian and Australian specimens of the morphospecies C. panesthiae and C. parapaesthiae did not form a distinct cluster each, but were mixed together or grouped with other Clevelandella species. Namely, C. panesthiae specimens from the Thai I population were classified in two clusters, one of which grouped with C. parapanesthiae isolates from the Thai II population. The Australian C. panesthiae KC139718 specimen clustered together with C. parapanesthiae specimens from the Cambodian, Thai I, and Vietnamese populations. The 18S rRNA gene sequences of all these C. parapanesthiae specimens were identical and differed from the C. panesthiae KC139718 exemplar in only 2 nucleotide positions. The remaining Australian C. panesthiae specimens (KC139715, KC139716, KC139717) formed a well-supported and structured clade, indicating that they belong to two distinct species. The Australian C. parapanesthiae KC139719 individual grouped with the Australian C. nipponensis KC139714 specimen with full statistical support in all analyses (Figs 17-18).
To summarize, monophylies of only three morphospecies, viz., C. constricta, N. galerus and P. brevis, were strongly statistically supported. Nevertheless, we cannot exclude that the C. constricta morphotype   Table 6. The scale bar denotes three substitutions per one hundred nucleotide positions.  Fig. 17. For specimen codes and further details, see Table 6. The scale bar denotes six substitutions per one thousand nucleotide positions.
comprises two distinct molecular species (genotypes), as suggested by the substructure within this morphotype in 18S rRNA gene trees inferred from the smaller dataset (Fig. 18). According to the present phylogenetic analyses, the C. panesthiae morphotype unites at least five distinct molecular species, the C. parapanesthiae morphospecies contains at least three molecular species, and the C. hastula morphospecies covers at least three molecular species (genotypes).

Morphological versus molecular species delimitation
Morphology has been the basic criterion to determine different species in protistology since the ancient times of Ehrenberg (1838) and Stein (1859). The morphospecies concept became more objectively analyzable with the application of numeric taxonomic methods, such as hierarchical clustering, principal component analysis, multidimensional scaling, and canonical discriminant analysis (for reviews, see Marhold 2011 andVďačný et al. 2014). These methods were utilized in ciliate species delimitation soon after their invention (e.g., Berger & Hatzidimitriou 1978;Foissner & Schubert 1983;Lynn & Malcolm 1983;Jones & Gates 1994), and rather recently have been also used to assess the taxonomic reliability of quantitative features and the consistency of classification to morphospecies, for instance, within the genera Semispathidium Foissner et al., 2002(Vďačný et al. 2014, Metopus Claparède & Lachmann, 1858 (Vďačný & Foissner 2017), Anoplophrya Stein, 1860and Metaradiophrya Jankowski, 2007(Obert & Vďačný 2019), or in Paraholosticha Wenzel, 1953(Jung & Berger 2019. The geometrical analyses, which serve to evaluate the shape information (for a review, see MacLeod & Forey 2002), commonly supplement the morphometric approach in zoology and botany. However, the shape analyses are still only rarely utilized in protistology (e.g., Healy-Williams & Williams 1981;Belyea & Thunell 1984), and geometrical techniques have not up to now been used to address the morphospecies concept in ciliates to our best knowledge. This method has very likely not been utilized for ciliate species delimitation hitherto due to some methodological concerns. Indeed, geometrical analyses are mostly used on rigid objects, such as shells or other skeletal structures. Therefore, a microphotograph depicting the actual shape of the ciliate is indispensable for this sort of analyses. A promising way out of this problem is examination of cells fixed with Bouin's solution, which preserves the body shape of clevelandellids in particular and of other ciliates in general (Foissner 2014) very well. The present study showed that eigenshapes, which had been calculated from the geometrical information of cells fixed with Bouin's solution and impregnated with protargol, correspond very well to the actual shape of clevelandellid morphospecies in vivo (Fig. 16B). Geometrical data thus might be successfully applicable in taxonomic studies of ciliates, when their shape is well-preserved after fixation and when cells are not misshapen after protargol impregnation. Moreover, the statistical assessment of the position of homologous cirri in hypotrichous ciliates with the aid of landmark analyses might be another useful tool in geometrical analyses of ciliates.
With the advent of molecular methods, the morphology-based delimitation of ciliate species has become more objectively testable also using various genetic markers (for a review, see Abraham et al. 2019).
There are multiple studies focused on integrated taxonomy at the specific and generic levels in various groups of ciliates, such as heterotrichs, euplotids, peritrichs, astomes or peniculines (e.g., Gentekaki & Lynn 2010;Boscaro et al. 2014;Przyboś et al. 2015;Zhao et al. 2018;Obert & Vďačný 2019, 2020Shazib et al. 2019). These studies suggest that molecular data tend to reveal a higher species diversity than the morphological approach per se. Such a discordant picture might be, however, also due to subtle morphological differences among ciliate species, which are often difficult to study.
In the present study, we integrated three approaches -the morphometric measurements, the cell geometrical information, and the 18S rRNA gene sequences -to assess the consistency of the morphological and molecular species delimitation with the example of the endozoic family Clevelandellidae. Our complex approach revealed that the assignment of clevelandellids to morphospecies, as defined by Kidder (1937), is fully consistent in both morphometric and geometrical analyses. In other words, each clevelandellid morphospecies forms a homogeneous cluster that is well isolated from all other morphospecies analyzed (Figs 15-16). Since both techniques are phenetic, their results refer only to the similarity of the analyzed species in the phenotypic space and not to their relatedness. To reveal the kinships within and among the clevelandellid morphospecies and to further independently test their homogeneity, we utilized the phylogenetic information contained in the 18S rRNA gene. However, monophylies of only two morphospecies, viz., C. constricta and P. brevis, were statistically strongly supported. The remaining morphospecies contained multiple more or less phylogenetically distant genotypes (Figs 17-18). The most dramatic situation was detected within the C. panesthiae and C. parapanesthiae morphospecies. More specifically, the former morphotype covers at least five distinct genotypes and the latter morphotype at least three genotypes. The taxonomic problem of the morphospecies C. panesthiae was encountered already by Lynn & Wright (2013), whose isolates from the Australian P. cribrata Saussure, 1864 were assignable even to three different genotypes. Lynn & Wright (2013) therefore suggested that this may be a cryptic species complex. Our data further indicate that C. hastula might also cover at least three more or less closely related genotypes, and two of them co-occur within the Thai population of Panesthia sp. The substructure of the C. constricta cluster is also suggestive of further diversification, since the Australian individuals isolated from P. cribrata are separated from the Asian specimens isolated from P. angustipennis (Figs 17-18).
The present study thus shows that morphology and genetics are rather inconsistent in classification of members of the family Clevelandellidae. Although the morphometric and geometrical data delimit the clevelandellid morphospecies consistently and unambiguously, their homogeneity is corroborated only in two out of the five testable morphospecies by the 18S rRna gene sequences and the remaining three morphospecies cover three to five distinct genotypes each. Consequently, a huge part of the real diversity of clevelandellids might be missed when only morphological data are used to identify species. Thus, the use of morphological identification alone is likely to miss distinct genotypes and can lead to the underestimation of diversity and overestimation of host range. Nevertheless, we cannot exclude that there are other characters, such as the resting cyst morphology, some ontogenetic peculiarities, or host specific symbiotic prokaryotes, which can help to delimit all clevelandellid genotypes also morphologically.

Phylogeny and molecular taxonomy of clevelandellids
The family Clevelandellidae was consistently recognized as a monophyletic group nested within the paraphyletic family Nyctotheridae (Lynn & Wright 2013;Li et al. 2018;Pecina & Vďačný 2020). In the present study, for the first time, we could test for the monophyly of the genus Clevelandella, because we obtained 18S rRNA gene sequences from another clevelandellid genus, Paraclevelandia. Kidder (1937) suggested that both Clevelandella and Paraclevelandia arose from a Nyctotherus-like progenitor. Furthermore, he speculated that Paraclevelandia could be considered as a basal lineage, because it lacks the peristomial projection, a peculiar structure that is not present in Nyctotherus Leidy, 1849 as well. Lynn & Wright (2013) suggested that the peristomial projection may have become longer towards the phylogenetically younger species of the genus Clevelandella. However, Paraclevelandia is not placed as a sister taxon of the genus Clevelandella but is depicted in a sister position to C. hastula, which has by contrast one of the longest peristomial projections. Interestingly, the evolutionary hypothesis of Lynn & Wright (2013) is partially corroborated in that C. constricta, which has the shortest and most inconspicuous peristomial projection among congeners, is classified as a sister taxon of a cluster uniting all other Clevelandella and Paraclevelandia species. Though the Fourier shape analysis is a phenetic method, it also placed Paraclevelandia within the genus Clevelandella (Fig. 15B), showing that the body shape of P. brevis is most similar to that of C. panesthiae (Fig. 16B). However, in the morphometric analyses, Paraclevelandia was placed outside the genus Clevelandella (Fig. 15A).
This classification was very likely caused by a conspicuous difference in the body size between these two genera. Although the geometrical analysis is a phenetic approach, it indicates along with molecular phylogenies that Paraclevelandia is not so dissimilar from Clevelandella as may one judge solely by the absence/presence of the peristomial projection. Indeed, there are also some other Clevelandella species (e.g., C. constricta, C. panesthiae, and C. parapanesthiae) that display an inconspicuous or a short peristomial projection.
The present phylogenetic analyses revealed that multiple Clevelandella species are complexes of more or less closely related taxa. The 18S rRNA gene sequences within the deep-branching C. constricta morphospecies differ by 0.3-0.5%, whereby the greatest divergence is between the Asian and Australian specimens. By contrast, all P. brevis sequences are identical and differ from the closest related sequences of the C. hastula morphospecies by up to 2.1%. The observed divergence within the C. hastula morphospecies was as high as 1.6% (Supplementary Table S2). Genotypes of each of those three clevelandellid morphotypes cluster together, which is indicative of a comparatively early stage of their diversification. On the other hand, genotypes of the C. panesthiae and C. parapanesthiae morphotypes do not cluster together. Since there are many more C. panesthiae genotypes scattered throughout the clevelandellid tree of life, the most parsimonious solution would be to consider the C. panesthiaelike morphotype as the ancestral phenotype from which the C. parapanesthiae-like lineages arose at least three times independently, as well as all other clevelandellid morphotypes. The only exception is represented by the C. constricta-like lineages, which branched off before the diversification of the C. panesthiae-like progenitor. As already suggested by Lynn & Wright (2013), C. panesthiae represents a cryptic species complex. Up to now, it contains at least five different genotypes, three being detected in the Australian P. cribrata (Lynn & Wright 2013) and two in the Asian P. angustipennis. The Australian genotypes differ from the asian ones by 1.0-4.1%. Divergences among the three australian genotypes range from 1.5% to 3.5%, and the two Asian genotypes differ by 1.5% (Supplementary Table S2). Three markedly different lineages were also detected within the C. parapanesthiae morphospecies. The australian genotype differs significantly from both asian genotypes by ~4.3%, and the asian genotypes exhibit a sequence divergence of about 0.8% (Supplementary Table S2).
although there is no interspecific divergence threshold for ciliates, there are Tetrahymena species that are morphologically indistinguishable and whose 18S rRNA gene sequences are identical. However, they are distinct taxa according to the cytochrome oxidase subunit I (COI) sequences. Some morphologically very distinct Tetrahymena species differ only by 0.3% in the 18S rRNA gene but can be over 10% divergent in the COI sequences (for details, see Rataj & Vďačný 2020 and references cited therein). A very similar picture has been detected, for instance, within the heterotrich genus Spirostomum Ehrenberg, 1834(Shazib et al. 2019. One must be, however, very cautious whether molecularly cloned or not cloned sequences are compared. More specifically, PCR products typically contain a mixture of gene copies, but some may be in distinctly lower quantities and could therefore be obscured by the predominant ribotype variant. Not cloned sequences need to therefore be considered as consensus sequences, and one should apply a much stricter interspecific divergence threshold to them. On the other hand, during molecular cloning also rare ribotypes might be captured and the within cell variability in the 18S rRna gene might be as high as 0.70% in litostomes (Vďačný et al. 2011), 0.92% in oligotrichs (Gong et al. 2013), and even 1.9% in astomes (Rataj & Vďačný 2019). Since the present clevelandellid sequences were not cloned, a much more conservative threshold should be applied. Considering the taxonomic practice in Tetrahymena Furgason, 1940 (Doerder 2019 and references cited therein), we find the 0.3% difference in the 18S rRNA gene as a strong indication of a different species. Nevertheless, we shall subject the clevelandellid genotypes delimited by the 18S rRNA gene sequences to further molecular analyses, using the hypervariable ITS-region and the first two barcoding domains of the 28S rRNA gene as well as the mitochondrial 16S rRNA gene in a following study. Likewise, we prefer to await further molecular data and analyses to solve the problem of paraphyly of the genus Clevelandella.
In their estimate of the free-living ciliate diversity, Foissner et al. (2008) speculated that the number of described morphospecies must be multiplied by at least 5 in the light of genetic and molecular data. Nevertheless, Foissner et al. (2008) remained conservative in their calculations and used multiplicators of only 2 and 3. Although we examined just four populations of panesthiine cockroaches, we detected two to three molecular species in four out of the five clevelandellid morphospecies investigated (a single sequence is available for C. lynni and therefore its monophyly could not yet be tested). When data from the Australian clevelandellids studied by Lynn & Wright (2013) are added, the number of genotypes within the C. panesthiae morphospecies rises to five. It is important to mention that molecular methods have been used to investigate clevelandellids isolated from only three panesthiine species (P. angustipennis, P. cribrata, and Panesthia sp.), and the genus Panesthia includes well over 50 taxa (Beccaloni 2020).
No other panesthiine genera have been examined for the presence of clevelandellid ciliates in their hindguts hitherto. This might be an important aspect in the diversity estimations, since our integrative taxonomic approach has revealed that ciliates are much less promiscuous to their hosts than previously anticipated. Indeed, hosts could constitute sharply isolated niches that might permit speciation of their symbiotic ciliates. For instance, four out of the five astome ciliates isolated from lumbricid earthworms and all three Tetrahymena species isolated from freshwater planarians were detected exclusively in the same host species (Obert & Vďačný 2019;Rataj & Vďačný 2019, 2020. The host specificity is indirectly also corroborated in that the specimens assigned to the morphospecies C. panesthiae and C. parapanesthiae isolated from the Asian P. angustipennis and the Australian P. cribrata did not group together in 18S rRNA gene phylogenies and obviously represent distinct evolutionary lineages (Figs 17-18). However, the extent of host specificity of clevelandellids needs to be analyzed on a much larger sample to confirm or reject the consistency of this eco-evolutionary trend.  estimated the net diversification rate of the order Clevelandellida to be on average 0.0282 lineages per one million years, with a 95% credibility interval ranging from 0.0055 to 0.0518. Under a birth-death model, the expected species (lineage) diversity is as follows: where E [N t ] is the expected number of species at time t, N 0 is the starting species diversity, r is the net diversification rate, and t is the age of a clade. According to our previous Bayesian relaxed molecular clock estimates, the family Clevelandellidae emerged on average about 50 Ma ago. Using a very conservative starting number of clevelandellid lineages to be only 12, we estimate the expected number of clevelandellids to be roughly on average 50 species, with an upper 95% credibility limit of about 160. Indeed, the present study shows that just in four panesthiine populations we have detected 10 distinct clevelandellid genotypes, and each very likely represents a separate species. Thus, the examination of only four panesthiine populations almost doubled the number of clevelandellid species known! This also indicates that we might still have underestimated the net diversification rate of clevelandellids in our previous study . Therefore, very likely a huge proportion of the clevelandellid diversity still waits to be discovered when further panesthiine cockroaches commence being studied.