Genome-wide characterization and expression profiling of NAC transcription factor genes under abiotic stresses in radish (Raphanus sativus L.)

Database search and detection of NAC transcription factors

The protein and nucleotide sequence datasets were obtained from the radish genome database for the detection of the RsNAC gene family (Mitsui et al., 2015). Radish protein sequences were retrieved using Hidden Markov Model profile (PF02365) generated from Pfam (Punta et al., 2012) with a cut-off E-value of 1e−3. To confirm the presence of conserved NAC domains, multiple alignment of full-length RsNAC protein sequences was carried out using the clustalW program with default pairwise and multiple alignment parameters. The data for Arabidopsis, B. rapa and rice NACs sequences were retrieved from the Arabidopsis Information Resource (Huala et al., 2001), Brassica Database (Cheng et al., 2011) and the Rice Genome Annotation Project (Kawahara et al., 2013), respectively. The sequences of all NAC TFs in other representative species were retrieved from the plant TF database (Jin et al., 2014). Additionally, TMHMM Server v. 2.0 (Krogh et al., 2001) was employed in discovering membrane-bound RsNAC proteins.

Phylogenetic analysis of RsNAC protein domains

Phylogenetic analysis was done using RsNAC domain sequences with previously grouped NAC domains from Arabidopsis and rice as reference sequences of eudicots and monocots, respectively (Abrouk et al., 2010). The phylogenetic tree was generated using MEGA6 program by the neighbor-joining method with 1,000 bootstrap replicates (Tamura et al., 2013).

Protein characterization and RsNAC sequence analyses

The Protparam program (Gasteiger et al., 2005) using ProtParam tool was employed to determine the molecular weight (MW) and isoelectric points (pI) using default settings. The online MEME program (Bailey et al., 2009) was used to search conserved motifs. The analysis was performed with default parameters except the following: number of repetitions, any; maximum number of motifs, 15; and optimum width of the motif, ≥50. The information for each RsNAC exon/intron organization was determined by mapping the CDS to DNA sequences using the Gene Structure Display Server 2.0 GSDS (Hu et al., 2014). To determine Arabidopsis orthologs, radish NAC proteins were used in BLASTP search against the Arabidopsis Information Resource, TAIR10 release (Huala et al., 2001). Reciprocal best BLAST hit blasting approach(Moreno-Hagelsieb & Latimer, 2008) was used and best hits were considered when bit-score of the compared proteins was more than 200, the sequence coverage value more than 50% and expected value E- value of ≤1e–3.

Chromosomal location of RsNAC genes

The sequences of RsNAC retrieved from the genomic sequences of the scaffolds were anchored to the integrated genetic map published byJeong et al. (2010). All the acquired potential RsNAC genes were searched against radish the chromosome sequence database Jeong et al. (2010) using BLASTN program. The sequences with a similarity of ≥99% and length coverage difference ≤5 base pairs (bp) were regarded as similar genes (Wang et al., 2015b). Subsequently, all the similar RsNAC gene locations (megabase, Mb) on the genetic map were annotated according to their corresponding pseudomolecules. All anchored genes were located on their specific genetic map against their location parameters using MapInspect software (http://www.plantbreeding.wur.nl/UK/software_mapinspect.html). Tandem repeats located within four predicted genes or within 50 kb from each other were identified manually and marked on the radish genetic map (Cannon et al., 2004).

Expression analysis of RsNAC genes based on the RNA-Seq data

The normalized RNA-Seq data were retrieved from Nodai Radish Genome Database (Mitsui et al., 2015), and then used for the expression profiling of the RsNACs in two tissues (roots and leaves) and six developmental stages (7, 14, 20, 40, 60 and 90 days after sowing, DAS) represented by RPKM (Reads Per Kilobase of transcript per Million mapped reads) value.

Plant growth, abiotic stress treatments and expression analysis

The ‘NAU-YH’ genotype seedlings were germinated on a wet filter paper in darkness for three days then transferred to the pots containing soil and peat media (1:1) and cultured in the growth chamber under 16-h (25 °C) and 8-h (18 °C) light and dark condition, respectively. Thirty days after transplanting, seedlings were subjected to the different treatments. For heat and salt treatments, seedlings were exposed to 42 °C and 200 mM NaCl conditions, respectively. For HMs treatments, seedlings were exposed to 20 mg L⁻¹ CdCl₂ ⋅2.5H₂O and 100 mg L⁻¹ Pb (NO₃)₂, respectively. For drought treatment, seedlings were treated with 20% polyethylene glycol 6000 (PEG) while 100 µM ABA was sprayed to the seedlings to induce ABA response. The samples from all the treatments were harvested after 24 h and frozen immediately in liquid nitrogen for RNA extraction. Total RNAs from triplicate samples were isolated using Trizol (Invitrogen, Carlsbad, CA, USA) according to manufacturer’s protocol and then digested with RNase-free DNase. The quality of total RNA was determined by agarose gel electrophoresis . The total RNA (2 µg) was used to synthesize the first-strand cDNA using PrimeScript™ RT Reagent Kit (TaKaRa, Dalian, Liaoning Sheng, China), and the obtained cDNA products were further diluted 1:9 with nuclease-free water for RT-qPCR. The transcript levels of 21 randomly selected RsNAC genes were investigated with RT-qPCR according to the method reported by Xu et al. (2013). The radish Actin gene was used as the internal control. The Beacon Designer 7.7 was used to generate gene-specific primers from the non-conserved region of NAC sequences (Table S1).

Identification of the RsNAC genes and comparative analysis in radish

In the determination of NAC genes in radish, a total of 181 NAC genes were obtained after searching the entire radish genome by the HMMER software. Furthermore, SMART and Pfam programs (Letunic, Doerks & Bork, 2012; Punta et al., 2012) were used to further screen the presence of the NAM domain in all the RsNAC sequences with a cut-off E-value of 1e–3. Finally, 172 RsNAC proteins were selected for subsequent analysis after exclusion of nine members with either no N-terminal NAM domain or with an E-value of >1e–3. For nomenclature, the family designation (NAC) was preceded by prefix ‘Rs’ for radish, and then followed by a number based on their generic order. The comprehensive list of RsNAC proteins analyzed in this study, the annotation of each sequence and their preatomic characteristics are shown in Table S2. The number of aa in RsNAC proteins ranged from 143 (RsNAC143) to 652 (RsNAC011) residues resulting in an average of ∼336 aa, whereas the pI ranged from 4.44 (RsNAC0730) to 9.63 (RsNAC173) with an average of ∼6.6. The results indicated that different NACs genes might function in diverse microenvironments (Kiraga et al., 2007). The RsNAC domains were aligned and observed by GeneDoc program. The RsNAC proteins possess an N-terminal NAM domain which is highly conserved compared to the diverse C-terminal (Figs. S1A–S1C). Moreover, the N-terminal NAM domains of RsNAC proteins were grouped into five distinct clusters (A–E) which were similar to that in A. thaliana and C. sativus (Tran et al., 2004; Zhang et al., 2017). However, among the 172 RsNACs, only the RsNAC155 lack the conserved subdomains A, B and C, while RsNAC161 and RsNAC127 do not contain subdomain A and C, respectively. Additionally, several proteins (RsNAC172, RsNAC171, RsNAC098, RsNAC100, RsNAC101, RsNAC156 and RsNAC143) lack subdomain E (Fig. S1C).

For comparative genomic analyses, we obtained NAC protein coding sequences for 21 representative plant species (Fig. 1), resulting to 2931 NAC TFs. Interestingly, all the vascular plants had NAC genes while in lower plants, two algae representatives had none. The phenomenon suggests that the NAC proteins might have expanded after the evolutionary split of the vascular plants from the lower plant species, and further indicate that this family is vascular plant-specific based on the missing NACs in algae. In conclusion, based on this study, the number of NAC in radish ranked among the members with high NAC representation in their genome. According to the density of NACs in the respective genomes, radish (0.33) ranged among the species with a high number of NAC in our analyses (Fig. 1). It was also noted that the densities of NAC TFs in most of the higher plant genomes were greater than those in all the lower plant genomes, revealing that NAC TFs increased considerably with the evolution of higher plants.

Phylogenetic analysis of NACs in radish

To gain further insights into the evolutionary relationship among NAC members from different plants, the NACs from radish were aligned against those in Arabidopsis and rice, and to plot an unrooted phylogenetic tree using MEGA 6 software (Fig. 2, Table S3). The RsNACs could be broadly categorized into two groups (group A and group B, respectively) based on their preferential homology to Arabidopsis and rice. As a result, group A and B showed high homology to Arabidopsis and rice, respectively. Additionally, 9 and 10 subgroups were obtained from group A (A1–A9) and B (B1–B10) in the phylogenetic analysis of the full-length NAC proteins of Arabidopsis, rice and radish based on similarity in NAC domain structures. Phylogenetic analysis revealed that most RsNAC members were closer to eudicot (AtNACs) than monocot (OsNACs) and preferentially clustered in group A. For example based on outer clades with a bootstrap value of >50, 12 RsNACs clustered with 10 AtNACs in A4 while 23 OsNACs clustered with only one RsNAC member (RsNAC085) in subgroup B7 (Fig. 2, Table S3).

Motif and gene structure analysis of the RsNAC genes

To explore the relationships among the RsNAC genes and predict their potential functions, the RsNAC proteins clustered into 9 subgroups based on the phylogenetic tree (Figs. 3A and 3C), named I to VIII. The gene number of each group ranged from 12 (subgroup I) to 28 (subgroup IV). RsNAC proteins belonging to a common cluster in the phylogenetic group (Figs. 3A and 3C) indicated a common motif distribution (Figs. 3B and 3D), suggesting functional similarities among members within the same subset. Fifteen (1–15)conserved motifs were analyzed in the RsNAC proteins (Figs. 3B and 3D). Notably, the majority of the discovered motifs were located in the N-terminal region of the NAC domain which could be an indication that these motifs are critical for the functioning of NACs. For example, subgroup I mainly contained six motifs (motifs 1, 2, 4, 5, 8 and 9), whereas subgroup II comprised of seven motifs (motifs 1, 2, 3, 4, 5, 7 and 13), which suggests that the RsNACs with similar functions tended to cluster into the same subgroups. Additionally, the structural features of the intron/exon distribution among RsNAC genes were analyzed based on their phylogenetic similarities (Figs. 4A–4D, Table S2). Exon/intron analysis depicted that introns distribution in RsNAC genes is diverse and varied from 0 to 6. Additionally, a total of 164 (∼95%) RsNAC genes had introns, while 8 (∼5%) RsNACs contained no intron (Figs. 4C and 4D). Among them, one RsNAC (RsNAC111) had six introns, two sets of 13 RsNAC members had five or four introns each and 109 RsNAC genes had two introns each. Particularly, the RsNACs with a maximum number of introns were located in subgroup VI with seven members having five introns, while the RsNACs without intron were located in subgroup I. These findings indicated a minimum diversity in gene structure among RsNAC genes belonging to a common group.

Membrane-associated RsNAC subfamily

Fifteen RsNAC proteins (∼8.72%) were identified as membrane-associated RsNACs proteins (Fig. 5, Table 1). Classification of NTL proteins among the subgroups was discrete. Arabidopsis orthologues analysis indicated that 8 RsNAC NTLs were orthologs to previously characterized NTL from Arabidopsis (Table 1). The subgroups IV, VI and VIII contained three NTL proteins, while groups I and VII contained 2 NTLs, indicating that the functions of radish NTLs were complex.

Chromosomal distribution of RsNAC genes in the radish genome

A total of 98 (∼57%) RsNAC genes were mapped onto 9 chromosomes (R1–R9) (Fig. 6, Table S4). For distribution and density of the RsNAC genes, R5 contained the maximum frequency (∼22%) of RsNAC genes followed by R6 (∼21%), while R3 had the lowest frequencies of RsNAC genes ∼3% (Fig. 6A). Thirteen RsNAC gene clusters containing 26 tandemly duplicated genes were identified in all chromosomes except chromosome R4 (Fig. 6B). Most tandem duplications were located on chromosome R5 with three gene clusters containing two genes each. The phenomenon indicated that some of the RsNAC genes have one or more paralogs, which could have been generated by whole-genome duplication (WGD) in radish. In general, majority of the RsNAC were located in the distal regions of the nine radish chromosomes while a few were found in proximal regions. Additionally, the RsNAC gene density in each chromosome ranged from 0.10 /Mb (R3) to 0.45 /Mb (R5). This indicated that the RsNAC genes were significantly non-random and unevenly distributed in the nine radish chromosomes.

RNA-Seq-based expression analysis of RsNAC genes

For further investigation on the roles of RsNAC gene, the RNA-Seq based analysis was utilized to explore the expression patterns of the 172 RsNACs in roots and leaves at six different developmental stages (Fig. S2A–S2D, Table S5). As shown in Figs. S2A and S2C, there was no detectable expression of four RsNACs genes (RsNAC002, 131, 039 and 008). Moreover, tissue-specific expression patterns were observed in several NAC genes. These genes could be subdivided into three expression categories, which contained genes preferentially expressed in roots, leaves and both. For example, RsNAC002 (group I), RsNAC160 (group III), RsNAC165 (group IV), RsNAC069 (group V), RsNAC008 (group VI), RsNAC126 (group VII), and RsNAC124 (group VIII) were preferentially expressed in roots, while RsNAC097 and RsNAC017 (group VII) showed biased expression in leaves. However, some RsNACs exhibited multiple expressions in more than one organ and developmental stages. For example, RsNAC001, RsNAC027, RsNAC080 and RsNAC040 (group VI), RsNAC007 (group VII) and RsNAC028 (group VIII) were detected in roots and leaves at the six stages of development.

RsNAC genes expression profiling under various abiotic stresses and phytohormone treatment

To explore the expression patterns of RsNAC genes under various abiotic stresses and hormone treatments, 21 representative genes were identified for expression analysis under Cd, Pb, salt, PEG and ABA treatments. The results showed differential expression patterns of these genes under one or more treatments(Fig. 7, Table S6). Some genes (RsNAC027, RsNAC038 and RsNAC062) showed stable up-regulation under Cd, Pb and salt treatments. However, the other genes showed a weak expression level under Cd and Pb stress. Most RsNAC genes were preferentially down-regulated under HM stress whereas Pb stress up-regulated only a few of them, including RsNAC038, RsNAC008 and RsNAC062 (Fig. 7).

Characterization of RsNAC genes

NAC TFs are one of the largest plant-specific transcriptional regulators that play critical roles in numerous stress and developmental processes. However, there is still no detailed information available on NAC family in radish. To further explore the NAC information in radish, 172 putative RsNAC genes were identified from the radish genome. Remarkably, the number of the RsNAC members ranges among the largest NAC families. This finding is probably due to the selective pressure conferred by environmental cues to facilitate the growth and development of plants (Hughes & Friedman, 2003). In comparative genomic analysis among different plant species, N. tabacum was revealed harboring the highest numbers of NAC TFs (280). However, their NAC protein density (N. tabacum, 0.055) was lower than those in radish because of their large genome sizes. Additionally, stress responses of several RsNAC genes were detected based on their corresponding Arabidopsis homolog, which could be of benefit for further studies in functional characterization of NACs genes in root vegetables.

According to multiple sequence alignment, all the identified RsNACs contain the expected NAM domain (Pfam PF02365) with approximately 160 aa apportioned into five distinct subdomains (A–E). However, several RsNACs have no typical NAC domain pattern. For example, RsNAC155 lacks the conserved domain A, B, C and D, while other 12 RsNACs lacked one or two sub-domains. Such variant NAC members are termed as NAC-like proteins as described in potato and rice NACs (Nuruzzaman et al., 2010; Singh et al., 2013). The C-terminal region of the RsNACs domains indicate a highly distinct pattern involved in biotic and abiotic stress processes (Mao et al., 2012). Previous studies revealed the involvement of NAC members in protein-binding activity which performs major functions in transcription regulation (Jensen et al., 2010; Nuruzzaman et al., 2010). For phylogenetic analysis, monocots specific phylogenetic clade was identified with only one RsNAC member. The results suggest the existence of evolutionary differences between NACs from monocots and dicots and a lineage-specific expansion after splitting from their common ancestor. Additionally, the majority of the monocots and dicots with diverse functions belonged to common clades indicating that greater expansion of the NAC family could have happened slightly before the splitting of monocots and dicots (Pereira-Santana et al., 2015). Based on the results of comparative genomic analysis, it was possible to predict the functions of several RsNAC genes in response to abiotic stress according to their Arabidopsis homologs, which could also be potentially utilized for further functional studies. It was found that ANAC019, ANAC055 and ANAC07 which are orthologs to RsNAC168, RsNAC140 and RsNAC105, respectively, have been implicated in positively responding to dehydration in Arabidopsis (Hu et al., 2006; Tran et al., 2004).

In this study, all the identified RsNAC NTLs proteins contained single TM at the C-terminal of their conserved domain except RsNAC096 with two TMs, which is similar to two soybean NAC proteins (GmNAC136 and GmNAC013) where each of them contained a pair of TMs (Le et al., 2011). Recently, NTLs were demonstrated to have critical functions in response to multiple stresses (Kim et al., 2007; Lee et al., 2014). Moreover, NTLs facilitates an efficient mechanism for gene regulation which is a survival feature that allows for timely responses to environmental cues (Kim et al., 2007). In Arabidopsis, it has been reported that four NTLs (At4g01540, At2g27300, At4g35580 and At3g49530) are involved in abiotic stress tolerance, which are membrane-related proteases located in the endoplasmic reticulum that initiates the activities of TFs in response to stresses (Chen, Slabaugh & Brandizzi, 2008; Kim et al., 2008; Seo, Kim & Park, 2008). It was found that the radish NTLs (RsNAC034, RsNAC023, RsNAC058, RsNAC158, RsNAC126, RsNAC014, RsNAC111 and RsNAC124) showed close phylogenetic relationship with AT3G49530.1, AT2G27300.1, AT5G22290.1, AT1G34190.1, AT1G34190.1, AT5G04410.1, AT1G33060.1 and AT1G32870.1, respectively. In Arabidopsis, these NAC were also identified as NAC NTLs (Kim et al., 2007). Thus, it is expected that some identified RsNAC NTLs could function critically in the adjustment of stress-responsive genes as a regulatory mechanism for the survival of radish under adverse growth conditions.

Chromosomal distribution and expression profiling of RsNAC genes

The vast number of genes in NAC family in radish confirms the occurrence of chromosome duplication events during evolution. In angiosperms, several genome duplications have occurred resulting in gene family enlargement (Cannon et al., 2004), this phenomenon occurred in radish about 5.1–8.4 and 12.8–21.4 million years ago (Shen et al., 2013). Genome duplication provides a critical evolutionary promoter in angiosperm genome enlargement (Van de Peer et al., 2009). Additionally, most studies have indicated that WGD is a key factor in angiosperm biology. During genome duplication processes, redundant genes are lost (fractionation), preserved, or fixed to perform a diverse task compared with the original genes. Therefore, gene family size could vary significantly among species due to WGD events that took place during their evolution (Abrouk et al., 2010; Cenci et al., 2014). In several TFs, genome duplication has been reported in plants including MYB and WRKY, which is the same in NAC (Karanja et al., 2017; Liu et al., 2014; Nuruzzaman et al., 2010). In this study, eleven clusters representing probable tandem duplication formed by RsNAC genes in the genome were identified to elucidate the mechanism behind a large number of the NAC family members in radish. In Arabidopsis, many gene families have enlarged due to the retention of original genes after WGD and transposition (Freeling et al., 2008; Rizzon, Ponger & Gaut, 2006). Each of these gene family expansions forms paralogs that potentially lead to genetic redundancy in addition to differential growth among different families (Pereira-Santana et al., 2015). This phenomenon could contribute to the expansion of the RsNAC TFs in radish, similar to B. rapa and M. domestica (Liu et al., 2014; Su et al., 2013). Tandem duplication also indicate the presence of redundant genes which according to existing theory cause loss or mutation of duplicate genes despite retention of several genes after WGD by functionalization (Lynch & Force, 2000).

In the present study, the publicly available RNA-Seq data was utilized in conducting in silico expression analysis of RsNAC genes. Transcript abundance of NAC genes in roots and leaves reveal a spatial, temporal and ubiquitous pattern of RsNAC genes of which RsNAC168 was preferentially expressed in leaves, which was similar to its corresponding orthologs in Arabidopsis (ANAC019) expressed in leaves for regulating drought and salinity responses (Jensen et al., 2010). Additionally, ANAC019 and ANAC055 (orthologs to RsNAC168 and RsNAC140, respectively) are involved in both ABA and JA-mediated regulation (Jensen et al., 2010). Notably, most expressed genes were responsive in more than one tissue and developmental stage, spatio-temporal expression pattern of the RsNACs could provide the specific usage of RsNACs in a transcriptional modification in any given tissue, while ubiquitously regulated RsNACs could control the transcription of a wide range genes. In rice, OsNAC10 is widely expressed in roots and flowers tissues in response to dehydration and salinization stresses (Jeong et al., 2010).

Differential expression of RsNAC genes under various abiotic stresses

Unfavorable environmental cues such as HMs, salinization, high temperature and dehydration could result in permanent and deleterious effects on the optimal survival of plants. In this study, the response of RsNAC genes to all the stress treatments indicate their integral role in regulating survival processes associated with stress response and signaling in radish and other related species. According to the previous report, NAC TFs were reported to coordinate stress by initiating stress-responsive genes (Jensen et al., 2010). Several RsNACs (including RsNAC008, RsNAC027, RsNAC038) were significantly upregulated in response to salt. In Arabidopsis, overexpression of ANAC072, ANAC055 (RsNAC039, RsNAC140) and ANAC019 (RsNAC168) showed improved stress tolerance under salt conditions than the wild-type (Hu et al., 2006; Jensen et al., 2010). The finding could suggest that up-regulating of these genes may have resulted in improved salt tolerance in radish, however, NAC genes were found to negatively control drought responses by up-regulating the expression of integral drought- responsive genes(Wang et al., 2013b). Some genes, such as RsNAC030, RsNAC041 and RsNAC145, showed diverse responses under different stresses, indicating that these may function as links between the various transduction pathways (Liu et al., 2014).

In recent years, NTLs have been extensively proven to play critical roles in cell division as well as abiotic stress responses (Lee et al., 2014). In selective breeding, genes with multiple stress responses possess high potential of stress tolerance manipulation in plants. However, an occasional overexpression of NAC genes may result in antagonistic responses such as delayed flowering and dwarfing in transgenic plants (Liu et al., 2011). Overexpression of ANAC019 and ANAC055 in Arabidopsis resulted in increased drought tolerance but also increased susceptibility to B. cinerea (Bu et al., 2008). Field evaluation and prudent choice of a suitable promoter can reduce the negative effects in the over-expression of some NAC genes (Shao, Wang & Tang, 2015). In this study, nine NTLs (RsNAC23, RsNAC34, RsNAC39, RsNAC041, RsNAC058, RsNAC080, RsNAC096, RsNAC124 and RsNAC145) were among the 21 responsive RsNAC genes that responded to multiple abiotic stresses. Among them, RsNAC023 and RsNAC080 responded to all stress treatments other than ABA, while RsNAC145 responded to salt, heat and drought only. Evidently, functional characterization of these genes could be of great interest in unraveling their potential for the development of radish cultivars tolerant to multiple abiotic stresses. Over-expression techniques will provide valuable insights into understanding the regulatory functions of NAC members in radish in response to abiotic stresses. Significantly, NTLs responded differently to various stress treatments, RsNAC058 was upregulated under Cd, heat and drought, but downregulated under Pb and ABA it was. The result implied that the NTL genes might be playing critical roles in plant response to abiotic stresses.