Research article | Open | Published:
Genome-wide expression profiles of Pyropia haitanensis in response to osmotic stress by using deep sequencing technology
BMC Genomicsvolume 16, Article number: 1012 (2015)
Pyropia haitanensis is an economically important marine crop grown in harsh intertidal habitats of southern China; it is also an excellent model system for studying mechanisms of stress tolerance. To understand the molecular mechanisms underlying osmotic tolerance and adaptation to intertidal environments, a comprehensive analysis of genome-wide gene expression profiles in response to dehydration and rehydration in Py. haitanensis was undertaken using digital gene expression profile (DGE) approaches combined with de novo transcriptome sequencing.
RNA-sequencing of the pooled RNA samples from different developmental phases and stress treatments was performed, which generated a total of 47.7 million clean reads. These reads were de novo assembled into 28,536 unigenes (≥200 bp), of which 18,217 unigenes (63.83 %) were annotated in at least one reference database. DGE analysis was performed on four treatments (two biological replicates per treatment), which included moderate dehydration, severe dehydration, rehydration, and normal conditions. The number of raw reads per sample ranged from 12.47 to 15.79 million, with an average of 14.69 million reads per sample. After quality filtering, the number of clean reads per sample ranged from 11.83 to 15.04 million. All distinct sequencing reads were annotated using the transcriptome of Py. haitanensis as reference. A total of 1,681 unigenes showed significant differential expression between moderate dehydration and normal conditions, in which 977 genes were upregulated, and 704 genes were downregulated. Between severe dehydration and normal conditions, 1,993 unigenes showed significantly altered expression, which included both upregulated (1,219) and downregulated genes (774). In addition, 1,086 differentially expressed genes were detected between rehydration and normal conditions, of which 720 genes were upregulated and 366 unigenes were downregulated. Most gene expression patterns in response to dehydration differed from that of rehydration, except for the synthesis of unsaturated fatty acids, several transcription factor families, and molecular chaperones that have been collectively implicated in the processes of dehydration and rehydration in Py. haitanensis.
Taken together, these data provide a global high-resolution analysis of gene expression changes during osmotic stress that could potentially serve as a key resource for understanding the biology of osmotic acclimation in intertidal red seaweed.
The marine red alga Pyropia is one of the most economically important mariculture crops. It has an annual production of at least 120,000 tons (dry weight), which is estimated to be worth at least US$1.3 billion per year [1, 2]. For hundreds of years, this crop has been cultivated in East Asian countries such as China, Korea, and Japan, of which it is currently recognized as one of its largest aquaculture industries . Pyropia haitanensis Chang et Zheng, an endemic species naturally distributed and widely cultivated along the coasts of South China, has comprised 75 %–80 % of the total production of cultivated Pyropia spp. in China .
As sessile organisms that inhabit in the intertidal zones of rocky coasts, the thallus of Py. haitanensis is totally submerged in water during high tide, but is exposed to the air during low tide. Therefore, this crop is constantly exposed to fluctuating and extreme abiotic conditions such as cyclic changes in light levels, abrupt temperature changes, and repeated desiccation/rehydration due to the turning tides. Consequently, Py. haitanensis experiences far more rapid and severe water loss during low tide, and unavoidably suffers from dramatic changes of osmotic potential, which were largely different from the environmental stresses experienced by desiccation-tolerant land plants. Py. haitanensis can tolerate extreme water loss (up to 70 %) and dramatic changes in osmotic potential, consistent with other Pyropia/Porphyra species [5–8]. Therefore, Py. haitanensis is considered an ideal model systerm for investigating the mechanisms of osmotic acclimation in intertidal seaweed.
Several studies have elucidated that anatomical properties, as well as physiological and biochemical changes allow Pyropia and Porphyra species to tolerate osmotic stress, which in turn allow them to thrive in the intertidal zone. For example, blade cells secrete a cell wall that is an agar-like sulfated galactan disaccharide (porphyran) with xylan microfibrils and proteoglycan elements ; this hydrophilic wall could reduce the rate of water and mineral loss from Pyropia and Porphyra species during emersion. Metabolic analyses showed that an increase in O-α-D-galactopyr-anosyl-(1→2)-glycerol (= floridoside) occurs under high salinity stress, suggesting that it acts as a compatible solute in reducing cytoplasmic and membrane damage . Abe et al.  reported that Pyropia dentata, a species that thrives at the highest intertidal level, can fully recover its photosynthetic activity after desiccation at a water potential of -158 MPa. Several protective proteins are secreted and accumulate during osmotic stress, which have been shown to play a role in cellular protection during stress. For example, Contreras-Porcia et al.  observed that in Porphyra columbina, the activity of a diverse range of antioxidant enzymes increased during desiccation, whereas this activity diminished to near basal levels during rehydration. These intracellular proteins provide protection against osmotic stress by stabilizing membranes and organelles, as well as counteracting oxidative stress. Furthermore, high levels of a dehydrin-like protein with a molecular weight of 17 kDa, which is extremely hydrophilic under dehydrative conditions and thermally stable, was detected in Porphyra umbilicalis, but not Pyropia yezoensis, which suggests that this protein plays a key role in the observed superior desiccation tolerance in the former species . In addition, expression of cyclophilin (PhCYP18) was shown to be dysregulated in the blades of Py. haitanensis under high salt stress, strong irradiance stress, and multifactorial stress compared to blades under normal conditions, thus suggesting that cyclophilin actively responded to stress situations and induced high stress tolerance .
In the advent of next-generation sequencing (NGS), whole genome-wide expression profiling have become a rapid and efficient method for identifying genes and specific metabolic processes involved in stress response, especially in organisms whose genomes have not been fully sequenced. Small-scale cDNA sequencing projects, cDNA microarrays, and de novo transcriptome sequencing have been performed in several stress-tolerant intertidal algae, including some Pyropia species [15–23]. Moreover, plastid and mitochondrial genomes of several Pyropia species have recently been sequenced [24–26]. In addition, the draft genome sequence of Py. yezoensis has been determined using NGS , whereas whole genome sequencing projects of other Pyropia species (i.e. Py. haitanensis) have not been completed. Although the Pyropia species may be a source of genetic determinants for osmotic acclimation, defense response-related gene discovery efforts in marine algae have been limited, and the molecular mechanisms involved in osmotic acclimation remain unknown to date.
To understand the molecular mechanisms underlying osmotic acclimation and adaptation to intertidal environments, a comprehensive analysis of genome-wide gene expression profiles in response to dehydration and rehydration in Py. haitanensis was performed. Because the genome sequence of Py. haitanensis is not yet available, deep expression profiling is feasible by large-scale sequencing of transcripts from eight hydrated and dehydrated Py. haitanensis algae samples using Illumina HiSeq 2000 platform. This sequencing dataset and analysis results will serve as a valuable resource for identifying the key genes and pathways involved in the response to osmotic stress in Py. haitanensis. The findings of the present study will lay the foundation for elucidating the molecular mechanisms of osmotic acclimation and provide useful information for the genetic breeding of Py. haitanensis.
Global transcriptome assembly
Due to the limited genome information on Py. haitanensis, to obtain a global transcriptome of Py. haitanensis, a cDNA library was constructed from an equal mixture of RNA isolated from different developmental phases and treatments (Additional file 1: Table S1). The library was sequenced using the Illumina HiSeq 2000 platform. A total of 53.46 million raw reads from the sequencing library were obtained. After quality filtering, a total of 47.76 million (about 89.35 % of the raw reads) clean reads were obtained, corresponding to 4.78 G bases. The GC content of the transcriptome was 64.93 %. An overview of the sequencing is presented in Additional file 1: Table S1.
Using the Trinity de novo assembly program , all clean reads were assembled into 31,838 transcripts, with an N50 length of 956 bp and an average length of 657 bp (Additional file 1: Table S1). The transcripts were subjected to cluster analysis. From these transcripts, a total of 28,536 unigenes were obtained, with an N50 length of 827 bp and an average length of 607 bp (Additional file 1: Table S1). The length statistics of assembled transcripts and unigenes revealed that 9,398 unigenes (32.9 %) were > 500 bp, and 3,854 unigenes (13.5 %) were > 1 kb (Additional file 1: Table S1). These results demonstrated the effectiveness of Illumina sequencing in rapidly capturing a large portion of the transcriptome. As expected, for a randomly fragmented transcriptome, there was a positive relationship between the length of a given unigene and the number of reads assembled into it (Additional file 2: Figure S1).
Function annotation of transcriptome unigenes
To assign accurate annotation information to all unigenes, multiple databases were interrogated, including NCBI non-redundant protein (Nr) database, NCBI non-redundant nucleotide sequence (Nt) database, the manually annotated and curated protein sequence (Swiss-Prot) database, Protein family (Pfam), euKaryotic Ortholog Groups (KOG), Gene Ontology (GO), and Kyoto Encyclopedia of Genes and Genomes (KEGG) database. The overall functional annotation is presented in Table 1. A total of 18,217 unigenes (63.83 %) were annotated in at least one database. Among these, 9,114 (31.93 %) had significant matches in the Nr database, whereas 7,562 (26.49 %) unigenes showed similarity to proteins in the Swiss-Prot database. In addition, 13,972 (48.96 %) and 4,952 (17.35 %) unigenes were successfully annotated in the GO and KEGG databases, respectively.
Sequencing and annotation for DGE libraries
To characterize the digital gene expression profiles (DGE) involved in Py. haitanensis response to dehydration/rehydration, a total of eight RNA samples generated from two biological replicates under the four treatments [normal conditions (CON), moderate water loss (MWL), severe water loss (SWL), and rehydration (REH)] were subjected to DGE sequencing analysis by using a Illumina HiSeq 2000 sequencing platform, respectively (Table 2). In total, 117.48 million raw reads were generated from the eight samples. The number of raw reads per sample ranged from 12.47 to 15.79 million, with an average of 14.69 million reads per sample. About 95 % of the raw reads passed the quality filters, resulting in a total of 111.22 million clean reads. The number of clean reads in each sample ranged from 11.83 to 15.04 million.
Clean reads from the eight samples were mapped to the global transcriptome of Py. haitanensis and obtained from the annotation information of each sample using RSEM software . The number of mapped clean reads in each sample ranged from 11.50 to 14.56 million, accounting for 96.78 %–97.74 % of the total clean reads (Additional file 3: Table S2). In addition, the number of read-mapped unigenes ranged from 19,843 to 21,954 (Additional file 3: Table S2). To evaluate the reproducibility of DGE library sequencing, Pearson correlation analysis was performed on every two replicates. A series of scatterplots comparing the normalized read counts (log transformed) for the replicate libraries in each treatment are shown in Additional file 4: Figure S2. The results showed that the square of the Pearson correlation coefficient (R2) between two replicate libraries in each treatment was > 0.91, indicating the reliability as well as operational stability of the experimental results.
Analysis of differentially expressed genes
Read count data was used to identify the differentially expressed genes (DEGs) among various treatments by using DESeq packages . A total of 1,681, 1,993, and 1086 differentially expressed genes were identified (adjusted P-value of < 0.05) between MWL and CON, SWL and CON, and REH and CON, respectively.
A hierarchical cluster analysis of normalized expression values for the 3,569 DEGs from pairwise comparisons among the four treatments (MWL vs. CON, SWL vs. CON, MWL vs. SWL, REH vs. CON, MWL vs. REH, and SWL vs. REH) was performed (Fig. 1). The analysis showed that the gene expression pattern in the MWL treatment was similar to that of the SWL treatment, whereas it was markedly distinct from the REH treatment. Furthermore, the 3,569 unigenes were divided into 10 subclusters based on their expression modulation, representing 10 different expression patterns (Fig. 2). In subcluster 1, the unigenes were significantly enriched in the pathway of aminoacyl-tRNA biosynthesis, DNA replication, mismatch repair, porphyrin and chlorophyll metabolism, and nucleotide excision repair. In subclusters 2 and 3, no significant enrichment was detected. Subcluster 4 was significantly enriched with the pathway of pyrimidine metabolism, folate biosynthesis, and RNA polymerase. In subcluster 5, the unigenes were significantly enriched in the pathway of carbon fixation in photosynthetic organisms, microbial metabolism in diverse environments, pentose phosphate pathway, fructose and mannose metabolism, glycolysis/gluconeogenesis, and biosynthesis of ansamycins. In subcluster 6, the unigenes were significantly enriched in the pathway of DNA replication, mismatch repair, base excision repair, nucleotide excision repair, cell cycle, and non-homologous end-joining. In subcluster 7, protein processing in endoplasmic reticulum was significantly enriched. In subclusters 8, 9, and 10, no significant enrichment was detected.
DEGs in response to dehydration treatment
A total of 1,681 unigenes were significantly (adjusted P-value of < 0.05) differentially expressed between CON and MWL. Among these, 977 genes were significantly upregulated, whereas 704 genes were significantly downregulated (Fig. 3, Additional file 5: Table S3). And 1,993 DEGs were found between CON and SWL, which included 1,219 upregulated unigenes and 774 downregulated unigenes (Fig. 3, Additional file 6: Table S4). With regard to DEGs between CON and MWL, GO enrichment tests showed that there were 39 and 19 significantly enriched GO terms in up- and downregulated unigenes, respectively (Additional file 7: Table S5). Similarly, between CON and SWL, there were 33 significantly enriched GO terms in upregulated unigenes (Additional file 8: Table S6 A). Interestingly, within the upregulated GO term associated with cell death, three unigenes (comp20313_c1, comp92397_c0, and comp19022_c0) functioned in the induction of apoptosis and regulation of apoptotic process under moderate dehydration stress (Additional file 5: Table S3). The enriched GO terms that were downregulated in MWL were mainly involved in tRNA aminoacylation and reproductive structure development. However, in contrast to that observed in MWL vs. CON, the downregulated unigenes in SWL were significantly enriched in 8 GO terms, which included DNA metabolic process, regulation of mitotic cell cycle, DNA-dependent DNA replication, and DNA repair (Additional file 8: Table S6 B).
Based on KEGG enrichment analysis, the upregulated unigenes in MWL were highly enriched in 7 pathways, including carbon fixation in photosynthetic organisms (PATHWAY: ko00710), pentose phosphate pathway (PATHWAY: ko00030), and biosynthesis of secondary metabolites (PATHWAY: ko01110) (Additional file 9: Table S7 A). Two categories of upregulated unigenes involved in the pathway of carbon fixation in photosynthetic organisms were found, which included C4-pathway and C3-pathway enzymes. Among these upregulated unigenes were 4 unigenes encoding C4-pathway enzymes, which included alanine transaminase (ALT), pyruvate kinase, and aspartate aminotransferase (AST), whereas 18 unigenes encoding C3-pathway enzymes were identified (Table 3). Furthermore, within the biosynthesis of secondary metabolites, 2 upregulated unigenes (comp7250_c0 and comp17143_c0) were found to be involved in trehalose biosynthesis (glucose-1-P → trehalose), which included starch phosphorylase and 1,4-alpha-glucan branching enzyme. In addition, the downregulated unigenes were highly enriched in 8 pathways such as DNA replication (PATHWAY: ko03030) and porphyrin and chlorophyll metabolism (PATHWAY: ko00860) (Additional file 9: Table S7 B). With respect to DEGs between CON and SWL, they were highly enriched in 11 pathways (Additional file 10: Table S8), nearly in accordance with that observed in MWL vs. CON. As noted in the heatmap of gene expression, SWL showed an expression pattern that was similar to that observed in MWL.
DEGs in the initial rapid response to rehydration treatment
Comparative analysis showed that there were 1,086 differentially expressed genes (adjusted P-value of < 0.05) between REH and CON, of which 720 unigenes were upregulated and 366 unigenes were downregulated (Fig. 3, Additional file 11: Table S9). GO enrichment testing showed that 28 GO terms, mainly including carbohydrate metabolism (such as carbon utilization, fructose metabolic process, and glycolysis), were significantly enriched among the upregulated genes, whereas no GO terms were enriched among the downregulated unigenes (Additional file 12: Table S10).
Based on KEGG analysis, the upregulated unigenes were highly enriched in carbon fixation in photosynthesis (PATHWAY: ko00710), fructose and mannose metabolism (PATHWAY: ko00051), pentose phosphate pathway (PATHWAY: ko00030), and protein processing in endoplasmic reticulum (PATHWAY: ko04141). Interestingly, within the carbon fixation pathway, only 18 unigenes putatively encoding C3-pathway enzymes were significantly upregulated in response to rehydration treatment, whereas none of genes encoding C4-pathway enzymes were detected (Table 3). Nevertheless, the downregulated unigenes were highly enriched in DNA replication (PATHWAY: ko03030), mismatch repair (PATHWAY: ko03430), nucleotide excision repair (PATHWAY: ko03420), ABC transporters (PATHWAY: ko02010), cell cycle (PATHWAY: ko04110), and base excision repair (PATHWAY: ko03410). Herein, there were 5 downregulated unigenes that were putatively encoding ABCB1, ABCB4, ABCG2, and mitochondrial ABC transporter ATM within the ABC transporters (Additional file 13: Table S11).
Comparison of DEGs responding to dehydration and rehydration treatments
The Venn diagram provides an illustration of the overlaps among differentially expressed genes in response to moderate dehydration, severe dehydration, and rehydration (Fig. 4). As shown in Fig. 4, a total of 355 differentially expressed genes were collectively involved in the response to moderate dehydration, severe dehydration, and rehydration. Of these DEGs, 7 unigenes encoding variety of fatty acid desaturases and fatty acid elongase that were involved in biosynthesis of unsaturated fatty acids (UFAs) were identified (Table 4). Meanwhile, 4 unigenes that encoded molecular chaperones, including HSPA4, endoplasmin, HSP70-2, and HSP81-1 were detected. In addition, 6 unigenes encoding various transcription factors, including basic leucine zipper (bZIP), Sigma-70 region 2/PDZ domain, MerR family regulatory protein, and zinc-finger protein, were collectively upregulated or downregulated in three pairwise comparisons (MWL vs. CON, SWL vs. CON, and REH vs. CON) (Table 5).
Additionally, 537 and 781 specific unigenes were only differentially regulated between MWL and CON, and SWL and CON, respectively. Moreover, a total of 489 specific unigenes showed significantly altered expression during rehydration. Most of these specific unigenes encoded a range of transporters such as aquaporins (comp15218_c0 and comp18064_c0), Na+/K+ ATPase alpha-subunit 1B (comp16722_c0), ZIP zinc transporters (comp14925_c0 and comp8630_c0), as well as ABC transporters (comp16736_c0, comp16736_c2, comp8118_c0, and comp21051_c0) (Additional file 13: Table S11). In addition, 27 unigenes encoding various ribosomal proteins were also detected, which were involved in ribosome biogenesis/translation, thus strongly suggesting a role for de novo translation during rehydration in Py. haitanensis.
Features of the gene encoding 1,4-alpha-glucan branching enzyme
With respect to the unigene (cDNA sequence) encoding 1,4-alpha-glucan branching enzyme, which was identified as PhSBE, it had 4,140 nucleotide residues in which was contained an 2,268 bp open reading frame, which was predicted to code for 755 amino acids with a molecular weight of 85.2 kDa (Additional file 14: Figure S3 A). Further alignment analysis revealed the gene was interrupted by the 291 bp and 208 bp introns, indicating that three exons and two introns were contained in the gene (Additional file 14: Figure S3 B). Analysis of amino acid sequence showed that it matched with Alpha-amylase_C domain, without signal peptide. Multiple sequence alignments showed that PhSBE was highly conserved with 1,4-alpha-glucan branching enzyme (starch branching enzyme) found in other species. It identified 70 % with Chondrus crispus (XP_005716136.1), 69 % with Gracilaria gracilis (AAB97471.1), 67 % with Cyanidioschyzon merolae strain 10D (XP_005536101.1), 59 % with Galdieria sulphuraria (XP_005703889.1), 56 % with Arabidopsis thaliana (NP_195985.3), 57 % with Zea mays (NP_001105316.1) and 56 % with Oryza sativa (BAA82828.1) (Additional file 15: Figure S4).
Validation of RNA-Seq-based gene expression
To validate the expression profiles obtained by RNA-Seq, qRT-PCR was performed on six genes selected at random with high or low expression levels. Expression comparisons were performed between MWL and CON, SWL and CON, and REH and CON by qRT-PCR. Generally, the expression profiles of the genes assayed show upregulated expression in response to MWL and SWL, and marginal increases in response to REH, nearly confirming the DGE expression data (Fig. 5). Two genes encoding omega-6 fatty acid desaturase (comp19838_c0) and pyruvate kinase (comp7242_c0) however showed downregulated expression in response to REH, which showed marginal upregulated expression in response to REH in DGE expression data. The relative expression of the gene encoding malate dehydrogenase (comp15628_c0) was inverse to the expected expression, with downregulated expression in response to MWL and SWL.
The mechanism of desiccation tolerance in Py. haitanensis
Carbon fixation is the most important biological process in all photosynthetic organisms, which can be divided into three general categories: C3, C4, and Crassulacean acid metabolism (CAM) . The Calvin cycle (C3 pathway) is the most basic and universal form of net carbon fixation in plants, algae, and cyanobacteria. However, the C4 pathway is an adjunct of the C3 pathway that developed novel and efficient CO2 concentration mechanisms to enhance ribulose bisphosphate carboxylase/oxygenase (Rubisco) performance even at limiting ambient CO2 levels . The C4 pathway has long been studied in higher plants with much higher photosynthetic rates than C3 plants [32, 33]. Nevertheless, the C4 pathway in algae has also been the subject of several reports in the past decade. The C4 pathway is important in carbon accumulation and photosynthetic carbon fixation in the marine diatom Thalassiosira weissflogii at conditions of low (atmospheric) CO2 . By genome analysis, all the required genes involved in C4 photosynthesis were presented in unicellular green algae Ostreococcus tauri, which is the smallest free-living eukaryote yet described . With respect to the red algae, Yang et al.  reported that, except for pyruvate phosphate dikinase, almost all genes involved in the C4 pathway have been identified in Py. yezoensis. Meanwhile, C. crispus genome possesses two NADP-malic enzymes (ME), three malate dehydrogenases (MDHs), and one phosphoenolpyruvate carboxylase (PEPC), and thus a malate-based C4-like carbon fixation pathway was considered to be present in the algal species . In Py. haitanensis, the expression level of Rubisco (key gene in the C3 pathway) was significantly lower in the conchocelis than that in the thallus, whereas the expression levels of PEPC and phosphoenolpyruvate carboxykinase (PEPCK), which are the key genes in the C4 pathway, were significantly higher in the conchocelis than that in the thallus . The C4-like pathway was thus considered to play an important role in the fixation of inorganic carbon in the conchocelis stage of Py. haitanensis . Interestingly, by DGE analysis, a large number of unigenes encoding C4 and C3 pathway-related enzymes were significantly upregulated in response to both MWL and SWL treatment, except for 2 unigenes (comp7349_c0 and comp17370_c0). Therefore, C4 pathway- and C3 pathway-related genes may collectively play an important role in the fixation of inorganic carbon when Py. haitanensis gametophytes were subjected to desiccation stress.
In most cases, the ability of the plant to survive desiccation correlates with the accumulation of carbohydrates. Trehalose (α-D-glucopyranosyl-(1,1)-α-D-glucopyranoside) is a non-reducing disaccharide of two glucose units that is predominantly present in desiccation-tolerant lower organisms, including some vascular plants such as Selaginella tamariscina, the moss, Tortula ruralis, and algae [37–39], and functions as a stress protection metabolite in the stabilization of biological structures under various abiotic stresses . In the present study, two upregulated unigenes encoding starch phosphorylase and 1,4-alpha-glucan branching enzyme were detected in dehydration treatment, which were predicted to involved in trehalose biosynthesis (glucose-1-P → trehalose). The two genes were also reported to be present in C. crispus genome, as well as G. gracilis, C. merolae, and G. sulphuraria genomes [36, 41–43]. In bacteria, fungi, plants and invertebrates, UDP-glucose is a precursor for making trehalose. The production of UDP-glucose from glucose 1-phosphate can be mediated by a classic UTP-glucose-1-phosphate uridylyltransferase (UGP2, EC 220.127.116.11) [44, 45]. And the release of the glucose-1-P (G1P) from the non-reducing ends of the outer chains of polysaccharides composed of a-1,4-linked glucose residues, such as amylopectin, is catalysed by starch phosphorylase. Furthermore, 1,4-alpha-glucan branching enzyme catalyses the biosynthesis of amylopectin from amylose. It is supposed that two upregulated unigenes may promote the release of the glucose-1-P (G1P) by degradation of amylopectin-like floridean starch, which further promote accumulation of trehalose, contributing to protecting the algae against osmotic stress in intertidal environment.
Chlorophyll biosynthesis in plants is subjected to modulation by various environmental factors. For example, in etiolated rice seedlings, ionic imbalance due to salinity stress resulted in additional downregulation (41 %–45 %) of seedling dry weight, as well as chlorophyll and carotenoid levels . Meanwhile, in Helianthus annuus L. plants, salt stress more drastically affects chlorophyll synthesis (decreasing α-linolenic acid synthesis) than chlorophyllase-mediated degradation . In the current study, 7 down-regulated unigenes encoding chlorophyll biosynthetic pathway enzymes in MWL, i.e. oxygen-dependent protoporphyrinogen oxidase and magnesium chelatase subunit H were detected, including one downregulated unigene that encoded for glutamate-1-semialdehyde 2,1-aminomutase in SWL. The downregulation of these enzymes resulted in a reduction in chlorophyll synthesis, which was similar to that observed in higher plants. Previous reports have shown that photosynthetic CO2 fixation in higher plants is inhibited by water-deficit stress, and the excitation energy by chlorophyll can greatly exceed the demand of the Calvin cycle for ATP and NADPH, thus resulting in the overreduction of the electron transport chain and enhanced generation of ROS, ultimately leading to inhibition of PS II reaction centers, damage to the ATP synthesizing machinery, and a decrease in photosynthetic rate [48–51]. Based on these findings, it is therefore assumed that in Py. haitanensis, the downregulation of chlorophyll biosynthesis contributes to the inhibition of the accumulation of highly photosensitive photodynamic tetrapyrroles that generates ROS under stress conditions. Thus, the downregulation of chlorophyll biosynthesis is advantageous in Py. haitanensis because the accumulation of toxic ROS is reduced.
Apoptosis is a basic biological process that functions in various aspects of animal and plant development and in their responses to stress [52, 53]. Cysteine proteases (caspases) are key components of animal apoptosis , and the activation of cysteine proteases constitutes the critical point in the apoptosis pathway of animal cells [55, 56]. However, while animals have caspase genes, plants such as Arabidopsis and rice do not harbor orthologous caspase sequences in their genomes, whereas caspase-like activities have also been detected in plant cells . Instead, plant genomes encode several related proteins, called metacaspases, which were found to mediate apoptosis . Meanwhile, there is evidence for a mitochondrial release of apoptosis activating molecules following a moderate level 55 °C heat shock . Recently, in red algae C. crispus, several homologues of apoptosis- and programmed cell death-related genes were also identified in the genome, including metacaspases, apoptosis inducing factor, Bax inhibitor, etc . Furthermore, in Py. yezoensis, the expression of genes encoding metacaspases was found to be upregulated under high-temperature stress, revealing their important roles in Py. yezoensis acclimation to heat stress . Similar to higher plants and other red algae, classical caspase genes or transcripts were not found in Py. haitanensis. It is worth noting that one unigene (comp18960_c0) encoding metacaspase was also found to be up-regulated under desiccation stress, but none was found in rehydration. And three upregulated unigenes putatively encoding for the mitochondria-derived activator of caspases were detected under desiccation stress, functioning in the induction of apoptosis and the regulation of apoptotic process (Additional file 5: Table S3). It is therefore likely that in response to desiccation stress, the apoptogenic proteins such as mitochondria-derived activator of caspases, were released from the mitochondria, which led to the activation of metacaspase molecules, and thereby resulting in apoptosis. Pyropia and other red algal species grown in harsh intertidal habitats posses an nearly overall apoptosis picture similar of the one found in plants.
The initial rapid response of Py. haitanensis to rehydration
Rehydration-induced revival of various metabolic processes in plants requires an efficient transport machinery. Water transport is of great importance to water homeostasis in plants subjected to dehydration and rehydration. Water channel proteins (known as aquaporins) often facilitate osmotic water flow across membranes, and are critical for osmotic regulation and glycerol transport in vascular plants, which in turn can mediate CO2 permeability and photosynthetic activity [60, 61]. The expression of several aquaporins is regulated in response to environmental factors such as dehydration stress, salinity stress, and rehydration [62–64]. Two recent studies clearly demonstrated the effective involvement of aquaporins in conferring stress tolerance in higher plants [65, 66]. The aquaporins in algae have also been the subject of several reports in recent decades. Henzler and Steudle  demonstrated that algal membranes can be modeled as composite structures in which water is transported through highly selective aquaporins, while solutes permeate via other routes through the membrane. Akai et al.  further reported that in the cyanobacterium Synechocystis sp. PCC 6803, the aquaporin (AqpZ) was involved in cell volume regulation and functioned as a water transport system that responds to daily oscillations of intracellular osmolarity. The published genome of Chlamydomonas contains three putative aquaporins: MIP1 (Cre12.g549300.t1.2), MIP2 (Cre17.g711250.t1.2), and MIP3 (Cre01.g038800.t1.2) . Moreover, MIP1 was found to be localized to the contractile vacuole in Chlamydomonas, functioning in vivo as a water channel . In the red alga Porphyra purpurea, some EST contigs putatively encoding aquaporins have been detected, although information on the specific functions of these channel proteins is currently limited . Similarly, in this work, two aquaporin unigenes were found to have altered expression in response to rehydration, of which one aquaporin unigene was upregulated and the other was downregulated, whereas none was observed in differentially expressed genes between MWL and CON, as well as SWL and CON. Plants aquaporins exhibit a high diversity with many different isoforms. Molecular analyses on regulation of the whole aquaporin family have often revealed complex transcriptional and post-translational response patterns, with sometimes opposite profiles between isoforms . It seems possible that above two unigenes code for different isoforms, and the increase in expression levels of one aquaporin gene might provide Py. haitanensis additional ability to increase the water permeability of the cells under rehydration. On the other hand, the Py. haitanensis may avoid water loss by downregulating another aquaporin gene during rehydration. Thus, in Py. haitanensis, these various aquaporins were actively involved in modulating transmembrane water transport during the initial phases following rehydration when subjected to dehydration stress. Although in evolutionary viewpoint, the genome architectures of red algae are distantly related with the vascular plants , the aquaporins could represent the ancient water transport system throughout evolution of plants. Further work is clearly required to completely understand the roles and involvement of the aquaporins in the intertidal algae response to rehydration.
Additionally, solute (osmolytes) transport in and out of the vacuoles of desiccation-tolerant angiosperms occurs during desiccation and rehydration . In addition, there is a significant representation within the transporter subcategory in the rehydration transcriptome of the desiccation-tolerant bryophyte T. ruralis . Similarly, several transport-associated unigenes were exclusively detected in response to rehydration in the present study, which encode ABC transporters, ion transporters, phosphate transporters, and aquaporins, thus suggesting that these may constitute the transport machinery of Py. haitanensis under rehydration conditions. In particular, the unigenes encoding ABC transporters were downregulated during rehydration (Additional file 13: Table S11). ABC transporters have been recognized to participate in a multitude of physiological processes that allow the plant to adapt to changing environments and cope with biotic and abiotic stresses, as well as detoxification processes . Chan et al. has also identified several ABC transporters that were related to multidrug resistance, bile salt pumps, and the transport of lipids into the plastid in P. purpurea and P. umbilicalis transcriptome . Therefore, it is likely that during rehydration, the downregulation of unigenes encoding ABC transporters instantly matched the decrease in the requirement for lipids and lipophilic compounds that protect the plant from biotic and abiotic stresses, with the revival of many biological metabolism.
Facultative C3-CAM species such as Guzmania monostachia (Bromeliaceae) and Talinum triangulare (Portulacaceae) can be induced by various environmental factors such as drought stress and salinity to utilize, and then returned to a typical C3 condition after a subsequent period of seven days of rehydration . Unlike that observed during dehydration treatment, many of the unigenes associated with the C4 pathway for photosynthetic carbon reduction were also not significantly altered, whereas many unigenes in Py. haitanensis that were involved in a putative C3 pathway for photosynthetic carbon assimilation were upregulated under rehydration, suggesting the return of a typical C3 condition after a subsequent period of one half hour of rehydration. These results revealed that Py. haitanensis employed different carbon-concentrating mechanisms based on the actual environmental factor.
Genes collectively responding to dehydration and rehydration in Py. haitanensis
Unsaturated fatty acids (UFAs) have profound effects on the fluidity and function of biological membranes. Plants, animals, and microorganisms regulate the synthesis of UFAs to remodel membrane fluidity during changing environmental conditions as well as in response to nutrients . Various fatty acid desaturases (FAD) are key enzymes in the synthesis of unsaturated fatty acids. In higher plants, drought tolerance is considered to be closely correlated with the level of unsaturated fatty acids [77–79]. For example, the increase in α-linolenic acid (18:3, LA) levels in Nicotiana tabacum cells due to the overexpression of the ω-3 fatty acid desaturases FAD3 and FAD7 enhances tolerance to drought stress , whereas non-tolerant plants decline their fraction of 18:3 . In addition, several studies have been performed to explore the role of unsaturated fatty acids in red algae. Sun et al.  reported that in Py. yezoensis cold stress caused an increase in the expression of FAD to improve the proportion of polyunsaturated fatty acids. Furthermore, the amounts of polyunsaturated fatty acids, especially those of 20:4 and 20:5, were decreased significantly in conchocelis of Py. haitanensis with the increase in temperature whereas thallus exhibited only minor differences in fatty acid levels between control and the heat stress treatments with a slight decrease of 20:5 and an increase of 20:4 . In the present study, different from that in response to temperature stress, all of differentially expressed unigenes encoding fatty acid desaturases and fatty acid elongase were collectively upregulated in response to dehydration and rehydration (Table 4). Thus, in Py. haitanensis, fatty acid desaturases and fatty acid elongase were involved in the increase in blade membrane unsaturation in response to dehydration and rehydration, all of which contributes to their osmotic acclimation to the intertidal habitat.
Families of transcription factors functioned downstream of signaling cascades that were related to biological and environmental stimuli. Several members of these families were previously identified to be responsive to various stresses  such as bZIP to drought and abscisic acid (ABA), MYB to dehydration, zinc finger to cold and drought, bHLH and NAC to drought, salinity, and ABA. Recently, transcription factors such as bZIP and zinc-finger protein were present both in stressed and rehydrated samples of poplar (Populus alba L.) . Similarly, of the 20 differentially expressed unigenes encoding transcription factors in MWL, SWL, and REH as compared to CON, 6 unigenes, which included bZIP and zinc-finger protein, were collectively upregulated or downregulated in dehydration/desiccation-stressed and rehydrated samples of Py. haitanensis (Table 5). In higher plants, transcription factors belonging to the class of DRE-binding protein (DREB)/C-repeat-binding factor (CBF), bZIP, MYC, and MYB homologs are involved in ABA-dependent or ABA-independent stress signaling pathways [83–85]. Interestingly, in the present study, one unigene encoding 9-cis-epoxycarotenoid dioxygenase (NCED), which is an important enzyme in the synthesis of the phytohormone ABA , was found to be upregulated in response to moderate dehydration, severe dehydration, and rehydration, thus revealing endogenous ABA levels may be induced as a result of dehydration and rehydration. Several of the drought-related genes can be induced by ABA. For example, the Arabidopsis rd22 BP1 gene that encodes a Myc homologue transcription factor has been shown to be induced by dehydration, high-salt conditions, and ABA . In Py. haitanensis, the upregulation of several unigenes encoding bZIPs by endogenous ABA in response to dehydration and rehydration still needs further verification. In addition, research studies that examine the roles and involvement of these transcription factors in stress signaling pathways of Py. haitanensis are also warranted.
The combination of de novo transcriptome sequencing and DGE analysis based on the NGS technology is a powerful method for identifying candidate genes and key metabolic processes involved in the response to osmotic stress in Pyripoa species. Using this method, genes associated with the C4 and C3 pathways, trehalose biosynthesis, porphyrin and chlorophyll metabolism, induction of apoptosis, reproductive structure development, and other carbohydrate metabolic process were determined to be involved in the desiccation response. On the other hand, we demonstrated the mechanisms involved in the initial response to rehydration based on the following perspectives: multiple transport machinery such as aquaporins and ABC transporters, and the return of the normal C3 condition. In addition, the synthesis of unsaturated fatty acids, various transcription factor families, and molecular chaperones have collectively been implicated in the process of dehydration and rehydration in Py. haitanensis. The transcriptome data generated in the present study can serve as an important resource for understanding tolerance mechanisms of Py. haitanensis as it thrives in its unique intertidal environment.
Algal materials and cultivation conditions
To eliminate the interference caused by genotypic differences, a lab-cultured pure line PH-38 of Py. haitanensis was used in the experiments. The sporophytes (conchocelis) of this pure line were directly developed from a single somatic cell that was enzymatically isolated from a farmed thallus and then cultured in bubbling sterilized seawater with Provasoli’s enrichment solution medium (PES) under 20 μmol photons · m–2 s–1 at 24 ± 1 °C and a 12:12 light:dark (L:D) photoperiod. Gametophytes (thalli) were formed from germination of the conchospores that were released from the mature sporophytes (conchosporangia), which were cultured in running sterilized seawater with PES under 50 μmol photons · m–2 s–1 at 20 ± 1 °C and a 12:12 light:dark (L:D) photoperiod.
Experimental design and sampling
To simulate conditions in the open sea, the gametophytes were allowed to acclimate for 2 weeks in running sterilized seawater (at 20 ± 1 °C, 12:12 light:dark (L:D) photoperiod, and 1,250 μmol photons · m-2 s-1) prior to performing the experiments. Following acclimation, the algae were exposed to various abiotic stress treatments (for whole transcriptome sequencing) and dehydration/rehydration treatments (for DGE sequencing).
Due to limited Py. haitanensis genome information, we first generated a global transcriptome for this species. The materials included samples from different developmental phases and treatments (Additional file 1: Table S1). After harvesting and weighing, the samples were immediately frozen in liquid nitrogen.
For DGE analysis, the gametophytes were subjected to dehydration and rehydration by exposing to air and transferring these back to seawater. The algal samples in normal conditions (CON) were harvested before the dehydration treatments. The algal samples of moderate water loss (MWL) were collected when the algae reached a water loss level of 30 ± 5 %, whereas the algal samples of severe water loss (SWL) were collected when the algae reached a water loss level of 80 ± 5 %. For rehydration (REH), severe dehydrated algae were transferred back to normal conditions and a representative sample was collected after 0.5 h. All treatments were performed at 20 ± 1 °C and 1,250 μmol photons · m-2 s-1. Two independent replicates were used for each treatment. Water loss was determined according to Kim et al. . Algal samples from each treatment were bulked separately prior to RNA isolation. All samples for stress treatments were immediately collected and frozen in liquid nitrogen.
Total RNAs were extracted from each sample using E.Z.N.A.™ Plant RNA Kit (Omega Bio-Tek, Norcross, GA, USA), and contaminating DNA was digested with RNase-Free DNase I (TIANGEN, Beijing, China), according to the manufacturer’s instructions. The status of RNA degradation and contamination were monitored on 1 % agarose gels. RNA purity was checked using a NanoPhotometer spectrophotometer (IMPLEN, CA, USA). RNA concentration was measured using a Qubit RNA Assay Kit and a Qubit 2.0 Fluorometer (Life Technologies, CA, USA). RNA integrity was assessed using the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, CA, USA). To obtain a comprehensive list of transcripts used in transcriptome sequencing, equal amounts of total RNA from samples of different developmental phases and treatments were pooled together.
Library preparation for transcriptome and DGE sequencing
For transcriptome and DGE sequencing, a total of 6 μg of RNA per sample was used as input material. All nine RNA samples (including a pooled RNA sample for de novo transcriptome sequencing and eight RNA samples for DGE sequencing) had RIN values > 7.0. Sequencing libraries were generated using an Illumina TruSeq RNA Sample Preparation Kit (Illumina, San Diego, USA) following the manufacturer’s recommendations, and nine index codes were added to attribute sequences to each sample. Briefly, mRNA was purified from total RNA using oligo (dT) magnetic beads. Following purification, mRNA was fragmented into smaller pieces using divalent cations under elevated temperature in an Illumina proprietary fragmentation buffer. The cleaved RNA fragments were used for first-strand cDNA synthesis using random oligonucleotides and SuperScript II. Second-strand cDNA synthesis was subsequently performed using DNA polymerase I and RNase H. Remaining overhangs were converted into blunt ends via exonuclease/polymerase activities and then the enzymes were removed. After adenylation of the 3' ends of DNA fragments, Illumina PE adapter oligonucleotides were ligated to prepare for hybridization. To select cDNA fragments of preferentially 200 bp in length, the library fragments were purified using the AMPure XP system (Beckman Coulter, Beverly, USA). DNA fragments with ligated adaptor molecules on both ends were selectively enriched using Illumina PCR Primer Cocktail in a 10-cycle PCR reaction. Products were purified (AMPure XP system) and quantified using the Agilent high-sensitivity DNA assay on the Agilent Bioanalyzer 2100 system. Clustering of the index-coded samples was performed on a cBot Cluster Generation System using the TruSeq PE Cluster Kit v3-cBot-HS (Illumina, San Diego, USA), according to the manufacturer’s instructions. After cluster generation, the library preparations were sequenced on an Illumina Hiseq 2000 platform. Together, the library from the pooled RNA sample was sequenced with the 100-bp pair-end reads for de novo transcriptome analysis. The other eight libraries were sequenced with the 100-bp single-end reads for DGE analysis.
Raw data (raw reads) in fastq format were first processed using in-house Perl scripts. In this step, clean data (clean reads) were obtained by removing reads containing adapter sequences, reads containing ploy-N, and low quality reads from the raw data. Simultaneously, the Q20, Q30, GC content, and sequence duplication level of the clean data were calculated. All downstream analyses were therefore based on clean, high-quality data.
De novo assembly and annotation
After quality filtering, transcriptome de novo assembly was performed using the short reads assembling program, Trinity , with min_kmer_cov set to 2 and all other parameters set to default. The optimal assembly results were chosen according to the assembly evaluation. Clustering analysis was then performed to generate a unigene database that comprised potential alternative splicing transcripts.
The unigenes were compared against the NCBI Nr and Nt databases, as well as Swiss-Prot database using BLAST 2.2.27+, with an E-value of 1e-10, 1e-5, and 1e-5, respectively. Gene names were assigned to each unigene based on the best BLAST hit (highest score). The unigene sequences were also aligned to the KOG database to predict and classify functions using BLAST 2.2.27+, with an E-value of 1e-3. The unigene sequences were searched against the Pfam database to predict functional domain and protein family using Hmmerscan (HMMER 3.0 package), using an inclusion E-value of 0.01.
To annotate the unigenes with GO terms that described biological processes, molecular functions, and cellular components, the Nr and Pfam annotation results were imported into the Blast2GO program , a software package that retrieves GO terms and subsequently identifies and compares gene functions.
To obtain an overview of gene pathways networks, KEGG pathways were assigned to the unigenes using the online KEGG Automatic Annotation Server (KAAS; http://www.genome.jp/kegg/kaas/). The bi-directional best hit (BBH) method was used to obtain
KEGG Orthology (KO) assignment . The output of KEGG analysis included KO assignments and KEGG pathways that were populated with the KO assignments.
Reads mapping to the reference genome
Prior to mapping reads to the reference database, raw reads were transformed into clean reads as earlier described. The de novo transcriptome database was selected as reference genome. The clean reads of the samples were mapped to the reference genome, and obtained from the annotation information of each sample using the RSEM (v1.2.0) software package .
Quantification of gene expression level
HTSeq software (www-huber.embl.de/users/anders/HTSeq/) was used to determine the number of reads that mapped to each gene. To compare expression abundance among samples, read counts were normalized to the reads per kilobase of exon model per million mapped reads (RPKM)  to obtain the relative levels of expression. An RPKM value of > 0.3 was defined as the threshold of significant gene expression.
Differential expression analysis
Differential expression analysis was performed using DESeq packages  for comparisons among samples with two biological replicates. Transcript abundances for each gene were expressed as the weighted mean of read counts from two replicates normalized to the overall library size (also known as ‘base mean’). The relationship between the relative abundance of each gene in parallel libraries was statistically assessed using Pearson’s correlation coefficient. P-values (adjusted for false discovery rate) were generated for each gene in pair-wise comparisons among samples. An adjusted P-value of 0.05 was set as the threshold for significant differential expression. Variance stabilized data obtained with DESeq was used to generate the heatmaps of differentially expressed genes.
GO and KEGG enrichment analysis of differentially expressed genes
To study the biological significance of differentially expressed genes, gene ontology (GO)-based enrichment tests were conducted using GOseq . GO terms with corrected P-values < 0.05 were considered significantly enriched by differentially expressed genes.
We used the KOBAS software  to test the statistical enrichment of differentially expressed genes in KEGG pathways. Pathways with corrected P-values < 0.05 were considered significantly enriched by differentially expressed genes.
Sequence analysis of the functional gene
The unigene encoding 1,4-alpha-glucan branching enzyme (PhSBE) was selected and further analyzed using the BLAST and ORF-finder at NCBI (http://www.ncbi.nlm.nih.gov/) for the homology study and putative ORF prediction. In order to identify the intron, alignment of the ORF and the gene sequence from draft genome (unpublished data) was performed using Spidey at NCBI. The analysis of deduced protein was carried out using ExPASy tool (http://www.expasy.org), which has hyperlinks to different prediction programs, including SingnalP for single peptide prediction. The domain prediction of PhSBE was assessed using the Pfam hidden Markov model (HMM) algorithm (http://pfam.xfam.org/search). The DNAMAN program was used for multiple sequence analysis.
Quantitative real-time PCR (qRT-PCR) validation
Total RNA was extracted as described for DGE library preparation and sequencing. For the first-strand cDNA synthesis experiment, a Transcriptor First Stand cDNA Synthesis Kit (Roche) was used following the manufacturer’s instructions. The expression levels of selected genes were analyzed by qRT-PCR with a Light-Cycle® 480 Real-Time PCR System. The reactions were performed in 20 μL volumes containing 10 μL of 2 × SYBR® Premix Ex Taq (TaKaRa Biotech Co., Dalian, China), 0.6 μL of each primer (0.3 μM final concentration of each primer), 2 μL of the diluted cDNA mix, and 6.8 μL of RNA-free water. The thermal profile for qRT-PCR was 95 °C for 30s, followed by 40 cycles of 95 °C for 5 s, 57 °C for 30s and 72 °C for 30s. Melting curves for each amplicon were then analyzed to verify the specificity of each amplification reaction. No-template controls were included for each primer pair and each PCR reaction was carried out in three biological replicates. The sequences of the primers used are given in Additional file 16: Table S12. The ubiquitin conjugating enzyme (UBC) and 18 s ribosomal RNA (18 s) genes were used as internal controls [80, 93]. The 2-ΔΔCt method was used to calculate relative gene expression values .
Availability of supporting data
The data sets supporting the results of this article are available in the Sequence Read Archive (SRA) database, accessible through NCBI BioProject ID PRJNA282903 for the transcriptome data (http://www.ncbi.nlm.nih.gov/bioproject/PRJNA282903) and PRJNA283027 for the DGE data (http://www.ncbi.nlm.nih.gov/bioproject/PRJNA283027).
Digital gene expression profile
NCBI non-redundant protein
NCBI non-redundant nucleotide sequence
The manually annotated and curated protein sequence database
euKaryotic Ortholog Groups
Kyoto Encyclopedia of Genes and Genomes
Moderate water loss
Severe water loss
Differentially expressed genes
Ribulose 1,5-bisphosphate carboxylase-oxygenase
Unsaturated fatty acids
Fatty acid desaturases
Basic leucine zipper
Reads per kilobase of exon model per million mapped reads
Ubiquitin conjugating enzyme
- 18 s:
18 s ribosomal RNA
Blouin NA, Brodie JA, Grossman AC, Xu P, Brawley SH. Porphyra: a marine crop shaped by stress. Trends Plant Sci. 2011;16:29–37.
Sutherland J, Lindstrom S, Nelson W, Brodie J, Lynch M, Hwang M, et al. A new look at an ancient order: generic revision of the Bangiales. J Phycol. 2011;47:1131–51.
Mumford TF, Miura A. Porphyra as food: cultivation and economics. In: Lemby CA, Walland JR, editors. Algae and human affairs. Cambridge: Cambridge University Press; 1988. p. 87–117.
Yan X, Li L, Aruga Y. Genetic analysis of the position of meiosis in Porphyra haitanensis Chang et Zheng (Bangiales, Rhodophyta). J Appl Phycol. 2005;17:467–73.
Smith CM, Satoh K, Fork DC. The effects of osmotic tissue dehydration and air drying on morphology and energy transfer in two species of Porphyra. Plant Physiol. 1986;80:843–7.
Sibbald PR, Vidaver W. Photosystem I-mediated regulation of water splitting in the red alga, Porphyra sanjuanensis. Plant Physiol. 1987;84:1373–7.
Cabello-Pasini A, Diaz M, Muñiz-Salazar R, Zertuche-Gonzalez JA, Pacheco-Ruiz I. Effect of temperature and desiccation on the photosynthetic performance of Porphyra perforata. J Phycol. 2000;36:10–4.
Ji Y, Tanaka J. Effect of desiccation on the photosynthesis of seaweeds from the intertidal zone in Honshu, Japan. Phycological Res. 2002;50:145–53.
Craigie J. Cell walls. In: Cole KM, Sheath RG, editors. Biology of the red algae. Cambridge: Cambridge University Press; 1990. p. 221–57.
Reed RH, CollinS JC, Russell G. The effects of salinity upon galactosyl-glycerol content and concentration of the marine red alga Porphyra purpurea (Roth) C.Ag. J Exp Bot. 1980;31:1539–54.
Abe S, Kurashima A, Yokohama Y, Tanaka J. The cellular ability of desiccation tolerance in Japanese intertidal seaweeds. Botanica Marina. 2001;44:125–31.
Contreras-Porcia L, Thomas D, Flores V, Correa JA. Tolerance to oxidative stress induced by desiccation in Porphyra columbina (Bangiales, Rhodophyta). J Exp Bot. 2011;62:1815–29.
Liu YC. Mechanism for diferential desiccation tolerance in Porphyra species. PhD thesis. Northeastern University, Department of Biology; 2009.
Jia Z, Niu J, Huan L, Wu X, Wang G, Hou Z. Cyclophilin participates in responding to stress situations in Porphyra haitanensis (Bangiales, Rhodophyta). J Phycol. 2013;49:194–201.
Fan X, Fang Y, Hu S, Wang G. Generation and analysis of 5318 expressed sequence tags from the filamentous sporophyte of Porphyra haitanensis (Rhodophyta). J Phycol. 2007;43:1287–94.
Collén J, Guisle-Marsollier I, Léger JJ, Boyen C. Response of the transcriptome of the intertidal red seaweed Chondrus crispus to controlled and natural stresses. New Phytol. 2007;176:45–55.
Dittami SM, Scornet D, Petit JL, Segurens B, Da Silva C, Corre E, et al. Global expression analysis of the brown alga Ectocarpus siliculosus (Phaeophyceae) reveals large-scale reprogramming of the transcriptome in response to abiotic stress. Genome Biol. 2009;10:R66.
Pearson GA, Hoarau G, Lago-Leston A, Coyer JA, Kube M, Reinhardt R, et al. An expressed sequence tag analysis of the intertidal brown seaweeds Fucus serratus (L.) and F. vesiculosus (L.) (Heterokontophyta, Phaeophyceae) in response to abiotic stressors. Mar Biotechnol. 2010;12:195–213.
Yang H, Mao Y, Kong F, Yang G, Ma F, Wang L. Profiling of the transcriptome of Porphyra yezoensis with Solexa sequencing technology. Chin Sci Bull. 2011;56:2119–30.
Chan CX, Zauner S, Wheeler G, Grossman AR, Prochnik SE, Blouin NA, et al. Analysis of Porphyra membrane transporters demonstrates gene transfer among photosynthetic eukaryotes and numerous sodium-coupled transport systems. Plant Physiol. 2012;158:2001–12.
Chan CX, Blouin NA, Zhuang Y, Zäuner S, Prochnik SE, Lindquist E, et al. Porphyra (Bangiophyceae) transcriptomes provide insights into red algal development and metabolism. J Phycol. 2012;48:1328–42.
Xie C, Li B, Xu Y, Ji D, Chen C. Characterization of the global transcriptome for Pyropia haitanensis (Bangiales, Rhodophyta) and development of cSSR markers. BMC Genomics. 2013;14:107.
Choi S, Hwang MS, Im S, Kim N, Jeong WJ, Park EJ, et al. Transcriptome sequencing and comparative analysis of the gametophyte thalli of Pyropia tenera under normal and high temperature conditions. J Appl Phycol. 2013;25:1237–46.
Wang L, Mao Y, Kong F, Li G, Ma F, Zhang B, et al. Complete sequence and analysis of plastid genomes of two economically important red algae: Pyropia haitanensis and Pyropia yezoensis. PLoS One. 2013;8, e65902.
Mao Y, Zhang B, Kong F, Wang L. The complete mitochondrial genome of Pyropia haitanensis Chang et Zheng. Mitochondrial DNA. 2012;23:344–6.
Kong F, Sun P, Cao M, Wang L, Mao Y. Complete mitochondrial genome of Pyropia yezoensis: reasserting the revision of genus Porphyra. Mitochondrial DNA. 2014;25:335–6.
Nakamura Y, Sasaki N, Kobayashi M, Ojima N, Yasuike M, Shigenobu Y, et al. The first symbiont-free genome sequence of marine red alga, Susabi-nori (Pyropia yezoensis). PLoS One. 2013;8, e57122.
Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, et al. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat Biotech. 2011;29:644–52.
Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics. 2011;12:323.
Anders S, Huber W. Differential expression analysis for sequence count data. Genome Biol. 2010;11:R106.
Taiz L, Zeiger E. Plant physiology. 5th ed. Sunderland: Sinauer Associates Inc; 2010.
Matsuoka M, Furbank RT, Fukayama H, Miyao M. Molecular engineering of C4 photosynthesis. Annu Rev Plant Physiol Plant Mol Biol. 2001;52:297–314.
Langdale JA. C4 cycles: past, present, and future research on C4 photosynthesis. Plant Cell. 2011;23:3879–92.
Reinfelder JR, Milligan AJ, Morel FM. The role of the C4 pathway in carbon accumulation and fixation in a marine diatom. Plant Physiol. 2004;135:2106–11.
Derelle E, Ferraz C, Rombauts S, Rouzé P, Worden AZ, Robbens S, et al. Genome analysis of the smallest free-living eukaryote Ostreococcus tauri unveils many unique features. Proc Natl Acad Sci U S A. 2006;103:11647–52.
Collen J, Porcel B, Carre W, Ball SG, Chaparro C, Tonon T, et al. Genome structure and metabolic features in the red seaweed Chondrus crispus shed light on evolution of the Archaeplastida. Proc Natl Acad Sci U S A. 2013;110:5247–52.
Liu MS, Chien CT, Lin TP. Constitutive components and induced gene expression are involved in the desiccation tolerance of Selaginella tamariscina. Plant Cell Physiol. 2008;49:653–63.
Smirnoff N. The carbohydrates of bryophytes in relation to desiccation tolerance. J Bryol. 1992;17:185–91.
Goddijn O, Smeekens S. Sensing trehalose biosynthesis in plants. Plant J. 1998;14:143–6.
Goddijn OJM, van Dun K. Trehalose metabolism in plants. Trends Plant Sci. 1999;4:315–9.
Lluisma AO, Ragan MA. Cloning and characterization of a nuclear gene encoding a starch-branching enzyme from the marine red alga Gracilaria gracilis. Curr Genet. 1998;34:105–11.
Nozaki H, Takano H, Misumi O, Terasawa K, Matsuzaki M, Maruyama S, et al. A 100 %-complete sequence reveals unusually simple genomic features in the hot-spring red alga Cyanidioschyzon merolae. BMC Biol. 2007;5:28.
Schönknecht G, Chen WH, Ternes CM, Barbier GG, Shrestha RP, Stanke M, et al. Extremophilic eukaryote gene transfer from bacteria and archaea facilitated evolution of an extremophilic eukaryote. Science. 2013;339:1207–10.
Kleczkowski LA, Geisler M, Ciereszko I, Johansson H. UDP-glucose pyrophosphorylase. An old protein with new tricks. Plant Physiol. 2004;134:912–8.
Meng M, Geisler M, Johansson H, Harholt J, Scheller HV, Mellerowicz EJ, et al. UDP-glucose pyrophosphorylase is not rate limiting, but is essential in Arabidopsis. Plant Cell Physiol. 2009;50:998–1011.
Turan S, Tripathy BC. Salt-stress induced modulation of chlorophyll biosynthesis during de-etiolation of rice seedlings. Physiol Plant. 2015;153:477–91.
Santos CV. Regulation of chlorophyll biosynthesis and degradation by salt stress in sunflower leaves. Sci Hortic. 2004;103:93–9.
Lawlor DW, Cornic G. Photosynthetic carbon assimilation and associated metabolism in relation to water deficits in higher plants. Plant Cell Environ. 2002;25:275–94.
Smirnoff N. The role of active oxygen in the response of plants to water deficit and desiccation. New Phytol. 1993;125:27–58.
Møller IM, Jensen PE, Hansson A. Oxidative modifications to cellular components in plants. Annu Rev Plant Biol. 2007;58:459–81.
Yuan S, Liu WJ, Zhang NH, Wang MB, Liang HG, Lin HH. Effects of water stress on major photosystem II gene expression and protein metabolism in barley leaves. Physiol Plant. 2005;125:464–73.
Wang H, Li J, Bostock RM, Gilchrist DG. Apoptosis: a functional paradigm for programmed plant cell death induced by a host-selective phytotoxin and invoked during development. Plant Cell. 1996;8:375–91.
Martins LM, Earnshaw WC. Apoptosis: alive and kicking in 1997. Trends Cell Biol. 1997;7:111–4.
Bonneau L, Ge Y, Drury GE, Gallois P. What happened to plant caspases? J Exp Bot. 2008;59:491–9.
Martin SJ, Green DR. Protease activation during apoptosis: death by a thousand cuts? Cell. 1995;82:349–52.
Martins LM, Kottke T, Mesner PW, Basi GS, Sinha S, Frigon N, et al. Activation of multiple interleukin-1β converting enzyme homologues in cytosol and nuclei of HL-60 cells during etoposide-induced apoptosis. J Biol Chem. 1997;272:7421–30.
Coll NS, Vercammen D, Smidler A, Clover C, Van Breusegem F, Dangl JL, et al. Arabidopsis type I metacaspases control cell death. Science. 2010;330:1393–7.
Balk J, Chew SK, Leaver CJ, McCabe PF. The inter membrane space of plant mitochondria contains a DNase activity that may be involved in programmed cell death. Plant J. 2003;34:573–83.
Sun P, Mao Y, Li G, Cao M, Kong F, Wang L, et al. Comparative transcriptome profiling of Pyropia yezoensis (Ueda) M.S. Hwang & H.G. Choi in response to temperature stresses. BMC Genomics. 2015;16:463.
Maurel C. Aquaporins and water permeability of plant membranes. Annu Rev Plant Physiol Plant Mol Biol. 1997;48:399–429.
Uehlein N, Lovisolo C, Siefritz F, Kaldenhoff R. The tobacco aquaporin NtAQP1 is a membrane CO2 pore with physiological functions. Nature. 2003;425:734–7.
Johansson I, Karlsson M, Johanson U, Larsson C, Kjellbom P. The role of aquaporins in cellular and whole plant water balance. Biochim Biophys Acta. 2000;1465:324–42.
Rodriguez MC, Edsgard D, Hussain SS, Alquezar D, Rasmussen M, Gilbert T, et al. Transcriptomes of the desiccation-tolerant resurrection plant Craterostigma plantagineum. Plant J. 2010;63:212–28.
Oliver MJ, Dowd SE, Zaragoza J, Mauget SA, Payton PR. The rehydration transcriptome of the desiccation-tolerant bryophyte Tortula ruralis: transcript classification and analysis. BMC Genomics. 2004;5:89.
Sade N, Gebretsadik M, Seligmann R, Schwartz A, Wallach R, Moshelion M. The role of tobacco aquaporin1 in improving water use efficiency, hydraulic conductivity, and yield production under salt stress. Plant Physiol. 2010;152:245–54.
Sade N, Vinocur BJ, Diber A, Shatil A, Ronen G, Nissan H, et al. Improving plant stress tolerance and yield production: is the tonoplast aquaporin SlTIP2;2 a key to isohydric to anisohydric conversion? New Phytol. 2009;181:651–61.
Henzler T, Steudle B. Reversible closing of water channels in Chara internodes provides evidence for a composite transport model of the plasma membrane. J Exp Bot. 1995;46:199–209.
Akai M, Onai K, Morishita M, Mino H, Shijuku T, Maruyama H, et al. Aquaporin AqpZ is involved in cell volume regulation and sensitivity to osmotic stress in Synechocystis sp. strain PCC 6803. J Bacteriol. 2012;194:6828–36.
Goodstein DM, Shu S, Howson R, Neupane R, Hayes RD, Fazo J, et al. Phytozome: a comparative platform for green plant genomics. Nucleic Acids Res. 2012;40:D1178–86.
Komsic-Buchmann K, Wöstehoff L, Becker B. The contractile vacuole as a key regulator of cellular water flow in Chlamydomonas reinhardtii. Eukaryot Cell. 2014;13:1421–30.
Li G, Santoni V, Maurel C. Plant aquaporins: Roles in plant physiology. Biochim Biophys Acta. 1840;2014:1574–82.
Smith DR, Hua J, Lee RW, Keeling PJ. Relative rates of evolution among the three genetic compartments of the red alga Porphyra differ from those of green plants and do not correlate with genome architecture. Mol Phylogenet Evol. 2012;65:339–44.
Ingram J, Bartels D. The molecular basis of dehydration tolerance in plants. Annu Rev Plant Physiol Plant Mol Biol. 1996;47:377–403.
Kretzschmar T, Burla B, Lee Y, Martinoia E, Nagy R. Functions of ABC transporters in plants. Essays Biochem. 2011;50:145–60.
Freschi L, Takahashi CA, Cambui CA, Semprebom TR, Cruz AB, Mioto PT, et al. Specific leaf areas of the tank bromeliad Guzmania monostachia perform distinct functions in response to water shortage. J Plant Physiol. 2010;167:526–33.
Aguilar PS, De Mendoza D. Control of fatty acid desaturation: a mechanism conserved from bacteria to humans. Mol Microbiol. 2006;62:1507–14.
Zhang M, Barg R, Yin M, Gueta-Dahan Y, Leikin-Frenkel A, Salts Y, et al. Modulated fatty acid desaturation via overexpression of two distinct ω-3 desaturases differentially alters tolerance to various abiotic stresses in transgenic tobacco cells and plants. Plant J. 2005;44:361–71.
Upchurch RG. Fatty acid unsaturation, mobilization, and regulation in the response of plants to stress. Biotechnol Lett. 2008;30:967–77.
Torres-Franklin ML, Repellin A, Huynh VB, d’Arcy-Lameta A, Zuily-Fodil Y, Pham-Thi AT. Omega-3 fatty acid desaturase (FAD3, FAD7, FAD8) gene expression and linolenic acid content in cowpea leaves submitted to drought and after rehydration. Environ Exp Bot. 2009;65:162–9.
Luo Q, Zhu Z, Zhu Z, Yang R, Qian F, Chen H, et al. Different responses to heat shock stress revealed heteromorphic adaptation strategy of Pyropia haitanensis (Bangiales, Rhodophyta). PLoS One. 2014;9, e94354.
Shameer K, Ambika S, Varghese SM, Karaba N, Udayakumar M, Sowdhamini R. STIFDB-Arabidopsis Stress Responsive Transcription Factor DataBase. Int J Plant Genomics. 2009;2009:583429.
Berta M, Giovannelli A, Sebastiani F, Camussi A, Racchi ML. Transcriptome changes in the cambial region of poplar (Populus alba L.) in response to water deficit. Plant Biol. 2010;12:341–54.
Abe H, Yamaguchi-Shinozaki K, Urao T, Iwasaki T, Hosokawa D, Shinozaki K. Role of Arabidopsis MYC and MYB homologs in drought- and abscisic acid-regulated gene expression. Plant Cell. 1997;9:1859–68.
Wang W, Vinocur B, Altman A. Plant responses to drought, salinity and extreme temperatures: towards genetic engineering for stress tolerance. Planta. 2003;218:1–14.
Sakuma Y, Maruyama K, Osakabe Y, Qin F, Seki M, Shinozaki K, et al. Functional analysis of an Arabidopsis transcription factor, DREB2A, involved in drought-responsive gene expression. Plant Cell. 2006;18:1292–309.
Chernys JT, Zeevaart JAD. Characterization of the 9-cis-epoxycarotenoid dioxygenase gene family and the regulation of abscisic acid biosynthesis in avocado. Plant Physiol. 2000;124:343–53.
Kim JK, Kraemer GP, Yarish C. Comparison of growth and nitrate uptake by New England Porphyra species from different tidal elevations in relation to desiccation. Phycological Res. 2009;57:152–7.
Gotz S, Garcia-Gomez JM, Terol J, Williams TD, Nagaraj SH, Nueda MJ, et al. High-throughput functional annotation and data mining with the Blast2GO suite. Nucleic Acids Res. 2008;36:3420–35.
Moriya Y, Itoh M, Okuda S, Yoshizawa AC, Kanehisa M. KAAS: an automatic genome annotation and pathway reconstruction server. Nucleic Acids Res. 2007;35:W182–5.
Mortazavi A, Williams BA, McCue K, Schaeffer L, Wold B. Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nat Meth. 2008;5:621–8.
Young MD, Wakefield MJ, Smyth GK, Oshlack A. Gene ontology analysis for RNA-seq: accounting for selection bias. Genome Biol. 2010;11:R14.
Mao X, Cai T, Olyarchuk JG, Wei L. Automated genome annotation and pathway identification using the KEGG Orthology (KO) as a controlled vocabulary. Bioinformatics. 2005;21:3787–93.
Li B, Chen C, Xu Y, Ji D, Xie C. Validation of housekeeping genes as internal controls for studying the gene expression in Pyropia haitanensis (Bangiales, Rhodophyta) by quantitative real-time PCR. Acta Oceanol Sin. 2014;33(9):152–9.
Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2−ΔΔCT method. Methods. 2001;25(4):402–8.
This work was supported by National Natural Science Foundation of China (Grant No. 31372517 and 30972247), National High Technology Research and Development Program of China (Grant No. 2012AA10A406, 2012AA10A411, and 2012AA10A401), Independent Innovation Foundation of Shandong Province (Grant No. 2013CXC80202), Public Science and Technology Research Funds Project of Ocean (Grant No. 201105021-8 and 201305030), National Infrastructure of Fishery Germplasm Resources, Fundamental Research Funds for the Central Universities (Grant No. DC201502070401), and Scientific Research Foundation for the Introduction of Talent of Dalian Nationalities University.
The authors declare that they have no competing interests.
LW and YM conceived and designed the experiment. LW performed the experiments and data analysis. YM and FK contributed by planning, supervising and financing the work. PS and MC helped to prepare the reagents and materials. LW and YM drafted and revised the manuscript. All authors read and approved the final manuscript.
Treatments and developmental stages for transcriptome sampling and summary of the Py. haitanensis transcriptome. (DOCX 17 kb)
Log-log plot showing the dependence of unigene lengths on the number of reads assembled into each unigene. (PDF 1390 kb)
Summary of clean reads aligned uniquely to the reference transcriptome for each library. (DOCX 13 kb)
Scatterplot showing the correlation of expression levels between two replicates. (PDF 2576 kb)
(A) The comparison of expressed gene lists between MWL and CON. (B) Differentially expressed genes between MWL and CON. Readcount: transcript abundances for each gene. Pval: P value. Padj: adjusted P-value. Adjusted P-value of < 0.05 were used as the threshold to judge the significance of gene expression difference. (XLSX 2517 kb)
(A) The comparison of expressed gene lists between SWL and CON. (B) Differentially expressed genes between SWL and CON. Readcount: transcript abundances for each gene. Pval: P value. Padj: adjusted P-value. Adjusted P-value of < 0.05 were used as the threshold to judge the significance of gene expression difference. (XLSX 2517 kb)
(A) GO enrichment analysis of the upregulated genes in MWL compared to CON. (B) GO enrichment analysis of the downregulated genes in MWL compared to CON. (XLSX 436 kb)
(A) GO enrichment analysis of the upregulated genes in SWL compared to CON. (B) GO enrichment analysis of the downregulated genes in SWL compared to CON. (XLSX 429 kb)
(A) KEGG pathway enrichment analysis of the upregulated genes in MWL compared to CON. (B) KEGG pathway enrichment analysis of the downregulated genes in MWL compared to CON. (XLSX 33 kb)
KEGG pathway enrichment analysis of the differentially expressed genes in SWL compared to CON. (XLSX 26 kb)
(A) The comparison of expressed gene lists between REH and CON. (B) Differentially expressed genes between REH and CON. Readcount: transcript abundances for each gene. Pval: P value. Padj: adjusted P-value. Adjusted P-value of < 0.05 were used as the threshold to judge the significance of gene expression difference. (XLSX 2323 kb)
(A) GO enrichment analysis of the upregulated genes in REH compared to CON. (B) GO enrichment analysis of the downregulated genes in REH compared to CON. (XLSX 334 kb)
Differentially expressed genes encoding ABC transporters between REH and CON based on KO database annotation. (DOCX 15 kb)
(A) Nucleotide and deduced amino acid sequence of cDNA encoding 1,4-alpha-glucan branching enzyme (PhSBE). Nucleotide residues are numbered in the 5’ to 3’ direction. Amino acid residues are numbered in the N- to C-terminal direction. The start codon is marked in red color font. An asterick denotes the stop codon. (B) Genomic sequence structure of the PhSBE ORF (5’ to 3’ direction). 1 (arrow) represents start site of translation, 2767 represents end site of translation. Exon 1 (1-28), exon 2 (320-401) and exon 3 (610-2767) are depicted as red boxes. Intron 1 (29-319) and intron 2 (402-609) are depicted as narrow blue boxes. (PDF 170 kb)
Multiple sequence alignment of 1,4-alpha-glucan branching enzymes from different species. 1,4-alpha-glucan branching enzymes were chosen from different species and aligned by DNAMAN. Sequence in frame means Alpha-amylase_C domain. The accession numbers of these sequences were showed as follows: Chondrus crispus (XP_005716136.1), Gracilaria gracilis (AAB97471.1), Cyanidioschyzon merolae strain 10D (XP_005536101.1), Galdieria sulphuraria (XP_005703889.1), Arabidopsis thaliana (NP_195985.3), Zea mays (NP_001105316.1) and Oryza sativa (BAA82828.1). (PDF 12907 kb)
Primers used in qRT-PCR for validating differentially expressed genes. (DOCX 14 kb)