E. amieti specimens were collected from a farmland pond in Qingpu District, Shanghai. The DNA extraction kit used was sourced from TaKaRa, while genome sequencing and annotation were conducted by Beijing Nuohe Zhiyuan Technology Co., Ltd. Transcriptome sequencing was carried out by Shanghai Pacino Biotechnology Co., Ltd.

To cultivate and collect the hypotrich ciliate E. amieti, the organisms were fed with Chlorogoium elongatum. After feeding, the E. amieti were starved for 5 days prior to DNA or RNA extraction. The starved cells were filtered using three layers of gauze to eliminate larger impurities. Subsequently, the cells were concentrated on qualitative filter paper to remove smaller impurities. The cells were then resuspended in pure water and centrifuged at 4,000 r/min for 5 min, preparing them for further use.

DNA extraction followed the protocol of the TaKaRa DNA Extraction Kit. Agarose gel electrophoresis was utilized to assess the purity and integrity of the DNA, while the purity (A₂₆₀/A₂₈₀ ratio) was measured using a Nanodrop spectrophotometer. Precise DNA concentrations were quantified with a Qubit fluorometer. Qualified DNA samples were fragmented into 350 bp lengths using Covaris technology. The DNA library was constructed using the NEB Next^® Ultra DNA Library Prep Kit (NEB, USA) and sequenced with the Illumina NovaSeq PE150 platform. Genome assembly was performed using the “careful” mode in SPAdes, followed by sequence fusion using CAP3. Sequences shorter than 500 bp were aligned to those longer than 500 bp using BLAST to filter out low-quality reads. This process removed sequences with a percentage of consistency ≥ 90% and coverage ≥ 80%, as well as bacterial, archaeal DNA sequences, non-telomeric sequences, mitochondrial genome sequences, and 11 sequences shorter than 100 bp. A total of 9,508 reads were excluded. The clean contigs were used for whole genome annotation.

RNA was extracted using the TaKaRa RNA extraction kit and fragmented into 300 bp fragments via ion interruption. A 450 bp library was then constructed and subjected to paired-end sequencing on the Illumina HiSeq platform. After removing sequences with 3’ end adaptors and reads with an average quality below Q20, the clean data was aligned with the genome using Trinity software. The Genome-guided pattern of Trinity was employed for de novo assembly and clustering based on alignment results. The longest transcript was selected as the unigene and used for gene functional annotation and gene structure prediction. Gene function annotation was conducted by comparing gene sets identified through gene structure annotation using InterProScan5 with established protein databases such as SwissProt, NR, Pfam, KEGG, and InterPro. The annotated species included Caenorhabditis elegans, Leishmania donovani, Reticulomyxa filosa, Paramecium duboscqui, Tetrahymena thermophila, Oxytricha trifallax, Stylonychina lemnae, Euplotes vannus, Euplotes octocarinatus, and Euplotes aediculatus.

Cilia-associated genes were identified by referencing the reported cilia and centrosomes genes database (Arnaiz et al., 2009).¹ Relevant sequences of cilia-associated genes were downloaded from the NCBI database using the BLAST similarity retrieval system. Utilizing sequence similarity alignment, the downloaded cilia and centrosomes genes sequences were employed as queries to identify homologous genes within the E. amieti genome by applying BLASTP software (E-value = 1e-5). The NR function and Pfam domain of these predicted protein sequences were annotated using BLASTP. For the genetic comparison among various hypotrich ciliates, genomic data for four additional hypotrich ciliates were retrieved from genome databases for Euplotes vannus,² Euplotes octocarinatus,³ and Oxytricha trifallax,⁴ and Stylonychia lemnae.⁵

Gene annotation, metabolic pathway annotation, gene function annotation, protein classification, and disease-associated annotation were conducted using various online tools, including NR,⁶ WebGestalt,⁷ METASCAPE,⁸ and GO.⁹ Conserved protein domains were analyzed using MEME¹⁰ online software, while Pfam¹¹ was used for protein domain prediction. Protein interactions were analyzed using the STRING database (PPI enrichment p-value < 1.0e-16). Genetic disease networks were visualized and analyzed using the DisGeNET database. Data visualization tools, including Venn diagrams, chord diagrams, and enrichment dot bubbles, were created using Origin Pro 2024b.

Screening and bioinformatics analysis of central genes The FASTA sequences of ciliopathy proteins obtained via BLAST software were imported into the STRING database¹² to create PPI networks. These network diagrams were further optimized and analyzed using Cytoscape software.¹³ The topological parameters of each protein in the PPI network were calculated using 12 topological analysis methods available in the CytoHubba plug-in of Cytoscape. The top 50 proteins, based on differential expression genes (DEGs) parameters from the 12 algorithms, were selected, and those appearing in the top 50 of at least six algorithms were identified as pivotal hub genes.

The identification of frameshifting genes was conducted following the methodology outlined by Wang (2017). Initially, all cilia-associated genes transcripts were aligned using BLASTX against the NCBI’s non-redundant protein database (E-value = 1e-5). Subsequently, transcripts that aligned to the same protein sequence but exhibited different open reading frames were identified and analyzed.

Based on the reported Euplotes PRF genes, the initial ORF consistently terminates with a 5′-AAA TAR-3′sequence. The Euplotes stop codon, either TAA or TAG was identified at the end of the starting reading frame. Subsequently, the “T” or “TA” was excised from the stop codon, contingent upon the potential type of frameshift. The resulting frameshift genes (fsgenes) were re-evaluated against the non-redundant protein database of NCBI. Once the C-terminal extended protein was produced as a results of the frameshift, the gene was classified as a PRF gene.

Cell division in E. amieti was induced by synchronized post-starvation feeding, with 80% of the cells entering the division phase 8 h after feeding on Chlorogoium elongatum. RNA extraction from E. amieti in both the division and G0 phases was performed using the RNeasy kit, following the manufacturer’s instructions. cDNA synthesis was conducted using the PrimeScrip™ RT Master Mix kit (Code No. RR036A, TaKaRa, Japan). Amplification was carried out using the TaKaRa PCR reaction kit. The reaction conditions included an initial denaturation at 94°C for 5 min, followed by 30–45 cycles of denaturation at 94°C for 45 s, annealing at 46°C for 45 s, and extension at 72°C for 1 min. A final extension at 72°C for 10 min was followed by storage at 4°C.

Quantitative PCR (qPCR) was performed using TB Green™ Premix Ex Taq™ II (Tli RNaseH Plus) (Takara Bio, Inc.) on samples from three independent experiments, following the manufacturer’s guidelines. Thermocycling conditions were pre-denaturation at 95°C for 30 s, followed by 35 cycles of denaturation at 95°C for 10 s, and annealing and extension at 60°C for 30 s. All experiments were conducted in duplicate, using the GAPDH gene as the housekeeping gene for normalization and quantification. The relative fold change in expression was calculated using the 2^−ΔΔC_q method [ΔC_q = C_q (target gene)−C_q (reference gene), ΔΔC_q = ΔC_q (target sample) −ΔC_q (control sample)], with C_q values below 30 considered valid. Data analysis followed the manufacturer’s guidelines. Primer sequences are listed in Supplementary Table S1.

Following the methodology outlined by Paschka et al. (2003), E. coli HT115 strains harboring recombinant expression vectors L4440-ARL2BP and L4440-DYNLRB2 were used to feed E. amieti. The L4440 vector is a widely used plasmid for RNAi studies, particularly in ciliates, due to its ability to express dsRNA under the control of the T7 promoter. The control group was fed E. coli HT115 containing the empty L4440 vector without target gene fragments to rule out non-specific effects of dsRNA expression. This control allowed us to rule out non-specific effects of dsRNA expression on cell viability and motility, ensuring that the observed phenotypic changes were specifically due to the silencing of ARL2BP and DYNLRB2. E. coli HT115 strains containing L4440, L4440-ARL2BP, and L4440-DYNLRB2 were diluted at a 1:100 ratio and inoculated into LB liquid culture medium, incubated at 37°C with shaking at 210 rpm for 180 min. Upon reaching an OD₆₀₀ of 0.4, isopropyl IPTG was added to a final concentration of 0.4 mM to induce dsRNA expression for 3 h. The bacteria were then washed three times with 1 mL of double-distilled water. Finally, 30 μL of E. coli HT115 was fed to every 2 × 10⁶ cells.

Cells were fixed in 1 mL of 4% paraformaldehyde at room temperature for 10 min and washed three times in ice-cold 0.1 M phosphate-buffered saline (PBS) (5 min per wash). Cells were then incubated in Flutax-2 (Oregon Green™ 488 conjugate, Invitrogen, United States) for 10 min at room temperature without permeabilization. After incubation, cells were washed briefly and mounted with 15 μL Slowfade Gold Antifade Mountant with DAPI (Thermo Fisher Scientific, United States). Images were acquired using an Olympus IX81 fluorescence microscope (Olympus, Japan).

According to the method of Arnaiz et al. (2009), 4–8 E. amieti were transferred onto a glass slide and observed under a microscope (IX81, Olympus). Images were captured at 0.3-s intervals under a 10 × objective lens and analyzed using Prism and ImageJ software.

Quantitative data with a normal distribution were presented as means ± standard error of the mean. Two independent sample t-tests were used to compare the two groups. Statistical analysis was performed using GraphPad Prism software, with a significance level set at p < 0.01 considered statistically significant.

The application of next-generation sequencing technologies to sequence and analyze the macronuclear genome of E. amieti has provided an in-depth understanding of its genetic composition. The gene annotation of genome sequencing of E. amiteti has been uploaded to the NCBI database (). The genome of E. amieti, which is 89 Mb in size with a GC content of 33%, aligns with the typical characteristics observed in ciliates, known for their elevated AT content. The identification of 105 rRNA genes and a substantial gene count of 27,650 highlights the intricate and diverse nature of the E. amieti genome. Particularly intriguing is the reassignment of the stop codon UGA to encode cysteine (Cys), a feature it shares with Euplotes octocarinatus and Euplotes vannus, highlighting a unique aspect of genetic code evolution within the Euplotes genus (Table 1).

To validate the sequencing results, a comparative analysis was conducted with other ciliate genomes, including Paramecium duboscqui, Tetrahymena thermophila, Oxytricha trifallax, Stylonychia lemnae, Euplotes vannus, and Euplotes octocarinatus. The findings indicate a higher number of successfully annotated genes in the E. amieti genome compared to other ciliates, suggesting potential unique gene features or the influence of increased sequencing depth. Additionally, these comparisons reveal variability in the evolution of the genetic code among ciliates, with E. amieti exhibiting genomic similarities to other Euplotes species, particularly in genome size, GC content, and the reassignment of UGA-encoded cysteines. These findings suggest a close evolutionary relationship and potentially similar biological functions and ecological roles in their respective environments. Consequently, the annotation and identification of the E. amieti genome align with existing knowledge of ciliates, providing additional genetic insights.

To elucidate the genetic components associated with cilia, 418 cilia-associated genes were identified in E. amieti by referencing human and reported cilia and centrosomes genes in the Cildb database. In hypotrich ciliates, dorsal-ventral differentiation influences the organization and clustering of cilia in the dorsal and ventral cortex, necessitating further investigation into cilia-associated genes to elucidate their unique ciliary assembly traits. This study undertakes a comparative analysis of cilia-associated genes in E. amieti with three other recognized hypotrich ciliates. Among the identified genes, 299 newly annotated cilia-associated genes are unique to E. amitei (Figure 1A; Supplementary Table S1). All analyzed cilia-associated genes in the four hypotrich ciliates include 52 gene families (Figure 1B; Supplementary Table S2). Notably, families such as tubulin, dyneins, kinesins, ARF (ADP-Ribosylation Factor), and STPK (Serine/Threonine Protein Kinase) have more than 50 members. The number of members within these protein families varies significantly across different species. Tubulin, a microtubule protein, is abundant in all four ciliate species, indicating its crucial role in the structure and function of cilia. Oxytricha trifallax expresses a higher number of dyneins and RAB proteins (Ras-related in brain proteins) compared to the other three species, suggesting a potential specialization or increased reliance on dynein-mediated processes and vesicular transport or signaling pathways.

Common cilia-associated genes across all four ciliates include the tubulin, dynein, and intraciliary transporter families. E. amieti shows higher expression of tubulin, kinesins, and dyneins, which are essential for ciliary movement and assembly. These proteins typically contribute to the structural integrity and motility functions of cilia. In contrast, Euplotes octocarinatus has higher expression levels of proteins such as centrin, CEP (centrosomal protein), and RSPH (radial spoke head protein). Centrin is involved in centriole stability, CEP proteins in centrosome and cilia function, and RSPH in ciliary beat regulation. Furthermore, 44 cilia-associated genes were found to be conserved across the four types of hypotrich ciliates, indicating the high conservation of these genes in hypotrich ciliates (Figure 1A). The subcellular localization of gene products in all four ciliates was enriched in five parts of the cilia: the ciliary basal body, transition zone, plasma, membrane, and tip (Figure 1C). These findings suggest that E. amieti can serve as a valuable model organism for investigating ciliopathies.

GO enrichment analysis of the predicted cilia-associated genes was performed using the MATESCAPE online analysis platform (Figure 2A). The results indicate that the majority of these genes encode enzymes involved in metabolic processes, such as those related to glucose metabolism, including phosphoglycerate mutase, malate dehydrogenase, and ATPase, which are crucial for energy metabolism, comprising 89% of the total. Additionally, 50% of these genes encode proteins containing specific domains, such as the protein kinase domain, tetratricopeptide repeat domain, and adenylate kinase domain. Furthermore, 32% of the genes belong to the family of flagellar proteins, including intra-flagellar transport proteins and flagellar/matrix proteins. Moreover, 20% of the genes are part of the dynamin family, encompassing dynein intermediate chain, dynein heavy chain, and dynein light chain. Another 15% of the genes are categorized within the binding protein family, which includes actin-binding proteins and ran-binding proteins, while 11% fall under the microtubule-related protein family (e.g., tubulin ε chain, α-tubulin N-acetyltransferase). A minority of the proteins are associated with transport or autophagy.

The top 20 enrichments from the GO analysis reveal the primary functional categories of the cilia-related proteins. Within the Biological Process (BP) category, the top three functions include plus-end-directed vesicle transport along microtubules, axonemal central apparatus assembly, and determination of axonemal central apparatus assembly. In the Cellular Component (CC) category, the predominant entries are the ciliary membrane, axoneme, and pericentriolar material. The top three functions in the Molecular Function (MF) category are tubulin binding, microtubule motor activity, and dynein light intermediate chain binding (Figure 2B).

In addition, KEGG pathway enrichment analysis of the predicted cilia-associated genes using the MATESCAPE platform revealed significant enrichment in pathways such as adenine ribonucleotide biosynthesis, oocyte meiosis, glycolysis/gluconeogenesis, vasopressin-regulated water reabsorption, Wnt signaling pathway, AMPK signaling pathway, and Huntington’s disease (Figure 2C).

The cilia-associated genes were categorized into 11 classes based on their roles with cilia (Figure 2D). The composition of these classes varies among the four ciliates studied. In E. amieti, the most prominent functional categories include ciliary structure, other cilia-associated proteins, and ubiquitin. Conversely, in Euplotes octocarinatus, the largest categories are ciliary structure, ciliary transition zone, and other cilia-associated proteins. For Stylonychia lemnae, the ciliary structure, ciliary transition zone total, and basal body/centrosome total are the most prominent categories. In Oxytricha trifallax, the dominant categories include ciliary structure, ciliary transition zone total, and other cilia-associated proteins.

The structural features of the cilia-associated genes in E. amieti were analyzed. Approximately 79.09% of these genes lacked introns, 13.4% contained one intron, and genes with multiple introns accounted for less than 10% of the total. All the genes utilized either TAA or TAG as the stop codon, with over three-quarters employing TAA and less than one-quarter utilizing TAG. Furthermore, around 30% of these genes undergo programmed frameshifting (Figure 3A). Localization analysis revealed that the majority of the gene products were found in the basal body, with 40 proteins localized in basal bodies encoded by PRF genes. The predominant + 1 programmed ribosomal frameshift (PRF) was observed, with AAATAA identified as the slippery sequence (Table 2).

Visual analysis of conserved domains of proteins encoded by cilia-associated genes was conducted using the MEME online software. Three conserved sequences were identified (Figure 3B). Among the 418 cilia-associated proteins analyzed, 37 were found to contain these conserved sequences (ARL3, ARL2, RAB11B, RABL2A, DNAL1, RAB6B, RAN, RAB8B, RAB1A, MAPK15, CSNK1E, AKT3, MDH1, STK38L, CAMK2D, PRKACB, ICK, MAP3K2, GSK3B, MOK, SLC47A1, PLK1, NEK1, RAB2A, PRPF4B, KIF9, NEK4, CDK1, KIF3B, AK7, STK36, CFAP65, PIK3R4, MARK4, PI4KA). These proteins are annotated as “cilia-associated” based on their gene origins, suggesting potential roles in ciliary functions. However, subcellular localization data from existing studies (e.g., UniProt) indicate that some members (e.g., RAN, AKT3) may also localize to the nucleus or cytoplasm, implying multifunctional roles beyond cilia. Further experimental validation is required to confirm their exclusive ciliary localization. The distribution of these sequences revealed that Motif-1 was the most prevalent, present in 26 proteins, followed by Motif-2 in 12 proteins, and Motif-3 in 10 proteins. Additionally, 11 proteins were identified to contain two distinct domains. It is hypothesized that proteins containing different domains may exhibit variations in functionality and evolution. For structural prediction, the BBS (Bardet-Biedl Syndrome) family proteins, known for their localization to basal bodies and close relation to cilia assembly, were selected. The predicted structures of BBS family proteins BBS2, BBS4, BBS7, and BBS9 are depicted in Figure 3C. Notably, BBS2, BBS7, and BBS9 share common conserved motifs, including leucine zipper and a basket structure formed by β-sheets.

An interaction network of 418 cilia-associated genes in E. amieti was constructed using the STRING database, illustrating direct or indirect interactions among these genes (Figure 4A). Within this network, 39 hub genes were identified, including WDR60, IFT88, IFT122, IFT140, IFT172, WDR35, IFT81, DYNLL2, DYNC2H1, DYNC2LI1, WDR19, IFT80, IFT22, IFT27, KIF3B, VCP, IFT74, BBS1, BBS2, BBS4, TTC26, IFT57, BBS5, BBS9, OFD1, TTC21B, IFT46, CEP290, CLUAP1, HSPA8, TRAF3IP1, DNAI2, CCDC65, CCT3, RSPH3, SPAG6, ACTB, CLTC, and DYNC1H1. To further investigate the function of these hub genes, their relationships with ciliopathies were explored using DisGeNET, a tool for analyzing genetic disease networks (Piñero et al., 2017). The analysis revealed that these hub genes are involved in various ciliopathies, with the top 10 including Jeune thoracic dystrophy, Kartagener syndrome, ciliary motility disorders, Meckel-Gruber syndrome, Bardet-Biedl syndrome, situs inversus, situs ambiguous, cone-shaped epiphyses, ulnar polydactyly of fingers, and retinitis pigmentosa (Figure 4B). These hub genes are linked to ciliopathic phenotypes involving the central nervous system, brain diseases, skeletal development disorders, retinal diseases, BBS syndrome, Joubert syndrome, Meckel syndrome, and primary ciliary dyskinesia. The relationship between hub genes and ciliopathic phenotypes was visualized using a GO chord diagram. The study found that SPAG6, TRAF3IP1, IFT22, and KIF3B showed lower correlations with ciliopathies, while IFT172, TTC21B, CEP290, ACTB, and DYNC1H1 were highly correlated. This suggests that certain genes are involved in multiple ciliopathies, while others are specific to individual types (Figure 4C).

Our analysis delved into the cellular components, molecular functions, and biological processes involved in these ciliopathies linked to the hub genes (Figure 4D). The GO enrichment revealed that the primary cellular components include cilia assembly, ciliary motility, centrosomes, and cilia basal bodies, as well as microtubule-based transport, microtubules, and cilia transition zones. The molecular functions and biological processes include the axial filament dynamic protein complex assembly, left/right symmetry determination, regulation of dynein complexes, regulation of microtubule processes, and protein localization to cilia. The correlation between ciliopathies and MF, BP, and CC suggests that ciliary assembly, motile cilia, ciliary basal bodies, and centrosomes are related to various types of ciliopathies. Conversely, factors such as abnormal ciliary beating patterns, extracellular fluid dynamics, and certain dynein-related processes (e.g., non-axonemal dynein assembly) are currently thought to play less direct roles in the pathogenesis of classical ciliopathies, though their contributions may vary across specific disease subtypes.

In ciliates, ciliogenesis exhibits a strong correlation with morphogenesis during the division phase. Extensive research has demonstrated that the dorsal and ventral cilia of the mother cells undergo replication prior to division and are systematically allocated to both daughter cells during division. Consequently, investigating the function of ciliary-associated proteins necessitates an initial examination of ciliary proteins that are regulated by the cell cycle. To verify the expression levels of cilia-associated genes in E. amieti, 48 genes were randomly chosen from the 52 cilia gene families for measurement using qPCR, comparing levels to those in the G0 phase. The findings revealed that the majority of the selected cilia-associated genes were downregulated during cell division. In contrast, six genes showed significant upregulation, including CKAP5, MAP1LC4, ARRDC, ADCY3, TULP3, and TTC, with ADCY3 (adenylate cyclase type 3) showing the highest level of upregulation (Figure 5).

To confirm the connection between the 418 cilia-associated genes and ciliary assembly, we selected ARL2BP and DYNLRB2 for functional validation. These genes were chosen based on their high expression levels in cilia-associated structures and their known roles in ciliary assembly and motility. We hypothesized that silencing these genes would disrupt ciliary function, leading to reduced motility and altered cytoskeletal organization in E. amieti. ARL2BP is a protein present in the basal body regions of photoreceptor cells and is needed for normal photoreceptor cilia doublets and outer segment structure. DYNLRB2 is a microtubule motor protein complex that is responsible for the retrograde transport of cellular cargoes, and is necessary for mitotic regulation. The amplified target fragments of ARL2BP and DYNLRB2 were inserted into L4440 to create the recombinant plasmids pL4440-ARL2BP and pL4440-DYNLRB2. The expression of dsRNA of the target genes in E. amieti was verified using qPCR. The findings indicated that the inserted gene fragments exhibited significant expression in pL4440-ARL2BP and pL4440-DYNLRB2 compared to the control group. Following the feeding of bacteria containing the interference expression vector, the relative expression levels of ARL2BP and DYNLRB2 mRNA were notably decreased in the interference group compared to the control group (Figure 6A). These results demonstrate that RNA interference targeting ARL2BP and DYNLRB2 significantly reduces their expression levels, leading to observable phenotypic changes such as increased mortality and reduced motility in E. amieti. This confirms the functional importance of these genes in ciliary assembly and motility.

The impact of RNAi on cell viability, motility, and cytoskeleton was further investigated. Analysis of the growth curve revealed a significant decrease in the number of surviving E. amieti over time following RNAi treatment compared to the control group. By the 10th day, the population of surviving E. amieti was notably diminished, and by day 15, all E. amieti had perished (Figure 6B).

The swimming trajectories and speeds of E. amieti following interference with ARL2BP or DYNLRB2 at 2, 4, and 6 days were also analyzed. In the control group, E. amieti exhibited typical linear swimming patterns. However, at 2, 4, and 6 days post-ARL2BP interference, E. amieti exhibited circular swimming trajectories with decreasing radii over time, eventually leading to a stationary spinning motion by days 4–6 (Figures 6C,D). Additionally, a significant decrease in swimming speed was observed with prolonged interference (Figures 6E,F).

The organization of cilia and the cytoskeleton was observed using Flutax. The results revealed that the depletion of ARL2BP or DYNLRB2 did not result in significant alterations in ciliary structure; however, notable changes were observed in the arrangement of cortical microtubules linked to the basal bodies following RNA interference of ARL2BP, exhibiting a curved and disordered pattern. Both the caudal cirri and the microtubules associated with frontoventral were observed to disappear (Figure 7).