Earthworms were collected at 36 localities in southwestern Slovakia, Central Europe (Supplementary Table S1). They were identified using a combination of external morphology (Pižl, 2002) and sequencing of the mitochondrial genes coding for NADH–ubiquinone oxidoreductase chain 1 (ND1) and cytochrome c oxidase subunit I (COI). Primers and PCR conditions used for the amplification of ND1 and COI are summarized in Supplementary Tables S2 and S3. Molecular identification of examined earthworms to species is shown in Supplementary Figures S1, S2. Earthworms were processed, dissected, and inspected for the presence of ciliates as described by Obert and Vďačný (2019, 2020a,b).

Living ciliates were manually isolated from the gut content of their earthworm hosts with Pasteur micropipettes and investigated in vivo at low (50–400 ×) and high (1,000 ×, oil immersion) magnifications with bright field and differential interference contrast under a Leica DM2500 optical microscope following Foissner (2014). The ciliary pattern and nuclear apparatus were revealed with protargol impregnation (Wilbert, 1975). Capturing of photomicrographs, preparation of illustrations, and measurements were conducted as described elsewhere (Zhang and Vďačný, 2021a,b). Terminology mostly follows Berger (1999, 2006).

Altogether, 21 quantitative features were scored on 52 protargol-impregnated specimens originating from epigeic, anecic, and endogeic lumbricid earthworms (Supplementary Table S4). The morphometric matrix was processed in Python ver. 3.6.6, using the libraries NumPy (Oliphant, 2015) and Pandas (McKinney, 2010). Pairwise similarities of Plagiotoma individuals were measured with Euclidean distance. Since this coefficient depends on the scale at which characters are measured, standardization by standard deviation was used to avoid unequal influences on the results (Marhold, 2011). Principal component analysis (PCA) was conducted on the pairwise standardized Euclidean similarity matrix, using the scikit-learn package (Pedregosa et al., 2011). Plotting of the PCA ordination diagram was done with Matplotlib (Hunter, 2007).

The same specimens that had been used for the morphometric multivariate analyses were captured by a Canon EOS 70D camera. Micrographs were imported into Inkscape ver. 0.92¹ and processed to generate cell outlines. These were statistically analyzed with the package Momocs,² as implemented in the environment R ver. 4.0.3 (R Development Core Team, 2020). First, the number of harmonics was determined, i.e., outlines were decomposed into a sum of trigonometric functions. Frequencies of these functions are integer multiples and hence harmonics of one another. The total cumulative harmonic power of 99.9% was obtained already with 13 harmonics. Consequently, the elliptical Fourier analysis was performed with the number of harmonics set to 13. Fourier coefficients were normalized and used for PCA and hierarchical clustering. Finally, multivariate analysis of variance (MANOVA) was performed on the PCA objects and served for testing of difference between body shapes of plagiotomids isolated from different ecological groups of lumbricid earthworms.

Single cells were placed in 180 μl of cell lysis buffer (Promega, Fitchburg, Wisconsin, United States), and their genomic DNA was isolated using the ReliaPrep^™ Blood gDNA Miniprep System (Promega, Fitchburg, Wisconsin, United States). Three nuclear genes (18S, 5.8S, and 28S), their internal transcribed spacers (ITS1 and ITS2), and the barcoding mitochondrial COI gene were PCR amplified. Primers and PCR conditions are provided in Supplementary Tables S2 and S3. PCR reactions were carried out with the GoTaq^® Long PCR Master Mix (Promega, Fitchburg, Wisconsin, United States), following the protocol described in our previous studies (Obert and Vďačný, 2019, 2020a,b; Obert et al., 2021). Bidirectional Sanger sequencing was conducted in Macrogen Europe B.V. (Amsterdam, Netherlands) on an ABI 3730 automatic sequencer.

Nuclear rRNA genes and their spacers were aligned according to the primary and the predicted secondary structures using the package 4SALE ver. 1.7.1 (Seibel et al., 2006). Secondary structures of 18S, 5.8S, and 28S rRNA molecules were constructed with R2DT (Sweeney et al., 2021), considering the models proposed for Tetrahymena thermophila (Lee and Gutell, 2012). Since helix 25 of the 28S rRNA molecule is hypervariable, its structure was adjusted with the free-energy minimization approach on the Mfold web server ver. 3.0 (http://www.unafold.org/) (Zuker, 2003). Modeling of ITS2 followed Obert and Vďačný (2020b). Secondary structures were plotted either in 4SALE (Seibel et al., 2006) or TRAVeLer (Elias and Hoksza, 2017). The helix number system of rRNA molecules was according to Lee and Gutell (2012) and Petrov et al. (2014). The mitochondrial protein-coding COI sequences were aligned with MEGA X (Kumar et al., 2018), using the protozoan mitochondrial genetic code (translation table 4) and the Muscle codon algorithm.

Nuclear rDNA and mitochondrial COI trees were constructed in a neighbor-joining (NJ), maximum likelihood (ML) and Bayesian framework. NJ analyses were conducted in MEGA X and included the maximum composite likelihood method, gamma-distributed rates among sites, a heterogeneous pattern among lineages, a pairwise deletion option to exclude alignment gaps, and 5,000 bootstrap pseudo-replicates. ML analyses were carried out in IQ-TREE ver. 1.6.10 (Nguyen et al., 2015), while Bayesian inferences (BI) were performed in MrBayes ver. 3.2.7 (Ronquist et al., 2012). Settings in ML analyses were as follows: (1) the best substitution model, as selected under the Bayesian information criterion, was assigned to each partition, (2) the edge-unlinked partition model that accounts for heterotachy (rate variation across sites and lineages) and allows each partition to have its own set of branch lengths, (3) thousand ultrafast bootstrap pseudo-replicates, and (4) the bnni algorithm to reduce overestimating bootstrap support. All other parameters were left default. Settings of Bayesian analyses were as follows: (1) prior parameters of evolutionary models as estimated with IQ-TREE, (2) model parameters were unlinked across partitions, (3) one million Markov chain Monte Carlo (MCMC) simulations, (4) a sampling frequency of trees and parameters at one hundred, and (5) a relative burn-in fraction of 25%. Convergence of MCMC analyses was confirmed with the in-built diagnostics of the program MrBayes. Trees were visualized in FigTree ver. 1.4.3.³

Species trees were calculated under the Bayesian multispecies coalescent model, using *BEAST ver. 1.8.3 (Heled and Drummond, 2010). Input files were prepared in BEAUti ver. 1.8.3 with the following settings: (1) best evolutionary substitution models as selected by IQ-TREE for each partition; (2) four categories for substitution rate heterogeneity; (3) uncorrelated lognormal clock; (4) ploidy scalars at 1.0 for the mitochondrial partition and 2.0 for the nuclear partition; (5) the Yule process model for the species tree prior; (6) piecewise constant population size, and (7) 100 million generations and a sampling frequency of 1,000 in Markov Chain Monte Carlo analyses. The convergence to stationary distribution (effective sample size > 200 for all parameters) was checked in Tracer ver. 1.6 (Rambaut et al., 2018). The maximum clade credibility tree was summarised in TreeAnnotator ver. 1.8.1 after discarding the first 10% of sampled trees.

Bayesian species delimitation was conducted in BP&P ver. 2.2 (Yang, 2015), with the same datasets as used in the construction of species trees. The maximum clade credibility tree obtained with *BEAST served as a guide tree for species delimitation. Each species delimitation model was assigned equal prior probability. Prior parameters for the ancestral population size θ and root age τ were estimated by running A00 analyses. The rjMCMC analyses were run for 100,000 generations with a sampling frequency of 2 and a burn-in of 10,000. A large fine-tuning parameter (ε = 15) was used to guarantee a good mixing in the reversible jump algorithm. All analyses were conducted twice to confirm consistency between runs.

To visualize conflict in the concatenated nuclear rDNA cistron and mitochondrial COI dataset, a neighbor-net analysis was undertaken with uncorrected distances in SplitsTree ver. 4 (Huson, 1998; Huson and Bryant, 2006). A species network was constructed in PhyloNet ver. 3.6.1, using the maximum pseudo-likelihood framework (Than et al., 2008; Wen et al., 2018). The network was computed from the nuclear rDNA and mitochondrial COI trees, allowing for a maximum of zero and ten reticulation nodes. Zero reticulations represented a null model corresponding to a species tree, while a maximum of 10 reticulations served to test how many reticulations could be present in the resulting phylogenetic networks. If there were exactly 10 reticulation nodes, another round of analyses with an increased number of reticulations would be needed. However, phylogenetic networks with a maximum of eight reticulations were recovered (see the ‘Results’ section) and hence no further analyses were needed. Each analysis was performed with ten runs and default settings, generating five optimal networks. The species networks were visualized with Dendroscope ver. 2.7.4 (Huson and Scornavacca, 2012).

To test whether hybridization could be the source of gene tree incongruence, the program JML ver. 1.3.1 was employed (Joly, 2012). This software uses a posterior distribution of species trees, population sizes, and branch lengths to simulate replicate sequence datasets under the coalescence with no migration. The minimum pairwise sequence distance between sequences of two species is evaluated on the simulated datasets and compared to the one estimated from the original data (i.e., from the rDNA cistron and COI datasets). JML analyses were run over 10,000 species trees sampled from the posterior distribution of *BEAST analyses. Settings were as follows: (1) relative mutation rate as estimated from the log file generated during the *BEAST analyses; (2) heredity scalars at 1.0 for the mitochondrial partition and 2.0 for the nuclear partition; and (3) best substitution models as selected under the Bayesian information criterion in IQ-TREE. All pairwise sequence distances with value of p < 0.05 were recorded as potential cases of hybridization.

The alpha-diversity of plagiotomids inhabiting the digestive tube of earthworms was studied at mesoscale, more specifically, at 36 localities in southwestern Slovakia, Central Europe (Figure 1). Altogether 880 specimens belonging to 19 earthworm species from three ecological (anecic, epigeic, and endogeic) groups were examined for the presence of plagiotomids (Supplementary Table S1). Only a single earthworm species from each ecological group was, however, colonized by plagiotomids. Namely, the anecic Lumbricus terrestris Linnaeus, 1758, the epigeic congener L. rubellus Hoffmeister, 1843, and the endogeic Aporrectodea tuberculata (Eisen, 1874).

According to the present phylogenetic analyses, as many as nine distinct mitochondrial lineages were recognized within the isolated plagiotomids. One lineage represents a new homogenous and molecularly well-delimited species, P. aporrectodeae sp. n., while the remaining eight lineages belong to the highly heterogenous P. lumbrici complex (see below). The former species occurred exclusively in the endogeic earthworm A. tuberculata and was detected at six out of the seven sites where this endogeic earthworm was sampled (Figure 1; Supplementary Table S1). Plagiotoma aporrectodeae sp. n. was never recorded in any Lumbricus species even when it co-occurred with A. tuberculata (Table 1).

For the sake of simplicity, mitochondrial lineages of the P. lumbrici complex were labeled by the first eight letters of the Greek alphabet (i.e., from alpha to epsilon). Lineages alpha, gamma, eta, and theta were found at nine localities but solely in L. terrestris; lineages beta, zeta, and epsilon were isolated only from L. rubellus collected just at a single locality; and lineage delta was detected in both Lumbricus species sampled at three localities. Theta is apparently the most widespread lineage, as it was recorded at seven out of the 10 sites where members of the P. lumbrici complex occurred. Just a single lineage of the P. lumbrici complex was found at six localities, while two lineages co-occurred at four localities. Interestingly, at a single spot (PUh), where both Lumbricus species carried plagiotomids, as many as seven lineages of the P. lumbrici complex were detected (Table 1).

To summarize, there is a correlation between the plagiotomid mitochondrial lineage and the earthworm species. Plagiotoma aporrectodeae sp. n. is present only in A. tuberculata, lineages alpha, gamma, eta, and theta exclusively in L. terrestris, and lineages beta, zeta, and epsilon solely in L. rubellus. On the other hand, there is no spatial relationship in the distribution of the plagiotomid lineages, causing the distribution of plagiotomid lineages to be mosaic-like.

Nuclear 18S rDNA (GenBank accession numbers OP538845–OP538959), ITS1-5.8S-ITS2 region and the first two domains of 28S rDNA (OP538960–OP539074) as well as the mitochondrial COI (OP562267–OP562381) sequences were obtained from 115 Plagiotoma specimens, which had been isolated from three earthworm species originating from 14 localities (Supplementary Table S5). Tree-building and network analyses of the nuclear rDNA cistron and the mitochondrial barcoding COI gene of all 115 Plagiotoma specimens revealed phylogenies with rather low resolution and multiple topological discrepancies (Figures 2A,B, 3A). Strong statistical supports and no conflicting signals were detected only in case of P. aporrectodeae sp. n. All specimens of this species shared identical rDNA cistron sequences and their COI variability ranged only from 0.00 to 0.58%. By contrast, much higher variability was detected within the P. lumbrici complex both in the rDNA cistron (0.00–1.88%, on average 0.73%) and in the COI gene (0.00–11.31%, on average 7.09%). The variability within the three rDNA regions was as follows: 0.00–1.86% (0.60% on average) in the 18S rRNA gene, 0.00–1.72% (0.89% on average) in the ITS1-5.8S-ITS2 region, and 0.00–2.31% (0.87% on average) in the first two barcoding domains of the 28S rRNA gene. Eight distinct mitochondrial lineages were robustly delimited within this species complex by all three tree-building methods (Figure 2B). Inter-lineage distances within the P. lumbrici complex ranged from 4.27 to 10.29% in COI sequences. The existence of these eight lineages was, however, either left statistically unsupported or they were paraphyletic in the nuclear rDNA cistron tree (Figure 2A).

The split decomposition diagram inferred from the concatenated rDNA cistron + COI dataset (115 specimens and 3,710 characters) by the neighbor-net algorithm displayed multiple reticulations across the eight lineages of the P. lumbrici complex. Individuals causing conflicting signals are placed close to the ‘core’ of lineages alpha, beta, and theta, and they are marked by blue arrows in Figure 3A. Conflicting positions of these specimens in the rDNA cistron and COI gene trees are marked by blue lines in Figures 2A,B. After exclusion of these specimens from the split decomposition analyses, the neighbor-net diagram became tree-like (Figure 3C), resembling the branching pattern suggested by the neighbor-joining rDNA cistron tree (Figure 2A). Gene trees built from the concatenated reduced dataset (99 specimens and 3,710 characters) were very well resolved and statistically supported (Figure 4A). Each lineage within the P. lumbrici complex received very high or full support. Lineage alpha grouped with P. aporrectodeae sp. n. with very strong support (100% NJ, 1.00 BI, 100% ML). Lineages beta, gamma, and delta formed a monophylum (96% NJ, 1.00 BI, 91% ML). Likewise, lineages eta, theta, zeta, and epsilon clustered together (100% NJ, 1.00 BI, 99% ML) but their interrelationships remain poorly resolved. A very similar picture was obtained also by coalescence-based analyses conducted in *BEAST (Figure 4B). Each plagiotomid lineage received a posterior probability of 1.00 in Bayesian species delimitation analyses conducted in BP&P. Individual lineages could not be, however, unambiguously delimited by barcoding analyses of the rDNA cistron per se, as intra- and inter-lineage distances partially overlapped (Figure 4C). On the other hand, a very distinct barcoding gap was revealed in the COI dataset. More specifically, the maximum intra-lineage distance was 2.33% while the minimum inter-lineage distance was 4.27% (Figure 4D).

Split decomposition (Figure 3A) and species network (Figure 3B) analyses of the 115 specimens dataset revealed the presence of non-historical signals, whose source might be incomplete lineage sorting, recombination, and/or hybridization. The exclusion of just 16 “problematic” specimens almost completely removed parallel splits from the neighbor-net diagram (see above and Figure 3C) and reduced the number of reticulations from eight to two in the species network (cp. Figures 3B with 3D). The two conflicting splits within the eta-theta-zeta-epsilon cluster are very likely due to the recombination within the nuclear rDNA cistron dataset (PhiTest, p = 0.0018). As expected, no recombination was detected within the mitochondrial COI dataset (PhiTest, p = 0.1261).

The presence of introgression during the evolutionary history of plagiotomids was tested also with the software JML ver. 1.3.1. This program uses posterior checking to test whether the minimum distance between sequences of two species is smaller than expected under a scenario that does not account for hybridization. Potential hybridization within the rDNA cistron dataset was detected between lineages beta and eta (p = 0.0150), beta and theta (p = 0.0224), beta and zeta (p = 0.0416), beta and epsilon (p = 0.0157), delta and theta (p = 0.0048), as well as between gamma and theta (p = 0.0048). As inferred from the COI dataset, lineage alpha was involved in potential hybridization with all other lineages of the P. lumbrici complex: alpha and beta (p = 0.0002), alpha and gamma (p = 0.0105), alpha and delta (p = 0.0002), alpha and eta (p = 0.0002), alpha and theta (p = 0.0187), alpha and zeta (p = 0.0263), alpha and epsilon (p = 0.0220).

To summarize, no traces of hybridization between P. aporrectodeae sp. n. and any lineage of the P. lumbrici complex were detected. On the other hand, all lineages of the P. lumbrici complex very likely experienced hybridization episodes during their evolutionary history. Despite that, only 16 out of the 115 specimens analyzed (13.91%) were identified to potentially have a hybrid origin.

As mentioned above, no variability was detected in the rDNA cistron of P. aporrectodeae sp. n., while unusually high variability (up to 1.88%) was observed within the P. lumbrici complex. Variable nucleotide positions were localized in the V2–V5, V7, and V9 regions of the 18S rRNA molecule (Figure 5). More specifically, out of the 42 polymorphic positions, 13 were situated in terminal loops of helices, six in bulges or adjacent positions, nine in single-stranded regions, and 14 in double-stranded regions. Mutations in the double-stranded regions were typically involved in compensatory base changes (CBCs) or retained helical structure when involved in non-canonical pairings (helix 30). Interestingly, two CBCs were found also in the V2 region mediating a long-range tertiary contact with helix 21es6d from the V4 region (Figure 5). Mutations tended to become accumulated in the single-stranded region of V2 and the double-stranded region of helix 44es12 of V9.

Only three variable positions were present in the 5.8S rRNA molecule (Figure 6, upper right panel). All were accumulated in the terminal loop of helix 9. Internal transcribed spacer 2 was similarly conservative, harboring only five polymorphic positions: one in the terminal loop of sub-helix B-3, one in the bulge of sub-helix B-4, and three in sub-helix B-2 (Figure 6, lower middle panel). There were as many as 31 mutations in the first two barcoding domains of the 28S rRNA molecule. All but one were situated in the first domain, whereby as many as 27 polymorphic sites were detected in the hypervariable helix 25. Interestingly, about 65% of these mutations were involved in CBCs (Figure 6).

The pronounced genetic variability provoked a question of whether plagiotomids isolated from the three ecological groups of lumbricid earthworms could be also delimited morphologically. However, any qualitative morphological feature that would enable their discrimination was not recognized. To assess the utility of morphometric features in delimitation of the three ecological groups of plagiotomids, a multidimensional statistical approach including principal component analysis (PCA) and elliptical Fourier analysis of body shape was utilized. Although some trends were recognizable in the PCA ordination diagram based on 21 morphometric characters, only plagiotomids living in the epigeic L. rubellus could be unambiguously separated from those occurring in the endogeic A. tuberculata (Figure 7A). The following morphometric features most contributed to the distinction of these two clusters along the first ordination axis: body length (loading 0.3295), body width (0.3047), and nuclear figure length (0.3168). Plagiotomids isolated from the anecic L. terrestris were partially mixed with those from L. rubellus and A. tuberculata. Nevertheless, L. terrestris-dwelling plagiotomids displayed a trend to separate from the other two groups along the second ordination axis, which reached the highest correlations with the length/width ratio of the anteriormost (0.3827) and posteriormost (0.4962) macronuclear nodule as well as with the ratio of the adoral zone to the body length (−0.3498). In other words, plagiotomids inhabiting the digestive tube of epigeic earthworms tend to have smaller and narrower bodies as well as a shorter nuclear figure and adoral zone of membranelles. On the other hand, plagiotomids colonizing endogeic earthworms are larger on average, possessing a higher number of adoral organelles and a longer nuclear figure. Plagiotomids living in anecic earthworms are morphometrically intermediate.

Shape analyses brought very similar results as did morphometric analyses. More specifically, plagiotomids isolated from the three ecological groups of earthworms were intermingled in the hierarchical cluster analysis based on Euclid distance and complete linkage algorithm (Figure 7B), and plagiotomids isolated from the epigeic L. rubellus could be discerned from those isolated from the endogeic A. tuberculata in the PCA diagram based on the Fourier coefficients (Figure 7C). Statistically significant differences among shapes of the three ecological groups of plagiotomids were also corroborated by MANOVA performed on the PCA objects (Hotelling-Lawley trace = 48.25, approximate F_{55, 382} = 67.03, p = 2.2e–16). Thus, plagiotomids isolated from epigeic earthworms tend to be narrower and more reniform, while those from endogeic earthworms are more stocky and oblong. The cell outline of plagiotomids living in anecic earthworms resembles that of plagiotomids inhabiting epigeic earthworms but is less reniform.

To summarize, the plagiotomid lineages as delimited by molecular data cannot be unambiguously recognized either by qualitative, morphometric (Figure 7A), or cell geometric data (Figures 7B,C). There is a continuum of variability among plagiotomids isolated from the three ecological groups of earthworms. Plagiotomids thus form a single highly variable cluster in the phenotypic space.

For over 180 years, all plagiotomids isolated from lumbricid earthworms were consistently identified as a single species, P. lumbrici. This obligate gut endosymbiont was reported from three biogeographic realms: the Palearctic (Dujardin, 1841; Pertzewa, 1929; Heidenreich, 1935; de Puytorac and Mauret, 1956; Dworakowska, 1966; Albaret, 1973, 1975), the Nearctic (Affaa et al., 2004), and the Neotropic (Cordero, 1928). Literature data thus suggests a wide distribution for P. lumbrici on one hand and low species differentiation within the genus Plagiotoma on the other hand. In accordance with this traditional view, the present multivariate statistical analyses of morphometric and cell geometric data did not reveal any well-delimited and discrete unites in the phenotypic space but a continuum of variability among plagiotomids isolated from the three ecological groups of lumbricid earthworms (Figures 7A–C). These findings, however, contrast with the high molecular diversity detected within the morphospecies P. lumbrici (Figures 5, 6) already at mesoscale, whereby as many as nine distinct mitochondrial lineages were identified in SW Slovakia (Figures 2B, 3C,D, 4A,B). We remain conservative and recognize, using Mayr’s (1995) biological species concept and Van Valen’s (1976) ecological species concept, only two discrete biological entities: the P. lumbrici complex comprising eight mitochondrial lineages and P. aporrectodeae sp. n. with just a single mitochondrial lineage. This, however, does not mean that the plagiotomid diversity is low. On the contrary, it is comparatively high but obscured by a complex evolutionary history shaped by secondary contacts of lineages whose speciation processes were not completed (see below).

Despite the hidden diversity, there are still pronounced differences (of at least one or two orders of magnitude) in alpha-diversities of plagiotomids and other endosymbiotic ciliate groups (nyctotherids, astomes, and hysterocinetids) living in the digestive tube of earthworms. A plausible explanation for this diversity paradox might be the differing mating-type system. Plagiotomids belong to the class Spirotrichea, which has a ‘high-multiple (open)’ type. This mating system strongly favors cross-mating phenomena, causing the genetic boundaries between species living in sympatry to be relaxed (e.g., Nobili et al., 1978; Valbonesi et al., 1988; Kuhlmann and Sato, 1993; Di Guiseppe et al., 2022). On the other hand, astomes and hysterocinetids belong to the class Oligohymenophorea whose mating system is ‘low-multiple (closed)’ (Orias et al., 2017). The closed system permits the establishment of strong species boundaries by strengthening incompatibilities among populations of different mating types. Mating incompatibilities fuel speciation processes by limiting gene flow, which supports evolutionary cohesion and ultimately leads to a comparatively low genetic variability within species. This hypothesis is corroborated by our observations on astomes (Obert and Vďačný, 2019, 2020a, 2021; Obert et al., 2021) and hysterocinetids (manuscript in preparation). By contrast to plagiotomids, individual species of astomes and hysterocinetids are homogenous in the nuclear as well as mitochondrial rDNA cistron and exhibit up to 1% variability in the barcoding COI gene. The ‘closed’ mating system thus might be responsible for the comparatively high alpha-diversity of astomes and hysterocinetids (cf. de Puytorac, 1954, 1969, 1972; Raabe, 1972; Obert et al., 2021), while the ‘open’ system might be accountable for the creation of highly variable complexes in plagiotomids.

The evolution of lumbricid earthworms has been largely driven by paleoclimatic and paleogeographic events (e.g., Pérez-Losada et al., 2011; Fernández et al., 2013; Domínguez et al., 2015; James et al., 2021). Central Europe is the ancestral area for a large portion of known lumbricid genera, whose diversification was very likely shaped by range expansion and retraction due to global climatic changes from the Paleocene–Eocene epochs to the Pleistocene (Marchán et al., 2021). Both epigeic and anecic earthworms evolved multiple times from endogeic ancestors, which suggests subsoil-dwelling, temporary and horizontally oriented burrow systems, and geophagy (feeding on soil) to represent the ground pattern of lumbricid earthworms (Domínguez et al., 2015). Interestingly, not only plagiotomids but also astomes inhabiting endogeic earthworms form deep-branching clades (Obert and Vďačný, 2019, 2020a, 2021; Obert et al., 2021). This suggests that endogeic earthworms might have served as ancestral hosts for intestinal ciliates and both anecic and epigeic earthworms were colonized later.

Plagiotomid mitochondrial lineages exhibit a distinct host-correlated pattern (Figures 4A,B,E), which also resembles the situation in astomes associated with lumbricid earthworms (Obert and Vďačný, 2019, 2020a, 2021; Obert et al., 2021). Diversification of plagiotomids was thus very likely governed by adaptive radiation in connection with the ecological group of their earthworm hosts. A similar diversification mode was suggested also for astome ciliates inhabiting lumbricid earthworms (Obert et al., 2021). The matter is, however, much more complex in plagiotomids than in astomes. The present phylogenetic trees indicate at least two independent transfers of plagiotomids from the anecic to the epigeic earthworms as well as multiple duplication events without host switching (Figures 4A,B). The latter coevolutionary event was very likely enabled by alternating glacial–interglacial periods, which led to the separation and reconnection of Lumbricus populations in the Palearctic. The approximately 41,000–100,000-year climatic cycles might have allowed the formation of distinct mitochondrial plagiotomid lineages in spatially isolated populations of earthworms but were not long enough to enable the evolution of strong species boundaries. After the secondary contact of earthworms during interglacial periods, hybridization of their plagiotomid endosymbionts might have taken place, as suggested by the present split decomposition (Figure 3A), species network (Figure 3B), and JML analyses. Introgression might have ultimately dissolved the emerging species boundaries and led to the collapse of multiple plagiotomid lineages into a single highly variable biological entity (cf. Bog et al., 2017; Kearns et al., 2018).

According to the time-calibrated phylogenies, Lumbricus and Aporrectodea s.s. are sister taxa that diverged about 65 Ma, whereby each genus further radiated about 20–40 Ma (Domínguez et al., 2015). According to the present phylogenetic and JML analyses, there are no signs of introgression between P. aporrectodeae sp. n. isolated from Aporrectodea and the P. lumbrici complex isolated from Lumbricus. Their separation is thus completed and very likely dates to the Aporrectodea–Lumbricus split in the Paleocene–Eocene. Plagiotoma aporrectodeae sp. n. differs from members of the P. lumbrici complex by 1.06–2.00% in the rDNA cistron, which corresponds to 1.58–2.99 × 10⁻⁴ substitutions per site per one million years. This clock rate matches quite well estimates for the 18S rRNA gene in various groups of ciliates, ranging from 1.24–3.96 × 10⁻⁴ substitutions per site per one million years (Wright and Lynn, 1997; Vďačný, 2015; Vďačný et al., 2019). On the other hand, the maximum distance in the rDNA cistron within the P. lumbrici complex is up to 1.03%, which corresponds to a divergence time of about 35 Ma. This estimate falls almost perfectly in the period of the Lumbricus radiation (Domínguez et al., 2015). Separation of the eight mitochondrial lineages within the P. lumbrici complex was, however, not completed during the Oligocene–Pleistocene epochs very likely due to the episodic reconnection of various Lumbricus populations/species during interglacial periods. Secondary contacts of earthworms might have enabled episodic hybridization of their plagiotomid symbionts, which in turn led to the weakening of the emerging species boundaries. Thus, the range dynamics of the earthworm hosts very likely strongly affected the diversification processes of their ciliate endosymbionts.

Hybridization is an important evolutionary process that affects diversification in various ways. It can accelerate diversification via formation of hybrid species (e.g., Stankowski and Streisfeld, 2015; Sefc et al., 2017; Ebersbach et al., 2020; White et al., 2020), but it can also slow down diversification via breakdown of species reproductive barriers, followed by fusion of previously separated lineages (e.g., Huxel, 1999; Gómez et al., 2015; Todesco et al., 2016). Organisms that were impacted by historical climate changes provide an excellent model system for studying how gene flow can drive but also prevent speciation. The evolution and distribution of lumbricid earthworms in the Palearctic were shaped by the Paleocene–Eocene Thermal Maximum and the Mid-Pleistocene Revolution, a fundamental change in the behavior of glacial cycles during the Quaternary glaciations (cf. Marchán et al., 2021). Repeated environmental oscillations thus might have enabled previously separated earthworm populations to come into secondary contact and hence allowed hybridization of their plagiotomid endosymbionts.

The present phylogenetic analyses revealed that eight out of the nine plagiotomid lineages show signs of hybridization (Figures 2A,B). All potential hybrids were detected only in two Lumbricus species, nevertheless, hybrids represented only 16 out of the 115 specimens analyzed (13.91%). The plagiotomid mitochondrial lineages differ by as much as 4.27–10.29% in the barcoding COI gene. Intra-species divergences in most ciliate species, however, typically range between 0 and 1% and only rarely are slightly higher than 1% (e.g., Kher et al., 2011, Doerder, 2014, 2019; Rataj and Vďačný, 2019, 2020, 2022; Obert et al., 2021; Zhang and Vďačný, 2021a,b, 2022). Doerder (2014) proposed a 4% interspecific divergence threshold for COI and Kher et al. (2011) even a more conservative 5% threshold. COI has been, indeed, recognized as a reliable barcoding marker in a variety of ciliates (e.g., Chantangsi et al., 2007; Strüder-Kypke and Lynn, 2010; Zhao et al., 2016; Park et al., 2019). In this light, all plagiotomid mitochondrial lineages could be classified as distinct molecular species. Despite this, we have found traces of hybridization and/or recombination among all lineages of the P. lumbrici complex (Figures 2A,B). These findings document that their speciation processes were not completed and introgression dissolved their emerging species boundaries. We, therefore, argue that gene exchange might have slowed down and ultimately hampered the diversification of Lumbricus-dwelling plagiotomids. In other words, hybridization led to their collapse into a single biological entity.

Hybridization is very likely responsible also for weakening the host structural specificity of plagiotomids by increasing their ecological amplitudes. This is documented by lineage delta, which was found both in L. terrestris and L. rubellus (Table 1). According to Albaret’s (1975) review, P. lumbrici has been so far isolated from four genera of the earthworm family Lumbricidae Rafinesque-Schmaltz, 1815: Lumbricus Linneus, 1758; Nicodrilus Bouché, 1972; Aporrectodea Orley, 1885; and Scherotheca Bouché, 1972. Specifically, it was recorded in four species of the genus Lumbricus [L. castaneus (Savigny, 1826); L. friendi Cognetti, 1904; L. rubellus Hoffmeister, 1843; and L. terrestris Linneus, 1758]; in two species of Nicodrilus [N. giardi (Ribaucourt, 1901) and N. longus (Ude, 1885), which is now considered a synonym of Aporrectodea longa (Ude, 1885)]; and in S. savignyi (Guerne and Horst, 1893). Although the molecular identity of plagiotomids reported in the literature is not known, all aforementioned reports are only from anecic and epigeic earthworms and hence might, indeed, refer to the P. lumbrici complex. This would suggest a rather broad host range and comparatively weak structural host specificity.

A further possible by-product of hybridization events might be the pronounced diversity (up to 1.88%) in the nuclear rDNA cistron within the P. lumbrici complex (Figures 5, 6). On the other hand, no signs of hybridization were noted in P. aporrectodeae sp. n. and no variability was detected in its rDNA cistron. Similarly, we did not observe any potential hybrids and any or very low variability in the highly conservative nuclear rDNA cistron or the hyper-variable COI gene in astomes and hysterocinetids (Obert et al., 2021). These observations indicate that introgression might increase molecular variability by overcoming concerted evolution, a process that converts copies of a gene in a multigene family into the same copy.