Transcriptome-Wide Identification and Expression Analysis of the NAC Gene Family in Tea Plant [Camellia sinensis (L.) O. Kuntze]

In plants, the NAC (NAM-ATAF1/2-CUC) family of proteins constitutes several transcription factors and plays vital roles in diverse biological processes, such as growth, development, and adaption to adverse factors. Tea, as a non-alcoholic drink, is known for its bioactive ingredients and health efficacy. Currently, knowledge about NAC gene family in tea plant remains very limited. In this study, a total of 45 CsNAC genes encoding NAC proteins including three membrane-bound members were identified in tea plant through transcriptome analysis. CsNAC factors and Arabidopsis counterparts were clustered into 17 subgroups after phylogenetic analysis. Conserved motif analysis revealed that CsNAC proteins with a close evolutionary relationship possessed uniform or similar motif compositions. The distribution of NAC family MTFs (membrane-associated transcription factors) among higher plants of whose genome-wide has been completed revealed that the existence of doubled TMs (transmembrane motifs) may be specific to fabids. Transcriptome analysis exhibited the expression profiles of CsNAC genes in different tea plant cultivars under non-stress conditions. Nine CsNAC genes, including the predicted stress-related and membrane-bound genes, were examined through qRT-PCR (quantitative real time polymerase chain reaction) in two tea plant cultivars, namely, ‘Huangjinya’ and ‘Yingshuang’. The expression patterns of these genes were investigated in different tissues (root, stem, mature leaf, young leaf and bud) and under diverse environmental stresses (drought, salt, heat, cold and abscisic acid). Several CsNAC genes, including CsNAC17 and CsNAC30 that are highly orthologous to known stress-responsive ANAC072/RD26 were identified as highly responsive to abiotic stress. This study provides a global survey of tea plant NAC proteins, and would be helpful for the improvement of stress resistance in tea plant via genetic engineering.


Introduction
Tea is a non-alcoholic drink known for its bioactive ingredients and health efficacy. Tea plant [Camellia sinensis (L.) O. Kuntze], which originated from the southwest region of China, had been cultivated and utilized for thousands of years, and is planted worldwide nowadays [1].
Surviving in the wild, tea plant often suffers from numerous detrimental environmental factors, such as extreme temperature, drought and salinity [2,3]. Enhanced environmental receptivity may significantly help improve the yield and quality of tea plant. The potential adversity resistance of plants is often determined by expressing stress-inducible genes that are regulated by specific transcription factors (TFs) [4]. TFs widely exiting in plants regulate the expression of target genes by binding directly or indirectly with specific cis-regulatory elements.
NAC (NAM-ATAF1/2-CUC) TF is a large TF family widely existing in plants. Most members of the NAC family contain a conserved NAC domain (~150 aa), which consists of five subdomains (A, B, C, D, and E), in the N-terminus. In general, subdomains A, C, and D are tightly conserved, whereas subdomains B and E are divergent [5]. Subdomains B and E may be closely related to the functional diversity of NAC genes, subdomain A may participate in dimer formation, and subdomains C and D are mainly involved in DNA binding (DB) [6][7][8]. The C-terminal regions of NAC TFs are highly variable, which might function as transcriptional activators or repressors that regulate the expression of downstream genes [9]. Moreover, some NAC TFs contain transmembrane motifs (TMs), which are working for the plasma membrane or endoplasmic anchoring, at the C-terminal end [10]. X-ray crystallography of Arabidopsis ANAC019 and rice SNAC1 revealed that the NAC domain monomer exhibits a novel TF folding pattern, consisting of a twisted antiparallel β-sheet, which is devoted to DB and encircled by an α-helix element on both sides [6,7].
Numerous NAC genes participate in the response to environmental stresses and hormone signaling. ANAC019, ANAC055 and ANAC072/RD26 are early-identified and well-characterized stress-related Arabidopsis NAC genes that are induced by drought, salinity, cold and ABA (abscisic acid) [33]. Microarray analysis of transgenic plants overexpressing either of these genes revealed that stress-inducible genes are upregulated and drought tolerance is significantly improved in these plants [33,34]. Similar to ANAC072/RD26, ATAF1 is another Arabidopsis NAC gene whose overexpression enhances plant tolerance to drought, ABA, salt, oxidative stress, and necrotrophic pathogen [35]. The overexpression of several rice NAC genes significantly enhances tolerance and maintains grain yield [36][37][38].
Genomic sequencing of tea plant has not been completed. The completed transcriptome sequencing of tea plant provides an opportunity to identify gene families at the transcriptome level [39]. Recent studies have identified and analyzed 89 AP2/ERF, 50 WRKY, 16 HSF, and 18 bZIP genes in tea plant; the roles of these genes in response to adversity stresses were also elucidated [40][41][42][43]. To date, detailed analysis of NAC family genes in tea plant remains to be conducted.
In the present study, the transcriptome data was just utilized to extract the NAC sequences of tea plant. A total of 45 CsNAC genes were identified in tea plant based on RNA-Seq data [39]. Subsequently, multiple sequence alignment, conserved motif analysis, phylogenetic construction and membrane-bound proteins identification were performed. On the basis of evolutionary relationships and structural analysis, nine NAC genes were selected for quantitative real-time PCR analysis. The expression levels of these genes were investigated in various tissues (bud, young leaf, mature leaf, young stem and root) and under diverse environmental stresses (drought, salt, heat, cold and ABA) in two cultivars of tea plant, namely, 'Huangjinya' and 'Yingshuang'. This study provides novel insights into the structures and functions of NAC genes in tea plant and serves as a valuable resource for the improvement of plant stress tolerance.

Identification of NAC gene family in tea plant
The key words "NAC" and "NAM" were used to retrieve potential unigenes from the summary of gene annotations of the tea plant transcriptome [39]. To identify the candidate NAC genes, BLASTp search (http://blast.ncbi.nlm.nih.gov) of the retrieved proteins were performed. Only those proteins with E-values less than 1e −10 and NAM domains above 100 were collected for further investigation. The accuracy of several CsNAC genes were validated by PCR paired-end sequencing, and completed by the GenScript Corporation (Nanjing, China). The primer designed from the 5' UTR and 3' UTR ( Table 1).
The Pfam database (http://pfam.sanger.ac.uk) was used, and hidden Markov models (HMMs) of the collected CsNACs were obtained. These models could serve as another proof of the confirmed CsNACs. To obtain the Arabidopsis orthologs of the NAC proteins of tea plant, all valid sequences were subjected to BLASTp search at Arabidopsis protein TAIR10 release (http://www.arabidopsis.org) by using default parameters. Furthermore, TMHHM server v. 2.0 (http://www.cbs.dtu.dk/services/TMHMM/) was used, and the membrane-bound proteins of CsNAC were predicted.
Multiple sequence alignment, phylogenetic analysis, and conserved motif analysis Full-length CsNAC protein sequences and three representative ANAC proteins (ANAC019, ANAC055 and ANAC072/RD26) were aligned using the ClustalX 1.83 program with default parameters. To survey the phylogenetic relationships of NAC genes in tea plant, domain sequences of NAC proteins in tea plant and Arabidopsis were used to construct a phylogenetic tree. After multiple alignments were carried out, phylogenetic trees were generated and displayed by the MEGA program (version 5.0) [44] with the following parameters: neighbor-joining (NJ) method, p-distance model, pairwise deletion, and 1000 bootstrap replicates. The Multiple Em for Motif Elicitation (MEME) suite (version 4.10.1) (http://meme-suite.org/tools/ meme) was used to identify the conserved motifs, and the default parameters were used, in addition to the maximum number of motifs (10).

Preparation of plant materials
The tea plant cultivars 'Yingshuang' and 'Huangjinya' are biennial cutting seedling and were grown in a greenhouse at the Nanjing Agricultural University (Nanjing, China). The artificial climate conditions were 23˚C temperature, 14/10 h light/dark, and 70% relative humidity. Tea plants were subjected to various stress treatments. For extreme temperature treatments, the seedlings were maintained at 4˚C and 38˚C, and some seedlings were irrigated with 20% PEG 6000 and 200 mM NaCl solution, respectively. Additionally, some seedlings were sprayed with 200 μM ABA solution.
The young leaves of all tested plants were harvested at 0, 2, 8, and 24 h after treatments, rapidly frozen in liquid nitrogen, and then stored at −80˚C.
To study the expression patterns of CsNAC genes in different tissues, including bud, young leaf, mature leaf, young stem and root, were collected from growing plants under normal conditions. All harvested samples were rapidly frozen in liquid nitrogen and then stored at −80˚C.

RNA isolation, cDNA synthesis, and qRT-PCR analysis
Total RNA from tea plant materials was extracted using the RNA Isolation Kit (Huyueyang, Beijing, China) in accordance with the manufacture's instruction. Before cDNA synthesis, RNA integrity was verified using 1.5% denaturing agarose gel, and purity and concentration were measured on a NanoDrop ND-2000 Spectrophotometer (Thermo Scientific, USA). Moreover, the PrimeScript RT reagent Kit (TaKaRa, Dalian, China) was used, and 1 μg of high-quality RNA samples were programed for reverse transcription into cDNA in accordance with the operation manual.
The specific primers used for calculating the relative expression were designed from the nonconserved region of genes using Primer Premier 5.0 software. Actin served as a reference gene [45]. The primer sequences used in this experiment are provided in Table 2. The SYBR Premix Ex Taq kit (TaKaRa, Dalian, China) was used for qRT-PCR on the Bio-rad IQ5 fluorescence quantitative PCR platform. The experiments were completed in a volume of 20 μL: 0.4 μL of each specific primer, 10 μL of 2 × SYBR Premix, 2 μL of diluted cDNA as a template, and 7.2 μL of ddH 2 O. The reactions were carried out under the following conditions: denaturation at 95˚C for 30 s, followed by 40 cycles of 95˚C for 10 s and 60˚C for 20 s. The experiments were performed with three independent biological replicates, and the relative expression level was calculated using the 2 −ΔΔCT method [46].

Identification of NAC family members in tea plant
According to the RNA-seq database, the annotated sequences with "NAC" and "NAM" were searched, and 120 NAC unigenes were obtained. Their corresponding amino acid sequences were subjected to the NCBI BLASTp program. Finally, 45 CsNAC genes, designated as CsNAC1-CsNAC45 (S1 Table), were identified and used for further investigation, although still some of these genes were still fragments. Eight CsNAC genes (CsNAC16, CsNAC17, CsNAC26, CsNAC29, CsNAC30, CsNAC33, CsNAC39 and CsNAC45) were examined by PCRdirect sequencing in tea plant cultivar 'Yingshuang', and the genes sequences showed very minor discrepancy (S1-S8 Figs). The remaining 75 sequences with an E-value more than 1e −10 or a NAM domain length less than 100 were excluded. HMMs were excellent implement in the bioinformatic analysis of biological sequences and the discovery of secondary structures [47]. In the present study, HMMs of the 45 predicted CsNAC proteins were determined using the Pfam database; the results showed that these CsNAC proteins indeed contained a conserved NAM domain (PF02365) at the N-terminus (S2 Table). This result could also be a powerful evidence for the correctness of the valid sequences.
We retrieved the functional annotations of the 45 CsNAC proteins predicted by BLASTX public databases GO, Nt, Nr, COG, Swiss-Prot, and TrEMBL (E-value 1E-5) (S3 Table). The predicted CsNAC proteins are involved in various abiotic stresses, such as osmotic stress, water deprivation, salt, heat, cold, and wounding, and biotic stresses, such as bacterium, fungus, nematode, and insect. CsNAC proteins are also involved in signal transduction, including various plant hormone signaling pathways, such as ethylene, abscisic acid, jasmonic acid, salicylic acid, and brassinosteroid. In addition, some CsNAC proteins may be implicated in the development of tissues and organs, such as leaf, seed, embryo, ovule, pollen, and lateral root. Moreover, four tea plant NAC proteins, namely, CsNAC2, CsNAC21, CsNAC26, and CsNAC33, are associated with flavonoid biosynthesis. Catechin (flavan-3-ols), a vital secondary metabolite of tea plant, is produced via flavonoid biosynthesis [48]. This finding indicates that NAC may be involved in catechin biosynthesis. The above annotations describe the potential functions of NAC proteins in tea plant, and the intensity of functions may be discrepant, which needs further research.

Phylogenetic analysis of the NAC TFs between tea plant and Arabidopsis
Several NAC genes have been elucidated in Arabidopsis, a typical model plant. To clarify the phylogenetic relationships of NAC family proteins in tea plant and Arabidopsis, an unrooted Table 2. Primers for qRT-PCR of CsNAC genes.

Target Gene Forward primer sequence (5'!3') Reverse primer sequence (5'!3')
tree was established on the basis of the aligned NAC domain (A-E) sequences. The results indicated that most of the obtained subgroups were consistent with previous phylogenetic analyses [5,20]. As shown in Fig 1, the phylogenetic tree clustered all of the NAC members into 17 subgroups. CsNAC proteins were diverse as ANAC proteins, and 45 CsNACs were unevenly distributed in 16 subgroups. By contrast, no member was found in subgroup ANAC001. Notably, the same phenomenon was observed in the phylogenetic analysis of rice and Arabidopsis NAC proteins. All subgroup ANAC001 members belong only to Arabidopsis [5], indicating that this subgroup may either be lost in tea plant and rice or acquired in Arabidopsis after the divergence of their last common ancestor. Furthermore, these NACs may have specialized roles in Arabidopsis. Therefore, the characteristics of this subgroup may be important in researching the genetic relationships among other plant species.
Genes showing a close evolutionary relationship generally exhibit similar functions. Phylogenetic analysis can be utilized to predict gene function, which is important in subsequent functional studies. Predicting stress-related functional genes on the basis of phylogenetic relationships is effective [14,49]. Subgroups AtNAC3, ATAF, and NAP share a close relationship, and most published stress-related NAC family members are included in these three subgroups. Arabidopsis contains three members (ANAC019, ANAC055 and ANAC072/RD26) of subgroup AtNAC3 and four members (ATAF1/ANAC002, ATAF2/ANAC081, ANAC032 and ANAC102) of subgroup ATAF, all of which are involved in multiple-stress response [9,33,50]. However, few reports refer to subgroup NAP for the stress resistance in Arabidopsis until now. The subdomain E of ATAF, AtNAC3 and NAP members is highly conserved and may determine the function of NAC proteins, indicating that subgroup NAP might also be involved in stress resistance [5]. ANAC029 is closely associated with leaf senescence in Arabidopsis [28]. Senescence is a common phenomenon when plants are subject to various serious adversity stresses. A minimum of three NAC genes belong to subgroup NAP in Chrysanthemum lavandulifolium, which could be induced by salt, drought, and salicylic acid [51]. These findings indicate that the subdomain NAP belong to the stress-related group. Therefore, we consider CsNAC17 and CsNAC30 from subgroup AtNAC3; CsNAC18, CsNAC31, and CsNAC32 from subgroup ATAF; and CsNAC26 from subgroup NAP as stress-related proteins.
To understand further the characteristics of NAC genes in tea plant, we confirmed that the NAC proteins of Arabidopsis are orthologous to 45 CsNACs (S4 Table). Using a Blastp search, the most satisfactory match of Arabidopsis protein was regarded as orthologous. The score and E-value were used to explain the satisfactory match of orthologous. As showed in S4 Table, CsNAC17, CsNAC29, and CsNAC30 were orthologous to Arabidopsis ANAC072/ RD26, of which CsNAC17 and CsNAC30 were highly orthologous to Arabidopsis ANAC072/ RD26 with more than 300 scores and powerful E-value support. CsNAC18, CsNAC31 and CsNAC32 were highly orthologous to Arabidopsis ANAC002/ATAF1 with at least 200 scores and powerful E-value support. The homologous relations of tea plant and Arabidopsis considerably matched the phylogenetic analysis. CsNAC29 exhibited a close evolutionary relationship to ANAC096, which functions synergistically with ABF2 and ABF4, and helps plants survive under dehydration and osmotic stresses [52]. Finally, CsNAC17, CsNAC18, CsNAC26, CsNAC29, CsNAC30, CsNAC31, and CsNAC32 were considered the potential stress-related genes.

Gene structure and conserved motif analysis
The diversity of plant protein sequence and structure generated in biological evolution is a possible mechanism for the formation of multigene families and functional diversities; the diversity characteristics of plants lead to the efficient use of natural resources or adapt to adverse environments [53]. To understand further the structural features of NAC proteins in tea plant, multiple sequence alignment of full-length CsNAC proteins was performed, and conserved motifs were predicted in concert with their phylogenetic relationships.
To examine the structure of CsNAC proteins, multiple sequence alignment of the fulllength CsNAC protein sequences, along with three representative ANAC proteins of Arabidopsis was conducted (S9 Fig). In addition, multi-sequence alignment of eight sequenced CsNAC proteins (CsNAC16, CsNAC17, CsNAC26, CsNAC29, CsNAC30, CsNAC33, CsNAC39 and CsNAC45) was performed (Fig 2). As expected, most of the CsNAC proteins shared a typical highly conserved N-terminal domain, which was divided into five subdomains (A, B, C, D, and E), and a highly variable C-terminal transcriptional regulation domain. The sequence alignments in S9 Fig showed that subdomains A, C, and D were more conserved than subdomains B and E, even though some domain regions were incomplete [5]. In subdomain D, a conserved bipartite nuclear localization signal was identified in most of the CsNAC proteins. This signal has been identified from many NAC proteins in other plant species, indicating these CsNACs are nuclearly localized [14,54,55].
To reveal the sequence features of CsNAC genes, the conserved motifs were predicted using the MEME program, and 10 conserved motifs were defined (Table 3 and S10 Fig). Members with a close evolutionary relationship displayed uniform or similar motif compositions (Fig 3). In the majority of the CsNAC proteins, motifs 2 plus 7, 5, 1 plus 10, 3 plus 4, and 8 corresponded to the conserved subdomains A, B, C, D and E, respectively. Almost all of the conserved motifs   [20,21,55]. Ooka et al. investigated the C-terminal region of the NAC proteins from Arabidopsis and rice; 13 motifs were found, despite the C-terminal region was highly divergent [5]. As a whole, the distribution of the main conserved motifs in tea plant and other species was similar; the different nature of growing habits, such as annual vs perennial and woody vs herbaceous, that may contribute to the differences. The distribution of the conserved motifs indicates that NAC protein functional characteristics may be mainly determined by the N-terminal conserved domain and that C-terminal region may also be involved in determining the functions of these proteins.

Membrane-associated CsNAC subfamily
Recently, membrane-associated TFs (MTFs) have received considerable attention, because of their critical roles in the gene regulatory network of plants [56]. Commonly, these TFs exist in a dormant state; once stimulated, the dormant form becomes activated through proteolytic cleavage; with the degradation of their cytoplasmic anchors, the activated TFs enter the nucleus and regulate the expression of target genes [56,57]. In the present study, using TMHHM server v. 2.0, we identified three CsNAC proteins (CsNAC2, CsNAC9, and CsNAC12) to contain a single α-helical TM at C-terminal ends (S5 Table). TM commonly exists in highly variable regions of the C-terminus, indicating that many potential MTFs are existing. Sufficient evidence has indicated that NAC MTFs play important roles in plant growth and development as well as in response to abrupt environmental changes [57,58]. All of the MTF NAC genes in Arabidopsis and rice could be induced by at least one type of environmental stress, namely, drought, salt, cold or heat [9,59]. The phylogenetic relationship of NAC MTFs among tea plant, Arabidopsis and rice showed that CsNAC MTFs have a close relationship with Arabidopsis MTFs. Thus, NAC MTFs in tea plant might have similar functions to those in Arabidopsis (Fig 4). The finding indicates that CsNAC MTFs may be involved in stress response. Several NAC MTFs were identified among various species through intensive genome-wide analysis. We investigated the NAC TFs whose MTFs have been analyzed in whole genomes and 12 species were collected (Fig 5), including A. thaliana [58], O. sativa [58], G. max [14], Solanum tuberosum [55], Sorghum bcolor [60], S. italic [13], B. rapa [61], Cicer arietinum [62], Zea mays [15], Brachypodium distachyon [63], M. esculenta [20], and J. curcas [21]. The NAC genes of some species were not used for analysis. For example, 104 NAC genes in S. lycopersicum [19], 145 in G. raimondii [53] and 74 in V. vinifera [17] were identified in genome-wide analysis, but the NAC MTFs were not analyzed. All of the identified NAC MTFs belonged to monocots and eudicots, although some NAC genes had been found in moss and lycophyte in previous studies. Obviously, the total number of NAC genes significantly differs in different species and may mainly undergo paleopolyploidy. For example, Chinese cabbage has undergone triplication since its divergence from Arabidopsis [64]. The proportion of NAC MTFs in monocots was obviously less than that in eudicots. Ha et al. [62] showed that only G. max and C. arietinum possess two TMs, suggesting that the existence of doubled TMs is specific to leguminous plants. In the present study, Euphorbiaceae (M. esculenta and J. curcas) also possessed doubled TMs. Coincidentally, Euphorbiaceae and leguminous plants could be classified in fabids, suggesting that the existence of doubled TMs is wide and specific to fabids plants. The finding indicates that NAC MTFs are closely related to plant evolution.

Expression profiles of CsNAC genes in four tea plant cultivar transcriptomes
Illumina RNA-seq data were obtained to assess the expression profiles among different tea plant cultivars of CsNAC genes under non-stress conditions [39]. The transcript abundance of 45 CsNAC genes was assessed in accordance with the RPKM values of four tea plant cultivars, namely, 'Yunnanshilixiang' (Tea_T1), 'Chawansanhao' (Tea_T2), 'Ruchengmaoyecha' (Tea_T3), and 'Anjibaicha' (Tea_T4), although some genes were lowly or barely expressed in some tea plant cultivars (S6 Table). A heat map was displayed on the basis of log2 transformed RPKM values (Fig 6). CsNAC2, CsNAC9, CsNAC12, CsNAC30, and CsNAC44 genes displayed high and stable expression levels in all of the tea plant cultivars. By contrast, CsNAC1 and CsNAC38 were barely expressed. Remarkably, CsNAC2, CsNAC9, and CsNAC12 belonged to NAC MTFs. The gene expression levels are highly associated with gene function, which indicates that NAC MTFs play important roles in the growth and development of tea plant. In addition, some genes showed obvious cultivar specificity. For example, CsNAC31 was highly expressed in Tea_T2, but weakly expressed in the three other tea plant cultivars. CsNAC4, CsNAC7, and CsNAC40 completely differed among the four tea plant cultivars. Twenty CsNAC genes obviously differed among three tea plant cultivars. In Tea_T1, CsNAC genes showed high abundance. Sixteen genes were highly expressed among at least one tea plant cultivar, and only one gene was not expressed in Tea_T1. In Tea_T3, CsNAC genes showed relatively low expression levels. Seventeen CsNAC genes were not detected in at least one tea plant cultivar, and all of them happened in Tea_T3; while Tea_T2 had four, and Tea_T1 and Tea_T4 only had one (S6 Table). In summary, the RPKM values revealed that the four tea plant cultivars differed in the expression of CsNAC genes. These tea plant cultivars are grown in different regions of China [39,40], and the differences may occur on the individual evolution of different cultivars.

Expression profiles of CsNAC genes in different tissues of tea plant
'Huangjinya' and 'Yingshuang' are two of the widely planted tea plant cultivars in China. 'Huangjinya', a light-sensitive albino tea plant cultivar, has a potential for processing high quality green tea even in summer, but it shows weak stress resistance [65]. 'Yingshuang', a tea plant cultivar with good resistance to abiotic stresses, especially cold stress, possesses considerable agronomic traits, such as appropriate phenol ammonia content and good cold resistance. Based on the above analysis of the transcriptome, four low expression genes (CsNAC1, CsNAC15, CsNAC35 and CsNAC38) and seven high expression genes (CsNAC2, CsNAC9, CsNAC12, CsNAC30, CsNAC32, CsNAC44 and CsNAC45) were selected for verification through qRT-PCR analysis in two tea plant cultivars, 'Huangjingya' and 'Yingshuang' (Fig 7). The four low expression genes were still barely expressed, and the seven high expression genes displayed high and stable expression levels in two tea plant cultivars, 'Huangjingya' and 'Yingshuang'.
NAC genes are involved in various developmental processes [17,66]. As natural properties of each gene, tissue-specific expression is closely associated with the growth and development of the particular tissue. In the present study, nine CsNAC genes, including six predicted stressrelated and three MTF CsNAC genes, were selected for subsequent qRT-PCR analysis. The expression patterns of these genes were assessed in five tissues, including root, young stem, mature leaf, young leaf and bud in two tea plant cultivars, namely, 'Huangjinya' and 'Yingshuang'.
These nine NAC genes widely exist in the five tissues of the two tea cultivars, and the expression patterns varied (Fig 8). Compared with other genes, three MTF CsNAC genes, namely, CsNAC2, CsNAC9, and CsNAC12, exhibited higher expression levels in all of the tissues. This result is consistent with the above expression analysis among different cultivars, suggesting the critical role of the genes in plant growth and development. In general, the same NAC genes exhibited similar tissue-specific expression profiles in the two tea cultivars. For example, in the two tea plant cultivars, CsNAC2 and CsNAC18 were highly expressed in young stem and bud, CsNAC12 in young stem and young leaf, and CsNAC26 in stem and mature leaf. However, CsNAC26 was lowly expressed in young leaf and bud. Discrepant expression patterns were also observed. For example, CsNAC2, CsNAC9, and CsNAC26 showed higher expression level in the root of 'Huangjinya' than in that of 'Yingshuang'; CsNAC17 showed the highest expression level in the root of 'Yingshuang' and in the bud of 'Huangjinya'; CsNAC30 and CsNAC32 showed higher expression levels in the root of 'Yingshuang' than in that of 'Huangjinya'. Tissue-specific genes facilitate the development of particular tissues. NAC1 and ANAC2, which are overexpressed in the roots of transgenic Arabidopsis, promote the development of lateral roots [25,67]. Arabidopsis ANAC036, which is strongly expressed in the leaves of transgenic plants, shows a semi-dwarf phenotype [68]. In general, the tissue expression profiles of NAC genes in different cultivars provide an evidence for further investigation on the development of tea plant.

Expression profiles of CsNAC genes under various stress treatments in tea plant
Environment stress adversely affect plant growth and productivity, and trigger a series of morphological, physiological, biochemical and molecular changes [69]. Considerable evidence has shown that NAC genes play vital roles in the response to biotic or abiotic stresses [9,70]. To explore the roles of NAC genes in tea plant under diverse environmental conditions, nine CsNAC genes, including six stress-related and three MTF CsNAC genes, were selected and subjected to qRT-PCR expression analysis in response to multiple stress treatments in 'Huangjinya' and 'Yingshuang'. Such treatments include drought (20% PEG 6000), salinity (200 mM NaCl), heat (38˚C), cold (4˚C) and ABA (200 μM).
Drought treatment. Under drought treatment (Fig 9), most of the CsNAC genes were upregulated in both tea plant cultivars. The relative transcript levels of seven CsNAC genes (CsNAC2, CsNAC9, CsNAC12, CsNAC17, CsNAC29, CsNAC30, and CsNAC32) gradually increased and reached the highest level at 24 h. CsNAC2 and CsNAC12 were initially downregulated. CsNAC18 reached a maximum of about two-fold despite the decrease at 2 and 8 h in 'Huangjinya' and 'Yingshuang', respectively. CsNAC26 was upregulated at all time points in 'Huangjinya' but was relatively stably expressed in 'Yingshuang' despite the decrease at 8 h. Notably, CsNAC17, CsNAC30 (ortholog of ANAC072/RD26), and CsNAC32 (ortholog of ANAC002/ATAF1) were highly induced (over 20-fold) in at least one tea plant cultivar, indicating their possible roles in drought stress responses.
Salt treatment. Under salt treatment (Fig 10), four CsNAC genes (CsNAC9, CsNAC17, CsNAC30, and CsNAC32) were mainly upregulated in both tea plant cultivars. By contrast, CsNAC12 exhibited a slight decline. CsNAC2 showed a relatively stable expression level. CsNAC26 in 'Huangjinya' and CsNAC29 in 'Yingshuang' were also relatively stably expressed. Meanwhile, CsNAC26 decreased after 8 h in 'Yingshuang', and CsNAC29 slightly increased in 'Huangjinya'. CsNAC18 was significantly upregulated in 'Yingshuang' at 2 h. CsNAC30 (ortholog of RD26) was the most induced by salt treatment with 107-and 28-fold induction in 'Huangjinya' and 'Yingshuang', respectively. However, this gene was declined at 8 h and reached another peak at 24 h in both tea plant cultivars, indicating the complexity of gene regulatory networks. High temperature treatment. Under high temperature treatment (Fig 11), seven CsNAC genes (CsNAC2, CsNAC17, CsNAC26, CsNAC29, CsNAC30, and CsNAC32) gradually increased and reached the highest level at 24 h. By contrast, CsNAC17, CsNAC26, and CsNAC30 were initially downregulated. CsNAC12 decreased in both tea plant cultivars. CsNAC9 also slightly decreased in 'Huangjinya' but showed a more stable expression in 'Yingshuang'. After 2 h, CsNAC18 was gradually upregulated in 'Huangjinya' but downregulated in 'Yingshuang'.
Low temperature treatment. Low temperature is a type of abiotic intimidation that significantly affects tea plant growth and productivity, particularly the cold of the late spring or late frost [2,69]. Under cold treatment (Fig 12), CsNAC32 was downregulated in both tea plant cultivars despite the high expression at 2 h in 'Huangjinya'. All of the other CsNAC genes gradually increased and reached the highest level at 12 or 24 h. Nevertheless, CsNAC17 and CsNAC26 were initially downregulated in 'Yingshuang'. CsNAC12 and CsNAC18 slightly increased with less than two-fold in both tea plant cultivars.
In general, most CsNAC genes were affected by drought, salinity, and high and cold stresses. Although the intensity of stress responses was varied, most CsNAC genes exhibited a similar tendency in the two tea plant cultivars. Whole-genome expression analysis in Arabidopsis showed that most of the NAC genes are responsive to salt and extreme temperatures [8,71].
SiNAC expression profiles showed that all the collected genes display varied expression patterns in response to one or more stresses and that cold stress induces relatively more dramatic changes in transcript abundance than dehydration or salinity [13]. The varied gene expression profiles suggest that CsNACs control a complex gene regulatory network and exert a regulatory effect on various physiologic functions for acclimatizing multiple challenges. Moreover, the expression levels of Arabidopsis RD26 orthologs StNAC072 and StNAC101 were highly induced by stress and ABA treatments [55]. Arabidopsis RD26 orthologs CaNAC06 and CaNAC67 are highly induced by dehydration [62], and MeNAC22 (RD26 orthologs) is strongly induced by osmotic and drought stresses [20]. Arabidopsis RD26 orthologs CsNAC17 and CsNAC30 were highly induced by at least one stress, which agrees with previous reports. The finding indicates that other species of Arabidopsis RD26 orthologs are most likely highly induced by stresses. Some research indicated that TF genes that are highly induced by stress could be preferentially utilized for further plant functional studies because of their potential in the development of improved stress-tolerant transgenic plants via overexpression. CsNAC17 and CsNAC30 were exhibited to be appropriate candidate genes for further plant research in tea plant.
ABA treatment. ABA is a plant hormone that not only participates in plant growth and development, but also plays a crucial role in the regulation of various stress responses [72]. NAC genes may be regulated by ABA-dependent or ABA-independent pathways because of the difference in their promoter elements [8,9,71,73]. Under ABA treatment (Fig 13), CsNAC17 initially decreased and then increased in both tea plant cultivars. CsNAC12, CsNAC26, and CsNAC30 were mainly downregulated, but CsNAC26 was upregulated in 'Huangjinya' at 24 h. CsNAC9 was upregulated in both tea plant cultivars, which reached a peak in 'Yingshuang' at 24 h and in 'Huangjinya' at 2 h. CsNAC32 was upregulated in 'Huangjinya' but was relatively stably expressed in 'Yingshuang'. The above six CsNAC genes responded to at least one stress, indicating that these genes may regulate stress responses in tea plant in an ABA-dependent manner. In addition, CsNAC2, CsNAC18, and CsNAC29 showed relatively stable expression levels in both tea plant cultivars. Therefore, CsNAC2, CsNAC18, and CsNAC29 may be regulated in ABA-independent pathways. Hu and his colleagues identified that 19 CaNAC genes are responsive to dehydration; of these genes, 7 are ABA independent and 12 are ABA dependent [62]. The response of CsNAC genes to ABA treatment suggests that these genes play roles in ABA signaling (Fig 14).

Conclusions
Genomic sequencing of tea plant has not been completed. Thus, transcriptome sequencing has a great potential to discover and identify novel genes and will accelerate the understanding of the regulatory functions of NAC TFs. On the basis of transcriptome sequences, 45 CsNAC genes were identified. Phylogenetic analysis, as well as the identification and analysis of motif and MTF members, provided an insight into the functional diversity of CsNAC proteins. The comparative analysis of CsNACs with their corresponding Arabidopsis orthologs helped predict the potential functions of CsNAC proteins. The cultivar and tissue expression profiles of CsNAC genes under normal growth conditions were surveyed. Furthermore, the expression analysis of certain CsNAC genes during various treatments, such as drought, salinity, heat, cold, and ABA, established substantial evidence to select candidate stress-resistant genes, which could be preferentially utilized to develop tea plants with improved resistance under stress conditions. Tea plant suffers from numerous abiotic stresses, such as heat, cold, drought and salinity (Fig 14). Through either ABA-dependent or -independent manner, CsNAC TFs could regulate downstream target genes, and accumulate stress-inducible genes or metabolites,