Transcriptome analysis reveals the encystment-related lncRNA expression profile and coexpressed mRNAs in Pseudourostyla cristata

Ciliated protozoans form dormant cysts for survival under adverse conditions. The molecular mechanisms regulating this process are critical for understanding how single-celled eukaryotes adapt to the environment. Despite the accumulated data on morphology and gene coding sequences, the molecular mechanism by which lncRNAs regulate ciliate encystment remains unknown. Here, we first detected and analyzed the lncRNA expression profile and coexpressed mRNAs in dormant cysts versus vegetative cells in the hypotrich ciliate Pseudourostyla cristata by high-throughput sequencing and qRT-PCR. A total of 853 differentially expressed lncRNAs were identified. Compared to vegetative cells, 439 and 414 lncRNAs were upregulated and downregulated, respectively, while 47 lncRNAs were specifically expressed in dormant cysts. A lncRNA-mRNA coexpression network was constructed, and the possible roles of lncRNAs were screened. Three of the identified lncRNAs, DN12058, DN20924 and DN30855, were found to play roles in fostering encystment via their coexpressed mRNAs. These lncRNAs can regulate a variety of physiological activities that are essential for encystment, including autophagy, protein degradation, the intracellular calcium concentration, microtubule-associated dynein and microtubule interactions, and cell proliferation inhibition. These findings provide the first insight into the potentially functional lncRNAs and their coexpressed mRNAs involved in the dormancy of ciliated protozoa and contribute new evidence for understanding the molecular mechanisms regulating encystment.

www.nature.com/scientificreports/ these DE proteins, 12 were specifically expressed proteins, and 14 were DE proteins 5 . Gao et al. found 6 specific proteins in dormant cysts of Pseudourostyla cristata (P. cristata) by using a shotgun LC-MS/MS approach 6 . More recently, two studies identified mRNA expression profiles related to the encystment of Colpoda aspera and P. cristata 4,8 . These two studies not only analyzed DE genes related to cyst formation but also identified the signaling pathways that may facilitate encystment, such as the cAMP, mTOR, PI3K/AKT, calcium and AMPK signaling pathways. The data from these molecular studies offer an in-depth view of the regulatory mechanism of cellular patterns during encystment and the function of organelles. LncRNAs are no less than 200 nucleotides in length and lack substantial protein-coding capability. In recent years, important regulatory functions have been established and extended to study various cellular processes 4,8 . Increasing evidence has demonstrated that lncRNAs play notable roles in a broad variety of biological processes, such as the modulation of epigenetic, transcriptional or posttranscriptional gene expression and the stress response [13][14][15] . However, to the best of our knowledge, the role of lncRNAs in ciliate encystment has not been studied.
In the present study, P. cristata, which can form massive cysts in culture, was used as the experimental model to study the potential function of lncRNAs in ciliate encystment. High-throughput sequencing (HiSeq), quantitative real-time PCR (qRT-PCR) and bioinformatic techniques were used to determine the expression profiles of lncRNAs in dormant cysts versus vegetative cells, and DE lncRNAs in the encystment process were analyzed. In addition, the roles of DE lncRNAs and their coexpressed genes in the encystment mechanism were addressed. The current research work shows, for the first time, lncRNAs that regulate ciliate encystment, contributing to a better understanding of the complicated molecular mechanisms regulating encystment stress resistance in ciliates.

Results of sequencing and characteristics of transcripts. The complete experimental design is
shown in Supplementary Fig. S1. A total of 98,557,186 and 99,463,640 raw reads were generated from the databases of vegetative cells and dormant cysts. After removing the low-quality reads, adaptor sequences, and heterogeneous sequences and performing quality control, 96,986,214 and 97,885,228 clean reads were obtained for the vegetative cell and dormant cyst libraries, respectively ( Table 1). The Q30 percentages of the reads were 95.46% in vegetative cells and 95.57% in cysts in our experiments; thus, the preprocessed read and base detection findings were reliable. Next, assessment of the read contamination indicated that the samples were not contaminated (Fig. 1).
To differentiate lncRNAs from protein-coding transcript units, three filtering processes mentioned in the system were used. We identified 11,959 accurately expressed lncRNA transcripts (Fig. 2). In addition, the results showed 11,959 lncRNAs, with an N50 of 919 bp and an average length of 768.49 bp. The maximum and minimum lengths of unigenes were 15,355 and 301, respectively (Table 2). Transcripts with a length of 301-400 bp accounted for the highest number of lncRNAs (3205), and transcripts containing 1901-2000 bp accounted for the lowest number (84) (Fig. 3a). We observed that the data were roughly symmetrical in the box plot of lncRNA expression levels (Fig. 3b). Due to the differences in the number of genes expressed and the distribution of gene expression values in the samples, we divided the expression values in the samples into different intervals and  S2). Next, we constructed an MA map and volcano plot to further understand the overall distribution of the DE lncRNAs ( Supplementary Fig. S3).
LncRNA-mRNA coexpression network. Most noncoding RNAs (ncRNAs) were not annotated and had an unknown function. We conducted correlation analysis between the ncRNAs and protein-coding mRNAs based on the expression of DE transcripts to predict the functions of the ncRNAs. To detect key lncRNAs and their potential functions in the encystment of P. cristata, we constructed a lncRNA-mRNA coexpression network and investigated the potential interactions between the lncRNA transcripts and mRNA transcripts. We constructed a coexpression network with the top 500 transcripts according to the p-value from low to high. There were more than 140 nodes in the network. Among these nodes, seven nodes with the highest correlation (COR) values were selected, and a coexpression network was constructed using R software (Fig. 4). The seven lncRNAs were DN12058_c0_g1_i1_2, DN25275_c0_g1_i1_2, DN28342_c0_g1_i1_2, DN16982_c0_g1_i1_2, DN15570_c0_g1_i1_2, DN17424_c0_g1_i2_2 and DN26140_c0_g1_i1_2. The mRNA transcripts coexpressed with these seven lncRNAs were involved in functions related to the ribosome and oxidative phosphorylation.
GO and KEGG analysis of mRNAs coexpressed with lncRNAs. In order to explore the potential functional implications of the DE lncRNAs, we performed Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) functional enrichment analyses for mRNAs in the lncRNA-associated network. We focused on three lncRNAs that were likely to be associated with cysts of P. cristata through bioinformatic analysis and comparative analysis of data from our previously published articles 8 . Among these three lncR-NAs, the upregulated lncRNA DN12058_c0_g1_i1_2 (DN12058) had 184 coexpressed mRNAs, and the top 30 coexpressed mRNAs were identified as upregulated after screening of the COR values. The top 30 coexpressed   Fig. S4. The term most enriched with the coexpressed mRNAs in the GO biological process (BP) category was autophagy. The most enriched term in the cellular component (CC) category was cytosolic small ribosomal subunit. The most enriched terms in the molecular function (MF) category was binding to mRNA. The mRNAs coexpressed with the lncRNA DN12058 were enriched in 81 signaling pathways. After screening the transcripts with an absolute fold change greater than 2, the coexpressed mRNAs were found to be significantly enriched in the proteasome, oxidative phosphorylation and ribosome pathways ( Table 4).
In addition, we focused on the upregulated lncRNA DN20924_c0_g1_i2_2 (DN20924) and downregulated lncRNA DN30855_c0_g1_i1_1 (DN30855). the lncRNA DN20924 had 376 coexpressed mRNAs, and the top 30   Fig. 5a. The term most enriched with coexpressed mRNAs in the GO BP category was de novo protein folding. The term protein binding involved in protein folding was the most enriched term in the CC category. The most enriched term in the MF category was chaperonin-containing T-complex. The mRNAs coexpressed with the lncRNA DN20924 were enriched in 139 signaling pathways. After screening the transcripts with an absolute fold change greater than 2, the coexpressed mRNAs were found to be significantly enriched in the following pathways: proteasome, endocytosis and protein processing in endoplasmic reticulum (Fig. 6a). The downregulated lncRNA DN30855 had 902 coexpressed mRNAs, and the top 30 coexpressed mRNAs were identified as downregulated after screening of the COR values. The top 30 coexpressed mRNAs included GGH, TTLL11, MELK, amyA, ICK, ATP2B, etc. (Supplementary Table S3). A total of 1287 GO terms were enriched with mRNAs coexpressed with the lncRNA DN30855. The top 30 GO terms are shown in Fig. 5b. The terms most enriched with coexpressed mRNAs in the GO BP category were pectin catabolic process, regulation of release of sequestered calcium ion into cytosol by sarcoplasmic reticulum, L-ascorbic acid metabolic process and positive regulation of ryanodine-sensitive calcium-release channel activity. The most enriched term in the CC category was extracellular region. The most enriched terms in the MF category were cellulose binding, glutathione dehydrogenase (ascorbate) activity and methylarsonate reductase activity. The mRNAs coexpressed with the lncRNA DN30855 were enriched in 150 signaling pathways. After screening the transcripts with an absolute fold change greater than 2, the coexpressed mRNAs were found to be significantly enriched in the lysosome, sphingolipid metabolism and peroxisome pathways (Fig. 6b). It was also found that the GO and KEGG annotations of the mRNAs coexpressed with the lncRNAs DN16982_c0_g1_i1_2 and DN15570_c0_g1_i1_2 were similar to those coexpressed with the lncRNA DN12058.
Validation of lncRNA expression using qRT-PCR. To validate the HiSeq results, 16 DE lncRNAs were randomly selected for qRT-PCR analysis. As shown in Fig. 7, the expression patterns of the 16 DE lncRNAs evaluated by RNA-Seq and qRT-PCR were significantly consistent, with similar trends. Indeed, the fold changes in the expression levels of the DE lncRNAs as evaluated by qRT-PCR were not precisely the same as those shown by RNA-Seq. For example, the lncRNA DN31340 showed an increased abundance of 1304-fold in dormant cysts by RNA-Seq, whereas the observed increase in abundance in dormant cysts was determined to be 4292-fold by qRT-PCR. This difference in transcript abundance can be attributed to the differential sensitivity of RNA-Seq and qRT-PCR regarding normalization of gene expression data. These findings revealed that the detection of lncRNAs in P. cristata was extremely credible. In addition, correlation analysis was conducted using the cor- www.nature.com/scientificreports/

Discussion
Research on the mechanism underlying cyst formation in ciliates has gradually advanced from the level of morphology to the level of molecular biology, but the research is still mainly focused on coding genes 4,7,8,16 . Studies on the molecular biological functions of lncRNAs in cyst formation in ciliates have not yet been reported.
LncRNAs are considered to be widely distributed in the genome, to be significant regulators of gene expression, and to participate in a wide range of important cellular processes. The regulatory roles of lncRNAs have been a focal point of research in the life sciences [13][14][15] . In this study, for the first time, we identified DE lncRNAs between dormant cysts and vegetative cells. As discussed below, these up-or downregulated lncRNAs can engage in various encystment-related activities through their coexpressed mRNAs. Although the results are preliminary, the large number of DE lncRNAs can form the basis of further studies on the molecular mechanisms of lncRNAs in ciliate encystment.

Regulation of lncRNAs and their coexpressed mRNAs. Among the 853 identified DE lncRNAs and
their coexpressed mRNAs, we mainly focused on analyzing three DE lncRNAs and their coexpressed mRNAs. GO term analysis showed that the mRNAs coexpressed with the upregulated lncRNA DN12058 were signifi- www.nature.com/scientificreports/ cantly enriched in the term autophagy, and KEGG pathway analysis showed that the mRNAs coexpressed with the upregulated lncRNA DN12058 were significantly enriched in the proteasome signaling pathway. These data suggested that two pathways, i.e., the ubiquitin-proteasome pathway and the autophagy pathway, were active in the degradation of proteins and organelles during cyst formation. These findings are consistent with the finding of Pan et al. 8 that the proteasome system was activated to degrade excess proteins during P. cristata encystment. Furthermore, these molecular events corresponded well with dramatic changes in cell morphology, including rapid cell volume contraction, autophagy and autophagosome formation during the cyst formation process 8,17 .
Our results showed that the mRNAs coexpressed with the upregulated lncRNA DN12058 included ATP1, ATP6, ND4, ND5, EEF1A, ND5, htpG (HSP90α), COX3, UBC, petF, acpP and PABP. www.nature.com/scientificreports/ ATP1 and ATP6 are ATP synthase subunits and key enzymes in mitochondrial energy metabolism in organisms, and they are involved in oxidative phosphorylation 18,19 . COX3 is a cytochrome C oxidase (COX) subunit and plays an important role in oxidative phosphorylation through the mitochondrial respiratory chain 20 . ND4 and ND5 are NADH dehydrogenase subunits and two of the seven subunits of respiratory complex I (NADH dehydrogenase) that span the inner membrane of mitochondria. They participate in the first step of mitochondrial oxidative phosphorylation through the electron transport chain 21 . Ferredoxin (petF) is an electron carrier www.nature.com/scientificreports/ active in the electron transport chain in mitochondria 22 . It was speculated that the lncRNA DN12058 is likely to facilitate the oxidative phosphorylation process by orchestrating the upregulation of these coexpressed mRNAs and then providing the ATP necessary for the rapid morphological changes during cyst formation. HtpG, called heat shock protein-90α (Hsp90α), is a member of the HSP90 family 23 . Hsp90α is secreted by a variety of cell types in response to extracellular stress signals and plays an important role in cell repair 24 . Hsp90α also participates in a variety of key cellular reactions, including signal transduction, cell cycle modulation and transcriptional regulation 25 . These observations suggest that lncRNA DN12058 probably participates in the abovementioned functions to facilitate the development of cysts by upregulating the coexpressed HSP90α. EEF1A is involved in cellular processes such as proteasomal protein degradation 26 and cytoskeletal rearrangement 27 , which suggests that upregulated expression of EEF1A in tandem with lncRNA DN12058 increases protein degradation and cytoskeletal rearrangement during encystment to facilitate rapid changes in cell morphology and the formation of compact rounded cysts.
Several studies have shown that polyadenylate-binding protein (PAB) binds directly to the poly(A) tail of mRNA and plays a critical role in the initiation of eukaryotic translation and stabilization/degradation of mRNA 28,29 . Acyl carrier protein (acpP) is a fundamental component of fatty acid biosynthesis and is a cofactor in many other biochemical pathways, such as those mediating the synthesis of polyketides, phospholipids, and cell wall polysaccharides 30 . The polyubiquitin gene ubiquitin C (UBC) is considered to be a stress-protective gene and is upregulated under various stress conditions, which is probably a consequence of an increased demand for ubiquitin to remove toxic misfolded proteins 31 . Therefore, we speculated that the lncRNA DN12058 can increase the levels of encystment-related substances (proteins, cyst wall polysaccharides) as well as the ability of cysts to defend against stress.
Our results showed that the mRNAs coexpressed with the upregulated lncRNA DN20924 included OSBPL5/8, ZIP9, TOP3, VAMP7, EIF1 and UBE1. UBE1 is a ubiquitin-activating enzyme, and its upregulated expression was consistent with our previous findings 8 . Under adverse conditions, the expression of UBE1 is often upregulated, and moreover, upregulated UBE1 enhances the activity of the ubiquitin-proteasome system, which in turn degrades various abnormal and damaged proteins during cyst formation to increase the ability of cysts to adapt to adverse conditions. VAMP7 is a v-SNARE protein in the longin family and is involved in various cell processes, including autophagy, cell membrane repair and secretion of lysosomes. It is also an important player in primary vesicle fusion 32 . Due to its different subcellular localization, VAMP7 is involved in several membrane transport steps, including late endosome-to-lysosome fusion, Golgi-to-plasma membrane exocytosis, and lysosome-toautophagosome vesicle transport [33][34][35][36] . Gu et al. observed that many membrane tubes were aggregated into tubular complexes in P. cristata cysts and speculated that their secretions are related to the formation of cyst wall precursors 8,16,37 . This observation implied that coexpression of VAMP7 with the upregulated lncRNA DN20924 is important to promote vesicle transport of cyst wall precursor-related substances during cyst formation and to break down excess components through phagocytosis and lysosomal degradation to facilitate the formation of compact rounded cysts.
Overexpression of OSBPL5/8 resulted in an increase in Ca 2+ in the mitochondrial matrix during histamine stimulation and an increase in the concentration of Ca 2+ in the endoplasmic reticulum-plasma vacuolar subregion 38 . Our previous publication also showed that an increase in the intracellular Ca 2+ concentration can stimulate or promote P. cristata cyst formation 8 . Therefore, coexpression of OSBPL5/8 with the upregulated lncRNA DN20924 probably promoted an increase in the intracellular Ca 2+ concentration as well as Ca 2+ signaling pathway activation and consequently stimulated cyst formation.
A large amount of evidence indicates that zinc is an essential metal that has many vital functions in eukaryotes and often exerts cell-specific effects on morphogenesis, differentiation and development. Zinc acts as a catalytic cofactor for more than 300 enzymes and also functions as an intracellular signaling molecule. Zinc deficiency can lead to various adverse effects on physiological processes. Zinc transporter member 9 (ZIP9) is a membrane-bound protein and an important member of the zinc transporter family that has essential roles in zinc homeostasis by regulating the transport of zinc across membranes into the cytoplasm. ZIP9 participates in the regulation of intracellular zinc homeostasis and is involved in the processes of morphogenesis and differentiation in eukaryotes 39,40 . The coexpression of ZIP9 with the upregulated lncRNA DN20924 may facilitate the transformation of vegetative cells to into cysts and maintain intracellular zinc homeostasis during encystment.
Recent studies have shown that DNA topoisomerases, including DNA topoisomerase III (TOP3), are ubiquitous enzymes that catalyze topological structural changes in DNA necessary for many steps in DNA processing, such as replication, transcription, repair and recombination. TOP3 plays roles in maintaining genomic stability and participates in various cellular processes. TOP3 can promote cyst formation by inducing gene expression in nonciliate protozoa 41,42 . Several studies have documented that nuclear changes, macronuclear chromatin reorganization and DNA modifications such as demethylation occur during ciliate encystment 43,44 . These results indicated that coexpression of TOP3 with the upregulated lncRNA DN20924 stimulated cyst formation by promoting DNA recombination and transcription and inducing cyst wall protein gene expression in P. cristata. Upregulated expression of TOP3 may also promote DNA repair and maintain genome stability to enhance the resistance of cysts.
It is worth noting that EEF1A, coexpressed with the upregulated lncRNA DN12058, and EIF1, coexpressed with the upregulated lncRNA DN20924, are a translation elongation factor and a translation initiation factor, respectively. Both of these factors are important players in protein synthesis. The synergistic upregulated expression of these two essential factors suggested that they stimulated the synthesis of encystment-related proteins such as cyst wall proteins and promoted encystment.
The lncRNA DN30855 and all of its coexpressed mRNAs, primarily TTLL11, MELK and ICK, were downregulated. Tubulin polyglutamylase (TTLL) can produce polyglutamine for an original posttranslational modification of tubulin. TTLL is widely found in protozoa 45  www.nature.com/scientificreports/ regulating microtubule structure or the interaction between microtubule-associated dynein and microtubules 46,47 .
Considering that microtubules in the ciliature are absorbed during the encystment of P. cristata 37 , it was proposed that coexpression of TTLL11 with the downregulated lncRNA DN30855 could result in poor microtubule transport functions and facilitate the entry of P. cristata into dormancy. Maternal embryonic leucine zipper kinase (MELK) is a cell cycle-dependent protein kinase, and its expression is specific to proliferating cells 48  It is worth noting that KEGG analysis of the mRNAs coexpressed with the lncRNA DN30855 showed enrichment in signaling pathways such as lysosome and peroxisome. Downregulated expression of mRNAs in these signaling pathways results in decreased signaling pathway activity, leading to slowing of metabolism during encystment.
Based on the above results, we generated a schematic diagram of a proposed hypothetical signaling network of the DE lncRNAs and coexpressed mRNAs that regulate or promote P. cristata encystment (Fig. 8).

Conclusions
In the current study, for the first time, we determined the lncRNA expression profiles of dormant cysts versus vegetative cells of P. cristata by HiSeq and revealed lncRNAs that were probably involved in the regulation of encystment. Three specific lncRNAs, DN12058, DN20924, DN30855, and their coexpressed mRNAs were studied. These lncRNAs may play roles in enhancing autophagy, protein degradation, cellular stress tolerance, synthesis of proteins required for encystment, intracellular calcium concentrations, and inhibition of cell proliferation via their coexpressed mRNAs during encystment. These three lncRNAs are likely to play synergistic regulatory roles through their coexpressed mRNAs to promote encystment. These findings provide new insights into the possible roles and detailed regulatory mechanism of lncRNAs in ciliate encystment and lead to a deeper understanding of the molecular mechanism underlying the resistance of ciliates to adverse conditions and dormancy in eukaryotes.

Methods
Cell culture and encystment induction. P. cristata (Ciliophora, Urostylida) was provided by Professor Fukang Gu from East China Normal University. P. cristata specimens were collected from paddy field water in southern Anhui, China, in May 2005. The ciliate cells were cultured in 10-cm Petri dishes with treated pond water and incubated at 25 °C. The pond water was filtered with filter paper produced from cotton linters and was then treated by autoclaving (LDZF-75L, Shanghai Shenan Co, Ltd., China) for 30 min at 121 °C. One or 2 www.nature.com/scientificreports/ boiled wheat grains were added to the culture to enrich bacteria, as the diet of ciliate cells. When the cultured vegetative cells reached a high density, the wheat grains were removed, and most of the vegetative cells formed cysts within 3-4 days. Both assays in our RNA-Seq study were performed with two biological replicates and technical replicates.
RNA isolation, library construction, and sequencing. Approximately  Preprocessing of raw reads. To avoid data error, quality control of the raw reads was completed with Trimmomatic 50 software by the following steps: (1) removal of adaptor sequences; (2) removal of low-quality reads; and (3) elimination of 3′ and 5′ low-quality bases via different methods. Then, the original amount of sequencing reads, effective quantity of sequencing reads, Q30 values and guanine and cytosine contents were determined and used to comprehensively evaluate the quality. Next, contamination of reads from the vegetative cells and the dormant cysts was assessed using BLAST and the NT library (fp://fp.ncbi.nih.gov/blast/db) with the criteria of an E value < 1e-10 and coverage > 80%, which confirmed that the samples were not contaminated.

Prediction of lncRNAs.
De novo transcriptome splicing was used to connect overlapping reads into longer sequences without relying on the reference genome, followed by splicing into transcripts through continuous extension. The assembled transcripts were selected for further analysis. Then, the retained transcripts were filtered according to the following process. (1) Transdecoder software was used to identify the coding region in the transcript. (2) The candidate sequence annotation information was combined to filter out the annotation sequences in the coding library. (3) 53 , and predictor of lncRNAs and messenger RNAs based on an improved k-mer scheme (PLEK, version 1.2, available online: https:// sourc eforge. net/ proje cts/ plek/ files/) 54 were used to filter transcripts with coding potential. Finally, the intersection of the CPC, CNCI, Pfam, and PLEK results was selected. The expression levels of the samples were calculated with the FPKM method 55 . FPKM, considers the impact of sequencing depth and unigene length on the fragment count and is currently the most commonly used expression level estimation method.
Coexpression correlations of the de lncRNAs and mRNAs. To reveal the correlated functions of these DE lncRNAs, we performed coexpression analysis of the DE lncRNAs and their coexpressed mRNAs.
The DE lncRNAs and their coexpressed mRNAs were used to construct lncRNA-mRNA pairs. The coexpressed mRNA in each lncRNA-mRNA pair can be selected to explore the main functional role of the lncRNA 56 . Therefore, by analyzing the functions of the mRNAs coexpressed with the DE lncRNAs, we can predict and judge their potential roles in cyst formation. Based on this approach, we analyzed and predicted the potential functions of the DE lncRNAs in the encystment of P. cristata.
To explore the functional roles of the lncRNAs, the coexpression correlations of the DE lncRNAs and DE mRNAs in vegetative cells and dormant cysts were analyzed according to the FPKM method. Then, the lncRNAcoexpressed mRNA pairs were analyzed using Pearson correlation analysis. A Pearson correlation coefficient ≥ 0.8 and p ≤ 0.05 were considered to indicate coregulation, and coregulated lncRNA-mRNA pairs were then screened for GO and KEGG pathway analyses.

GO and KEGG enrichment analyses.
Based on the coexpression analysis results, mRNAs coexpressed with lncRNAs were subjected to GO and KEGG enrichment analyses. LncRNA gene functions were predicted based on the GO and KEGG functional annotations of the coexpressed mRNAs. GO provides controlled and structured vocabularies that model BPs, CCs and MFs 57 . KEGG is a database containing 16 main databases roughly divided into systems information, chemical information and genomic information 58 . The most enriched GO terms might reflect potential functions of the lncRNAs. KEGG pathway analysis was also performed to understand the functions of the DE lncRNAs and their interactions with their coexpressed genes. These analyses revealed the kinds of roles that these DE lncRNAs play in the encystment of P. cristata.

qRT-PCR validation of lncRNAs.
To verify the reliability of the HiSeq results, we randomly selected and validated the expression of 16 lncRNAs in vegetative cells and dormant cysts by qRT-PCR analysis. Total RNA was extracted using a MiniBEST Universal RNA Extraction Kit (Takara, RR9767) according to the manufacturer's instructions. Total RNA was transcribed into cDNA by using a PrimeScript RT Reagent Kit (Takara, RR036A) following the manufacturer's guidelines. The primers for qRT-PCR verification were designed by Sangon Biotech www.nature.com/scientificreports/ (Shanghai, China). qRT-PCR analysis was performed using a TB Green Premix Ex Taq II Kit (Takara, RR820A). Briefly, the 20-μl qRT-PCR mixtures consisted of 10 μl of TB Green Premix Ex Taq II (Tli RNase H Plus), 0.8 μl of PCR forward primer (10 μM), 0.8 μl of PCR reverse primer (10 µM), 0.4 μl of ROX Reference Dye II, 2 μl of cDNA and 6 μl of nuclease-free water. The amplification reactions were incubated at 95 °C for 30 s, followed by 40 cycles at 95 °C for 5 s and 60 °C for 30 s. Glyceraldehyde-3-phosphate dehydrogenase (GAPDH) was used as the reference gene, and all of the primers used are listed in Supplementary Table S4. Quantification of lncRNA expression was performed using the comparative Ct method, and the specificity of the amplification products was evaluated by melting curve analysis. For each sample, reactions were performed in three independent wells. The relative gene expression levels were calculated by using the 2 −ΔΔCt method 59 . Student's t-test was performed, and the results were considered to be significant at p < 0.05.

Data availability
Statistical analysis of the lncRNA data is described in the corresponding sections. Group comparisons were performed with independent t-tests. p < 0.05 was considered to be statistically significant. The raw Illumina sequences have been deposited in the Sequence Read Archive under accession numbers SRR11413602 and SRR11413601. www.nature.com/scientificreports/ Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.