Paramecium Diversity and a New Member of the Paramecium aurelia Species Complex Described from Mexico

Paramecium (Ciliophora) is an ideal model organism to study the biogeography of protists. However, many regions of the world, such as Central America, are still neglected in understanding Paramecium diversity. We combined morphological and molecular approaches to identify paramecia isolated from more than 130 samples collected from different waterbodies in several states of Mexico. We found representatives of six Paramecium morphospecies, including the rare species Paramecium jenningsi, and Paramecium putrinum, which is the first report of this species in tropical regions. We also retrieved five species of the Paramecium aurelia complex, and describe one new member of the complex, Paramecium quindecaurelia n. sp., which appears to be a sister species of Paramecium biaurelia. We discuss criteria currently applied for differentiating between sibling species in Paramecium. Additionally, we detected diverse bacterial symbionts in some of the collected ciliates.


Introduction
Paramecium is one of the most studied genera of ciliates. Currently, at least fourteen morphological species of Paramecium are recognized as valid, and several more require reinvestigation [1][2][3]. Most of the morphological species include a number of genetically isolated groups, referred to as syngens. In some cases, syngens have been elevated to full species (the P. aurelia species complex, [4]) or could be described at least as genetic species due to the well-proved absolute reproductive barrier separating them, such as syngens in P. bursaria [5]. Additionally, several cryptic species of Paramecium have been reported [3], yet they are disputed as true species due to the failure to establish their laboratory culturing and to collect more specimens in nature. Finally, some contested Paramecium species were documented only once from specific localities and were not subjected to thorough morphological, physiological or molecular description [1,6,7].
Most of the Paramecium morphospecies have recognizable morphological traits, allowing easy and fast identification of new representatives isolated from nature [1]. Numerous molecular phylogenetic analyses of inter-and intraspecific diversity of Paramecium have been performed, and molecular barcoding of different genes has allowed researchers to discriminate between morphospecies when morphological features were blurred [8,9], or even between sibling species or syngens within morphological species [5,[10][11][12]. While the 18S rRNA gene has proven to be too conservative to disclose the intraspecies relationships in Paramecium, the cytochrome C oxidase subunit I (COI) gene is routinely used as a barcoding sequence in Paramecium molecular phylogenetic analysis [3,9,13,14]. This gene is sufficiently divergent and permits the inference of reliable genus level tree configuration and existence of intraspecific groups within each Paramecium morphospecies [14][15][16].
Extensive sampling, even in well-studied territories, allows researchers to find new Paramecium species or to confirm and validate species previously described [2,3,17,18]. Moreover, as many geographic regions remain poorly studied for Paramecium species occurrence, it is very possible that knowledge of the diversity of this genus is incomplete. Central America has not been surveyed systematically for Paramecium diversity. In this study, we provide new data on Paramecium occurrence and distribution in Mexico, including description of a new species of the P. aurelia complex found in two remote localities. Additionally, we report the presence of bacterial symbionts in some collected Paramecium strains since their occurrence has been considered as a possible criterion for Paramecium species [19].

Sampling and Maintenance of Paramecium Strains
About 130 freshwater samples were taken from more than 40 different waterbodies (natural and artificial lakes and ponds, canals, streams, drains, wetlands, and even a water reservoir in a roof-top garden) in several localities of seven states of Mexico: Ciudad de México, Estado de México, Hidalgo, Querétaro, Veracruz, Quintana Roo, and Yucatán in January-March 2019. The volume of each water sample collected was 10-30 mL. Samples from the same waterbody were taken at a distance of at least 20 m from each other and, thus, were considered as representing separate populations of ciliates. The samples were quickly transported to the laboratory and screened for paramecia within 24 h after sampling. Ciliates were detected under a Nikon SMZ 800 (Nikon Corporation, Tokyo, Japan) stereomicroscope. Then, all samples were kept on rice grains for 10-14 days and monitored every three days for previously unnoticed paramecia to show up. Several cells from each "positive" sample were isolated separately into depression slides; when possible, we isolated up to 10 cells from each sample, aiming to represent paramecia of different sizes and cell shapes. The ciliates introduced in culture were maintained on lettuce medium bacterized the day before use with Enterobacter cloacae, and supplemented with 0.8 mg/L of β-sitosterol (Merck, Darmstadt, Germany), as described earlier [19]. All currently alive strains used in the study are available upon request from RC CCM collection (World Data Centre for Microorganisms, RN 1171), Saint Petersburg State University, Saint Petersburg, Russia.

DIC Microscopy and Stainings
Live cells observations were made with differential interference contrast (DIC) microscopy with a Nikon Labophot-2 microscope equipped with a Nikon Digital Sight DS2Mv (Nikon Corporation) camera. We observed the cytological features important for quick species identification in Paramecium, namely cell size and shape, size, number and structure of micronuclei, structure of contractile vacuoles, and presence of algal symbionts [1]. Several staining techniques were employed, including the Feulgen procedure in De Lamater protocol, Harris hematoxylin, silver nitrate impregnation after Champy's fixation, silver carbonate and protargol [20]. Morphometric measurements were taken from stained cells.

Molecular Identification of Paramecium Strains and Bacterial Symbionts
Paramecium strains isolated from all samples and attributed to different morphospecies were subjected to sequencing of the mitochondrial COI gene. Additionally, in order to reconstruct a complete molecular phylogenetic tree of the P. aurelia species complex, COI gene sequences were obtained for the P. primaurelia strains Ir 4-2 (Russia) and FT11 (Pakistan), P. pentaurelia strains NR-2 (USA) and Nr1-9 (Russia), P. septaurelia strains 227 and 38 (USA). The total cell DNA was extracted from 100 to 200 cells of each strain using the GenElute Mammalian Genomic DNA Purification Kit (Sigma, Germany), according to the protocol "Genomic DNA from tissue". The PCR was performed using Encyclo Taq polymerase (Evrogen, Russia). The 767 bp-long partial COI gene sequences were amplified using the primers F388 and R1184, which are suitable for the majority of Paramecium species as described by Strüder-Kypke et al. [9]. For P. primaurelia, P. triaurelia, P. pentaurelia and P. septaurelia, the primers COI-long F (GATAAGGCTTGAGATGGCATACCCAGGAAG), and COI-longR (CAAAACCCATGTAAGCCATAACGTAGACAG) were designed, and 35 cycles of PCR were performed with an annealing temperature of 60 • C. Additionally, the partial sequences were obtained for the mitochondrial cytochrome C oxidase subunit II (COII) gene of some strains of P. biaurelia and presumably new species of the P. aurelia complex (see below; GenBank accession numbers MT318927-MT318930). In all cases, the same primers were used for PCR and sequencing. The partial 16S rRNA gene sequence for bacterial symbionts inhabiting P. putrinum strain K8 was amplified by PCR using the primers 16S alfa F19 and 16S alfa R1517, and sequenced with the primer 16S F343 [21]. All oligonucleotides were synthesized by Eurofins DNA (Germany). The PCR products were directly purified and sequenced at the Core Facility Center "Molecular and Cell Technologies" (St Petersburg State University, Saint Petersburg, Russia).

Molecular Phylogenetic Analysis
The COI gene sequences obtained in this study (GenBank accession numbers MT078136-MT078152, MT318931-MT318935) were aligned with COI gene sequences manually selected from GenBank or retrieved from ParameciumDB [22] for the strains with sequenced mitochondrial genomes [23]. Manual selection of entries was performed, since many COI gene sequences in GenBank are either too short, incorrectly assigned to a species or simply missing for some species in GenBank. The longer sequences were trimmed manually to obtain the 767 bp-long COI gene fragment, and the incomplete COI gene sequences from Genbank were chosen so that all sequences in the final set were at least 600 bp long. MUSCLE algorithm was used for alignment (online multiple alignment program [24]). Phylogenetic trees were constructed to infer the relationships of all isolated Paramecium strains using Phylogeny.fr [25,26]. The trees were computed by the bootstrapping procedure (500 bootstraps) and approximate likelihood ratio test method PhyML 3.1/3.0 aLRT [27]. The Maximum Likelihood analysis was performed with the HKY model.

Mating Tests
For the strains of the presumably new species of the P. aurelia complex, preparation of cell lines for mating tests and studies of autogamy were carried out by method of daily re-isolations [28]. The sexually reactive cultures were mixed with each other and also with P. biaurelia tester strains IST, Rieff, and Ts from the RC CCM collection. Strains Rieff and IST belong to the same mating type compatible with the Ts strain mating type. The conjugation was observed only between +18 • C and +21 • C. The conjugating couples were picked with the Pasteur pipette, then exconjugant cells were isolated into separate microaquariums, and F1 clones were established. The fragment of the nuclear mtB gene involved in mating type control [29] was amplified for F1 clones in 32 cycles of PCR using the primers bi-mtBF (GCACACCCTCTTAAAATAAGT) and bi-mtBR (AAATCTCGCAAACAACTACTG) with an annealing temperature of 55 • C. This fragment was sequenced using the same primers to confirm the heterozygosity of F1 clones (i.e., to confirm the exchange of pronuclei in conjugation), as allelic single nucleotide polymorphisms were visible on chromatograms as double peaks (data not shown). F2 progeny were obtained after autogamy of F1 clones, and survival rates of the F2 generation were counted. The F2 clones were considered as viable if they completed more than 6 divisions after autogamy, since old macronuclei remain functional during the first 5-6 vegetative divisions in Paramecium [19], and confer otherwise inviable cells to survive during that period.

Diversity of Paramecium and Its Bacterial Symbionts Revealed by Extensive Sampling in Several Regions of Mexico
From all sampled Mexican localities, paramecia were found in 19 waterbodies. Strains from 30 populations of six morphospecies, namely representatives of the P. aurelia species complex, P. jenningsi, P. caudatum, P. multimicronucleatum, P. bursaria, and P. putrinum were introduced in the laboratory cultures. The most frequently recorded species was P. multimicronucleatum (15 populations from 9 waterbodies). Six populations from different waterbodies contained representatives of the P. aurelia species complex, while other morphospecies were relatively rare ( Table 1). It was not uncommon to find several Paramecium species in the same community. The greatest diversity was detected in the lakes, ponds, and streams of the Cantera Oriente reserve (Mexico City), where we isolated a number of strains of P. caudatum, P. multimicronucleatum, P. bursaria, and P. putrinum. In several populations, some Paramecium specimens were inhabited by bacterial endosymbionts (Table 1, Figure 1). Most of these symbiotic bacteria still have to be identified and are currently being studied. Among the most interesting findings were presumably Trichorickettsia sp. abundantly present in host cytoplasm in several P. putrinum strains from Cantera Oriente ( Figure 1A), as well as unknown cytoplasmic ( Figure 1C,D) and intranuclear ( Figure 1E,F) symbionts in different strains of P. multimicronucleatum. Tiny bacteria able to produce R-bodies resembling Caedibacter sp. or Caedimonas sp. [30] were found in cytoplasm of P. tetraurelia cells from Xochimilco Lake (Mexico City), while other cytoplasmic symbionts were detected in the P. octaurelia strain from Cenote Azul (Quintana Roo).

Phylogenetic Analysis of the Collected Strains
The strain attribution to certain morphospecies by DIC microscopy of the living cells was confirmed by COI gene sequencing and further positioning of a strain within the Paramecium

Phylogenetic Analysis of the Collected Strains
The strain attribution to certain morphospecies by DIC microscopy of the living cells was confirmed by COI gene sequencing and further positioning of a strain within the Paramecium phylogenetic tree. Molecular characterization by COI and COII genes sequencing is the fastest and most reliable way to discriminate between different species of the P. aurelia complex [9,10]. All sibling species form separate branches on trees inferred from these molecular markers [10,31], and COI and COII gene barcodes make it possible to identify each species of the P. aurelia complex, eliminating the need for laborious round-robin mating tests with representatives of all species of the complex. COI gene sequencing can also reveal different haplotypes that cluster into intraspecific groups within P. multimicronucleatum [16,32] and into reproductively isolated syngens in P. bursaria [5]. However, the haplotypes revealed within P. caudatum do not form pronounced and well-supported branches [3,14]. Thus, analysis of COI gene sequences of Mexican Paramecium strains ( Figure 2) showed that P. bursaria isolated from the lakes in Cantera Oriente (Mexico City) and Amealco (Querétaro) belonged to syngen R3, which is known to be widespread in Far East Russia, China, Japan, and South America [5]. Paramecium multimicronucleatum strains sorted into three branches within this morphospecies cluster. Paramecium jenningsi from Mexico City grouped together with strains of the species found in Asia and Africa ( Figure 3). Finally, we succeeded in recovering four known species of the P. aurelia complex: P. primaurelia, P. triaurelia, P. tetraurelia, and P. octaurelia. Strains from two populations (Dinamos, Mexico City and Amealco, Querétaro) clustered in a separate branch as a sister species for P. biaurelia ( Figure 3). This suggested to us that we may have discovered a novel member of the complex, so, these strains were thoroughly studied in order to figure out if they actually represented a new species of the P. aurelia complex.
Diversity 2020, 12, x FOR PEER REVIEW 7 of 20 phylogenetic tree. Molecular characterization by COI and COII genes sequencing is the fastest and most reliable way to discriminate between different species of the P. aurelia complex [9,10]. All sibling species form separate branches on trees inferred from these molecular markers [10,31], and COI and COII gene barcodes make it possible to identify each species of the P. aurelia complex, eliminating the need for laborious round-robin mating tests with representatives of all species of the complex. COI gene sequencing can also reveal different haplotypes that cluster into intraspecific groups within P. multimicronucleatum [16,32] and into reproductively isolated syngens in P. bursaria [5]. However, the haplotypes revealed within P. caudatum do not form pronounced and well-supported branches [3,14]. Thus, analysis of COI gene sequences of Mexican Paramecium strains (Figure 2) showed that P. bursaria isolated from the lakes in Cantera Oriente (Mexico City) and Amealco (Querétaro) belonged to syngen R3, which is known to be widespread in Far East Russia, China, Japan, and South America [5]. Paramecium multimicronucleatum strains sorted into three branches within this morphospecies cluster. Paramecium jenningsi from Mexico City grouped together with strains of the species found in Asia and Africa ( Figure 3). Finally, we succeeded in recovering four known species of the P. aurelia complex: P. primaurelia, P. triaurelia, P. tetraurelia, and P. octaurelia. Strains from two populations (Dinamos, Mexico City and Amealco, Querétaro) clustered in a separate branch as a sister species for P. biaurelia ( Figure 3). This suggested to us that we may have discovered a novel member of the complex, so, these strains were thoroughly studied in order to figure out if they actually represented a new species of the P. aurelia complex. bursaria are included as outgroups. Groups I, II and III within P. multimicronucleatum (see Table 1) are indicated. The tree was computed by the bootstrapping procedure (500 bootstraps) and approximate likelihood ratio test method PhyML 3.1/3.0 aLRT. Numbers at nodes represent posterior probabilities higher than 0.4. The scale bar represents the branch length, corresponding to 0.3 substitution per site. The COI gene sequence accession numbers of the strains collected in this study are shown in blue.

Figure 2.
Phylogenetic position of Paramecium caudatum and Paramecium multimicronucleatum strains collected in this study on the mitochondrial COI gene tree. The sequences of P. chlorelligerum and P. bursaria are included as outgroups. Groups I, II and III within P. multimicronucleatum (see Table 1) are indicated. The tree was computed by the bootstrapping procedure (500 bootstraps) and approximate likelihood ratio test method PhyML 3.1/3.0 aLRT. Numbers at nodes represent posterior probabilities higher than 0.4. The scale bar represents the branch length, corresponding to 0.3 substitution per site. The COI gene sequence accession numbers of the strains collected in this study are shown in blue.

New Species of the Paramecium aurelia Complex, P. quindecaurelia n. sp.
The ciliates were collected from a small lake near the highway in Amealco, Querétaro (population A65) and from the lake drain in Los Dinamos, the National Park in the mountains on Mexico City territory (population D88). Several cells from both populations were isolated, and the obtained strains were introduced into laboratory cultures. They can be maintained in a variety of temperatures but demonstrated the best growth rate, three divisions per day, at 21-23 °C. The ciliates look like typical representatives of the P. aurelia complex (Figure 4). The cells are rather large, and maximum body length of 80 fixed specimens ranged from 104 to 168 µm with a mean length of 125 µm (16.9 µm standard error). Maximum body width of 80 specimens ranged from 19 to 39 µm with Figure 3. The phylogenetic tree of the Paramecium aurelia species complex inferred from mitochondrial COI gene sequences. Paramecium multimicronucleatum was used as an outgroup. The red arrow indicates position of P. quindecaurelia n. sp. The red brace shows the pair of sister species P. biaurelia and P. quindecaurelia n. sp.; the green brace shows the pair of sister species P. primaurelia and P. pentaurelia; the blue brace shows the pair of sister species P. tetraurelia and P. octaurelia. The tree was computed by the bootstrapping procedure (500 bootstraps) and approximate likelihood ratio test method PhyML 3.1/3.0 aLRT. Numbers at nodes represent posterior probabilities higher than 0.3. The scale bar represents the branch length, corresponding to 0.5 substitution per site. The COI gene sequence accession numbers of the strains collected in this study are shown in blue.

New Species of the Paramecium Aurelia Complex, P. Quindecaurelia n. sp.
The ciliates were collected from a small lake near the highway in Amealco, Querétaro (population A65) and from the lake drain in Los Dinamos, the National Park in the mountains on Mexico City territory (population D88). Several cells from both populations were isolated, and the obtained strains were introduced into laboratory cultures. They can be maintained in a variety of temperatures but demonstrated the best growth rate, three divisions per day, at 21-23 • C. The ciliates look like typical representatives of the P. aurelia complex (Figure 4). The cells are rather large, and maximum body length of 80 fixed specimens ranged from 104 to 168 µm with a mean length of 125 µm (16.9 µm standard error). Maximum body width of 80 specimens ranged from 19 to 39 µm with a mean width of 28 µm (8.0 µm standard error). The average kineties number was 55. The length of the infundibulum ranged from 15 to 36 µm with an average of 26 µm. All cells had one roundish macronucleus that was longer than wide ( Figure 4I), as well as two micronuclei of the vesicular type ( Figure 4G,H) typical for P. aurelia [1,4]. The average length of the macronucleus was 32 µm, and average width was 12 µm. Two contractile vacuoles were present per cell, and each vacuole had one pore and six or seven canals in different cells ( Figure 4E,F). The macronucleus is fragmented into 35-40 pieces during autogamy, and two macronuclear anlagen are formed in the cell ( Figure 4J). All morphometric data confirm that these strains are very similar to P. biaurelia ( [33], Supplementary Table S1).
Diversity 2020, 12, x FOR PEER REVIEW 9 of 20 a mean width of 28 µm (8.0 µm standard error). The average kineties number was 55. The length of the infundibulum ranged from 15 to 36 µm with an average of 26 µm. All cells had one roundish macronucleus that was longer than wide ( Figure 4I), as well as two micronuclei of the vesicular type ( Figure 4G,H) typical for P. aurelia [1,4]. The average length of the macronucleus was 32 µm, and average width was 12 µm. Two contractile vacuoles were present per cell, and each vacuole had one pore and six or seven canals in different cells ( Figure 4E,F). The macronucleus is fragmented into 35-40 pieces during autogamy, and two macronuclear anlagen are formed in the cell ( Figure 4J). All morphometric data confirm that these strains are very similar to P. biaurelia ([33], Supplementary  Table S1).  The COI gene sequencing was performed to assign A65 and D88 strains to a certain species of the P. aurelia complex. The strains from the same populations had identical sequences (Table 2), while between A65 and D88 strains the sequence was 99.3% similar (5 single nucleotide polymorphisms (SNPs) per 760 bp-long sequence). The best match in GenBank for both strains was KX008305.1, belonging to P. biaurelia. This sequence was 94.9% similar to A65 (39 SNPs per 760 bp-long sequence) and 94.5% similar to D88 (42 SNPs per 760 bp-long sequence). In the P. aurelia phylogenetic tree inferred from the COI gene sequences, strains A65 and D88 formed a separate branch sister to P. biaurelia (Figure 3, red brace), mirroring another pair of closely related species, P. primaurelia and P. pentaurelia (Figure 3, green brace). The COI gene sequence similarity between and within the latter two species is of the same range ( Table 3).
The first cells were entering autogamy only after 20 divisions in daily re-isolations cycle, while synchronous autogamy of the culture was observed after 25-27 vegetative divisions. Two macronuclear anlagen were formed after the sexual process, also typical for the P. aurelia species. We never observed selfing (intrastrain conjugation) in A65 and D88 strains, thus, the mating type determination is not stochastic as in some species of the complex [31]. Furthermore, we never obtained conjugation by mixing mature cells of different strains from both populations, so, presumably, all isolated strains were of the same mating type. In closely related P. biaurelia, the mating types are inherited maternally [19,29], and it is very likely that cytoplasmic inheritance is characteristic for these strains also.
Since we were unable to obtain conjugation between A65 and D88 strains, we also tested the possibility of conjugation between these strains and tester strains of P. biaurelia. Several conjugating couples were observed in crosses between IST and D88 and between Rieff and D88 strains after incubation of slightly starved reactive cells at + 18 • C, while only single couples were observed in mating tests of A65 with Rieff and IST strains. No conjugation was observed between Ts strain and either A65 or D88. In the control mating tests between P. biaurelia strains, conjugation was well pronounced at the same conditions. None of the few exconjugants in crosses utilizing A65 survived.
Exconjugant cells were isolated and F1 clones were established for the crosses IST × D88 and Rieff × D88. All F1 clones fit well and grew like healthy paramecia. F2 progeny were obtained by autogamy from all F1 clones. Survival rate of F2 clones was 3% (2 of 60) in the IST × D88 cross, and 8% (5 of 60) in the Rieff × D88 cross. At the same time, survival rate of F2 clones obtained in the IST intrastrain cross was 95%, and in the cross of P. biaurelia strains Rieff × Ts, it was 77%. Moreover, during the following month of maintenance, no F2 clones from the IST × D88 cross survived, and only two F2 clones remained viable from the Rieff × D88 cross. This rate of survival is very low and could be due to a technical error of the experiments. For example, some cells could have remained vegetative and were occasionally selected from autogamous culture of F1 heterozygous clones or some F1 cells might regenerated old macronuclei in autogamy.

Overview of Paramecium Diversity Revealed in Mexico
Screening of samples collected by us in Mexico showed unsurprisingly that Paramecium is widely represented in fresh waterbodies in this part of the world. The territory of Mexico, and in general, of Central America, remains unexplored by ciliatologists compared to Europe or some parts of the USA. Several Paramecium species, which we found in the current studies, were reported for the first time from these vast territories. To compensate the lack of knowledge about the Paramecium biogeography of Central America, we summarized the retrospective records of Paramecium according to [34][35][36][37][38][39][40][41][42][43][44][45][46][47]; all Paramecium records available from Mexico as well as Central and South America are shown in Table 4.
Besides the most important discovery of a new member of the P. aurelia species complex, it is worth noting several new facts on Paramecium biogeography. First is the occurrence of P. putrinum in the waterbodies of Mexico City. This species is not very frequent but can be considered common for temperate climate zones [1,11]. To our best knowledge, it has never been documented from tropical climates, but our observation is the southern-most collection of P. putrinum and may represent the southern extent of the species range. The same concerns P. triaurelia collected in Chapultepec Lake, Mexico City, as this species was not previously isolated as far south or in a tropical environment [4]. Record of P. jenningsi is notable, as this species is rare, and, unlike P. putrinum, usually inhabits waterbodies of tropical or subtropical zones [48]. This is only the third finding of P. jenningsi in the region, after Panama and Florida, USA [49]. Paramecium jenningsi strain DK from Mexico belonged to the same genotypic group as P. jenningsi strains collected in Japan and Madagascar [50]. Interestingly, P. bursaria strains found in two localities in Mexico belonged to syngen R3, common in Japan, China and also known from South America [5]. As for P. caudatum, all Mexican strains appeared to be related to each other and grouped in molecular phylogenetic trees with strains from Brazil, China and the European part of Russia (Figure 2), thus, confirming that the genotypes within these species do not show any special geographical pattern [3,14,32]. Paramecium multimicronucleatum strains found in Mexico belonged to three branches within this morphospecies. All branches include strains from all over the world.

Paramecium Quindecaurelia n. sp. and a Species Concept in Paramecium
Since Sonneborn [4] gave species rank to fourteen syngens of the P. aurelia complex, only one more member of this complex has been discovered: P. sonneborni [52]. Actually, it would be amazing if there were no more currently unknown sibling species of this complex in nature. Here, we report the sixteenth species of the P. aurelia complex, suggesting the name Paramecium quindecaurelia n. sp., as numerical tradition was interrupted with the description of P. sonneborni.
Morphologically, all species of the P. aurelia complex are indistinguishable [4], except P. sonneborni, which has a unique micronuclei morphology [52]. Morphometric characteristics do not show any significant variation among the sibling species of the complex [33]. The number of kineties is not considered important to differentiate between the P. aurelia species, while features of oral cortex are very conservative in Paramecium and cannot be used as a species characteristic within this genus [53]. Thus, morphological analysis of P. quindecaurelia n. sp., expectedly, revealed that it has no discriminating features and in particular, is extremely similar to its closest relative, P. biaurelia.
We reconstructed the complete phylogenetic tree for the P. aurelia species complex inferred from the COI gene sequence (Figure 3). Paramecium quindecaurelia n. sp. branches in the same cluster with P. biaurelia, but the COI gene identity between these two species does not exceed 95.0%, while within each of these species, it is not less than 99.3% ( Table 2). The comparison of COII gene sequences of P. quindecaurelia n. sp. and P. biaurelia strains gives very similar values: 94.8-95.1% between two species and 99.6-99.9% within each species (Supplementary Table S2). The phylogenetic distance between P. quindecaurelia n. sp. and P. biaurelia inferred from COI gene sequence analysis is of the same range as between the most closely related sister species P. primaurelia and P. pentaurelia (Figure 3, Table 3) where identity of the COI gene sequences is 94.0-94.4%. According to Sonneborn's data [19], P. primaurelia and P. pentaurelia were genetically isolated from each other but had identical isozyme patterns while all other species of the P. aurelia complex were characterized by unique zymograms. The only additional differences between these two species were unstable mating type O in P. pentaurelia, while in the clonal life of P. primaurelia the mating types never changed [19], and "a very weak and unreliable mating reaction, not leading to conjugation, occurred between type E of P. pentaurelia and type O of P. septaurelia", but the same was not true for P. primaurelia [4]. Thus, P. primaurelia and P. pentaurelia are also very similar physiologically to each other.
It is not surprising that P. quindecaurelia n. sp. is able to mate with P. biaurelia, and that P. primaurelia can conjugate with P. pentaurelia [19]. However, existence of the pronounced reproductive barrier between P. quindecaurelia n. sp. and P. biaurelia, despite the isolation is not absolute, further confirms that we found a novel member of the P. aurelia complex. The survival of F2 progeny in crosses between P. biaurelia and P. quindecaurelia n. sp. is slightly greater than zero. It is assumed that the reproductive species criterion is reliable but not absolute, as in many zoological species interspecific hybrids are known, and sometimes they also are considered as separate species (for example, water frogs from the Pelophylax esculentus complex- [54]). The F1 interspecies hybrids of P. biaurelia and P. quindecaurelia n. sp. cross are perfectly viable, but lethality in the F2 generation is very high, so the hybrids can be considered effectively sterile. Unfortunately, there are no data in the literature on F1 and F2 survival in crosses between P. primaurelia and P. pentaurelia. In the third pair of closely related species, P. tetraurelia and P. octaurelia (Figure 3), F1 hybrids had problems with survival and growth, and F2 post-autogamous progeny never survived [55]. Still, their COI gene sequence similarity is just 83.9%, and the phylogenetic distance between the latter two species is bigger than between P. primaurelia and P. pentaurelia or between P. biaurelia and P. quindecaurelia n. sp. (see Figure 3 and Tables 2  and 3). Different molecular markers have different resolution in Paramecium. For example, comparison of 18S rRNA gene provides good overview of the whole genus phylogeny [8,15], while it becomes useless if applied to the P. aurelia complex [8]. Mitochondrial COI and COII genes are rather conserved within P. aurelia sibling species (polymorphism in these genes among strains of the same species does not exceed 1-2%), while even the closest sibling species differ at least in 5% of these gene sequences ( [9,13] and this work). At the same time, in related to Paramecium genus Tetrahymena, the divergence for more than 1% in the COI gene sequence is considered reliable interspecific difference [56]. Obviously, only molecular phylogenomic analysis would allow us to understand the diagnostic level of molecular divergence between closely related species of ciliates.
In our opinion, at present, our data are sufficient to claim that P. quindecaurelia n. sp. is a separate species and not just a divergent group in P. biaurelia. Interestingly, P. biaurelia, the most common species of the P. aurelia complex in cold and moderate climate zones, at least in Europe [57], is not known from tropical environments [4,58], which may shelter its twin species, P. quindecaurelia n. sp. Further remarks. The strains of this species are able to conjugate with P. biaurelia strains. No endosymbionts have been detected in the species so far.

Conclusions
In this work, we confirmed that sampling in poorly studied regions, such as Central America, may reveal broad diversity of Paramecium, making it possible also to find new species. We found representatives of six Paramecium morphological species, and the collected strains belonged to different groups within some of these species (as in the case of the P. aurelia species complex and P. multimicronucleatum). Numerous molecular phylogenetic data demonstrate that each Paramecium morphospecies includes a number of intraspecific groups, which in some cases are known to correspond to syngens [5,11]. Genomic analyses showed that the P. aurelia complex emerged as a result of three whole genome duplications followed by species radiation [59], and speciation has been actively going on in this group [60,61]. It is known that the survival rate of F2 progeny in crosses of strains belonging to the same species, for example, in P. sexaurelia [62], can vary significantly. Inability to form conjugating couples seems to be one of the most limiting components of the reproductive barrier between sibling species of the P. aurelia complex. If couples can be formed, the possibility that a small proportion of the progeny from interspecies crosses may survive cannot be completely ruled out. Comparative genomics can potentially elucidate if all sibling species of the P. aurelia complex are absolutely isolated from each other. Pairs of twin species, like P. primaurelia / P. pentaurelia and P. biaurelia / P. quindecaurelia n. sp. may serve as the best models to study genetic isolation and gene flow between the P. aurelia species.