Pseudokeronopsis carnea and P. flava cells were isolated from a freshwater pond in Baihuayuan Park (36°04′N, 120°22′E), Qingdao, China. Species were initially determined by morphological features and later confirmed by sequencing their SSU-rRNA genes. A single cell was picked, washed, and cultured in flasks using filtered and autoclaved pond water. Cells were incubated with rice grains at 23°C for 21 days, then collected using a glass micropipette under a stereomicroscope. Genomic DNA was extracted using the MagAttract HMW DNA kit (QIAGEN, #67563, Germany). A DNA library was constructed with NEBNext DNA Library Prep Master Mix Set for Illumina (NEB, United States) following the manufacturer’s instructions. RNA extraction was performed with the RNeasy Plus Mini Kit (Qiagen, Germany) following the manufacturer’s instructions. The RNA libraries were generated using NEBNext Ultra RNA Library Prep Kit for Illumina (NEB, United States) following the manufacturer’s instructions.

Pair-end 150 bp sequencing reads were performed on the Illumina Hiseq 2500 platform, producing 20 Gb and 10 Gb of clean data for the DNA and RNA libraries, respectively. Genomes were assembled using SPAdes v3.12 (Bankevich et al., 2012) (parameters: -k 21,33,55,77 –careful). Contigs with low coverage (< 2×) or small size (< 200 bp) were removed. QUAST v5.0.2 (Gurevich et al., 2013) was used to measure genomic statistics including GC content and N50. RSEM v1.3.3 (Li and Dewey, 2011) was used to calculate the sequencing depth of each contig. For homologous gene annotation, genomic contigs were aligned with protein sequences from the SWISS-PROT database using BLASTX v2.3.0 (Camacho et al., 2009) (parameters: evalue = 1e-5, querygenecode = 6). Gene model annotation information was extracted using a custom Perl script and used to train AUGUSTUS v2.5.5 (Stanke et al., 2006) (parameters: –species = pseudokeronopsis –min_intron_len 15,39) for gene model prediction. RNA-seq reads were assembled into transcripts using rnaSPAdes v3.11.1 (Bushmanova et al., 2019) and aligned with the genome assembly by BLAT v3.6 (Kent, 2002) to optimize the gene models. Predicted genes without start and stop codons were filtered out using a custom Perl script. The RNA-seq reads were mapped to genome contigs using Tophat2 v2.0.10 (Kim et al., 2013). The mapped read count of each gene was measured by featureCounts v1.6.1 (Liao et al., 2013). Potential cis-regulatory sequence motifs were searched within subtelomeric regions using MEME v5.3.3 (Bailey et al., 2015). The sequence motifs identified were visualized using WebLogo 3 (Crooks et al., 2004). Frequency of stop codon usage (TAA, TGA, and TAG) was measured from the homolog sequence alignment between the CDS or transcript sequences of each species and the ciliate protein library using BLASTX v2.3.0 (parameter: evalue = 1e-5), as previously described (Pan et al., 2019).

A total of 238 orthogroups were identified among 31 ciliates (see Supplementary Table 2) using OrthoFinder and were aligned using mafft (Emms and Kelly, 2019). The concatenated ortholog sequence alignment dataset was used for phylogenomic analysis on CIPRES Science Gateway server v3.3 (Miller et al., 2010). RAxML-HPC2 v8.2.9 (Stamatakis, 2014) under LG model of amino acid substitution (Γ distribution + F, four rate categories, 1,000 bootstrap replicates) was used to perform maximum likelihood (ML) analysis. PhyloBayes MPI 1.5a (Lartillot et al., 2009) (CAT-GTR model + Γ distribution, four independent chains, 4,000 generations with 10% burn-in, convergence Maxdiff < 0.3) was used to perform Bayesian inference (BI) analysis. The phylogenetic tree was visualized using MEGA v7.0.20 (Kumar et al., 2016). The time of speciation was estimated using r8s (Sanderson, 2003) and corrected using calibration times obtained from the TimeTree database (Hedges et al., 2006). Computational analysis of gene family evolution (CAFE) (De Bie et al., 2006) was performed to identify gene families that have undergone significant expansion or contraction. Gene families were annotated against the Gene Ontology (GO) database using InterProscan (Jones et al., 2014). R package clusterprofiler (Yu et al., 2012) was used to conduct an enrichment analysis of expanded and contracted gene families.

Using high-throughput sequencing data, we assembled the MAC genomes of two Pseudokeronopsis species (Table 1). The genome assemblies of P. carnea and P. flava are 76.8 Mb and 52.5 Mb in size, respectively, in line with other known ciliate genomes (Figure 1). Most of the contigs bear telomeric repeats (C₄A₄) on at least one end (P. carnea, 81.1%; P. flava, 82.8%). The average size of contigs capped with telomeric repeats on both ends is 1.7 kb (P. carnea) and 1.1 kb (P. flava), and the vast majority (P. carnea, 97%; P. flava, 96%) appear to be gene-sized chromosomes (Figures 2A,B), consistent with the nanochromosome architecture found in other spirotrich genomes. Predicted gene numbers in P. carnea and P. flava are 12,734 and 9,520, respectively. 84% of P. carnea and 78% of P. flava genes are annotated using SWISS-PROT or NR databases. Introns in these species are tiny (only 36 nucleotides on average for both P. carnea and P. flava), a feature also observed in several other ciliate genomes (Figure 2C). The subtelomeric regions between the transcription start site (TSS) or transcription end site (TES) and the adjacent 5′ telomere repeat are also short (Figure 2D), similar to previous observations in the nanochromosomes of Oxytricha, Euplotes and Strombidium (Kim et al., 2013; Swart et al., 2013; Chen et al., 2019). Overall, the combined genomic features of Pseudokeronopsis represent an extremely compact eukaryotic genome architecture, with single genes containing minimal introns, nested between short subtelomeric regions (Figures 2D,E).

To perform the phylogenomic analysis, we collected public genomic/transcriptomic datasets available for 31 ciliates (Supplementary Table 2), and identified 238 orthologs among P. carnea, P. flava, and these species. The system assignment of species we describe here based on maximum likelihood (ML) and Bayesian inference (BI) methods generally agrees with previous studies (Gentekaki et al., 2017; Chen et al., 2018). Phylogenomic analysis of the P. carnea and P. flava sequenced in the current study shows full support for their cluster with two previously reported Pseudokeronopsis species (Figure 2F). The analysis also supports a larger cluster that includes the spirotrichs Oxytricha trifallax and Stylonychia lemnae. Standard stop codons in ciliates are frequently reassigned to code for amino acids (Swart et al., 2016). Similar to these two species, the stop codons “TAA” and “TAG” are consistently reassigned in all Pseudokeronopsis species, leaving “TGA” as the only stop codon (Figure 2F and Supplementary Table 2), reflecting the close evolutionary relationship between these groups. Interestingly, compared to previous phylogenomic analysis result in which the litostomatean Entodinium is not included (Chen et al., 2018), we find that the clustering of class Litostomatea and class Spirotrichea is poorly supported in the current BI tree, and not supported at all by the ML tree (Litostomatea clusters with Colpodea first after including Entodinium). This opens the possibility that the assignment of litostomateans could require further review in the future by expanding sampling.

Although Pseudokeronopsis nanochromosomes have limited space, we found an AT-rich region exists in the 5′ subtelomeric regions of most complete nanochromosomes that carry single genes (92.7 and 98.5% for P. carnea and P. flava, respectively) (Figure 3A). Further analysis 5′ subtelomeric regions reveals a 15 nt putative CRE in the sequence motif W_nAWTW_n which is positioned between the transcription start site (TSS) and translation start codon in both species (Figure 3B). Similar CREs were identified in the 5′ subtelomeric regions on the nanochromosomes of Strombidium (Li et al., 2021) and Nyctotherus (McGrath et al., 2007). We identified an additional AT-rich element in the Pseudokeronopsis species, this time in the 3′ subtelomeric regions, with the sequence motif “TTNATTTCNTTAA.” This sequence is found between the translation stop codon and the transcription end site (TES) (Supplementary Figure 1A), and is distinct from the 5′ subtelomeric sequence at the other end.

To demonstrate that these potential CREs are conserved evolutionarily, we investigated the base substitution rate of nucleotides in the 5′ subtelomeric regions. The base substitution rate in the putative CRE is 50% lower than that in the entire subtelomeric region (48.2 and 55.2% for P. carnea and P. flava, respectively), indicating that nucleotides in the CRE are under stronger selection. The substitution pattern also shows a distinct difference between the 5′ CRE and surrounding nucleotides (Figure 3C and Supplementary Figure 1B). Compared with the entire subtelomeric region, the substitution rate of A-to-G and T-to-C is greatly reduced in the CRE, but dramatically increased in C-to-A substitutions, indicating that G/C nucleotides in the CRE are consistently being replaced by A/T nucleotides.

Although not positioned upstream of the TSS, as is seen for TATA-box-like elements in other eukaryotic organisms like yeast (Lin et al., 2010), these sequences may still act as non-canonical regulatory CREs that bind transcriptional trans-activating factors. To determine whether this motif is related to transcription initiation, we compared the expression of genes on nanochromosomes that either possess or lack this CRE in their 5′ subtelomeric regions. We observed that genes with this CRE have significantly higher transcription levels in both species (Figure 3D). The majority of genes without the 5′ CRE are silent and the putative CRE is more associated with active genes, which could be the source of this transcription activity difference (Figure 3E). On the contrary, most of the genes with the adjacent CRE are actively transcribed (94.2 and 90.4% for P. carnea and P. flava, respectively). These observations suggest that transcription initiation on Pseudokeronopsis nanochromosomes depends on this CRE near the transcription start site, though the nature of this CRE is not clear. The sequence may act as a promoter by binding directly to a transcription factor, or it may contribute a necessary structural feature to the DNA in this region. Future studies should further test whether chromatin accessibility is greater at this location, and if active chromatin marks are enriched (Sheng et al., 2021). A similar association between the 3′ CRE and gene transcription was also identified (Supplementary Figures 1C,D). Together with the conserved base substitution patterns, our results reveal a strong evolutionary selection pressure upon the AT-rich CRE, and suggest it plays a functionally important role in transcription regulation. Considering the compact nanochromosome architecture, the inclusion of these sequences in the primary transcript UTRs, and reassignment of stop codons in these species, it is also possible that these CREs assist in translation initiation/termination at the 5′ and 3′ ends, respectively.

To help understand the origins of spirotrich nanochromosomes, we investigated the expansion and contraction of gene families in 31 ciliates based on 238 orthogroups. Nanochromosome architecture has been reported in species of several disparate ciliate clades (classes Spirotrichea, Armophorea, and Litostomatea), suggesting multiple origins of extensive fragmentation within ciliates (Riley and Katz, 2001; McGrath et al., 2007; Huang and Katz, 2014; Špaková et al., 2014; Park et al., 2021). The spirotrich clade, which features gene-sized chromosomes, originated approximately 900 million years ago, accompanied by a dramatic expansion/contraction of several gene families (Figure 4). Within spirotrichs, Pseudokeronopsis species originated 220 million years ago, which in turn separated from the clade containing Oxytricha and Stylonychia 430 million years ago. As gene family expansion reflects phenotypic diversity and genetic adaptions during evolution (Harris and Hofmann, 2015; Yan et al., 2019), we identified 1079 and 768 expanded gene families (p < 0.05) in P. carnea and P. flava, respectively. These expanded gene families are enriched in a variety of pathways (Supplementary Figure 2). Compared with three representative species that do not carry nanochromosomes (Supplementary Figure 3), the expanded gene families in both Pseudokeronopsis species contribute to transcription factor binding and sequence-specific DNA binding. Although given the incomplete nature of the current ciliate genomic datasets as a limitation, our analyses provide a baseline about the transcription regulation pathways rewiring in species with nanochromosomal organization for the future studies.