- Research article
- Open Access
Comparative transcriptome profiling of Pyropia yezoensis (Ueda) M.S. Hwang & H.G. Choi in response to temperature stresses
© Sun et al. 2015
- Received: 28 August 2014
- Accepted: 27 April 2015
- Published: 17 June 2015
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).
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.
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.
- Pyropia yezoensis
- High temperature stress
- Chilling stress
- Freezing stress
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 . Stressful temperature conditions can damage the enzymes needed for photosynthesis, respiration, and protein synthesis and so affect the growth and development of plants . 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 . 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 . 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 . In Vitis vinifera, it was used to describe how plant transcriptomes change during three developmental stages, post setting, veraison, and ripening .
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 . Gracilaria cornea showed photosynthetic and respiratory responses adapted to oceanic salinities and subtropical to tropical water temperatures . 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 . Choi et al. reported that most of the transcripts produced by Pyropia tenera under high-temperature conditions were from the heat shock protein family ; 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 .
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.
Global transcriptome assembly and annotation
Summary of assembly and annotation results for P. yezoensis using Trinity
Number and ratio of annotated unigenes
Unigene with annotation
Annotated in Nr
Annotated in KEGG
Annotated in SwissProt
Annotated in PFAM
Annotated in GO
Annotated in KOG
Annotated in all databases
Total number of unigenes
Sequencing and annotation of eight DGE libraries
Summary of sequencing and mapping results for eight DGE libraries
Error rate (%)
GC content (%)
Global analysis of the unigenes in the four temperature treatment groups
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
Gene Ontology functional analysis of DEGs
DEGs in response to chilling treatment
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
Role of unsaturated fatty acids under low temperature stress conditions
Differentially expressed unigenes encoding fatty acid desaturase in CS and FS relative to NT
delta-5 fatty acid desaturase
delta-5 fatty acid desaturase
delta-12 fatty acid desaturase
delta-15 fatty acid desaturase
delta-6 fatty acid desaturase
delta-9 fatty acid desaturase
fatty acid desaturase
fatty acid desaturase
fatty acid desaturase
delta-4 fatty acid desaturase
delta-8 fatty acid desaturase
delta-4 fatty acid desaturase
HSP in response to temperature stresses
Variation of HSP unigene expression under different treatments compared to NT
HSP20 family protein
HSP20 family protein
HSP20 family protein
Mitochondrial chaperonin, Cpn60/Hsp60p
Mitochondrial chaperonin, Cpn60/Hsp60p
Mitochondrial chaperonin, Cpn60/Hsp60p
Mitochondrial chaperonin, Cpn60/Hsp60p
Heat shock protein 90
Heat shock protein 90
Histones in response to temperature stresses
Variation in the expression of unigenes related to histone and acetylation in HT, CS, and FS relative to NT
Histone acetyltransferase complex
Histone deacetylase complex
Sperm histone P2
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 . 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 . 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 . 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 . MYB proteins are key to regulatory networks controlling development, metabolism, and responses to biotic and abiotic stresses . 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 .
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 . 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 . 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 . Many studies have been performed to explore the role of metacaspases in PCD. In Arabidopsis, metacaspases AtMC1 and AtMC2 were found to mediate PCD . In yeast cells, metacaspases are involved in PCD when the yeast cells suffer different extracellular stresses such as viral infections and heat shock . 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 . 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 . 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 . 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 . 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) . It was here assumed that P. yezoensis resists high-temperature stress by inhibiting photosynthesis, as has been found in P. haitanensis .
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 . 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 . In Pyropia, abscisic acid has been detected using liquid chromatography–tandem mass spectrometry . 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 . For this reason, it was inferred that the direct pathway was active in P. yezoensis, as in the green algae Dunaliella  and in fungi, Cercospora cruenta . This is not the preferred pathway of category higher plants, however . 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 . The overexpression of ABSCISIC ACID INSENSITIVE 3 (ABI3) confers increased freezing tolerance in Arabidopsis .
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.
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
Samples used to construct the global transcriptome
Free-living conchocelis filamentous
Natural seawater with Provasoli’s enrichment solution medium (PES),18°C, 24 μmol photos m−2 s−1, light:dark=12:12
Natural seawater with Provasoli’s enrichment solution medium (PES), 24°C, 24 μmol photos m−2 s−1, light:dark=12:12
Natural seawater, 9°C, 24 μmol photos m−2 s−1, 4 h
Natural seawater added with deionized water of equal volume, 9°C, 24 μmol photos m−2 s−1, 4 h
Natural seawater added with 33 g NaCl per liter, 9°C, 24 μmol photos m−2 s−1, 4 h
Natural seawater, 9°C, 1500 μmol photos m−2 s−1, 4 h
Natural seawater, 9°C, dark, 4 h
--(drought), 4°C, dark, 4 h
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.
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.
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 . 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 .
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 . 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 .
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  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 . 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 . 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 .
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).
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.
- 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.PubMedView ArticleGoogle Scholar
- Larcher W. Physiological plant ecology: ecophysiology and stress physiology of functional groups: Springer Science and Business Media 2003.Google Scholar
- 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.PubMed CentralPubMedView ArticleGoogle Scholar
- 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.PubMed CentralPubMedView ArticleGoogle Scholar
- 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
- 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.PubMed CentralPubMedView ArticleGoogle Scholar
- 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.PubMedView ArticleGoogle Scholar
- 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.View ArticleGoogle Scholar
- 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.Google Scholar
- 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.Google Scholar
- 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.Google Scholar
- Storey JD. The positive false discovery rate: A Bayesian interpretation and the q-value. Ann Stat 2003,31(6):2013-2035.Google Scholar
- Murata N, Los DA. Membrane fluidity and temperature perception. Plant Physiol. 1997;115(3):875.PubMed CentralPubMedGoogle Scholar
- 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.View ArticleGoogle Scholar
- 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.PubMedView ArticleGoogle Scholar
- Bossie MA, Martin CE. Nutritional regulation of yeast delta-9 fatty acid desaturase activity. J Bacteriol. 1989;171(12):6409–13.PubMed CentralPubMedGoogle Scholar
- 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.PubMedView ArticleGoogle Scholar
- Bray EA, Bailey-Serres J, Weretilnyk E. Responses to abiotic stresses. Biochemistry and molecular biology of plants 2000, 1158: e1203.Google Scholar
- 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.PubMedView ArticleGoogle Scholar
- 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.PubMedView ArticleGoogle Scholar
- Schroda M, Vallon O. Chaperones and proteases. Chlamydomonas Sourcebook, Ed. 2008;2:671–729.Google Scholar
- Mayer M, Bukau B. Hsp70 chaperones: cellular functions and molecular mechanism. Cell Mol Life Sci. 2005;62(6):670–84.PubMed CentralPubMedView ArticleGoogle Scholar
- 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
- 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.Google Scholar
- 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.PubMed CentralPubMedView ArticleGoogle Scholar
- Franklin KA. Plant chromatin feels the heat. Cell. 2010;140(1):26–8.PubMedView ArticleGoogle Scholar
- 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.PubMed CentralPubMedView ArticleGoogle Scholar
- 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.PubMedView ArticleGoogle Scholar
- 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.PubMedView ArticleGoogle Scholar
- Grover A, Mittal D, Negi M, Lavania D. Generating high temperature tolerant transgenic plants: achievements and challenges. Plant Sci. 2013;205:38–47.PubMedView ArticleGoogle Scholar
- 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.View ArticleGoogle Scholar
- 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.PubMed CentralPubMedView ArticleGoogle Scholar
- 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.PubMedView ArticleGoogle Scholar
- 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.PubMedView ArticleGoogle Scholar
- 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.PubMedView ArticleGoogle Scholar
- Koonin E, Aravind L. Origin and evolution of eukaryotic apoptosis: the bacterial connection. Cell Death Differ. 2002;9(4):394–404.PubMedView ArticleGoogle Scholar
- 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.Google Scholar
- 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.PubMedView ArticleGoogle Scholar
- 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.PubMedView ArticleGoogle Scholar
- 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.PubMedView ArticleGoogle Scholar
- 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.View ArticleGoogle Scholar
- 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.Google Scholar
- 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.View ArticleGoogle Scholar
- Mathur S, Agrawal D, Jajoo A. Photosynthesis: Response to high temperature stress. J Photoch Photobio B 2014,137:116-126.Google Scholar
- 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.Google Scholar
- Penfield S. Temperature perception and signal transduction in plants. New Phytol. 2008;179(3):615–28.PubMedView ArticleGoogle Scholar
- Zeevaart J, Creelman R. Metabolism and physiology of abscisic acid. Annu Rev Plant Physiol Plant Mol Biol. 1988;39(1):439–73.View ArticleGoogle Scholar
- 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.Google Scholar
- 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.PubMed CentralPubMedView ArticleGoogle Scholar
- Bopp-Buhler M, Wabra P, Hartung W, Gimmler H. Evidence for direct ABA synthesis in Dunaliella (Volvocales). Cryptogamic botany 1991,2(2):192-200.Google Scholar
- 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.PubMedView ArticleGoogle Scholar
- 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.PubMedView ArticleGoogle Scholar
- 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.PubMedView ArticleGoogle Scholar
- 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.PubMedView ArticleGoogle Scholar
- 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.View ArticleGoogle Scholar
- 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.PubMed CentralPubMedView ArticleGoogle Scholar
- Conesa A, Götz S. Blast2GO. A comprehensive suite for functional analysis in plant genomics. International journal of plant genomics 2008, 619832.Google Scholar
- 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.PubMed CentralPubMedView ArticleGoogle Scholar
- Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics. 2011;12(1):323.PubMed CentralPubMedView ArticleGoogle Scholar
- 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.PubMedView ArticleGoogle Scholar
- Anders S, Huber W. Differential expression analysis for sequence count data. Genome Biol. 2010;11(10):R106.PubMed CentralPubMedView ArticleGoogle Scholar
- 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.Google Scholar
- 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.PubMedView ArticleGoogle Scholar
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.