Phylogenetic study of several representative ciliate species in the Changbai Mountain region of Northeast China

Ciliates represent an exceptionally diverse group of unicellular eukaryotes. They are ubiquitously distributed across freshwater, marine, and soil habitats, where they fulfill critical ecological functions within microbial food webs^[1]. As model organisms, molecular phylogenetic studies of ciliates have provided key insights into fundamental biological questions, including microbial species delimitation, mechanisms of population differentiation, and evolutionary histories^[2,3]. Accumulating evidence indicates that ciliate population genetic structures exhibit two contrasting patterns. Some species display pronounced phylogeographic structuring, suggesting genetic divergence driven by geographic isolation, whereas others show genetic homogeneity across broad geographic ranges, reflecting either efficient gene flow or recent common ancestry^[4,5]. Nevertheless, most of these investigations have focused predominantly on temperate regions of Europe and North America, leaving the phylogenetic relationships of ciliate populations inhabiting Northeast China, particularly those occupying habitats with complex geological backgrounds, largely unexplored.

The Changbai Mountain region, situated on the border between China and North Korea, is home to the largest active volcano in China^[6]. Its volcanic landscapes were shaped by multiple episodes of volcanic activity spanning approximately the last five million years, making it the most geologically distinctive area in Northeast China. The region also serves as the headwater source of the Songhua River^[7]. Centered around the Tianchi crater lake, this area encompasses a diverse array of freshwater microhabitats, including volcanic barrier lakes, hot springs, and streams, exhibiting pronounced environmental heterogeneity, including extreme pH gradients, marked temperature variations, and a wide range of nutrient conditions^[8,9]. Such unique geological settings provide a natural laboratory for investigating the molecular evolution of ciliates under extreme environmental conditions. Current research on ciliates in the Songhua River Basin has largely focused on faunal surveys in the middle and lower reaches. Recent taxonomic and phylogenetic studies based on samples collected from these areas have revealed that the ciliate communities are dominated by three major classes: Oligohymenophorea (particularly the genera Tetrahymena, Glaucoma and Paramecium), Spirotrichea (represented by Euplotes, Stylonychia, and Uroleptus), and Colpodea (represented by Colpoda and Paracolpoda)^[10−13]. In contrast, the headwater region, the Changbai Mountain area, remains a blank slate with respect to molecular phylogenetic studies of key ciliate lineages. Although a recent morphological survey has identified a Bryometopus sp. in the Changbai Mountain region, hinting at an underappreciated reservoir of ciliate diversity, systematic investigations at the molecular level concerning the phylogenetic relationships and genetic structure of ciliates in this area are still lacking^[14].

Notably, the approximately 500 km of geographic distance between the headwater region of the Changbai Mountains and the Harbin section of the Songhua River Basin, together with the starkly contrasting environmental conditions between the volcanic headwater zone and the Songnen Plain, provide an ideal study system for examining whether volcanic landscapes are sufficient to drive population genetic differentiation in ciliates^[15]. Studies on other taxonomic groups have likewise demonstrated that volcanic activity can shape population genetic structure through multiple mechanisms. Volcanic eruptions may create physical barriers that directly impede population migration and gene flow. The habitat heterogeneity arising from volcanic landscapes may generate strong selective pressures, driving populations toward local adaptive evolution. Meanwhile, microhabitats that have remained spatially isolated over extended periods, such as crater lakes, create favorable conditions for the occurrence of genetic drift and allopatric speciation^[16,17]. Nevertheless, whether the aforementioned mechanisms apply to ciliates remains unclear. Although ciliates possess potential dispersal strategies such as resting cysts, their responses to extreme environmental heterogeneity and their capacity for long-distance gene flow remain uncertain^[18,19]. Against this background, the present study proposes the following hypothesis. The Changbai Mountains, as an active volcanic region with a long and dynamic geological history, would act as a 'species pump' driving population differentiation in ciliates. Consequently, ciliate populations in the volcanic headwater region of the Changbai Mountains and those in the Harbin section of the Songhua River Basin should exhibit significant genetic differentiation, forming distinct geographic lineages in the phylogenetic tree, with genetic distances between the two populations being substantially greater than the genetic variation observed within each population.

This study focuses on the Changbai Mountain region with the aim of systematically obtaining molecular systematics information on ciliates in this area. Samples were collected from the Changbai Mountain headwater section and the Harbin section of the Songhua River, encompassing 31 species across three classes: Oligohymenophorea, Spirotrichea, and Colpodea. Through amplification and analysis of SSU rRNA and mtSSU rRNA gene sequences, this study seeks to reconstruct the higher-level phylogenetic relationships of the target lineages and evaluate their patterns of intraspecific differentiation. For the representative species shared between the two regions, population-level phylogenetic analyses based on SSU rRNA were conducted for six species: Paramecium caudatum, P. primaurelia, Euplotes octocarinatus, E. daidaleos, Colpoda inflata, and C. maupasi. Haplotype network analyses based on mtSSU rRNA were performed for P. caudatum, P. primaurelia, E. chongmingensis, E. daidaleos, C. grandis, and C. inflata (the latter two replacing E. octocarinatus and C. maupasi due to data availability). By comparing genetic distances between populations from the two regions, as well as between these populations and those from other regions worldwide, this study provided new data supporting the role of volcanic landscapes in the molecular evolution of ciliates and laid a foundation for understanding the mechanisms underlying the formation of the ciliate fauna in the Songhua River Basin.

A total of 50 populations representing 31 species belonging to 14 genera and three classes were collected in this study. The sampled taxa included the genera Paramecium, Anteglaucoma, Glaucoma, Uronema, Colpidium, and Lembadion (class Oligohymenophorea); Euplotes, Holostichides, Hemiamphisiella, Hemiurosomoida, Stylonychia, and Uroleptus (class Spirotrichea); together with Colpoda and Tillina (class Colpodea). Among these, specimens of Tillina and some Colpoda populations were isolated from soil samples collected in Erdaobaihe Town, Antu County, Jilin Province, China. The remaining Colpoda populations together with all other species from the other classes were collected from freshwater habitats in Baishan City (Jilin Province), Anda City (Heilongjiang Province), and Harbin City (Heilongjiang Province), China. Detailed sampling information is provided in Supplementary Table S1, and sampling locations are shown in Supplementary Fig. S1. Soil ciliates were enriched using the non-flooded Petri dish method at room temperature (ca. 25 °C), with wheat grains added to promote bacterial growth. Other species were maintained in their native habitat water and temporarily kept in Petri dishes at 25 °C, with wheat grains also added as a nutrient source to stimulate bacterial proliferation. The infraciliature was revealed by silver carbonate staining and protargol staining^[20,21]. Terminology and classification mainly followed Gao et al.^[22] and Lynn^[23].

DNA extraction, amplification, and sequencing

For each species, one to five cells from clonal culture were isolated under a stereomicroscope using a micropipette, washed with distilled water at least thrice to remove potential contaminants, and then incubated in non-nutrient distilled water for 6 h to ensure complete removal of residual food particles. Cells were then transferred to an Eppendorf tube in a volume of no more than 5 μL distilled water. Genomic DNA was extracted using the DNeasy Blood & Tissue Kit (QIAGEN, Germany, supplied by Shanghai) following the manufacturer's instructions.

SSU rRNA and mtSSU rRNA genes were amplified by the polymerase chain reaction (primers see Supplementary Table S2). Taq polymerase (Takara Ex Taq; Takara Bio Inc., Otsu, Japan) was used for PCR amplification to minimize experimental errors. PCR conditions for SSU rRNA gene amplification were as follows: initial denaturation at 94 °C for 5 min; followed by five cycles of denaturation at 94 °C for 30 s, annealing at 56 °C for 1 min 45 s, and extension at 72 °C for 2 min; then another 25 cycles of denaturation at 94 °C for 45 s, annealing at 60 °C for 1 min 45 s, and extension at 72 °C for 2 min; with a final extension at 72 °C for 8 min. The mtSSU rRNA gene was amplified under the following conditions: initial denaturation at 94 °C for 5 min; five cycles of denaturation at 94 °C for 45 s, annealing at 58 °C for 1 min 45 s, and extension at 72 °C for 2 min; then another 25 cycles of denaturation at 94 °C for 45 s, annealing at 60 °C for 1 min 45 s, and extension at 72 °C for 2 min; with a final extension at 72 °C for 10 min. The PCR product purification was performed using the TIANgel Midi Purification Kit (Beijing, TIANGEN BIOTECH, China), cloned using the PMD 18-T vector cloning kit (Takara Biomedicals, China), and a randomly selected clone was sequenced bidirectionally in Shanghai Sangon Biological Engineering and Technical Service Company (Shanghai, China). Sequenced contigs were assembled with SeqMan v.7.1.0 (DNAStar).

Datasets and alignments

The sequences used in this study included newly determined sequences as well as other related ciliate sequences obtained from the National Center for Biotechnology Information (NCBI) database (accession numbers are shown in Figs. 1−3). To ensure the reliability of subsequent phylogenetic analyses, downloaded sequences that were ambiguous in morphology, unknown in origin, or excessively short were excluded from further analyses. The SSU rRNA gene sequences, including those from the three class-level datasets used for higher-taxon phylogenetic reconstruction, as well as the population sequences of six representative species used for intraspecific genetic differentiation analyses, were aligned using the Clustal W program with default parameters implemented in BioEdit 7.0.1^[24]. To eliminate ambiguous alignment sites, all datasets were manually checked and edited using BioEdit 7.0.1. The sequences used for phylogenetic analysis were compiled into four datasets: (1) the Oligohymenophorea SSU rRNA dataset, containing 105 taxa and 1,781 sites; (2) the Spirotrichea SSU rRNA dataset, containing 102 taxa and 1,960 sites; (3) the Colpodea SSU rRNA dataset, containing 83 taxa and 1,749 sites; and (4) six population-level SSU rRNA datasets for Paramecium caudatum, P. primaurelia, Euplotes octocarinatus, E. daidaleos, Colpoda inflata, and C. maupasi for intraspecific genetic differentiation analyses. Phylogenetic trees for Oligohymenophorea and Spirotrichea were rooted with three Colpodea species as outgroups, while the phylogenetic tree for Colpodea was rooted with three species from Oligohymenophorea as outgroups.

In addition to the SSU rRNA datasets, three independent datasets were compiled for haplotype network analyses based on the mtSSU rRNA gene. These datasets included: (1) a dataset for P. caudatum and P. primaurelia, comprising eight taxa and 990 sites; (2) a dataset for E. chongmingensis and E. daidaleos, comprising seven taxa and 740 sites; and (3) a dataset for C. grandis and C. inflata, comprising eight taxa and 1,036 sites. Sequences from both the Changbai Mountain and Harbin populations were included for all six species. All sequences were aligned using Clustal W implemented in BioEdit 7.0.1, followed by manual adjustment. As the analyses focused on intraspecific haplotype relationships, no outgroup sequences were included in these datasets.

Phylogenetic analyses

Maximum likelihood (ML) analyses were conducted based on the three aforementioned alignments with 1,000 bootstrap replicates, using RAxML-HPC2 on XSEDE 8.2.12 via the CIPRES Science Gateway under the GTRGAMMA substitution model^[25,26]. Bayesian inference (BI) was employed for the identical alignments with MrBayes 3.2.7 on XSEDE on the CIPRES Science Gateway, under the GTR + I + G model determined by the Akaike Information Criterion (AIC) in MrModeltest 2.2^[27,28]. Markov chain Monte Carlo (MCMC) simulations were executed for 10,000,000 generations, with tree sampling conducted every 100 generations, and the initial 10,000 trees were discarded as burn-in. To investigate intraspecific genetic differentiation among six representative species, phylogenetic trees were reconstructed based on SSU rRNA gene sequences from different populations of each species. Neighbor-joining (NJ), maximum parsimony (MP), and maximum likelihood (ML) analyses were conducted using MEGA12 with default settings, and nodal support was assessed using 1,000 bootstrap replicates^[29].

All resulting phylogenetic trees were visualized using MEGA12. The bootstrap values interpretation follows Vďačný & Foissner^[30], that is, considering values ≥ 95 as high, from 71 to 94 as moderate, from 50 to 70 as low, and < 50 as no support^[31]. Bayesian posterior probability < 0.95 is considered low, and values ≥ 0.95 are considered high^[32].

Haplotype networks

Because mtSSU rRNA sequences were not available for all species, the haplotype network analyses included a slightly different set of species from the intraspecific phylogenetic analyses based on SSU rRNA. Specifically, haplotype networks based on mtSSU rRNA gene sequences were constructed for six representative species shared between the Changbai Mountain and Harbin regions: Paramecium caudatum, P. primaurelia, Euplotes chongmingensis, E. daidaleos, Colpoda grandis, and C. inflata. E. octocarinatus and C. maupasi were not included due to the lack of available mtSSU rRNA sequences. The median-joining method implemented in PopART 1.7 was used, and mutations are represented as line segments in the networks^[33,34].

Due to the limited number of mtSSU rRNA sequences available per species (typically three to four haplotypes per species), conventional haplotype network analyses for individual species would provide limited information. Therefore, we constructed combined networks for closely related species pairs to visually compare the magnitude of intraspecific haplotype variation against interspecific divergence. This approach allows a direct assessment of whether intraspecific genetic variation approaches or exceeds the level of interspecific differentiation.

A total of 62 new sequences have been deposited in GenBank, including 15 SSU rRNA and four mtSSU rRNA gene sequences from the class Oligohymenophorea, 26 SSU rRNA and three mtSSU rRNA gene sequences from the class Spirotrichea, and nine SSU rRNA and five mtSSU rRNA gene sequences from the class Colpodea. Accession numbers, length, and GC content of each sequence are provided in Supplementary Table S1.

Phylogenetic analysis of the class Oligohymenophorea based on SSU rRNA gene

The ML and BI trees exhibited largely congruent topologies; only the ML tree is presented, with support values from both methods indicated at the nodes (Fig. 1). Phylogenetic analyses revealed that the class Oligohymenophorea is divided into two major monophyletic lineages: clade I comprises Hymenostomatia and Peritrichia (86% ML, 0.99 BI), whereas clade II consists of Scuticociliatia and Apostomatia (64% ML, 1.00 BI).

Figure 1. 

The maximum-likelihood (ML) tree based on the SSU rRNA gene of major members of the class Oligohymenophorea. Newly added sequences in this study are bolded in red type. Node support is shown as: ML bootstraps/BI posterior probability. Dashes (−) indicate the discordance in branching patterns between ML and BI. Fully supported (100%/1.00) branches are marked with solid circles. The scale bar corresponds to ten substitutions per 100 nucleotide positions. Pentagon reflect populations obtained from the upper reaches of the Songhua River, while triangles reflect populations obtained from the lower reaches of the Songhua River.

Within Hymenostomatia, four families were recovered: Tetrahymenidae, Turaniellidae, Glaucomidae, and Ichthyophthiriidae. Turaniellidae and Ichthyophthiriidae form a sister group nested between three Tetrahymenidae assemblages (borealis group, australis group, and paravorax group). The newly sequenced Colpidium sp. clusters with full support with the Changbai Mountain population of C. striatum (100% ML, 1.00 BI), and together with another Colpidium sp. sequence (87% ML, 0.90 BI), they group with other C. striatum populations (76% ML, 0.67 BI). Two Anteglaucoma harbinensis populations from Changbai Mountain form a highly supported sister group with the Harbin population (99% ML, 1.00 BI), and this clade is sister to a clade containing two Glaucoma sp. sequences (86% ML, 1.00 BI). Within Scuticociliatia, Pleuronematida and Philasterida are monophyletic and sister to each other; the newly sequenced Uronema nigricans clusters with full support with previously sequenced conspecific populations (100% ML, 1.00 BI). Peniculia comprises five families (Parameciidae, Frontoniidae, Paranassulidae, Stokesiidae, and Lembadionidae) forming a well-supported monophyletic group. For P. caudatum, the Changbai Mountain population (pop1) is sister to the Harbin population (87% ML, 1.00 BI), then to pop3 (88% ML, 1.00 BI), and together with pop2 (96% ML, 1.00 BI) they cluster with populations from Japan and the USA (100% ML, 1.00 BI). P. primaurelia from Changbai Mountain clusters with the Harbin population (100% ML, 1.00 BI), and further with other populations (81% ML, 0.98 BI). P. bursaria from Changbai Mountain clusters with populations from Harbin, Japan, and elsewhere (100% ML, 1.00 BI). Three newly sequenced Paramecium species (P. multimicronucleatum, P. fokini, and P. putrinum) and Lembadion lucens each cluster with their respective previously sequenced counterparts.

Phylogenetic analysis of the class Spirotrichea based on the SSU rRNA gene

The ML and BI trees yielded congruent topologies; only the ML tree is presented, with support values from both methods indicated at the nodes (Fig. 2). Within Spirotrichea, Oligotrichia and Hypotrichia form a sister group that clusters closely with Euplotia. Branching successively, Phacodiniidia and Protohypotrichia constitute a sister clade outside the remaining assemblage. Diverging basally, Licnophoria, represented by Licnophora cf. lyngbycola, is recovered as sister to all other subclasses.

Figure 2. 

The maximum-likelihood (ML) tree based on the SSU rRNA gene of major members of the class Spirotrichea. Newly added sequences in this study are bolded in red type. Node support is shown as: ML bootstraps/BI posterior probability. Dashes (−) indicate the discordance in branching patterns between ML and BI. Fully supported (100%/1.00) branches are marked with solid circles. The scale bar corresponds to ten substitutions per 100 nucleotide positions. Pentagon reflect populations obtained from the upper reaches of the Songhua River, while triangles reflect populations obtained from the lower reaches of the Songhua River.

Euplotia comprises two monophyletic orders, Euplotida and Discocephalida. Within Euplotida, Euplotidae groups with Certesiidae, forming a tight cluster with Aspidiscidae, whereas Gastrocirrhidae is positioned outside and Uronychidae is sister to these clades.

All newly sequenced Euplotes populations clustered with their conspecific or closely related congeners. Clustering with full support (100% ML, 1.00 BI) was E. parkei from the upper and lower Songhua River basin. E. daqingensis pop2 from Harbin first formed a sister group with a Daqing population (61% ML, 0.63 BI), subsequently clustering with E. daqingensis pop1 from Changbai Mountain (100% ML, 1.00 BI). Also recovered with full support was the grouping of E. heilongjiangensis from Changbai Mountain with the Harbin population (100% ML, 1.00 BI). Two E. euryhalinus populations (Changbai Mountain and Harbin) clustered together (100% ML, 1.00 BI), forming a sister group to other E. euryhalinus sequences (99% ML, 1.00 BI). E. chongmingensis pop1–2 from Changbai Mountain clustered with two Harbin populations (95% ML, 1.00 BI), a clade that was sister to a clade comprising E. chongmingensis pop3 from Anda and the Shanghai population (99% ML, 1.00 BI). Full support (100% ML, 1.00 BI) was also observed for the clustering of three E. octocarinatus populations from Changbai Mountain with E. octocarinatus pop4 from Harbin, with this clade further grouping with other previously sequenced conspecific populations (98% ML, 1.00 BI). Positioned outside this cluster were four E. daidaleos sequences, among which E. daidaleos pop1–2 from Changbai Mountain formed a highly supported sister pair (97% ML, 1.00 BI) and further clustered with the Harbin population (73% ML, 0.76 BI). E. encysticus pop1 from Harbin clustered with the Shenzhen population and an Euplotes cf. encysticus isolate (96% ML, 1.00 BI), a clade that subsequently grouped with E. encysticus pop2–3 from Changbai Mountain (100% ML, 1.00 BI). Finally, E. aediculatus clustered with full support with its conspecific populations (100% ML, 1.00 BI).

Within Hypotrichia, the monophyly of Stichotrichida, Sporadotrichida, and Urostylida was not supported. The Stylonychia ammermanni isolate from Changbai Mountain formed a highly supported sister group with other known strains of this species (95% ML, 1.00 BI), a clade that was recovered as sister to S. lemnae (86% ML, 1.00 BI) and subsequently clustered with S. harbinensis from Harbin (57% ML, 0.89 BI). Hemiamphisiella sp. clustered with H. terricola with moderate support (64% ML, 0.90 BI). The Changbai Mountain population of Hemiurosomoida longa clustered with full support with its previously sequenced conspecific counterpart (100% ML, 1.00 BI). Holostichides sp. formed a moderately supported sister group with three other species within the family Urostylidae (84% ML, 1.00 BI). The newly sequenced Changbai Mountain population of Uroleptus stueberi clustered with high support with previously reported U. stueberi sequences (99% ML, 1.00 BI), whereas a Changbai Mountain population of U. longicaudatus grouped with its published counterpart with moderate support (86% ML, 0.96 BI).

Phylogenetic analysis of the class Colpodea based on SSU rRNA gene

Maximum likelihood (ML) and Bayesian inference (BI) analyses were performed on the class Colpodea based on SSU rRNA gene sequence data, yielding nearly identical topologies. The ML tree is presented in Fig. 3, with support values from both methods indicated at the nodes. All four orders within Colpodea were recovered as monophyletic. A clade was formed by Colpodida and Cyrtolophosidida, with Bursariomorphida as its sister, whereas Platyophryida occupied the basal position within the class. All nine newly sequenced species were nested within the core of the Colpodea clade. For C. inflata, the two newly sequenced Changbai Mountain populations (pop1–2) were found to form a sister clade with the Harbin population (55% ML, 0.70 BI), and were further clustered with C. inflata pop3 from Changbai Mountain with high support (99% ML, 1.00 BI). Subsequently, this clade was recovered as sister to a clade consisting of the Harbin C. inflata population and C. cf. inflata (89% ML, 1.00 BI). The newly sequenced C. grandis from Changbai Mountain and the Harbin population were observed to form a sister group (59% ML, 0.89 BI). This clade was later grouped with full support (100% ML, 1.00 BI) into a clade composed of C. reniformis and C. henneguyi. Notably, C. maupasi pop4 from Harbin was found to cluster with C. ecaudata, albeit with low bootstrap support (31% ML). This clade was recovered as sister to C. maupasi pop1–3 from Changbai Mountain (63% ML, 0.73 BI), whereas the Japanese population of C. maupasi was placed outside this cluster with full support (100% ML, 1.00 BI). The newly sequenced Tillina sp. was clustered with T. minima with full support (100% ML, 1.00 BI).

Figure 3. 

The maximum-likelihood (ML) tree based on the SSU rRNA gene of major members of the class Colpodea. Newly added sequences in this study are bolded in red type. Node support is shown as: ML bootstraps/BI posterior probability. Dashes (−) indicate the discordance in branching patterns between ML and BI. Fully supported (100%/1.00) branches are marked with solid circles. The scale bar corresponds to ten substitutions per 100 nucleotide positions. Pentagon reflect populations obtained from the upper reaches of the Songhua River, while triangles reflect populations obtained from the lower reaches of the Songhua River.

Intraspecific phylogenetic relationships and genetic distances of six representative species

Phylogenetic trees were constructed for six representative species (Paramecium caudatum, P. primaurelia, Euplotes octocarinatus, E. daidaleos, Colpoda inflata, and C. maupasi) using ML, NJ, and MP methods. The three methods yielded largely congruent topologies; only the ML trees are presented, with support values from all three methods indicated at the nodes (Fig. 4). For five of the six species (P. caudatum, P. primaurelia, E. octocarinatus, E. daidaleos, and C. maupasi), populations from Changbai Mountain and Harbin consistently clustered together, forming a sister lineage without clustering with other geographic populations (Fig. 4a–f). This clustering pattern received robust support for P. caudatum (87/93/87) and E. octocarinatus (100/100/100), while moderate support was observed for P. primaurelia (78/98/65), E. daidaleos (88/55/83), and C. maupasi (78/67/58).

Figure 4. 

Maximum likelihood phylogenetic trees based on SSU rRNA gene sequences for six representative species from the classes Oligohymenophorea, Spirotrichea, and Colpodea, showing intraspecific relationships among geographic populations of (a) P. caudatum, (b) P. primaurelia, (c) E. octocarinatus, (d) E. daidaleos, (e) C. inflata, and (f) C. maupasi. Different geographic populations are color-coded. Nodal supports are shown in the order ML/NJ/MP; dashes (−) indicate nodes where the topology of the ML tree is inconsistent with those inferred by other methods. The scale bar corresponds to one substitution per 100 nucleotide positions.

In contrast, C. inflata exhibited a distinctive phylogenetic pattern. Unlike the other five species, its Guangzhou population did not form an independent clade but was nested within the clade comprising the Changbai Mountain and Harbin populations. The Canadian population clustered with another Harbin population, forming a sister clade to the aforementioned clade, whereas the Russian population was positioned at the base of this clade.

BLAST searches revealed that the newly obtained SSU rRNA gene sequences of the Changbai Mountain population of P. caudatum pop1–3 exhibited the highest similarity (99.7%–99.8%) to the Harbin population (PP236958). For E. octocarinatus, the three newly obtained sequences from the Changbai Mountain population pop1–3 showed the highest similarity with the Harbin population pop4, reaching 99.6%–99.8%. In contrast, when compared with other published E. octocarinatus sequences in GenBank, the highest similarity was only 99.0%, which was with the Mexico population. Notably, the three SSU rRNA gene sequences of the C. inflata Changbai Mountain population pop1–3 showed 98.7%–99.6% similarity with the Harbin populations (MZ557833, MZ557838, and MZ557804), whereas the similarity with the Guangzhou population reached as high as 99.9%.

To further elucidate the degree of genetic differentiation between upstream and downstream populations in the Songhua River basin, three species with high support values in the phylogenetic tree—P. caudatum, E. octocarinatus, and C. inflata—were selected for in-depth intraspecific genetic distance analysis. To assess whether the unique volcanic landscape of the Changbai Mountains has led to genetic differentiation between ciliate populations in this region and those from the Harbin section of the Songhua River basin, we compared SSU rRNA gene genetic distances between the Changbai Mountain and Harbin populations, as well as between these Chinese populations and conspecific populations from other countries (Brazil, Japan, and Italy). The genetic distance matrices (Tables 1−3) revealed that the genetic distances between the Changbai Mountain and Harbin populations of P. caudatum and E. octocarinatus were extremely small (P. caudatum: 0.001–0.004; E. octocarinatus: 0.005–0.007), markedly lower than those between these populations and foreign populations from Brazil, Japan, and Italy (the latter ranging from 0.004–0.019 and 0.011–0.020, respectively). In contrast, the genetic distance between the Changbai Mountain and Harbin populations of C. inflata was larger, ranging from 0.009 to 0.019, which is higher than that observed between its populations from other regions.

Haplotype network analysis

To investigate population structure, haplotype networks were constructed for six representative species—P. caudatum, P. primaurelia, E. chongmingensis, E. daidaleos, C. grandis, and C. inflata—based on their mtSSU rRNA gene sequences (Fig. 5a–c). No haplotype sharing was observed, and only mild genetic differentiation was detected among geographic populations, with no strong geographic structuring. As shown in Fig. 5a, P. caudatum and P. primaurelia were highly divergent (131 substitutions). For P. caudatum, the Harbin and Changbai Mountain populations differed by only two substitutions, and the Japanese population differed by one to three substitutions. For P. primaurelia, the Changbai Mountain and Harbin populations differed by three substitutions, while the British population differed by three to six substitutions, indicating slightly greater but still low divergence. A similar pattern was observed in Fig. 5b. E. chongmingensis and E. daidaleos exhibited high interspecific divergence (185 substitutions). Within E. chongmingensis, the Changbai Mountain and Harbin populations differed by two to eight substitutions, suggesting slight differentiation. In contrast, E. daidaleos showed greater divergence between the two regions (32–43 substitutions). In Fig. 5c, C. grandis and C. inflata were separated by 122 substitutions, while C. cf. inflata differed from C. inflata haplotypes by 54–61 substitutions. For C. grandis, the Changbai Mountain and Harbin populations differed by only two to four substitutions. For C. inflata, a more pronounced divergence was observed (eight–11 substitutions), though no strong geographic isolation was evident overall.

Figure 5. 

Haplotype networks based on mtSSU rRNA gene sequences for three ciliate classes: (a) Oligohymenophorea (P. caudatum and P. primaurelia), (b) Spirotrichea (E. chongmingensis and E. daidaleos), and (c) Colpodea (C. inflata and C. grandis). All sequences represent unique haplotypes (no shared haplotypes) and are labeled as 'species abbreviation_gene_sequence number'. Colored circles indicate different geographical origins. Red nodes represent unsampled intermediate haplotypes (median vectors). Numbers on branches indicate nucleotide substitutions (only shown for branches with > 10 substitutions).

In this study, we conducted SSU rRNA gene sequence analyses of 31 ciliate species collected from the Changbai Mountains–Songhua River basin in Northeast China, providing the first molecular data for representatives of the classes Oligohymenophorea, Spirotrichea, and Colpodea from this region. The topologies recovered from our newly generated molecular data are highly congruent with those of previous studies, further corroborating the robustness of the higher-level phylogenetic framework^[10−12,35].

Within the Oligohymenophorea, the newly sequenced species—Colpidium sp., Anteglaucoma harbinensis, Paramecium multimicronucleatum, P. fokini, and P. putrinum—all nest within their respective lineages with strong nodal support. This result further supports the division of the Oligohymenophorea into two major lineages: one comprising the Hymenostomatia and Peritrichia, and the other encompassing the Scuticociliatia and Peniculia^[36]. Within the Hymenostomatia, the family Glaucomidae is nested within the Tetrahymenidae, rendering the latter paraphyletic^[37]. Moreover, the monophyly of the Peniculia—comprising the Parameciidae, Frontoniidae, Paranassulidae, Stokesiidae, and Lembadionidae—is strongly supported^[13,22,38].

In the Spirotrichea, the nine newly sequenced Euplotes species, along with six representative species from the genera Stylonychia, Hemiamphisiella, Hemiurosomoida, Holostichides, and Uroleptus, all cluster with their conspecific or congeneric counterparts with high support, providing molecular evidence for the distribution of these species in the Songhua River basin. The phylogenetic topology shows that Licnophoria forms an early-diverging independent lineage at the base of the Spirotrichea, whereas Oligotrichia and Hypotrichia emerge as sister groups, a topology consistent with the most recent reconstructions^{[3,11,22,39,40]}. The monophyly of the Stichotrichida, Sporadotrichida, and Urostylida is not supported, reinforcing the consensus that these groups require taxonomic revision^[11,22,40].

Within the Colpodea, the newly sequenced Colpoda inflata, C. grandis, C. maupasi, and Tillina sp. each cluster with high support with conspecific individuals, attesting to the reliability of the sequence data. Phylogenetic analyses support the monophyly of the orders Colpodida, Cyrtolophosidida, Bursariomorphida, and Platyophryida, with the Platyophryida occupying a basal position^[12,14,41]. Bardeliella pulchra represents the earliest-diverging lineage within the Colpodida, a topology consistent with previous studies^[42]. C. maupasi and C. ecaudata cluster with species of Exocolpoda and Ropoma, while being distantly placed from other congeners on the phylogenetic tree, a pattern also reported in earlier work^[12]. Notably, the phylogenetic position of the Harbin population C. maupasi pop4 could not be stably resolved; it groups with C. ecaudata with only 31% bootstrap support. Previous studies have suggested that species of the genus Colpoda may have undergone rapid radiation diversification, which may partly explain the low support at certain nodes and the topological instability observed in this study^[41]. This low support also indicates that SSU rRNA gene sequences have limited resolution for distinguishing closely related species within the genus Colpoda, and further hints at the possible existence of cryptic species diversity within this genus.

Intraspecific genetic differentiation and haplotype network analysis of upstream and downstream populations along the Songhua River

Through intraspecific phylogenetic analyses, genetic distance estimation, and haplotype network construction for six selected species, including Paramecium caudatum, P. primaurelia, Euplotes octocarinatus, E. daidaleos, Colpoda inflata, and C. maupasi, this study revealed varying degrees of genetic differentiation between populations from the upper reaches of the Songhua River, specifically the Changbai Mountains, and the lower reaches, namely Harbin. Phylogenetic and haplotype network analyses indicated that the extent of genetic differentiation between Changbai Mountains and Harbin populations differs among species, as shown in Figs. 4 and 5. Genetic distance estimates further corroborated these observations: the genetic distances between Changbai Mountains and Harbin populations of P. caudatum and E. octocarinatus were remarkably low, with values of 0.001–0.004 and 0.005–0.007, respectively, markedly lower than those between these populations and conspecifics from Brazil, Japan, and Italy. These results reveal that the genetic differentiation between the Changbai Mountain and Harbin populations is species-specific. For P. caudatum and E. octocarinatus, despite the considerable geographic distance of approximately 500 km and the starkly contrasting environmental conditions between the two regions, no substantial genetic differentiation was observed between upstream and downstream populations. This phenomenon likely reflects the high dispersal capacity of these species, particularly their planktonic or epiphytic lifestyles, as well as the effective transport of resting cysts via water currents, wind, or animal activities within the river basin, which may reflect potential dispersal-mediated genetic connectivity between upstream and downstream populations^[43,44]. Consistent with this interpretation, extensive gene flow across broad geographical ranges has been documented in oligotrich ciliates^[45]. As a key survival strategy enabling ciliates to withstand unfavorable environmental conditions, resting cysts facilitate long-distance dispersal. In contrast, ciliate species incapable of forming resting cysts, such as Tetrahymena thermophila, exhibit restricted distribution and limited gene flow^[46]. Among the species examined in this study, C. inflata exhibited relatively higher genetic divergence within the Songhua River basin, suggesting a comparatively lower dispersal ability relative to the other five species or more pronounced population structure. In summary, for species with high dispersal capacity like P. caudatum, the species pump effect may be counteracted by efficient gene flow, leading to genetic homogenization; for species with low dispersal capacity like C. inflata, the same volcanic landscape may promote genetic differentiation by creating habitat heterogeneity and dispersal barriers. Thus, the Changbai Mountains may act as a species pump for ciliates in the Songhua River basin, with its effects manifested in a species-specific manner depending on the dispersal capacity of each species.

Among the six species examined, the upstream-downstream genetic distances of P. caudatum and E. octocarinatus, which ranged from 0.001 to 0.004 and from 0.005 to 0.007, respectively, were lower than that of C. inflata, which ranged from 0.009 to 0.019. This interspecific variation may be linked to differences in ecological habits: P. caudatum and E. octocarinatus primarily inhabit freshwater planktonic or epiphytic environments, and their resting cysts may be more readily dispersed passively by water currents^[47−49]; by contrast, C. inflata occupies both soil and freshwater habitats, and its dispersal may rely more on soil particles or animal vectors, potentially resulting in lower dispersal efficiency^[30]. Nevertheless, given that only two species per class were included in this study and that intraspecific sample sizes are limited, these interpretations remain preliminary. The observed patterns of genetic differentiation should not be generalized to the class level without further evidence.

Widespread genetic homogeneity and phylogeographic pattern of C. inflata

This study showed that Colpoda inflata had a relatively higher genetic differentiation between upstream and downstream populations compared to the other five species within the Songhua River basin. However, this does not preclude its ability for long-distance dispersal. As detailed below, C. inflata exhibits widespread genetic homogeneity. Among the six representative species examined, C. inflata exhibited a distinctive phylogeographic pattern. Phylogenetic analyses revealed that the Guangzhou population nests within the Changbai Mountains–Harbin clade, the Russian population occupies a basal position within this clade, and the Canadian population forms a sister group to the Harbin population (Fig. 4). In the haplotype network, the Changbai Mountains and Harbin populations differ by eight–11 mtSSU substitutions, with an SSU rRNA gene genetic distance of 0.009–0.019, which is higher than those observed for P. caudatum and E. octocarinatus.

These results are consistent with previous studies. A PCR-RFLP analysis based on 15 globally distributed isolates revealed extremely low genetic differentiation in C. inflata, with 0.0038–0.033 substitutions per nucleotide site, and haplotypes randomly distributed across the globe^[50]. An expanded study comprising 114 isolates further confirmed the absence of geographical clustering of ribotypes^[51]. Physiologically, resting cysts of C. inflata exhibit strong tolerance to desiccation for up to six weeks and to freezing conditions^[52]; more broadly, cysts of the genus Colpoda are capable of withstanding a wide range of extreme stresses^[53].

In the present study, the basal position of the Russian population and the clustering of the Canadian population with the Harbin population indicate that geographical distance is not the sole determinant of genetic relatedness. Instead, these patterns likely reflect stochastic long-distance dispersal events or post-glacial range expansion-mediated gene flow^[54]. This widespread genetic homogeneity is thus best explained as a consequence of random dispersal driven by highly resistant resting cysts^[55,56]. By contrast, the congeneric species C. maupasi did not exhibit a similar pattern in this study; its northeastern populations form an independent lineage distinct from populations in Japan. This contrast may reflect differences in dispersal capacity or population history between the two Colpoda species, although it could also stem from limited global sampling of C. maupasi, warranting further investigation with broader geographic coverage.

Impact of Changbai Mountains volcanic activity on genetic differentiation of ciliates

We hypothesized that the Changbai Mountains would act as a 'species pump' driving population differentiation in ciliates. Accordingly, significant genetic differentiation was expected between upstream (Changbai Mountains) and downstream (Harbin) populations. However, our results did not support this hypothesis: the six representative species exhibited extremely low genetic distances between the two regions (Paramecium caudatum: 0.001–0.004; Euplote octocarinatus: 0.005–0.007; Tables 1 and 2). Even for Colpoda inflata, which showed relatively higher divergence (0.009–0.019; Table 3), its genetic distance was far smaller than the expected value corresponding to a geographical distance of approximately 500 km, indicating that geographical distance did not lead to significant genetic differentiation.

Several speculative explanations may account for this finding. First, the strong dispersal capacity of ciliate resting cysts may counteract the isolating effects of volcanic topography^[43,44]. Second, the most recent large-scale eruption occurred approximately 1,000 years ago (946–947 CE)^[56]; ongoing gene flow since then may have counteracted the isolating effects of volcanic activity. Furthermore, the Songhua River system may facilitate gene flow between upstream and downstream populations.

Notably, the 'species pump' hypothesis has been documented in the freshwater amphipod Gammarus nekkensis^[57]. In contrast, the ciliates in the present study did not exhibit a similar pattern, likely reflecting a fundamental difference in dispersal capacity between microorganisms and macroorganisms. Macroorganisms are constrained by geographical barriers, whereas ciliates, equipped with highly resistant resting cysts, may achieve stochastic long-distance dispersal, allowing gene flow to counteract volcanic isolation^[55,56]. For most ciliate species, the Songhua River system may function as an integrated ecological unit, with ongoing gene flow potentially maintaining genetic homogeneity across the basin, although this inference awaits validation through more direct estimates of gene flow.

Limitations and future perspectives

Several limitations of this study should be acknowledged. First, only two species per class were included, which limits the generalizability of interspecific comparisons and constrains the extension of our conclusions to the class level. Second, the molecular markers employed in this study, namely SSU rRNA and mtSSU rRNA gene, are both conserved markers with limited mutation rates, which may impede the detection of recent population differentiation. Third, certain species, including Paramecium primaurelia, Euplotes daidaleos, and Colpoda maupasi, exhibited low nodal support in phylogenetic analyses, and therefore the relevant conclusions require further validation. Fourth, global sampling coverage is limited, particularly with the absence of samples from regions such as Africa and South Asia, which restricts a comprehensive understanding of global phylogeographic patterns.

Future research should integrate multi-gene data, including the LSU rRNA gene, which provides higher phylogenetic resolution than SSU rRNA, and COI, which serves as an effective DNA barcode for species-level identification, thereby achieving finer-scale resolution of population genetic structure. In addition, cross-continental systematic sampling is recommended, along with comparative physiological experiments, to thoroughly investigate the relationship between dispersal capacity and genetic differentiation in ciliates. Furthermore, targeted sampling of extreme microhabitats, such as volcanic hot springs in the Changbai Mountains, may reveal specialized lineages adapted to extreme environments and provide critical insights into the impact of volcanic activity on microbial evolution.

In this study, we analyzed SSU rRNA gene sequences of 31 ciliate species from the Changbai Mountains-Songhua River basin, as well as mtSSU rRNA gene sequences of six representative species. The phylogenetic topologies of Oligohymenophorea, Spirotrichea, and Colpodea were highly consistent with previous studies, validating the current higher-level taxonomic framework. For most species, no significant genetic differentiation was detected between the Changbai Mountain populations and the Harbin populations in the Songhua River basin, with only Colpoda inflata showing relatively higher genetic differentiation. Colpoda inflata exhibited widespread genetic homogeneity, supporting stochastic long-distance dispersal via resting cysts. Unlike macroorganisms, the ciliates in this study do not support the 'species pump' hypothesis for the Changbai Mountains, likely because the high dispersal capacity of their resting cysts counteracts the isolating effects of volcanic topography. This study provides fundamental molecular data for ciliates in this region and highlights the importance of integrating microbial life history traits into biogeographical research.