Two to three-year old Chinese yew (T. chinensis) seedlings were grown in a green house with the settings at 12 light/12 h dark, and a temperature of 23/15°C. The relative humidity was 70% and photosynthetically active radiation at the leaf level was 250–350 μmol m⁻² s⁻¹. Before stress, stems of the control-grown seedlings were wrapped by cotton and cling film and the seedlings were transferred into the growth chamber. The stress treatment was carried out in a growth chamber, with a three-step incubation of plants at 10°C for the first 48 h, followed by 5°C for the next 48 h, and then 0°C until the samples harvest. After plants had been stressed for 12 h at 0°C, plant tissues were harvested for RNA extraction. Tissues from control seedlings and stressed seedlings were sampled from the longest leaf (leaf 5–8). Each treatment (control and cold stress) was carried out three times, resulting in three biological replicates. Collected samples were frozen in liquid nitrogen immediately and stored at −80°C for further molecular analysis.

Activity of three enzymes, superoxide dismutase (SOD), peroxidase (POD), and catalase (CAT) was determined in leaf tissue samples on the Epoch Microplate Reader (BioTek, USA) using an ELISA enzyme activity Kit (TianGen Inc., Beijing, China), following the manufacturers' instructions. The leaf tissue samples were harvested and enzyme activities were measured after 1, 2, 4, 8, 12, 24, and 48 h of the 0°C treatment. Enzyme activities of the control plants (control) were also measured by directly sampling leaf tissues from seedlings that were grown at normal temperature (23°C). Three biological replicates were analyzed.

Frozen tissue samples were ground in liquid nitrogen into fine powder, which was used for RNA extraction. Total RNA was extracted using the Column Plant RNAout Kit (Tiandz, Inc. Beijing, China), and genomic DNA was digested using Desoxyribonuclease I Amplification Grade (Invitrogen, San Diego, California, USA) following the manufacturer's instructions. RNA concentration and purity was determined Nanodrop Spectrophotometer (ND-1000 Spectrophotometer, Nanodrop products, Wilmington, USA), and Qubit 2.0 (Thermo Fisher Scientific, USA). Poly(A) mRNA was enriched using magnetic beads with oligo(dT). The mRNA library was constructed using the Illumina TruSeq RNA sample prep kit (Illumina, San Diego, CA) following the manufacturers' instructions. Paired-end sequencing was performed using Miseq 150 bp kit (Illumina, San Diego, CA) on Illumina Miseq platform (Illumina, San Diego, CA). Paired-end sequencing was performed using the following Illumina adapters: adapter-R1 5′-AGATCGGAAGAGCACACGTCTGAACTCCAGTCA-3′ and adapter-R2 5′-AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGT-3′. The raw data of Miseq paired-end sequencing was the fastq format.

The barcodes and adapters in reads were removed before assembly, using cutadapt (Martin, 2011), and sequencing quality was checked using Fastqc. A total of 4.16 Gb clean data was obtained for all sample reads. As none genome data was available for Taxus, a de novo RNA-seq assembly was performed toward transcriptome analysis, using Trinity platform (Grabherr et al., 2011; Haas et al., 2013) following Trinity users' guild. The transcriptome was assembled from all six samples (three control samples and three cold stressed samples). The quality of the transcriptome was assessed by TrinityStats scripts and blasting the transcriptome to Swissprot database (Gasteiger et al., 2001). Before building a transcript expression matrix, transcript abundance in each sample was estimated by aligning each sample reads to the transcriptome using Trinity scripts and combined RSEM (Wang et al., 2009) and bowtie packages. To identify transcripts expressed differently, the transcript expression abundance was normalized by TMM (Trimmed Mean of M-values) method (Robinson and Oshlack, 2010; Dillies et al., 2012). The transcript matrix and normalized TMM matrix were used for downstream analysis such as quality check of samples and replicates and different expression (DE) analysis.

To assess the quality of replicates across all samples, the trinity PtR script was used to generate a correlation matrix (heat map) for all sample replicates using the counts-per-million (CPM) and log2 transformed data. Principal Component Analysis was also performed to explore any relationship among the sample replicates. For DE analysis, the Bioconductor R package, edgeR, was used to perform pairwise comparison between control and cold-stressed samples (Robinson et al., 2010). MA and volcano plots were generated to show transcripts that were expressed differently. Subsequently, all isoforms that were at least 2^∧2 fold differently expressed and had p values at most 0.001, 0.01, and 0.05 respectively, were extracted. Heat maps were made to (i) describe the sample correlations using DE isoforms and (ii) show all isoform correlations in each sample.

Open Reading Frames (ORFs) of transcriptome and DE transcripts (both up- and down-regulated isoforms) at p = 0.05 were predicted by Transdecoder (https://github.com/TransDecoder/TransDecoder/releases). ORFs were used for function annotations. KEEG pathway annotation was performed on the KAAS (KEEG Automatic Annotation Serve) website (http://www.genome.jp/tools/kaas/). Gene Ontology (GO) annotation analyses of differently expressed isoforms was performed using BinGo (https://www.psb.ugent.be/cbd/papers/BiNGO/Home.html) by assigning the annotation to the Arabidopsis. Plant transcription factors (TFs) were predicted using PlantTFDB database (http://planttfdb.cbi.pku.edu.cn/).

The reads data has been submitted to NCBI SRA database. The accession numbers of reads data in were SRR5167062, SRR5167061, SRR5167060, SRR5167059, SRR5167058, and SRR5167057 in Bioproject PRJNA360896. Other statistical analyses were carried out using the Minitab 16.0 software (Minitab Inc.). A p-value of smaller than 0.05 was considered as significant.

RT-qPCR was carried out to validate the RNA-seq data. Digestion of genomic DNA and synthesis of cDNA were similar as in section RNA extraction, cDNA synthesis, and sequencing. Quantitative PCR experiments were performed on Bio-Rad iQ5 Optical System (BIO-RAD Laboratories Inc., USA) using SYBR Premix, and relative expression was calculated as described previously (Meng et al., 2016; Meng and Fricke, 2017). Primers used in the present study were designed using PrimeQuest Tools, and primer sequences are given in Table S1. A linear correlation analysis of qPCR and RNA-seq (Figure S1) showed that expression data obtained from different approaches correlated well with each other (Pearson's correlation = 0.712, and p < 0.001).

Cold stress caused a general increase in the activity of all three enzymes in leaf tissue (Figure 1). Hydroperoxidase (CAT) and peroxidase (POD) activity increased rapidly within 1–2 h when plants were exposed to 0°C and remained at a relative high level. Compare to CAT and POD, superoxide dismutase (SOD) activity increased initially slower, but increased rapidly after 2 h and remained at a high level.

The de novo assembled transcriptome had a GC content of 41.83% and included 83,963 transcripts of which 62,654 genes were identified. The median and N50 length of the transcriptome was 598 and 1,497 bp, respectively. Of all transcripts assembled, 55,135 transcripts could be annotated in reviewed Uniprot database at expectation value (e) of 1e-5, which account for 65.7% of the total transcripts. A total of 44,368 ORFs were predicted for the transcriptome by Transdecoder, of which 39,672 sequences (89.4%) could be annotated in reviewed Uniprot database at an expectation value of 1e-5.

Principal components analysis (PCA, Figure S2) of CPM showed that (i) replicates of stressed samples were well-grouped, (ii) while replicates of control samples were not as well-grouped, yet they clearly grouped different from stressed samples. Similar results were obtained by correlation matrix heat map analysis.

Different expression (DE) analyses revealed that 37,934 transcripts expressed differently between control and cold stressed samples at a fold-change (FC) ratio of more than 4 (|log₂FC| > 2), as indicated by red dots in Volcano and MA plots (Figures 2A,B). Of all transcripts that expressed differently, 29,783 transcripts (~78.5%) were up-regulated and 8,151 (~21.5%) were down-regulated under cold stress. A fold-change of more than four (|log₂FC| > 2) was shown for 417 DE transcripts at p < 0.001, 980 transcripts at p < 0.01 and 2,025 transcripts at p < 0.05 (Figure 2C).

Analysis of sample correlation (Figure S3A) using differently expressed isoforms revealed a close correlation between replicates in each treatment, whereas the relationship between treatments were distant from each other. Isoform features correlation heatmap (Figure S3B) showed that DE isoforms could be partitioned into clusters with similar expression patterns.

Annotations were made for 1,437 up-regulated and 588 down-regulated transcripts under cold stress, which had a |log₂FC| > 2 and a significant level of p value < 0.05. Following the prediction by Transdecoder, 1,424 ORFs were predicted from 1,437 up-regulated transcripts and 521 were predicted from 588 down-regulated transcripts.

The plasma membrane is the first major transport barrier of a cell which experiences environmental changes, and a sensitive cold signal receptor system at the plasma membrane would help plants to adjust to stressed conditions. GO enrichment analysis indicated significant changes in the expression of membrane related genes. Genes such as plasma membrane-associated cation-binding protein and plasma membrane proton pump were identified in the up regulated ORFs which were assigned to “response to stress” term (Table S2), indicating significant changes in plasma membrane function in response to cold stress. This caused changes in the following aspects: (a) membrane fluidity/rigidity, (b) membrane permeability of signal receptors e.g., Ca²⁺, and (c) membrane intrinsic protein family (Alonso et al., 1997). In the present study, the CaM4 gene of T. chinensis, was up-regulated in response to cold stress, pointing to an increased intracellular Ca²⁺ concentration. It suggested that membrane permeability of Ca²⁺ receptors was enhanced in response to cold stress to allow more Ca²⁺ entry into the plasma. Ca²⁺ influx was a very important change (initial event) in membrane during cold and (Chinnusamy et al., 2007) suggested the increase of Ca²⁺ influx ability may be induced by accumulation of ROS caused by cold stimulation. The intracellular Ca²⁺ sensors such as calmodulin (CaM) or calcium-dependent protein kinases could bind the Ca²⁺ and then interact with target proteins to activate a series of regulation evens (e.g., MAPK signaling path) for regulating cold-responsive gene expression (DeFalco et al., 2010). Nevertheless, genes such as Calmodulin-like protein and Calcium-dependent protein kinase were up-regulated in response to cold stress (Table S2, ORFs were annotated as AT5G66210, AT1G24620, AT2G38910, and else).

Although the secondary messengers of ROS signaling are unknown at present, it has been suggested that ROS (i) inhibits phosphatases directly, change the membrane structure to enhance Ca²⁺ influx and also act as the receptors of Ca²⁺, (3) the unknown messengers were accepted by redox-sensitive transcription factors, e.g., MAPK signaling pathway (Miller et al., 2008). What's more, (Moller and Sweetlove, 2010) proposed oxidized peptides as the possible secondary ROS messengers that are also involved in MAPK signaling. It is commonly accepted that the plant MAPK signaling pathway plays an important role in plant abiotic stress responses (Pitzschke et al., 2009; Huang et al., 2012; Danquah et al., 2014). The secondary ROS messenger of the MAPK cascade caused indirect activation/up-regulation of cold-responsive TFs to enhance the tolerance capability to cold. One of the most effective defense response to stress is plant signal regulation by TFs, which are components of complex regulatory network in plants.

Eleven TFs in plant-MAPK signaling pathway were up-regulated by cold stresses as predicted by KEGG pathway annotation, particularly TFs that located upstream of the pathway. It has been reported that plants with over-expressed AtNDPK2 had reduced ROS levels (Moon et al., 2003), indicating the role of NDPK2 in regulating ROS scavenging-related genes. This is consistent with the present finding that the increased ROS scavenger activity in T. chinensis could be well-explained by the up-regulated NDPK2 in the MAPK signaling pathway. Two of the TFs in the MAPK pathway, ETR/ERS and CaM4 are secondary signal receptors of ethylene and Ca²⁺, respectively, both substances being messengers induced by stresses. It has been reported for Arabidopsis that the MAPKKK (MAP kinase kinase kinases)-MAPKK (MAP kinase kinases)-MAPK cascade could be activated by ROS, and assist plants in cold acclimation (Teige et al., 2004). Cold stress might have also activated this cascade in T. chinensis to regulate gene expression, as the MEKK1 and MKK9 members of the cascade showed significant up-regulation by cold stress. MYC2 has been reported to be an osmotic responsive gene and to be involved in ABA signaling (Abe et al., 2003) and enhance osmotic stress tolerance of transgenic plants (Gao et al., 2011). The up-regulated MYC2 in T. chinensis might have been due to the osmotic stress component of cold stress and enhanced ABA signaling caused by low temperature. Some studies have supported the idea that PP2C is an important member of the ABA signal transduction pathway and that PP2C functions as regulator of various signal transduction pathways (Rodriguez, 1998). Our observation that PP2C is up-regulated by cold stress supports the ideas that PP2C is activated by accumulated ABA.

Many of the downstream TFs in the MAPK signaling pathway were not up-regulated by cold stress by a factor of four. This may be because many other TFs fulfill “complementary” functions of those TFs. Nevertheless, many cold responsive TF families were found to be up-regulated by cold stress by a factor of more than four. Five TF families including AP2-EREBP (APETALA2/ET-Responsive Element Binding Protein), MYB (Myeloblastosis), NAC (NAM, ATAF1/2, CUC2), bHLH (basic Helix-Loop-Helix), and WRKY (named after the WRKY amino acid motif) have been shown to be involved in the cold stress response in Arabidopsis (Changsong and Diqiu, 2010). All these five families were affected by cold stress in T. chinensis, and therefore, they support a role in the Taxus cold stress response. This applies in particular to NAC, WRKY, and MYB which were the three TF families most-affected in expression by low temperature. Two of the other most numerous TF families up-regulated by cold stress in Taxus were ERF and bZIP. They have also been reported to be involved in the cold stress response in other species (Heidarvand and Maali Amiri, 2010; Mizoi et al., 2012; Calzadilla et al., 2016).

Almost twice as many genes (62,654 vs. 36,493) were assembled in the present study as previously reported (Hao et al., 2011). This may be due to the existence of cold activated genes, as overall, 29,783 transcripts up-regulated under cold stressed condition by a factor of 4 (|log₂FC|> 2). At all significant levels, more isoforms were up-regulated than down-regulated in response to cold stresses, similar to the findings in wheat (Winfield et al., 2010), Lotus japonicas (Calzadilla et al., 2016), and Arabidopisis (Fowler, 2002).

Reactive oxygen species (ROS), for example superoxide (O2-) and H₂O₂, play important roles in regulating plant gene expression, particularly under abiotic stress conditions (Chinnusamy et al., 2007; Heidarvand and Maali Amiri, 2010; Tamás et al., 2010; Aroca et al., 2012; Li et al., 2014). It has also been reported that ROS which accumulate under cold stress activate the membrane Ca²⁺ channels to increase Ca²⁺ flux for gene regulation signaling. The roles which ROS play are distinct (del Río et al., 2006; Huang et al., 2012). On the one hand, they are thought to play a key role in regulating stress-induced signal transduction events. On the other aspect, ROS caused many harmful effects as they inhibit proteins such as aquaporins (Knipfer et al., 2011), induce damages to DNA and lipids, and can ultimately lead to cell death. Therefore, a common problem that plants grown under stressed condition face is how to remove or detoxify ROS effectively. Genes that were assigned in the present study to GO terms “Response to abiotic stimulus,” “Response to stress,” and “Response to endogenous” are listed in Table S2. Genes of ROS scavenging enzymes, such as peroxidase (AT5G05340 and else) and phospholipid hydroperoxide glutathione peroxidase (PHGPx, AT4G11600) were transcriptionally up-regulated in response to cold stress. It has been reported that the up-regulation of these genes inhibits cell death induced by ROS (Chen et al., 2004). Also, enzyme activity of ROS scavengers (including SOD, CAT, and POD) increased under cold stress in the present study. It has been suggested that the scavenging system responds positive to cold stress by removing the toxic ROS that are, initially, necessary for inducing cold-responsive gene expression (Bailey-Serres and Mittler, 2006). In the present study, the increase in CAT activity was the highest of the three ROS-related enzymes studied, might due to that H₂O₂ as a signaling factor might be the most abundant ROS and other ROS should firstly converse to H₂O₂ for further deoxidation. Taken together, the present data suggest that that up-regulation of the genes encoding ROS scavengers and activation of enzymes enhanced the ROS detoxification ability of T. chinensis in leaf tissues.

“Metabolism” was the most abundant category in the KEGG classification and the category most affected by cold stress. The most numerous four pathways in this classification were purine metabolism, thiamine metabolism, starch and sucrose metabolism, and amino and nucleotide sugar metabolism.

Purine metabolism is considered as a housekeeping function, however, a recent study by Watanabe et al. (2014) and Takagi et al. (2016) indicated that the purine metabolite, allantoin is involved in stress response such as ABA metabolism and jasmonic signaling, and therefore enhances plants cold tolerance. Our results showed that purine metabolism in Taxus was up-regulated. An increase in purine nucleotides in response to cold has been reported previously for plants and animals (Desautels and Himmshagen, 1979; Frischknecht and Baumann, 1985).

Up regulation of thiamine metabolism and starch and sucrose metabolism and amino and nucleotide sugar metabolism in response to cold stress may have played a role in ROS scavenging. It has been reported that thiamine (vitamin B1) is induced by the oxidative stress which accompanies cold stress rather than cold stress per se (Rapalakozik, 2011) and that thiamin inhibits lipid peroxidation and free radical oxidation of oleic acid to protect cells from ROS. Some of the carbohydrates are important components of the ROS scavenging system, function together with other components such as phenylpropanoids (Van den Ende and Valluru, 2009). Therefore, changes in starch and sucrose metabolism in response to cold stress may contribute to scavenging of excessive ROS.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.