Identification and characterization of NAC genes in response to abiotic stress conditions in Picea wilsonii using transcriptome sequencing

Abstract

Picea wilsonii is widely used for a forestation and landscaping in China. However, the molecular mechanisms underlying abiotic stress tolerance in this species remain elusive. In this study, we examined the transcriptomic changes induced by drought and salt stress in 8-week-old seedlings of P. wilsonii via RNA-Seq analysis; seedlings grown under normal conditions (without abiotic stress) served as a control. Approximately 36.26 Gb of transcriptome data were obtained from the three treatments. A total of 4171 (2579 up-regulated and 1732 down-regulated) and 3438 (1592 up-regulated and 1706 down-regulated) differentially expressed genes were identified under drought and salt treatments, respectively. Both RNA-Seq and qRT-PCR analyses showed that the expression of PwNAC genes, including PwNAC11, 24, 36, 38, 45, 47, 51 and 52, was greatly induced not only by salt stress but also by drought treatment, suggesting that these genes are involved in both salt and drought stress resistance in P. wilsonii. This is the first report of the transcriptome analysis of P. wilsonii. Our data will serve as a valuable genomic resource for the further characterization of the molecular mechanisms underlying abiotic stress tolerance in P. wilsonii.

Keywords:

• Picea wilsonii
• transcriptome
• abiotic stress
• NAC transcription factors

Introduction

Picea wilsonii (Pinaceae) is a dominant coniferous forest tree species in northern China characterised by cold and shade tolerance. Given its broad environmental adaptability, it could serve as a good raw material for the timber industry. However, compared to other congeneric species, the natural distribution range of P. wilsonii is declining because of climate change and anthropogenic interference.

Transcription factors play an important role in plant growth and development and response to external stimuli [1]. A large number of transcription factor families are widely implicated in plant development and stress response, including NAC (short for NAM, ATAF and CUC), Arabidopsis transcription activation factor (ATAF), cup-shaped cotyledon (CUC), v-myb avian myeloblastosis viral oncogene homolog (MYB), basicregion-leucine zipper (bZIP) and dehydration responsive element binding factor. NAC is a large transcription factor family that plays a pivotal role in plant growth, development, stress resistance and hormone signaling [2–9]. To date, 117 NAC genes have been identified in Arabidopsis thaliana, 151 in rice (Oryza sativa), 79 in grapes (Vitis vinifera), 32 in white spruce (Picea glauca), 163 in poplar (Populus spp.), 152 in soybean (Glycine max) and 152 in tobacco (Nicotiana tabacum) [8,10–13]. In Arabidopsis, overexpression of NAC genes increases resistance to abiotic stresses [14,15]. For example, the ATAF-1 overexpressing plants showed greater resistance to drought stress than ataf-1 and ataf-2 mutant plants [3]. The expression of ANAC019, ANAC055 and ANAC072 (RD26) genes is induced by salt, drought and abscisic acid (ABA) treatments, and overexpression of these genes enhances drought resistance in plants [11]. In poplar, overexpression of the NAC13 gene significantly enhances salt tolerance, whereas suppression of NAC13 expression by silencing significantly reduces salt tolerance [16]. In addition to their role in abiotic stress resistance, NAC transcription factors also promote plant growth and development. In wheat (Triticum aestivum), transgenic plants overexpressing the TaNAC6 gene showed not only stronger resistance to drought stress but also higher yield than the wild type [17]. In white birch (Betula papyrifera), overexpression of BpNAC012 enhanced lignin accumulation in transgenic lines and improved resistance to salt and osmotic stresses [18]. However, Arabidopsis thaliana ANAC069 was reported to negatively regulate salt and osmotic stress tolerance [19], suggesting functional divergence of NACs. Despite the numerous studies on NAC transcription factors in model plant species and a variety of crops, the information available on NAC function in woody plants, and especially spruce, is limited.

In this study, we aimed to determine the molecular basis of abiotic stress response in P. wilsonii. Deep transcriptome sequencing of P. wilsonii seedlings was performed using the Illumina technology to investigate the effect of salt and drought stresses on gene expression. Based on the analysis of differentially expressed genes (DEGs), we further examined the expression of 13 candidate NAC genes in P. wilsonii plants treated with low temperature, hydrogen peroxide (H₂O₂), phytohormones and wounding. The results of this study will provide the basis for the improvement of agronomic traits in P. wilsonii.

Materials and methods

Plant materials and treatments

Eight-week-old seedlings of P. wilsonii were used for transcriptome sequencing. Seeds were germinated in the greenhouse, and seedlings were grown at 19 °C and 16 h light/8 h dark photoperiod [20]. To carry out the salt treatment, seedlings were grown on Murashige and Skoog (MS) liquid medium supplemented with 300 mM NaCl for 6 h [20,21]. To conduct the drought treatment, seedlings were placed on dry filter paper with 35% relative humidity at 20 °C for 6 h [22,23]. Seedlings grown on MS liquid medium for 6 h were used as a control.

Additionally, qRT-PCR experiments were performed using 8-week-old P. wilsonii seedlings treated with various abiotic stresses. Briefly, seedlings grown on solid matrices were supplied with water (control) or water containing 300 mM NaCl, 100 μM ABA, 20 mM H₂O₂, 100 μM gibberellin (GA₃), 100 μM indole-3-acetic acid (IAA), 100 μM methyljasmonate (Me-JA), 200 μM salicylic acid (SA) and ethylene (ET) released from 5 mM ethephon (ETH) solution for 0, 3, 6 and 12 h [20,24–26]. To perform the drought treatment, watering was withheld for 0, 7, 14 and 21 days. To determine the effect of low temperature, seedlings were subjected to 4 °C for 0, 3, 6 and 12 h [27]. To conduct the wounding treatment, seedlings were wounded with a sharp needle and harvested at 0, 0.5, 1, 3 and 6 h. All treatments were performed in three biological replications, with 100 seedlings per replicate.

RNA-Seq library construction and transcriptome sequencing

Total RNA was extracted from P. wilsonii seedlings (Aidlab, Beijing) and analyzed using the Nanodrop Spectrophotometer. Then, mRNA was purified from total RNA samples using magnetic oligo (dT) beads. First-strand cDNA was synthesized from RNA fragments using random hexamers, and the second strand of cDNA was synthesized using buffer, dNTPs, DNA polymerase I and RNase H. The double-stranded cDNA was purified using AMPure XP beads, end-repaired and ligated to sequencing adaptors. Subsequently, the cDNA fragments (200 ± 25 bp) were separated on 2% TAE-agarose gels. RNA-Seq libraries were constructed using the Illumina’s mRNA-seq Sample Preparation Kit (Illumina Inc., San Diego, CA, USA) and sequenced using the Illumina HiSeq™ 2500 platform.

Sequence assembly and annotation

Clean sequence reads were assembled using the Trinity software program [28]. These fragments were then extended to generate longer fragments (contigs), which were clustered and further assembled according to the paired-end information. The longest transcripts in the cluster were defined as unigenes. The online BLAST search tool [29] was used to compare unigene sequences with non-redundant (Nr) [30], Swiss-Prot [31], Gene Ontology (GO) [32], Clusters of Orthologous Groups (COG) [33], euKaryotic Orthologous Groups (KOG) [34] and Kyoto Encyclopedia of Genes and Genomes (KEGG) [35]databases. Aminoacid sequences encoded by the unigenes were annotated with the HMMER [36] software using the Pfam database [37].

qRT-PCR

To validate RNA-Seq data, the expression of selected unigenes was examined by qRT-PCR using gene-specific primers designed with the Primer Premier version 3.0 software; the PwEF-1α gene was used as an internal reference (Supplementary material, Table S1) [38]. The qRT-PCR was performed on a Step One Plus Real-time PCR system (ABI, USA) with the Power Up SYBR Green Kit (ABI) using the following conditions: 5 °C for 2 min; 95 °C for 3 min; 40 cycles of 95 °C for 15 s and 60 °C for 1 min. Melting curve analysis was performed from 60 °C to 95 °C by increasing the temperature at a rate of 1.6 °C s⁻¹. All qRT-PCR experiments were performed in three biological replicates, each comprising three technical replications. Relative transcript levels were calculated using the 2^−ΔΔCT method [39,40].

Table 1. Length distribution of the contigs, transcripts, unigenes from the Trinity de novo assembly.

Length range	Contig	Transcript	Unigene
200–300	19,417,751 (99.53%)	163,556 (50.32%)	144,876 (64.70%)
301–500	46,655 (0.24%)	56,916 (17.51%)	39,152 (17.48%)
501–1000	23,872 (0.12%)	42,324 (13.02%)	19,962 (8.91%)
1001–2000	12,845 (0.07%)	36,053 (11.09%)	11,963 (5.34%)
2000+	8062 (0.04%)	26,196 (8.06%)	7972 (3.56%)
Total number	19,509,185	325,045	223,925
Total length	908,041,432	225,061,179	104,765,699
N50 length	47	1406	572
Mean length	46.54	692.40	467.86

Results and discussion

RNA-Seq and transcriptome assembly

After removing adaptor sequences and low-quality reads, 36.26 Gb of clean reads were obtained from three treatments, including salt, drought and control, each comprising 42,328,519, 38,526,704 and 40,026,038 reads, respectively (Supplementary material, Table S2).

Table 2. Functional annotation of Picea wilsonii transcriptome.

Annotated databases	Unigene	≥300 nt	≥1000 nt
COG	23,348	13,045	6145
GO	38,980	21,465	8792
KEGG	15,416	8915	3932
KOG	38,392	21,713	9542
Pfam	42,586	27,008	13,443
Swiss-Prot	31,787	21,144	10,183
nr	68,973	38,917	16,144
All	73,940	40,159	16,274

Shortread sequences from the three treatments were assembled into 19,509,185 contigs greater than 200 bp in length (Table 1). These contigs were further connected into 104,537 transcripts (>500 bp), with N50 value of 1046 bp and a mean length of 692.40 bp. Cluster analysis and assembly of these transcripts revealed 39,897 unigenes (>500 bp), which accounted for 17.81% of the total unigene number (223,925). The N50 value and mean length of unigenes were 572 bp and 467.86 bp, respectively.

Functional annotation and classification of unigenes

In the NCBI Nr database, homologs of P. wilsonii unigenes were identified in a wide range of plant species. Approximately 30.01% of the annotated unigenes, with the best score, showed matches in Picea sitchensis, Meyeozyma guilliermondii and Amborella trichopoda (Figure 1).

Figure 1. Species distribution in the top BLAST hits in the NCBI non-redundant (Nr) database.

Diverse protein databases (COG, GO, KEGG, KOG, Pfam, Swiss-Prot and Nr) were used to annotate the unigenes by BLAST, with an E-value of less than 10⁻⁵ (Table 2). According to GO annotation results, 89,968 (43.27%), 48,655 (23.40%) and 69,278 (33.32%) unigenes were assigned to the biological process, molecular function and cellular component categories, respectively (Figure 2, Supplementary material, Table S3).

Figure 2. Functional annotation of the assembled sequences based on gene ontology (GO)categories.

According to COG annotations, 223,925 unigenes were assigned to functional prediction and classification (Figure 3, Supplemental Table S4). The largest group was the general function prediction (5187 unigenes, 16.83%), followed by translation, ribosomal structure and biogenesis (3363 unigenes, 10.91%), replication, recombination and repair (2446 unigenes, 7.29%), transcription (2194 unigenes, 7.12%), posttranslational modification, protein turnover and chaperones (2049 unigenes, 6.79%), signal transduction mechanisms (1544 unigenes, 5.01%) and carbohydrate transport and metabolism (1805, unigenes 5.86%). Notably, no unigene was assigned to extracellular structures.

Figure 3. Classification of unigenes based on the Clusters of Orthologous Groups (COG) database.

The KEGG pathway database was used to predict the pathway annotations of unigenes. A total of 18,983 unigenes were assigned to 120 biological pathways. The enriched pathways included cellular processes (13 unigenes), gene information processing (143 unigenes), metabolism (437 unigenes), organismal systems (8 unigenes) and environmental information processing (23 unigenes) (Figure 4). Since the characteristics of many genes have not yet been elucidated, the present gene annotation information was incomplete, and only 33.02% of all the assembled unigenes could be annotated. Nevertheless, the data set generated in this study undoubtedly serves as a valuable resource for the investigation of specific processes, functions and pathways in P. wilsonii.

Figure 4. Pathway assignment of the unigenes based on the Kyoto Encyclopedia for Genes and Genomes (KEGG) database.

Expression and phylogenetic analyses of 13 candidate NAC genes

A total of 4171 (2579 up-regulated and 1732 down-regulated) and 3438 (1592 up-regulated and 1706 down-regulated) DEGs were identified under drought and salt treatments, respectively. Commonly regulated by drought and salt stress of these were 1364 genes.

NAC transcription factors play an important role in plant growth, development, hormone signal transduction and abiotic stress response [41,42]. Transcriptome analysis of Arabidopsis and rice showed that 20–25% of the NAC genes are induced by at least one kind of abiotic stress such as salt, drought, cold or phytohormones [9,42]. In this study, to validate RNA-Seq data and perform sequence alignments, we selected 13 candidate NAC genes probably involved in salt or drought stress tolerance. In plants subjected to salt stress, expression levels of PwNAC11, 24, 36, 38, 41, 45, 47, 51 and 52 were highly up-regulated throughout the treatment. The peak expression time points were at 3 h for PwNAC36 and 51 (>7-fold higher than the control), 6 h for PwNAC11 and 38 (∼8-fold higher than the control) and 12 h for PwNAC24, 45, 47, 52 (>3-fold higher than the control). Unlike other NAC genes, which showed a gradual increase in expression, PwNAC41 expression peaked at 3 h during the salt treatment and remained unchanged thereafter (Figure 5). Overall, these results showed that PwNACs potentially function at different stages to regulate salt stress resistance in the plant.

Figure 5. Expression of 13 NAC genes under salt and drought treatments.

In plants treated with drought stress, PwNAC11 and 18 showed the highest expression at day 7, which was >5-fold greater than that in the control. The expression of PwNAC52 and 55 increased rapidly; compared to the control, both these genes showed 20-fold higher expression at 7 days and >66-fold higher expression (peak) at 21 days. On the contrary, the expression of PwNAC34 declined during the drought treatment, reaching the lowest level at 14 days (Figure 5), suggesting that PwNAC34 plays a negative role in plant response to drought stress.

A phylogenetic tree of 13 P. wilsonii NACs and 21 Arabidopsis homologs was constructed using the neighbor-joining method with the MEGA software (Figure 6). PwNAC4 and PwNAC11 were close to ANAC019, ANAC055, ANAC072 (RD26), ANAC002, ANAC032 and ANAC081 (ATAF2). PwNAC52 was similar to ANAC010 and ANAC073. PwNAC41 and PwNAC51 showed high homology to ANAC090 and ANAC036. PwNAC34, PwNAC38 and PwNAC45 were similar to ANAC025, ANAC029 and ANAC047.

Figure 6. Phylogenetic analysis of Picea wilsonii and Arabidopsis thaliana NAC genes.

Expression profiles of 13 candidate PwNAC genes under low temperature and H₂O₂ treatments

In plants subjected to low temperature (4 °C), expression levels of PwNAC4, 24, 36, 41 and 45 were up-regulated, whereas those of PwNAC18, 34, 38, 47 and 52 were repressed; three genes including PwNAC11, 51 and 55 showed almost unchanged expression in response to low temperature (Figure 7). Notably, the expression level of PwNAC4 in cold-treated plants was 3.63-fold greater than that in the control at 12 h (peak expression). In plants treated with H₂O₂, the expression of PwNAC4,11, 24, 34, 36, 41, 45, 47 and 55 was induced rapidly; the expression level of these genes in H₂O₂-treated plants was >1.3-fold higher at 3 h and that of PwNAC36 was 16.12-fold higher than that in the control (Figure 7). In contrast to the substantial changes in gene expression levels observed under salt and drought treatments, the low temperature treatment resulted in minimal changes in gene expression, suggesting that the NAC genes investigated in this study play more important roles in drought and salt stresses than in cold stress.

Figure 7. Expression profiles of 13 NAC genes under low temperature and hydrogen peroxide (H₂O₂) treatments.

Expression profiles of 13 candidate PwNAC genes in response to phytohormones and wounding

We examined the effect of five phytohormones including GA₃, IAA, SA, ETH, JA and ABA on the expression of PwNAC genes. GA₃ treatment induced the expression of PwNAC4, 11, 34, 36, 38, 41, 45, 51 and 55, while IAA treatment up-regulated the expression of PwNAC4, 18, 34, 36, 38, 41, 45, 52 and 55. Previously, SA, JA and ETH have been reported as three key signaling molecules involved in plant defense responses [43]. Here, we found that PwNAC4, 36, 38, 41, 45, 47 and 55 were induced by Me-JA, SA and ETH (Figure 8), suggesting these genes are potentially involved in plant defense against pathogens. NAC transcription factors are reported to participate in drought stress via the ABA-dependent or ABA-independent pathway [42]. Our results showed that ABA enhanced the expression of PwNAC4, 11, 18, 24, 38, 41, 47, 51, 52 and 55, and the expression patterns of these genes under ABA treatment were similar to those in the drought treatment, suggesting the presence of ABA-dependent pathway in P. wilsonii, which is involved in the response to drought stress (Figure 8).

Figure 8. Expression profiles of 13 NAC genes under phytohormone treatments.

Additionally, we detected the expression pattern of NAC genes in plants subjected to wounding stress (Figure 9). The expression of PwNAC4, 11, 18, 34, 36, 38, 41, 45, 52 and 55 was rapidly induced by wounding, suggesting the possibility that these genes act as early regulators of the response to mechanical damage.

Figure 9. Expression profiles of 13 NAC genes under the wounding treatment.

Conclusions

In this study, we investigated the transcriptome of P. wilsonii plants under salt and drought treatments using the Illumina sequencing technology. Approximately 36.26 Gb of data were generated from three treatments (control, salt and drought). A total of 4171 and 3438 DEGs were identified in plants exposed to drought and salt treatments, respectively. The RNA-Seq data of PwNAC genes were validated using the qRT-PCR assay. The expression of PwNAC11, 18, 24, 36, 38, 45, 47, 51, 52 and 55 was greatly induced by drought stress, whereas that of PwNAC11, 24, 36, 38, 41, 45, 47, 51 and 52 was induced by salt stress, suggesting the involvement of these genes in abiotic stress resistance in P. wilsonii. Overall, this study provides a database for further functional genomic research on woody plants. Data generated in this study represent a rich genetic resource for breeding.

Acknowledgment

We thank Biobreeding Biotechnology Corporation (Shaanxi, China) for conducting RNA-Seq and data analysis.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This work was supported by the Agricultural Ministry of China (under grant number 2016ZX08009003-002).