Skip to main content
Download PDF

Comparative transcriptome profiling of Pyropia yezoensis (Ueda) M.S. Hwang & H.G. Choi in response to temperature stresses

BMC Genomics volume 16, Article number: 463 (2015) Cite this article

Abstract

Background

Pyropia yezoensis is a model organism often used to investigate the mechanisms underlying stress tolerance in intertidal zones. The digital gene expression (DGE) approach was used to characterize a genome-wide comparative analysis of differentially expressed genes (DEGs) that influence the physiological, developmental or biochemical processes in samples subjected to 4 treatments: high-temperature stress (HT), chilling stress (CS), freezing stress (FS) and normal temperature (NT).

Results

Equal amounts of total RNAs collected from 8 samples (two biological replicates per treatment) were sequenced using the Illumina/Solexa platform. Compared with NT, a total of 2202, 1334 and 592 differentially expressed unigenes were detected in HT, CS and FS respectively. Clustering analysis suggested P. yezoensis acclimates to low and high-temperature stress condition using different mechanisms: In heat stress, the unigenes related to replication and repair of DNA and protein processing in endoplasmic reticulum were active; however at low temperature stresses, unigenes related to carbohydrate metabolism and energy metabolism were active. Analysis of gene differential expression showed that four categories of DEGs functioning as temperature sensors were found, including heat shock proteins, H2A, histone deacetylase complex and transcription factors. Heat stress caused chloroplast genes down-regulated and unigenes encoding metacaspases up-regulated, which is an important regulator of PCD. Cold stress caused an increase in the expression of FAD to improve the proportion of polyunsaturated fatty acids. An up-regulated unigene encoding farnesyl pyrophosphate synthase was found in cold stress, indicating that the plant hormone ABA also played an important role in responding to temperature stress in P. yezoensis.

Conclusion

The variation of amount of unigenes and different gene expression pattern under different temperature stresses indicated the complicated and diverse regulation mechanism in response to temperature stress in P. yezoensis. Several common metabolism pathways were found both in P. yezoensis and in higher plants, such as FAD in low-temperature stress and HSP in heat stress. Meanwhile, many chloroplast genes and unigene related to the synthesis of abscisic acid were detected, revealing its unique temperature-regulation mechanism in this intertidal species. This sequencing dataset and analysis may serve as a valuable resource to study the mechanisms involved in abiotic stress tolerance in intertidal seaweeds.

Background

Because plants lack the ability of locomotion, they are exposed to various environmental stresses. Temperature stress is the most common type of stress to which plants are subjected. Temperature-related stress can occur at (a) temperatures below freezing, (b) low temperatures above freezing, and (c) high temperatures [1]. Stressful temperature conditions can damage the enzymes needed for photosynthesis, respiration, and protein synthesis and so affect the growth and development of plants [2]. Consequently, plants have evolved mechanisms to monitor their environments and to respond with cellular, physiological, and developmental changes to optimize growth and reproductive success.

Intertidal seaweeds inhabit an inherently stressful environment with rapidly changing physical conditions. Seaweed in this zone undergoes extreme environmental changes that include desiccation, osmotic shock, exposure to intense sunlight, and high and freezing temperatures. A better understanding of the mechanisms involved in abiotic stress tolerance in seaweeds could help elucidate their successful survival, reproduction, and distribution in the intertidal region. Transcriptome analysis is an efficient way of achieving this goal. A transcriptome is a complete set of the transcripts in a cell that becomes active at a specific developmental stage or under certain physiological conditions. Transcriptomes provide information that can be used to identify the functional elements of the genome and the molecular constituents of cells and tissues [3]. Digital gene expression (DGE), which is based on the sequencing of genome-wide expression profiles, is an efficient method that can be used to analyze transcriptome data and so identify, quantify, and annotate expressed genes at the genome level even without prior sequence knowledge. This allows for higher confidence in target discovery and pathway studies. Currently, this technique is widely used in higher plant research. In the Chinese bayberry (Myrica rubra), it has been used to examine gene expression in developing bayberry fruit. Results showed energy-related metabolism to be enhanced and all genes involved in anthocyanin biosynthesis to be up-regulated during the fruit ripening processes [4]. In Lycoris sprengeri, DGE was performed to evaluate differential gene expression between bulbs and bulblets and determine the biological and molecular mechanisms underlying bulb development [5]. In Vitis vinifera, it was used to describe how plant transcriptomes change during three developmental stages, post setting, veraison, and ripening [6].

Pyropia yezoensis (Ueda), M.S. Hwang & H.G. Choi, one of the most economically important marine crops, is widely cultivated in China, Japan, and Korea, with an annual harvest of more than 1 million tons (fresh weight) and a value of over U.S. $1.5 billion per year (http://www.fao.org/fishery/statistics/en). P. yezoensis is naturally distributed in the intertidal zone of the temperate region in the northern hemisphere. In this region, the temperature may change dramatically between seawater and air, especially in the changeover seasons between autumn and winter as well as winter and spring. The thallus of P. yezoensis is totally submerged in the water during high tide but exposed to air during low tide. Consequently, the thallus may suffer from the stress of high or low temperature and abrupt temperature changes. This makes this alga an ideal research model for investigations of the mechanisms underlying temperature stress tolerance in intertidal seaweed.

The effect of temperature on algae has been studied. In Chondrus crispus it has been reported that plants subjected to high natural stress have more differentially expressed genes and more potential marker genes, and express more antioxidative genes and HSPs. Many of these up-regulated genes are stress genes [7]. Gracilaria cornea showed photosynthetic and respiratory responses adapted to oceanic salinities and subtropical to tropical water temperatures [8]. Studies of Pyropia have been conducted at the physiological, genomic, and proteomic level. For example, Kayama et al. reported that water temperature can affect the fatty acid composition of Pyropia [9]. Choi et al. reported that most of the transcripts produced by Pyropia tenera under high-temperature conditions were from the heat shock protein family [10]; Xu et al. conducted a comparative proteomic analysis of Pyropia haitanensis in response to high-temperature stress using liquid chromatography-tandem mass spectroscopy and database searching, indicating that the algal blades resisted high-temperature stress by inhibiting photosynthesis and other nonessential metabolic processes [11].

In this paper, gene expression and regulation in response to temperature stresses were examined by performing genome-wide high-throughput transcriptomic sequencing for P. yezoensis. In view of the unique intertidal environment, four treatments were used: high-temperature stress, chilling stress, freezing stress and normal temperature. The DGE approach was used to comparatively analyze the expression patterns and characterize the genes that influenced physiological, developmental, and biochemical processes on a genome-wide basis when P. yezoensis was subjected to extreme temperature fluctuations. This sequencing dataset and analysis may serve as a valuable resource for identifying the key genes and pathways involved in responses to temperature stresses in P. yezoensis. This will lay the foundation for identifying the mechanisms underlying extreme temperature tolerance and may provide information needed for the genetic breeding of high- and low-temperature tolerant varieties.

Results

Global transcriptome assembly and annotation

A global transcriptome was generated using samples from different developmental phases and treatments to obtain as many functional gene transcripts as possible. A total of 1.34 × 107 quality paired-end reads were obtained after filtering out low-quality data (tags containing unknown base N and only adaptor tags). The GC content of the transcriptome was 63.2%. After assembly and annotation, a total of 18,640 unigenes with a mean length of 527 bp and an N50 value of 641 bp were obtained. Unigenes with lengths between 200 and 500 bp were overrepresented, comprising approximately 69.4% of the total number of unigenes; the next most abundant size class was 500–1000 bp, constituting approximately 19.7% of the total unigenes. There were fewer unigenes with lengths of 1–2 kbp and more than 2 kbp, i.e., approximately 8.6% and 2.3% of the total, respectively. The longest unigene was 8372 bp (Figure 1).

Figure 1
figure 1

Size distribution of the assembled unigenes.

Full size image

BLAST analysis was performed on all 18,640 unigenes using the following databases: NCBI non-redundant protein sequences (Nr), NCBI non-redundant nucleotide sequence (Nt), Swiss-Prot, Kyoto Encyclopedia of Genes and Genomes database (KEGG), Protein family (PFAM), and Clusters of Orthologous Groups of proteins (KOG). There were 10,597 (56.9%) unigenes with homologous sequences in at least one of these databases. Among them, 5937 (31.85%), 4845 (26.0%), and 3225 (17.3%) unigenes were found in NR, Swiss-Prot, and KEGG, respectively (Table 1).

Table 1 Summary of assembly and annotation results for P. yezoensis using Trinity
Full size table

Sequencing and annotation of eight DGE libraries

A total of eight RNA samples, generated from two biological replicates of P. yezoensis under the four treatments were subjected to RNA-seq (Table 2). Approximately 11.1–12.6 million raw reads were sequenced per sample. After filtering out low-quality data (tags containing unknown base N and only adaptor tags), approximately 10.7–12.1 million clean reads remained in each library. Approximately 73.1–84.5% reads in the eight DGE libraries were mapped to the global transcriptome of P. yezoensis (Table 2), suggesting that the transcriptome was a reliable reference. To evaluate the reproducibility of DEG library sequencing, a Pearson correlation analysis was performed for every two replicates. The square of the Pearson correlation coefficient (R2) was greater than 0.92, indicating both operational stability and the reliability of the experimental results (Figure 2).

Table 2 Summary of sequencing and mapping results for eight DGE libraries
Full size table
Figure 2
figure 2

Correlation tests for the replicates. The abscissa represents the value log10 (RPKM + 1) of one duplicate; the ordinate represents the value log10 (RPKM + 1) of the other duplicate. R2 is the square of Pearson Correlation Coefficient.

Full size image

Global analysis of the unigenes in the four temperature treatment groups

DESeq analysis was performed to identify the genes that had been differentially expressed between the different temperature treatment groups. A total of 17,399 expressed unigenes were obtained, of which 16,296 were found in all samples (Figure 3). The comparison showed 171 unigenes to be unique to the high-temperature treatment group, 51 to the chilling group, and 32 to the freezing group.

Figure 3
figure 3

Venn diagram showing the number of expressed unigenes in the NT, HT, CS, and FS treatments. This Venn diagram was drawn using the number of genes expressed (RPKM > 0.3) in the four treatments (FS, HT, CS, and NT).

Full size image

A hierarchical cluster was used to determine the profiles of the differentially expressed unigenes among the four temperature stress treatments (Figure 4). Hierarchical clustering showed that the differentially expressed unigenes could be divided into 8 subclusters based on the modulation of their expression, representing 8 different expression models (Figure 4). Unigenes with the same or similar expression patterns were gathered in the same cluster. The overall pattern of transcript changes in the freezing temperature treatment at −8°C and the chilling temperature at 0°C were similar to that of the normal temperature at 8°C. The pattern under the high temperature treatment at 24°C was markedly different from that of the other three treatment groups (Figure 4A). The results clearly indicate that P. yezoensis used different mechanisms in response to high- and low-temperature stresses, but it used similar mechanisms in response to the chilling and freezing treatments.

Figure 4
figure 4

Clustering analyses of DEGs in the four temperature treatments. Clusters were obtained using the hierarchical-means method based on the 5592 DEGs (the intersection set of the DEGs in the four treatments). (A) Heatmap of hierarchical clustering. Each column represents a treatment, and each row represents a unigene. Differences in expression were shown in different colors. Negative numbers indicate down-regulated and positive number means up-regulated. (B) Hierarchical clustering figure of differentially expressed genes. The “x” axis represents the treatments. The “y” axis represents the value of the relative expression level (log2 (HT, CS, FS/ NT)). In the subcluster, the gray line represents the value of the relative expression levels, the blue line represents the mean value of the relative expression level, and the red line is the reference. The subclusters in Figure B represent the models in Figure A from the top down.

Full size image

In subcluster 1, the unigenes were found to be closely related to translation of genetic information. These included genes related to the pathway of ribosome and RNA transport. In subcluster 2, the pathway associated with carbohydrate metabolism and energy metabolism were significantly enriched, including genes for the pentose phosphate pathway, glycolysis and carbon fixation in photosynthetic organisms. In subcluster 3, the unigenes related to the synthesis and decomposition of carbohydrates were affected in a manner similar to those in subcluster 2, except that there were more unigenes related to amino acid metabolism. In subcluster 4, subcluster 5, and subcluster 8, the unigenes were not concentrated but rather related to many processes (see Additional file 1). In subcluster 6, the unigenes were up-regulated in response to high-temperature stress, and the pathway associated with replication and repair of DNA were significantly enriched, showing the same expression trends as in subcluster 7. These unigenes were closely related to protein processing in the endoplasmic reticulum.

Number of differential expression unigenes

Differentially expressed unigenes with a q-value ≤ 0.05 [12] between pairwise combinations of the treatments were identified. A total of 2202, 1334, and 592 significantly altered unigenes were detected between NT and HT, NT and CS, and NT and FS, respectively. A total of 74 unigenes showed expression different from NT in all three treatment groups (Figure 5).

Figure 5
figure 5

Venn diagram of DEGs under the control and temperature-stressed conditions. This Venn diagram was drawn using the number of different expression genes in the three stress treatments (FS, HT, and CS) relative to the control group (NT).

Full size image

Gene Ontology functional analysis of DEGs

Gene Ontology functional analysis of DEGs was performed to describe the properties of genes and their products in P. yezoensis subjected to different temperature stresses. For the biological process category, the two mostly highly represented terms among the 19 level-2 categories were cellular process and metabolic process; for the cellular component category, the mostly highly represented terms among the 13 level-2 categories were cell, cell part, macromolecular complex and organelle; for the molecular function category, the two mostly highly represented terms among the 8 level-2 categories were binding and catalytic activity (Figure 6).

Figure 6
figure 6

Comparative distribution of the GO terms in the four treatment groups.

Full size image

DEGs in response to chilling treatment

A total of 1334 differentially expressed unigenes were identified in CS, relative to NT. These included 762 up-regulated unigenes and 572 down-regulated unigenes (Figure 7). The top 100 down- and up-regulated genes are shown in Additional files 2 and 3. Based on KEGG analysis, the up-regulated unigenes were highly enriched in ribosome biogenesis in eukaryotes, alanine, aspartate and glutamate metabolism, pyrimidine metabolism, and biosynthesis of unsaturated fatty acids. The down-regulated unigenes were highly enriched in aminoacyl-tRNA biosynthesis and porphyrin and chlorophyll metabolism. GO enrichment showed there to be 1 and 7 significantly enriched GO terms in the up- and down-regulated unigenes, respectively (see Additional file 4).

Figure 7
figure 7

Number of DEGs in the different temperature stress groups.

Full size image

In higher plants, increased synthesis of polyunsaturated fatty acids (predominantly trienoic) occurs in response to chilling stress. DEG analysis revealed 12 fatty acid desaturase (FAD) unigenes among the 762 up-regulated unigenes, but no FAD unigenes were found among the down-regulated unigenes. Among the 1334 DEGs, unigenes encoding glutathione S-transferase, oxidoreductase, and thioredoxin were found. These promote chilling tolerance by maintaining cell redox homeostasis. Two categories of DEGs functioning as temperature sensors were analyzed, both histone deacetylase complex and transcription factor. Results showed an up-regulated unigene encoding the histone acetyltransferase complex and a down-regulated unigene encoding the histone deacetylase complex. 13 unigenes encoding transcription factor that included basic leucine zipper (bZIP), upstream activation factor (UAF), and myeloblastosis oncogene were differentially regulated.

DEGs in response to freezing treatment

A total of 592 differentially expressed unigenes were found in FS, relative to NT; 383 unigenes were up-regulated and 209 unigenes were down-regulated (Figure 7). The top 100 down- and up-regulated genes are shown in Additional files 5 and 6. There were fewer DEGs in FS than in CS and HT. It can be inferred that, under freezing conditions, many physiological functions are suspended and the expression levels of the unigenes were maintained at a normal level. Among the up-regulated unigenes, there were 7 FAD unigenes and one unigene encoding farnesyl pyrophosphate synthase. There were also many unigenes encoding glutathione S-transferase, oxidoreductase, thioredoxin, and the histone acetyltransferase complex, which promote chilling tolerance, like the unigenes found in CS. As noted in the global analysis of the gene expression component, FS has a pattern of expression similar to that of CS. GO enrichment analysis showed there to be 3 significantly enriched GO terms among the up-regulated unigenes and 26 significantly enriched GO terms among the down-regulated unigenes (see Additional file 7).

DEGs in response to heat treatment

A total of 2202 differentially expressed unigenes were identified in HT, relative to NT; 1103 unigenes were up-regulated and 1099 unigenes were down-regulated (Figure 7). The top 100 down- and up-regulated genes are shown in Additional files 8 and 9. KEGG analysis showed that the up-regulated DEGs were highly enriched in protein processing in endoplasmic reticulum and the down-regulated unigenes were enriched in ribosome biogenesis in eukaryotes, one carbon pool by folate, alanine, aspartate, and glutamate metabolism, and RNA transport. GO enrichment analysis showed that 116 GO terms were significantly enriched in the down-regulated unigenes, but among the up-regulated unigenes, no GO terms were enriched (see Additional file 10).

In this study, three categories of DEGs functioning as temperature sensors were found, including heat shock proteins (HSPs), H2A, and transcription factors. Among the up-regulated unigene pools, there were 18 well-characterized HSPs (i.e., 3 Hsp20, 2 Hsp90, 4 Hsp60 and 9 Hsp70). In addition, 4 up-regulated unigenes encoding H2A, which is believed to be a temperature sensor, were found. An additional 22 transcription factor (TF) unigenes were differentially regulated in HT and NT, including bZIP, C-repeat binding factor (CBF), TFIIB, and UAF. A multiprotein bridging factor 1 (MBF1) was up-regulated. Among the top 20 down-regulated unigenes, 15 were chloroplast genes. Four unigenes encoding metacaspase were found to be up-regulated under high temperature stress which was consider to play an important role in programmed cell death.

Validation of RNA-Seq-based gene expression

To validate the expression profiles obtained by RNA-Seq, RT-PCR was performed on seven genes selected at random with high or low expression levels. Expression comparisons were performed between HT and NT, CS and NT, and FS and NT by RT-PCR. For all of the genes, the trend in RT-PCR expression was in agreement with the RNA-Seq data except for comp8549_c0 (Figure 8).

Figure 8
figure 8

RT-PCR analysis of 7 randomly selected unigenes. ((A) Gene expression data for DGE analysis. The fold changes of the genes were calculated as the log2 vaule of HT/NT, CS/NT and FS/NT and shown on the y-axis. (B) The qRT-PCR analysis of gene expression data. Expression ratios of these genes in HT, CS and FS were compared to NT, respectively).

Full size image

Discussion

Role of unsaturated fatty acids under low temperature stress conditions

In general, plants regulate the proportion of polyunsaturated fatty acids to change membrane fluidity in response to temperature stress. In both yeast and cyanobacteria, it has been hypothesized that temperature-mediated alterations of membrane fluidity may themselves be the primary temperature sensing event. It has been speculated that the same might be true in higher plants [13]. In this study, a total of 30 unigenes encoding fatty acid desaturases were found, in which 12 unigenes belonging to 7 diverse types were differentially expressed (Table 3). In higher plants, it has been reported that cold tolerance is closely correlated to the level of unsaturated fatty acids in phosphatidylglycerol from the chloroplast membrane, especially in the sn-1 position [14]. Plants synthesize 18:2 n-6 and 18:3 n-3 from saturated fatty acid 18:0; these are catalyzed by soluble delta-9 acyl-ACP, membrane-bound delta-12 fatty acid and ω3 desaturases, respectively [15]. Delta-9 fatty acid desaturase catalyzes the formation of the initial double bond between the 9th and 10th carbons of palmitoyl (16:0) and stearoyl (18:0) coenzyme A (CoA) substrates to produce 16:1 and 18:1 fatty acids, respectively [16]. Delta-15 fatty acid desaturase is located in plastids and catalyzes the introduction of a double bond into 18:2 n-6 esterified in the sn-1 position of glycerolipids [17]. In the CS and FS treatments, all 10 FAD were up-regulated, promoting the synthesis of fatty acid desaturases, thereby increasing membrane fluidity, which results in acclimation to chilling stress in the intertidal zone. In the high-temperature treatment group, 2 FAD were down-regulated. In the KEGG database, 12 unigenes were successfully annotated to the pathway of unsaturated fatty acid biosynthesis (Figure 9). In the high-temperature treatment group, in addition to the 2 down-regulated FADs, the unigene comp704_c0 encodding 3-hydroxy acyl-CoA dehydratase was up-regulated which takes part in fatty acid elongation.

Table 3 Differentially expressed unigenes encoding fatty acid desaturase in CS and FS relative to NT
Full size table
Figure 9
figure 9

Differential gene expression in the biosynthesis pathway of unsaturated fatty acids in P. yezoensis. (Green represents up-regulated unigenes; red represents down-regulated unigenes. ★ represents unigenes differentially expressed between NT and CS; ◆ represents unigenes differentially expressed between NT and FS. ● represents unigenes differentially expressed between NT and HT).

Full size image

HSP in response to temperature stresses

It is well established that HSP contributes greatly to cell protection through the folding and translocation of nascent proteins, the refolding of denatured proteins, and the disassembly of already formed protein aggregates under both stress and non-stress conditions [18-21]. Based on subunit sizes, HSP can be divided into 5 groups, HSP100, HSP90, HSP70, HSP60, and small HSP. In this study, a total of 56 unigenes were found to encode HSP. The expression level of the most of the unigenes was not significantly different, while their expression levels were high (see Additional file 11). In addition, HSP70, HSP90, HSP60, and HSP20 were found to significantly differently expressed as a result of the temperature treatment (Table 4). Recently, HSP70 has attracted widespread attention; this protein plays a vital role in the transport of nascent proteins across membranes into organelles, the folding of newly translated proteins, the repair of misfolded proteins, and targeting damaged proteins for degradation [20]. When plants suffer from stress, HSP70s are up-regulated, participate in the refolding of denatured proteins, maintain cell homeostasis, and protect organisms from damage [22]. Nine unigenes encoded HSP70; six were up-regulated and three were down-regulated (Table 4). HSP60 has been shown to function as a chaperonin in the assembly of mitochondrial enzyme complexes composed of proteins encoded by nuclear genes and imported from the cytosol [23]. All 4 HSP60 genes were down-regulated in both in HT and FS. In general, the HSP unigenes did not exhibit regular changes. It is here assumed that the HSP unigenes maintain high expression levels prior to frequent habitat changes to allow the organism to acclimate to the intertidal environment.

Table 4 Variation of HSP unigene expression under different treatments compared to NT
Full size table

Histones in response to temperature stresses

Histones are basic proteins that package DNA into nucleosomes, and histone gene expression is closely correlated with the cell cycle and cell proliferation [24]. It has been reported that histone genes are repressed by cold stress [25]. Some researchers believe that, at cooler temperatures, H2A.Z occupancy represses gene expression by creating a physical block to transcription or by preventing the binding of complexes that activate transcription [26]. CS treatment showed less histone expression than NT, but FS treatment showed no differently expressed unigenes that encoded histone. A unigene encoding histone acetyltransferase complex was found to be up-regulated and a unigene encoding histone deacetylase complex that was found to be down-regulated (Table 5). Histone acetylation and deacetylation processes occur at specific lysine residues within the N-terminal. One theory suggests that ethanoyl is electronegative and can neutralize the electropositivity of histones. Consequently, ethanoyl can reduce the affinity of histones and DNA, which is electronegative. In this way, the acetylation of histone can lose nucleosomes in order to activate transcription. In the heat treatment, 4 histone unigenes including H2A, H2B, H3, and H4 were found to be up-regulated.

Table 5 Variation in the expression of unigenes related to histone and acetylation in HT, CS, and FS relative to NT
Full size table

Multi-protein bridging factor and transcription factors

Multiprotein bridging factor 1 (MBF1) is a highly conserved transcriptional coactivator involved in the regulation of diverse processes, such as endothelial cell differentiation, histidine metabolism, hormone-regulated lipid metabolism, and central nervous system development [27-29]. MBF1 functions up-stream of salicylic acid, trehalose, and ethylene during high-temperature stress by causing the accumulation of numerous transcription factors and transcription products of signal transduction genes [30]. In the current study, the expression of MBF was higher in HT but remained unchanged under cold stress (both FS and CS), which suggests that cold stress had little effect on the expression patterns of PyMBF1. A similar result was reported in a previous study which found that PyMBF1 transcripts are up-regulated in P. yezoensis cells during exposure to oxidative and heat stresses, and that the heat activation of PyMBF1 requires membrane fluidization [31]. In addition, 5 Zn-finger protein unigenes were up-regulated. These are believed to be controlled by MBF1 at the transcriptional level during high-temperature stress. In Arabidopsis thaliana, MBF1c acts as a transcriptional regulator. It binds DNA and controls the expression of 36 different transcripts during heat stress, including the important transcriptional regulator DRE-binding protein 2A (DREB2A), two heat shock transcription factors (HSFs), and several zinc finger proteins [32]. The current work suggests that PyMBF1 plays a role in the oxidative and heat stress response pathways in P. yezoensis.

Transcription factors, as key regulators of gene expression can be found in any organism. These factors act downstream of signaling cascades related to biological and environmental stimuli. In this study, transcription factors played an important role in response to temperature stresses. Five up-regulated bZIP unigenes were detected in HT and two down-regulated bZIP unigenes were detected in CS. The expression of MYB was higher in CS. Two TF families, bZIP and MYB, are involved in ABA signaling and its gene activation, which plays a vital role in plant stress responses [33]. MYB proteins are key to regulatory networks controlling development, metabolism, and responses to biotic and abiotic stresses [34]. It has been suggested that MYB15 is part of a complex network of transcription factors that control the expression of CBFs and other genes in response to cold stress [35].

Programmed cell death in response to high-temperature stress

In unicellular organisms, programmed cell death (PCD) is a protective action that keeps their cellular best adapted to survival under stress conditions. It may also serve as a possible mechanism for managing cell cycles and differentiation [36]. Caspases are cysteine proteases that are important regulators of PCD in animals. Plant genomes do not contain structural homologs of caspases. Instead, they encode several related proteins, called metacaspases. These are also present in other organisms such as fungi, parasitic protozoa, and some bacteria [37]. Metacaspases are suggested to be the ancestors of metazoan caspases, and plant metacaspases have previously been shown to be genuine cysteine proteases that automatically process their substrates in a manner similar to that of caspases [38]. Many studies have been performed to explore the role of metacaspases in PCD. In Arabidopsis, metacaspases AtMC1 and AtMC2 were found to mediate PCD [39]. In yeast cells, metacaspases are involved in PCD when the yeast cells suffer different extracellular stresses such as viral infections and heat shock [40]. Up until now, there have been only a few studies on the relationship between metacaspases and PCD in algae in response to temperature, For example, shifting Symbiodinium microadriaticum from 27 to 32°C resulted in an increase in mortality, an increase in caspase 3-like activity, and an increase in nitric oxide (NO) production [41]. In the current work, four unigenes encoding metacaspase were found to be up-regulated under high-temperature stress, but none were found in other treatment groups, from which it can be inferred that the overexpression of the genes encoding metacaspases played an important role in P. yezoensis acclimation to heat stress. It has been reported that PCD is a central theme during plant reproductive development, and precise control of PCD execution, or its prevention, are intimately linked with successful plant reproduction [42]. Sexual reproduction could be accelerated in an elevated temperature and light climates; in Sargassum horneri it occurred at least 3 months earlier than in the wild [43]. This inspired the present investigation of the expression of genes encoding metacaspases activated under high-temperature stress conditions to regulate the progress of development, accelerating the formation of generative cell in P. yezoensis.

Chloroplast genes in response to high-temperature stress

Under high-temperature treatment conditions, the chloroplast genes encoding 3-oxoacyl-acyl-carrier-protein synthase 3, photosystem I subunit XI, photosystem II 47 kDa protein, ribosomal protein L19, and ribulose-1,5-bisphosphate carboxylase were significantly down-regulated relative to the 8°C group, indicating that less photosynthesis was taking place. Photosystem II (PSII) is one of the most thermosensitive components of photosynthesis and PSI activity is much more heat stable than PSII. The two factors that make PSII electron transport most susceptible to heat stress are (i) increases in the fluidity of thylakoid membranes at high temperature which causes PSII light harvesting complexes to become dislodged from the thylakoid membrane and (ii) dependence of PSII integrity on electron dynamics [44]. Research on P. tenera may provide support for the current finding that the gross photosynthesis for the gametophytes of P. tenera were determined over a temperature range of 8–34°C. This showed that the gross photosynthetic rates of 46.3 μmol O2 mgchl-a −1 min−1 was highest at 9.3 (95% Bayesian credible interval (BCI): 2.3–14.5°C) [45]. It was here assumed that P. yezoensis resists high-temperature stress by inhibiting photosynthesis, as has been found in P. haitanensis [11].

Role of hormones in response to temperature stress in P. yezoensis

A growing body of evidence suggests that plants use endogenous hormones to couple their growth rate to temperature rather than subject themselves to temperature-induced inhibition of growth [46]. Although abscisic acid (ABA) contains 15 carbon atoms, it is formed in two different ways, directly from the C15 sesquiterpene precursor farnesyl diphosphate (FDP), which is the direct pathway, and formed by cleavage of C40 carotenoids originating from the MEP pathway, which is considered the indirect pathway [47]. In Pyropia, abscisic acid has been detected using liquid chromatography–tandem mass spectrometry [48]. Here, transcriptome data were searched for unigenes related to the biosynthesis of abscisic acid to determine which pathway P. yezoensis took to synthesize abscisic acid. There was one unigene encoding farnesyl pyrophosphate synthase, which acts in the direct pathway, but no unigenes related to the indirect pathway were found, not even 9-cis-epoxycarotenoid dioxygenase (NCED), which cleaves carotenoid precursors to produce xanthoxin, which can subsequently be converted into ABA [49]. For this reason, it was inferred that the direct pathway was active in P. yezoensis, as in the green algae Dunaliella [50] and in fungi, Cercospora cruenta [51]. This is not the preferred pathway of category higher plants, however [52]. In addition, the unigene encoding farnesyl pyrophosphate synthase was up-regulated in CS. From this, it was inferred that the concentration of abscisic acid increased to acclimate to temperature in FS. In many studies, abscisic acid biosynthesis was found to be an essential requirement for cold acclimatization and acquired thermotolerance [53]. The overexpression of ABSCISIC ACID INSENSITIVE 3 (ABI3) confers increased freezing tolerance in Arabidopsis [54].

Conclusions

This study has demonstrated the usefulness of the digital gene expression approach to identify the differentially expressed genes of P. yezoensis under different temperature treatments. First, in response to high temperature stress, P. yezoensis showed a constitutive gene expression profile different from that of other treatments, but the two low temperature stresses (i.e., chilling and freezing stress) were similar. Under heat stress, unigenes related to replication and repair of DNA and protein processing in endoplasmic reticulum showed an up-regulated trend. More unigenes related to carbohydrate metabolism and energy metabolism were up-regulated under low-temperature stress conditions than at normal temperatures. The mechanisms involved in temperature response were demonstrated from the following perspectives: membrane fluidity, histones, PCD, photosynthesis, and transcription factors. A large number of genes of unknown function showed significant variations in expression under temperature stress, providing a new and valuable clue for future research into mechanisms by which P. yezoensis tolerates the unique environment of the intertidal zone.

Methods

Plant materials

To eliminate the interference caused by genotypic differences, a lab-cultured pure line PY4-7 of P. yezoensis was used in these experiments. The gametophytes of this pure line were generated from the clonal cultivation of a single somatic cell that originated from a farmed thallus, which was cultured in bubbling natural seawater with Provasoli’s enrichment solution medium (PES) under 50 μmol photons m−2 s−1 at 8 ± 1°C and a 12:12 light:dark (L:D) photoperiod. Sporophytes were generated from the self-fertility of the pure line and were cultured in bubbling natural seawater with PES under 20 μmol photons m−2 s−1 at 20 ± 1°C and a 12:12 light:dark (L:D) photoperiod.

Experimental design and sampling

The material included the samples from different developmental phases and treatments (Table 6). After harvesting and weighing, the samples were immediately frozen in liquid nitrogen [55].

Table 6 Samples used to construct the global transcriptome
Full size table

For DGE analysis, four temperature treatments were set up: normal temperature (NT, 8°C), high temperature (HT, 24°C), chilling stress (CS, 0°C), and freezing stress (FS, −8°C). Two biological replicates were used in each treatment. After 6 h of treatment, all samples were collected and stored in liquid nitrogen.

RNA isolation

Total RNA was extracted from the thallus using the Plant RNA Kit (Omega, U.S.) in accordance with the manufacturer’s instructions. RNA degradation and contamination were monitored on 1% agarose gels. The purity was evaluated using a NanoPhotometer spectrophotometer (IMPLEN, CA, U.S.). RNA concentration was measured using a Qubit RNA Assay Kit and Qubit 2.0 Fluorometer (Life Technologies, CA, U.S.). RNA integrity was assessed using an RNA Nano 6000 Assay Kit and the Bioanalyzer 2100 system (Agilent Technologies, CA, U.S.).

Library preparation for transcriptome and DGE sequencing

A total of 6 μg RNA per sample was used as input material for the RNA sample preparation. All nine RNA samples (including a pooled RNA sample for transcriptome sequencing and eight RNA samples for DGE sequencing) had RIN values above 7.0. Sequencing libraries were generated using an Illumina TruSeq RNA Sample Preparation Kit (Illumina, San Diego, CA, U.S.) in accordance with the manufacturer’s recommendations. Eight index codes were added to attribute sequences to each sample. Briefly, mRNA was purified from total RNA using poly-T oligo-attached magnetic beads. Fragmentation was carried out using divalent cations at elevated temperatures in Illumina proprietary fragmentation buffer. First-strand cDNA was synthesized using random oligonucleotides and SuperScript II (Invitrogen). 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 enzymes were removed. After adenylation of 3′ ends of DNA fragments, Illumina PE adapter oligonucleotides were ligated to prepare for hybridization. To select cDNA fragments of 200 bp in length, the library fragments were purified using an AMPure XP system (Beckman Coulter, Beverly, U.S.). 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 an Agilent high-sensitivity DNA assay and Agilent Bioanalyzer 2100 system. Clustering of the index-coded samples was performed on a cBot Cluster Generation System using a TruSeq PE Cluster Kit v3-cBot-HS (Illumia) according to the manufacturer’s instructions. After cluster generation, the libraries were sequenced on an Illumina Hiseq 2000 platform. The sample for transcriptome was sequenced with 100 bp pair end reads. The other eight libraries for DGE were sequenced with 100 bp single end reads for DGE analysis.

Quality control

Raw data (raw reads) of fastq format were first processed through in-house Perl scripts. In this step, clean data (clean reads) were obtained by removing reads containing adapter, reads in which more than 10% of the bases were unknown,and low quality reads (where more than 50% of bases in a read had a quality value Q ≤5) from the raw data. The Q20, Q30, and GC content and sequence duplication level of the clean data were calculated. All the downstream analyses were based on clean, high-quality data.

De novo assembly and annotation

De novo transcriptome assembly was performed using the short-reads assembly program, Trinity (v2012-10-05) with min_kmer_cov set to 2 and all other parameters set default [56]. Using pair-end reads, we detected contigs from the same transcript as well as the distances between these contigs. Next, we used Trinity to connect the contigs and obtain sequences that could not be extended on either end, known as unigenes. The optimal assembly results were chosen according to the assembly evaluation. Then the clustering analysis was performed to achieve a unigene database which comprised the potential alternative splicing transcripts.

After clustering, the unigenes were divided into two classes: clusters and singletons. Finally, BLASTx alignment (E-value < 0.00001) was performed between unigenes and the protein databases, including Nr, Swiss-Prot, KEGG, and KOG. The best alignment results were used to decide the sequence direction of unigenes.The unigenes were compared against those in the NCBI Nr and Nt database and 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 unigenes sequences were searched against PFAM database to predict functional domain and protein family using Hmmerscan (HMMER 3.0 package) with an inclusion E value of 0.01.

To annotate the unigenes with GO terms describing biological processes, molecular functions and cellular components, the Nr and PFAM annotation results were imported into Blast2GO program, a software package that retrieves GO terms, allowing gene functions to be determined and compared [57].

In order to gain an overview of gene pathways networks, KEGG pathways were assigned to the unigenes using the online KEGG Automatic Annotation Server, (http://www.genome.jp/kegg/kaas/). The bi-directional best hit (BBH) method was used for KEGG Orthology (KO) assignment [58]. The output of KEGG analysis here includes KO assignments and KEGG pathways populated with the KO assignments.

DGE data analysis

Reads mapping to the reference genome

Transcriptome data were selected as the reference. Reads were aligned to the reference genome using RSEM (v1.2.0) software package [59].

Quantification of gene expression level and differential expression analysis

HTSeq software (www-huber.embl.de/users/anders/HTSeq/) was used to count the reads mapped to each orthologous genomic region. For all comparisons, read counts were normalized to the aligned RPKM [60] to obtain the relative levels of expression. RPKM > 0.3 was defined as the threshold of significant gene expression. Differential expression analysis was performed using the R packages of DESeq for comparisons among samples with two biological replicates [61]. The correlation of the detected number of counts between parallel libraries was assessed statistically by calculating the Pearson correlation. P-values (adjusted for false discovery rate) were generated for each gene in pair-wise comparisons among samples. The P-values were adjusted using the method described by Benjamini and Hochberg [62]. A corrected P-value of 0.05 was set as the threshold for significant differential expression.

Variance-stabilized data obtained using DESeq was used to generate the heatmaps of differentially expressed genes. Clustering analysis was performed using the command hclust in R, and the heatmap was drown using the R packages of ggplot2 and pheatmap.

Gene ontology (GO) and KEGG enrichment analysis of differentially expressed genes

GO enrichment analysis of differentially expressed genes was implemented using the GOseq R package, in which gene length bias was corrected. GO terms with corrected P-values below 0.05 were considered significantly enriched by differentially expressed genes.

KEGG is a database resource meant to facilitate understanding of the high-level functions and utilities of biological systems, such as the cells, organisms, and ecosystems, using molecular-level information, especially large-scale molecular datasets generated by genome sequencing and other high-throughput experimental technologies (http://www.genome.jp/kegg/). KOBAS software was here used to test the statistical enrichment of differentially expressed genes in KEGG pathways. Pathways with corrected P-values below than 0.05 were considered significantly enriched by differentially expressed genes.

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. Quantitative real-time PCR was conducted using SYBR Green dye (LightCycle® 480 SYBR Green I Master). The selected genes were verified using the LightCycle® 480 Real-Time PCR System with the following cycling conditions: 95°C for 5min, followed by 45 cycles of 95°C for 10 s, 60°C for 10 s and 72°C for 20s. Specificity of the qPCR product was analyzed by melting curve analysis. The sequences of the primers used are given in Additional file 12. β-actin (ACT3) and translation initiation factor 4A (eIF4a) served as internal controls. The 2-△△Ct method was used to calculate relative gene expression values [63].

Availability of supporting data

The data sets supporting the results of this article are available in the Sequence Read Archive (SRA), accessible through NCBI BioProject ID PRJNA235353 for the transcriptome data (https://www.ncbi.nlm.nih.gov/sra/?term=PRJNA235353) and PRJNA236059 for the DGE data (https://www.ncbi.nlm.nih.gov/sra/?term=PRJNA236059).

References

  1. Iba K. Acclimative response to temperature stress in higher plants: approaches of gene engineering for temperature tolerance. Annu Rev Plant Biol. 2002;53(1):225–45.

    Article  CAS  PubMed  Google Scholar 

  2. Larcher W. Physiological plant ecology: ecophysiology and stress physiology of functional groups: Springer Science and Business Media 2003.

  3. Wei W, Qi X, Wang L, Zhang Y, Hua W, Li D, et al. Characterization of the sesame (Sesamum indicum L.) global transcriptome using Illumina paired-end sequencing and development of EST-SSR markers. BMC Genomics. 2011;12(1):451.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  4. Feng C, Chen M, Xu C-j, Bai L, Yin X-r, Li X, et al. Transcriptomic analysis of Chinese bayberry (Myrica rubra) fruit development and ripening using RNA-Seq. BMC Genomics. 2012;13(1):19.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  5. Chang L, Xia Y-p, Chen J-j, Xiao Y-m. Application of digital gene expression tag profiling on differential gene expression of two developmental stages in bulbs of Lycoris sprengeri. Indian J Horticulture. 2011;68(4):522–8.

    Google Scholar 

  6. Zenoni S, Ferrarini A, Giacomelli E, Xumerle L, Fasoli M, Malerba G, et al. Characterization of transcriptional complexity during berry development in Vitis vinifera using RNA-Seq. Plant Physiol. 2010;152(4):1787–95.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  7. 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(1):45–55.

    Article  PubMed  Google Scholar 

  8. Dawes C, Orduna-Rojas J, Robledo D. Response of the tropical red seaweed Gracilaria cornea to temperature, salinity and irradiance. J Appl Phycol. 1998;10(5):419–25.

    Article  Google Scholar 

  9. Kayama M, Iijima N, Kuwahara M, Sado T, Araki S, Sakurai T. Effect of water temperature on the fatty acid composition of Porphyra. Bulletin of the Japanese Society of Scientific Fisheries 1985, 51(4):687-691.

  10. Choi S, Hwang MS, Im S, Kim N, Jeong W-J, Park E-J, Gong Y-G, Choi D-W. Transcriptome sequencing and comparative analysis of the gametophyte thalli of Pyropia tenera under normal and high temperature conditions. J Appl Phycol 2013,25(4):1237-1246.

  11. Xu Y, Chen C, Ji D, Hang N, Xie C. Proteomic profile analysis of Pyropia haitanensis in response to high-temperature stress. J Appl Phycol 2014,26(1):607-618.

  12. Storey JD. The positive false discovery rate: A Bayesian interpretation and the q-value. Ann Stat 2003,31(6):2013-2035.

  13. Murata N, Los DA. Membrane fluidity and temperature perception. Plant Physiol. 1997;115(3):875.

    CAS  PubMed Central  PubMed  Google Scholar 

  14. Nishida I, Murata N. Chilling sensitivity in plants and cyanobacteria: the crucial contribution of membrane lipids. Annu Rev Plant Biol. 1996;47(1):541–68.

    Article  CAS  Google Scholar 

  15. Tocher D, Leaver M, Hodgson P. Recent advances in the biochemistry and molecular biology of fatty acyl desaturases. Prog Lipid Res. 1998;37(2–3):73–117.

    Article  CAS  PubMed  Google Scholar 

  16. Bossie MA, Martin CE. Nutritional regulation of yeast delta-9 fatty acid desaturase activity. J Bacteriol. 1989;171(12):6409–13.

    CAS  PubMed Central  PubMed  Google Scholar 

  17. Sakamoto T, Los DA, Higashi S, Wada H, Nishida I, Ohmori M, et al. Cloning of ω3 desaturase from cyanobacteria and its use in altering the degree of membrane-lipid unsaturation. Plant Mol Biol. 1994;26(1):249–63.

    Article  CAS  PubMed  Google Scholar 

  18. Bray EA, Bailey-Serres J, Weretilnyk E. Responses to abiotic stresses. Biochemistry and molecular biology of plants 2000, 1158: e1203.

  19. Feder ME, Hofmann GE. Heat-shock proteins, molecular chaperones, and the stress response: evolutionary and ecological physiology. Annu Rev Physiol. 1999;61(1):243–82.

    Article  CAS  PubMed  Google Scholar 

  20. Nelson RJ, Ziegelhoffer T, Nicolet C, Werner-Washburne M, Craig EA. The translation machinery and 70 kd heat shock protein cooperate in protein synthesis. Cell. 1992;71(1):97–105.

    Article  CAS  PubMed  Google Scholar 

  21. Schroda M, Vallon O. Chaperones and proteases. Chlamydomonas Sourcebook, Ed. 2008;2:671–729.

    Google Scholar 

  22. Mayer M, Bukau B. Hsp70 chaperones: cellular functions and molecular mechanism. Cell Mol Life Sci. 2005;62(6):670–84.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  23. Cheng MY, Hartl F-U, Martin J, Pollock RA, Kalousek F, Neupert W, et al. Mitochondrial heat-shock protein hsp60 is essential for assembly of proteins imported into yeast mitochondria. Nature. 1989;6208:585–674.

    Google Scholar 

  24. Sanchez MdlP, Caro E, Desvoyes B, Ramirez-Parra E, Gutierrez C. Chromatin dynamics during the plant cell cycle. In: Semin Cell Dev Biol 2008,19(6): 537-546.

  25. Steward N, Kusano T, Sano H. Expression of ZmMET1, a gene encoding a DNA methyltransferase from maize, is associated not only with DNA replication in actively proliferating cells, but also with altered DNA methylation status in cold-stressed quiescent cells. Nucleic Acids Res. 2000;28(17):3250–9.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  26. Franklin KA. Plant chromatin feels the heat. Cell. 2010;140(1):26–8.

    Article  CAS  PubMed  Google Scholar 

  27. Takemaru K-I, Li F-Q, Ueda H, Hirose S. Multiprotein bridging factor 1 (MBF1) is an evolutionarily conserved transcriptional coactivator that connects a regulatory factor and TATA element-binding protein. Proc Natl Acad Sci. 1997;94(14):7251–6.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  28. Brendel C, Gelman L, Auwerx J. Multiprotein bridging factor-1 (MBF-1) is a cofactor for nuclear receptors that regulate lipid metabolism. Mol Endocrinol. 2002;16(6):1367–77.

    Article  CAS  PubMed  Google Scholar 

  29. Liu Q-X, Jindra M, Ueda H, Hiromi Y, Hirose S. Drosophila MBF1 is a co-activator for Tracheae Defective and contributes to the formation of tracheal and nervous systems. Development. 2003;130(4):719–28.

    Article  CAS  PubMed  Google Scholar 

  30. Grover A, Mittal D, Negi M, Lavania D. Generating high temperature tolerant transgenic plants: achievements and challenges. Plant Sci. 2013;205:38–47.

    Article  PubMed  Google Scholar 

  31. Uji T, Sato R, Mizuta H, Saga N. Changes in membrane fluidity and phospholipase D activity are required for heat activation of PyMBF1 in Pyropia yezoensis (Rhodophyta). J Appl Phycol. 2013;25(6):1887–93.

    Article  CAS  Google Scholar 

  32. Suzuki N, Sejima H, Tam R, Schlauch K, Mittler R. Identification of the MBF1 heat‐response regulon of Arabidopsis thaliana. Plant J. 2011;66(5):844–51.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  33. Wang W, Vinocur B, Altman A. Plant responses to drought, salinity and extreme temperatures: towards genetic engineering for stress tolerance. Planta. 2003;218(1):1–14.

    Article  CAS  PubMed  Google Scholar 

  34. Dubos C, Stracke R, Grotewold E, Weisshaar B, Martin C, Lepiniec L. MYB transcription factors in < i> Arabidopsis</i> Trends Plant Sci. 2010;15(10):573–81.

    Article  CAS  PubMed  Google Scholar 

  35. Agarwal M, Hao Y, Kapoor A, Dong C-H, Fujii H, Zheng X, et al. A R2R3 type MYB transcription factor is involved in the cold regulation of CBF genes and in acquired freezing tolerance. J Biol Chem. 2006;281(49):37636–45.

    Article  CAS  PubMed  Google Scholar 

  36. Koonin E, Aravind L. Origin and evolution of eukaryotic apoptosis: the bacterial connection. Cell Death Differ. 2002;9(4):394–404.

    Article  CAS  PubMed  Google Scholar 

  37. Saheb E, Trzyna W, Bush J. Caspase-like proteins: Acanthamoeba castellanii metacaspase and Dictyostelium discoideum paracaspase, what are their functions? J Biosciences 2014,39(5):909-916.

  38. Belenghi B, Romero-Puertas MC, Vercammen D, Brackenier A, Inzé D, Delledonne M, et al. Metacaspase activity of Arabidopsis thaliana is regulated by S-nitrosylation of a critical cysteine residue. J Biol Chem. 2007;282(2):1352–8.

    Article  CAS  PubMed  Google Scholar 

  39. 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(6009):1393–7.

    Article  CAS  PubMed  Google Scholar 

  40. Flower TR, Chesnokova LS, Froelich CA, Dixon C, Witt SN. Heat shock prevents alpha-synuclein-induced apoptosis in a yeast model of Parkinson's disease. J Mol Biol. 2005;351(5):1081–100.

    Article  CAS  PubMed  Google Scholar 

  41. Bouchard JN, Yamasaki H. Implication of nitric oxide in the heat-stress-induced cell death of the symbiotic alga Symbiodinium microadriaticum. Mar Biol. 2009;156(11):2209–20.

    Article  CAS  Google Scholar 

  42. Olvera-Carrillo Y, Salanenka Y, Nowack MK. Control of programmed cell death during plant reproductive development. In: Biocommunication of Plants. Springer 2012,14:171-196.

  43. Pang SJ, Liu F, Shan TF, Gao SQ, Zhang ZH. Cultivation of the brown alga Sargassum horneri: sexual reproduction and seedling production in tank culture under reduced solar irradiance in ambient temperature. J Appl Phycol. 2009;21(4):413–22.

    Article  Google Scholar 

  44. Mathur S, Agrawal D, Jajoo A. Photosynthesis: Response to high temperature stress. J Photoch Photobio B 2014,137:116-126.

  45. Watanabe Y, Nishihara GN, Tokunaga S, Terada R. Effect of irradiance and temperature on the photosynthesis of a cultivated red alga, Pyropia tenera (= Porphyra tenera), at the southern limit of distribution in Japan. Phycol Res 2014,62(3):187-196.

  46. Penfield S. Temperature perception and signal transduction in plants. New Phytol. 2008;179(3):615–28.

    Article  CAS  PubMed  Google Scholar 

  47. Zeevaart J, Creelman R. Metabolism and physiology of abscisic acid. Annu Rev Plant Physiol Plant Mol Biol. 1988;39(1):439–73.

    Article  CAS  Google Scholar 

  48. Wang X, Zhao P, Liu X, Chen J, Xu J, Chen H, et al. Quantitative profiling method for phytohormones and betaines in algae by liquid chromatography electrospray ionization tandem mass spectrometry. Biomedical Chromatography. 2013.

  49. Chernys JT, Zeevaart JA. Characterization of the 9-cis-epoxycarotenoid dioxygenase gene family and the regulation of abscisic acid biosynthesis in avocado. Plant Physiol. 2000;124(1):343–54.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  50. Bopp-Buhler M, Wabra P, Hartung W, Gimmler H. Evidence for direct ABA synthesis in Dunaliella (Volvocales). Cryptogamic botany 1991,2(2):192-200.

  51. Yamamoto H, Inomata M, Tsuchiya S, Nakamura M, Uchiyama T, Oritani T. Early biosynthetic pathway to abscisic acid in Cercospora cruenta. Biosci Biotechnol Biochem. 2000;64(10):2075–82.

    Article  CAS  PubMed  Google Scholar 

  52. Hirai N, Yoshida R, Todoroki Y, OHIGASHI H. Biosynthesis of abscisic acid by the non-mevalonate pathway in plants, and by the mevalonate pathway in fungi. Biosci Biotechnol Biochem. 2000;64(7):1448–58.

    Article  CAS  PubMed  Google Scholar 

  53. Gilmour SJ, Thomashow MF. Cold acclimation and cold-regulated gene expression in ABA mutants of Arabidopsis thaliana. Plant Mol Biol. 1991;17(6):1233–40.

    Article  CAS  PubMed  Google Scholar 

  54. Tamminen I, Mäkelä P, Heino P, Palva ET. Ectopic expression of ABI3 gene enhances freezing tolerance in response to abscisic acid and low temperature in Arabidopsis thaliana. Plant J. 2001;25(1):1–8.

    Article  CAS  PubMed  Google Scholar 

  55. 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(20):2119–30.

    Article  CAS  Google Scholar 

  56. 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 Biotechnol. 2011;29(7):644–52.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  57. Conesa A, Götz S. Blast2GO. A comprehensive suite for functional analysis in plant genomics. International journal of plant genomics 2008, 619832.

  58. Moriya Y, Itoh M, Okuda S, Yoshizawa AC, Kanehisa M. KAAS: an automatic genome annotation and pathway reconstruction server. Nucleic Acids Res. 2007;35 suppl 2:W182–5.

    Article  PubMed Central  PubMed  Google Scholar 

  59. Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics. 2011;12(1):323.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  60. Mortazavi A, Williams BA, McCue K, Schaeffer L, Wold B. Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nat Methods. 2008;5(7):621–8.

    Article  CAS  PubMed  Google Scholar 

  61. Anders S, Huber W. Differential expression analysis for sequence count data. Genome Biol. 2010;11(10):R106.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  62. Benjamini Y, Hochberg Y: Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc B 1995,57(1):289-300.

  63. 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.

    Article  CAS  PubMed  Google Scholar 

Download references

Acknowledgements

This work was supported by National Natural Science Foundation of China (Grant No. 31372517), National High Technology Research and Development Program of China (Grant No. 2012AA10A401 and No.2012AA10A406), Independent Innovation Foundation of Shandong Province (Grant No. 2013CXC80202), National Infrastructure of Fishery Germplasm Resources and the Key Project of Chinese Ministry of Education (Grant No. 107070). We would like to thank Zeeshan Niaz for polishing the language of the manuscript.

Author information

Authors and Affiliations

  1. Key Laboratory of Marine Genetics and Breeding (MOE), College of Marine Life Sciences, Ocean University of China, Qingdao, 266003, China

    Peipei Sun, Yunxiang Mao, Min Cao, Fanna Kong & Guiqi Bi

  2. Key Laboratory for Sustainable Utilization of Marine Fisheries Resources, Ministry of Agriculture, Yellow Sea Fisheries Research Institute, Chinese Academy of Fishery Sciences, Qingdao, 266071, China

    Guiyang Li

  3. Institute of Plant Resources, Dalian Nationalities University, Dalian, 116600, China

    Li Wang

Authors
  1. Peipei Sun
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Yunxiang Mao
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Guiyang Li
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Min Cao
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Fanna Kong
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Li Wang
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. Guiqi Bi
    View author publications

    You can also search for this author in PubMed Google Scholar

Corresponding author

Correspondence to Yunxiang Mao.

Additional information

Competing interests

The authors declare that they have no competing interests.

Authors' contributions

YM conceived and designed the project. PS and MC performed the experiments. YM, PS, FK, GL and GB analyzed and interpreted the data. YM, PS, and LW drafted and revised the manuscript. All authors read and approved the final manuscript.

Additional files

Additional file 1: Table S1.

KEGG enrichment of the 8 subclusters).

Additional file 2: Table S2.

The top 100 up-regulated unigenes (annotated) in CS compared with NT.

Additional file 3: Table S3.

The top 100 down-regulated unigenes (annotated) in CS compared with NT.

Additional file 4: Table S4.

GO enrichment analysis of (down and up)-regulated genes in CS compared with NT.

Additional file 5: Table S5.

The top 100 up-regulated unigenes (annotated) in FS compared with NT.

Additional file 6: Table S6.

The top 100 down-regulated unigenes (annotated) in FS compared with NT.

Additional file 7: Table S7.

GO enrichment analysis of (down an up)-regulated genes in FS compared with NT.

Additional file 8: Table S8.

The top 100 up-regulated unigenes (annotated) in HT compared with NT.

Additional file 9: Table S9.

The top 100 down-regulated unigenes (annotated) in HT compared with NT.

Additional file 10: Table S10.

GO enrichment analysis of (down and up)-regulated genes in HT compared with NT.

Additional file 11: Table S11.

Expression abundance of HSP genes in the eight samples.

Additional file 12: Table S12.

Primers used in qRT-PCR for validating differentially expressed genes.

Rights and permissions

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 https://creativecommons.org/licenses/by/4.0/.

The Creative Commons Public Domain Dedication waiver (https://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Sun, P., Mao, Y., Li, G. et al. Comparative transcriptome profiling of Pyropia yezoensis (Ueda) M.S. Hwang & H.G. Choi in response to temperature stresses. BMC Genomics 16, 463 (2015). https://doi.org/10.1186/s12864-015-1586-1

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12864-015-1586-1

Keywords

BMC Genomics

ISSN: 1471-2164

Contact us