Decryption of the survival "black box": gene family expansion promotes the encystment in ciliated protists

Abstract

Background: Encystment is an important survival strategy extensively employed by microbial organisms to survive unfavorable conditions. Single-celled ciliated protists (ciliates) are popular model eukaryotes for studying encystment, whereby these cells degenerate their ciliary structures and develop cyst walls, then reverse the process under more favorable conditions. However, to date, the evolutionary basis and mechanism for encystment in ciliates is largely unknown. With the rapid development of high-throughput sequencing technologies, genome sequencing and comparative genomics of ciliates have become effective methods to provide insights into above questions.

Results: Here, we profiled the MAC genome of Pseudourostyla cristata, a model hypotrich ciliate for encystment studies. Like other hypotrich MAC genomes, the P. cristata MAC genome is extremely fragmented with a single gene on most chromosomes, and encodes introns that are generally small and lack a conserved branch point for pre-mRNA splicing. Gene family expansion analyses indicate that multiple gene families involved in the encystment are expanded during the evolution of P. cristata. Furthermore, genomic comparisons with other five representative hypotrichs indicate that gene families of phosphorelay sensor kinase, which play a role in the two-component signal transduction system that is related to encystment, show significant expansion among all six hypotrichs. Additionally, cyst wall-related chitin synthase genes have experienced structural changes that increase them from single-exon to multi-exon genes during evolution. These genomic features potentially promote the encystment in hypotrichs and enhance their ability to survive in adverse environments during evolution.

Conclusions: We systematically investigated the genomic structure of hypotrichs and key evolutionary phenomenon, gene family expansion, for encystment promotion in ciliates. In summary, our results provided insights into the evolutionary mechanism of encystment in ciliates.

Keywords: Chitin synthase; Encystment; Gene family expansion; Hypotrich; Pseudourostyla cristata.

Fig. 1

Sequencing, assembly, and features of the macronuclear genome of Pseudourostyla cristata. A Characteristics of all contigs and morphology of P. cristata. Ventral view in vivo is shown in a. Tracks b to f represent the distribution of gene density in sense strand (+), the distribution of gene density in antisense strand (-), GC content, genome coverage of reads, gene number of complete chromosomes, with densities calculated in 100-kb, 100-kb, 1-kb, 100-kb, 2-kb windows, respectively. All contigs are arranged from small to large in size (the outermost circle). B Statistics on assembly and annotation information of the macronuclear genome of P. cristata. Alveolate database (alveolata_odb10) was used for BUSCO analysis. C Distribution pattern of sequencing depth of all contigs. D The cumulative distribution of contig length. E The length distribution of contigs with different telomeres (0, 1 and 2). F The usage of three typical stop codons in P. cristata (Q, Glutamine; *, stop). G The length distribution of coding sequence (CDS). H The length distribution of intron sequence. Motifs of 31 and 34 bp introns are listed on the top right

Fig. 2

Analyses of heterozygosity of MAC genome (A), intron structure (B) in Pseudourostyla cristata and nine representative ciliates. A Log plot of a Kmer spectral genome composition analysis of P. cristata. Heterozygosity was estimated by jellyfish and GenomeScope2, based on 21-mers in Illumina sequence reads of the MAC genome. The len = inferred haploid genome length, uniq = percentage non-repetitive sequence, aa = Homozygous, ab = Heterozygous, kcov = mean kmer coverage for heterozygous bases, err = error rate of the reads, dup = average rate of read duplications, k = Kmer, p = ploid. The observed 21-mers frequency distribution is depicted in purple. The black lines represent the modeled distribution of 21-mers in the full genome. The orange lines represent the modeled distribution of the unique fraction of the genome. We find ~ 104× and ~ 208× coverage for heterozygous and homozygous peaks in our dataset, respectively. B Motif sequences of intron with most abundant size category in ten representative ciliates. The purple shadow and polyline show the percentage of base A at each position in the second half of intron sequences. Asterisks represent the conserved A nucleotide (The percentage of base A at the corresponding position is greater than 50%) which most likely represents a branch site

Fig. 3

Comparison of chromosomes and genes length of MAC genome in six hypotrichous ciliates. A Comparison of the distribution of chromosome length. B Comparison of the distribution of gene length. The gray solid dots represent the average length of chromosomes in each species

Fig. 4

Characteristics of MAC genomes among six hypotrichous ciliates. A–D Comparison of the distribution of GC content, gene number per complete chromosome, subtelomeric regions length, and exon number per gene. The numbers in (C) represent the p-value. The number in the legend of (D) represents the exon number per gene

Fig. 5

Analyses of orthogroups and gene family expansion/contraction within ciliates. A Upset plot of intersecting sets of orthogroups in six hypotrichous ciliates. B Phylogenomic tree with divergence time, and gene family expansion/contraction for Pseudourostyla cristata and 16 other species. The numbers at nodes indicate the number of expanded (blue) and contracted (orange) gene families at different evolutionary time points. Numbers following species names represent expanded and contracted gene families for species alone. MRCA, most recent common ancestor

Fig. 6

GO and KEGG pathway enrichment analyses of significantly expanded gene families. A Level 2 GO terms associated with significantly expanded gene families in Pseudourostyla cristata. B, C GO (B) and KEGG (C) enrichment bubble plots of Pseudourostyla cristata. The top 20 GO terms and KEGG pathways with the smallest Q-values are shown, with the ordinate as the GO term or KEGG pathway and the abscissa as the RichFactor. GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes. RichFactor = the ratio of enriched gene number to all gene number in this pathway term

Fig. 7

Consensus tree based on eighty-six chitin synthases (CHS) proteins from nine ciliates and six CHS proteins from two fungi, with corresponding information of gene expression, motif, conserved domain, and gene structure. A Normalized gene expression level of different CHS genes. Data are represented as log10-transformed ratios of TPM values [log10(TPM value)]. Grey dots at the vertical axis indicate corresponding genes not supported by RNA-seq reads. B Consensus tree inferred from 86 CHS proteins of nine ciliates and six CHS proteins of two fungi. Sc, Stentor coeruleus; Bs, Blepharisma stoltei; Ev, Euplotes vannus; Sl, Stylonychia lemnae; Ot, Oxytricha trifallax; Pc, Pseudourostyla cristata; Tt, Tetrahymena thermophila; Pp, Pseudocohnilembus persalinus; Pt, Paramecium tetraurelia. Elliptical and rectangular profiles represent that the pellicle of ciliates is soft and rigid, respectively. Black circles indicate that node support values greater than 95% (high confidence value). C MEME motif distribution of each protein. D Conserved domain of each protein. E Corresponding gene structure of each CHS protein. The green boxes in (E) represent exons, lines represent introns. The four boxes on the black background below the figure are the sum of TPM for CHS genes in each species and the legends of (C), (D), and (E), respectively