Genome-Wide Identification of the NAC Gene Family in Zanthoxylum bungeanum and Their Transcriptional Responses to Drought Stress

NAC (NAM, ATAF1/2, and CUC2) transcription factors (TFs) are one of the largest plant-specific TF families and play a pivotal role in adaptation to abiotic stresses. The genome-wide analysis of NAC TFs is still absent in Zanthoxylum bungeanum. Here, 109 ZbNAC proteins were identified from the Z. bungeanum genome and were classified into four groups with Arabidopsis NAC proteins. The 109 ZbNAC genes were unevenly distributed on 46 chromosomes and included 4 tandem duplication events and 17 segmental duplication events. Synteny analysis of six species pairs revealed the closely phylogenetic relationship between Z. bungeanum and C. sinensis. Twenty-four types of cis-elements were identified in the ZbNAC promoters and were classified into three types: abiotic stress, plant growth and development, and response to phytohormones. Co-expression network analysis of the ZbNACs revealed 10 hub genes, and their expression levels were validated by real-time quantitative polymerase chain reaction (qRT-PCR). Finally, ZbNAC007, ZbNAC018, ZbNAC047, ZbNAC072, and ZbNAC079 were considered the pivotal NAC genes for drought tolerance in Z. bungeanum. This study represented the first genome-wide analysis of the NAC family in Z. bungeanum, improving our understanding of NAC proteins and providing useful information for molecular breeding of Z. bungeanum.


Introduction
Transcription factors (TFs) are proteins that regulate gene expression by binding to specific cis-acting promoter regions of target genes [1]. The NAC (NAM, ATAF1/2, and CUC2) family is one of the largest plant-specific TF families; its name is derived from the three super families-no apical meristem (NAM), Arabidopsis transcription activation factor (ATAF), and cup-shaped cotyledon (CUC) [2]. NAC proteins contain a highly conserved N-terminal region (NAC domain) and a variable C-terminal region. The Nterminal domain, which consists of 150 amino acids, is further divided into 5 subdomains with diverse functions: the A subdomain is related to functional dimerization, the B and E subdomains are responsible for protein functional diversification, and the C and D subdomains are associated with DNA binding [2,3]. The variable C-terminal region of NAC TFs may also participate in transcriptional regulation through the transcriptional activation/repression region (TAR or TRR) [3]. listed in Table S1. The coding sequences (CDSs) ranged from 579 to 2307 bp, and the predicted proteins varied from 192 to 768 amino acids in length, with molecular weights between 48.18 and 193.36 kDa. Their theoretical isoelectric points (pIs) varied from 4.90 to 5.19, with an average of 5.07, indicating that most ZbNAC proteins were weakly acidic. Subcellular location prediction suggested that most ZbNAC proteins were located in the nucleus (81), some in the plasma membrane (31) and chloroplast (15), and a few in the cytoplasm (5), peroxisome (4), and Golgi apparatus (1). These subcellular location results suggest that ZbNAC TFs may participate in the regulation of nuclear genes and have multiple functions in adaptation to diverse environments.  Duplication events were analyzed in the 109 ZbNAC genes (Table 1, Figure 2). Five pairs of tandem duplicates were distributed on two chromosomes, Chr 1 and Chr 8. Interestingly, four pairs of tandem duplications were located on Chr 8 and consisted of five adjacent ZbNAC genes, indicating that these five ZbNAC genes probably share a common ancestor. In addition, 42 pairs of segmental duplications were identified among the 109 ZbNAC genes. These results suggest that segmental duplication may have been the main driving force during ZbNAC family evolution. Ka and Ks can be used to measure the course of To further explore the evolutionary relationships between Z. bungeanum and other plants, we constructed six collinearity maps of NAC genes for Z. bungeanum, A. thaliana, Camellia sinensis, Ziziphus jujube, Oryza sativa, Triticum aestivum, and Zea mays ( Figure 3). There were 41, 72, and 39 orthologous gene pairs between Z. bungeanum and A. thaliana, C. sinensis, and Z. jujube. However, there were fewer orthologous pairs between Z. bungeanum and O. sativa (14), T. aestivum (5), and Z. mays (8). To further explore the evolutionary relationships between Z. bungeanum and other plants, we constructed six collinearity maps of NAC genes for Z. bungeanum, A. thaliana, Camellia sinensis, Ziziphus jujube, Oryza sativa, Triticum aestivum, and Zea mays ( Figure 3). There were 41, 72, and 39 orthologous gene pairs between Z. bungeanum and A. thaliana, C. sinensis, and Z. jujube. However, there were fewer orthologous pairs between Z. bungeanum and O. sativa (14), T. aestivum (5), and Z. mays (8).

Phylogenetic Analysis of ZbNAC Proteins
To explore the evolutionary relationships among NAC proteins, a phylogenetic reconstruction was performed based on a multiple sequence alignment using MEGA 7.0 ( Figure 4). An unrooted phylogenetic tree was constructed with 109 ZbNAC TFs and 88 AtNAC TFs based on the conserved NAC domain, and the proteins were clustered into four distinct groups (I-IV). Group I contains the most AtNAC members (62) and the fewest ZbNAC members (6), suggesting that there may be a divergence from the last common ancestor in this group. There were fewer NAC members in A. thaliana than in Z. bungeanum in groups II-IV: group II contains 9 AtNAC proteins and 24 ZbNAC proteins, group III contains 9 AtNAC proteins and 38 ZbNAC proteins, and group IV contains 8 AtNAC proteins and 41 ZbNAC proteins. ZbNAC and AtNAC proteins with close phylogenetic relationships may have similar gene structures and biological functions.

Phylogenetic Analysis of ZbNAC Proteins
To explore the evolutionary relationships among NAC proteins, a phylogenetic reconstruction was performed based on a multiple sequence alignment using MEGA 7.0 ( Figure 4). An unrooted phylogenetic tree was constructed with 109 ZbNAC TFs and 88 AtNAC TFs based on the conserved NAC domain, and the proteins were clustered into four distinct groups (I-IV). Group I contains the most AtNAC members (62) and the fewest ZbNAC members (6), suggesting that there may be a divergence from the last common ancestor in this group. There were fewer NAC members in A. thaliana than in Z. bungeanum in groups II-IV: group II contains 9 AtNAC proteins and 24 ZbNAC proteins, group III contains 9 AtNAC proteins and 38 ZbNAC proteins, and group IV contains 8 AtNAC proteins and 41 ZbNAC proteins. ZbNAC and AtNAC proteins with close phylogenetic relationships may have similar gene structures and biological functions.

Gene Structures and Conserved Motifs in ZbNAC Proteins
To further explore the relationships among ZbNAC family members, the phylogeny, gene structures, and conserved motifs of the ZbNAC TFs were analyzed ( Figure 5). Eight conserved motifs were identified in the 109 ZbNAC proteins, and their amino acid lengths varied from 15 to 50 ( Figure 5B,D). The motif numbers in the ZbNAC proteins varied from 2 (ZbNAC104) to 9 (ZbNAC080). Motifs 1-6 were found in the NAC domain regions of 79 ZbNAC proteins. Motif 5 may be the most important motif in the ZbNAC TFs, as it was present in almost all proteins except for ZbNAC109, ZbNAC054, ZbNAC071, and ZbNAC095. Conversely, Motif 8 was least common and was present in only 10 ZbNAC proteins. Although motif types and numbers differed among ZbNACs from different groups, most members with a close phylogenetic distance shared the same motifs and similar motif locations.  Gene structure analysis showed that the ZbNAC encoding genes contained 1 to 8 introns ( Figure 5C). Most of the ZbNAC genes (67) contained two introns, two (ZbNAC073, ZbNAC047) contained one intron, and one (ZbNAC080) contained eight introns. As expected, genes that were closely phylogenetically related shared common motifs and similar motif locations. For example, ZbNAC068, ZbNAC069, and ZbNAC070 are clustered together in the phylogeny tree; all have five exons with similar locations and lengths.

Cis-element Analysis of ZbNAC Genes
Twenty-three types of cis-elements were distributed unevenly in the 2000-bp upstream regions of the 109 ZbNAC genes ( Figure 6 and Figure S1). These cis-elements were classified into three categories: abiotic stress (6 types), plant growth and development (7 types), and response to phytohormones (10 types). Judging from the numbers of these elements, they are mainly involved in phytohormone pathways and responses to abiotic stress, thus revealing the important role of ZbNACs in the regulation of hormone synthesis and adaptations to environmental stress. With respect to abiotic stress, TC-rich repeats (75) were the most abundant elements; they were located in 45 ZbNAC promoters and were annotated as drought response elements ( Figure S2). With respect to hormone response, the cis-elements were mainly involved in abscisic acid, auxin, gibberellin, methyl jasmonate (MeJA), and salicylic acid pathways. Among these elements, ABRE elements were the most abundant (369, abscisic acid responsive) followed by CGTCA motifs (206, MeJA responsive) and TGACG motifs (206, MeJA responsive). The numbers and types of cis-elements were similar in ZbNAC genes that had a close phylogenetic relationship. In group I, there were fewer elements, and they were seldom assigned to the plant growth and development category. In group II, the number of cis-elements was highest, and TC-rich repeats, ABRE elements, CGTCA motifs, and TGACG motifs were the most abundant. For example, there were 18 ABRE elements in ZbNAC092, 14 ABRE elements in ZbNAC065 and ZbNAC046, and 8 CGTCA motifs in ZbNAC004 and ZbNAC102, and these elements were annotated as ABA and MeJA responsive. The ZbNAC TFs in group II may therefore be important for drought response and ABA and MeJA regulation.

RNA-seq Analysis of 109 ZbNAC Genes under Drought Stress
A hierarchical clustering heatmap was used to visualize the expression patterns of 109 ZbNAC genes of two Z. bungeanum cultivars under drought stress. These genes were divided into four clusters (clusters 1-4) based on their FPKM values ( Figure 7A,B). In cluster 1, the total expression level of 30 ZbNACs decreased with progressive drought stress in both FJ and HJ, and the F4/F1 and H4/H1 ratios were less than one, with the exception of ZbNAC085 (3.59) and ZbNAC102 (1.97) in HJ. Thirty-three ZbNACs were placed into cluster 2, and most of them increased in expression with progressive drought stress, with higher expression in HJ than in FJ. The F4/F1 and H4/H1 ratios were greater than one in almost all members of cluster 2 except for ZbNAC081 (0.95) and ZbNAC094 (0.84). Interestingly, the ratios of six ZbNACs from cluster 2 (ZbNAC007, ZbNAC046, ZbNAC047, ZbNAC056, ZbNAC065, and ZbNAC079) were greater than ten in both cultivars, suggesting that these six genes may play pivotal roles in response to drought stress. There were 19 ZbNACs in cluster 3; most were significantly enhanced in FJ during early drought, but the changes in HJ were less marked, indicating that the genes in cluster 3 were specifically involved in adaptation to drought stress in FJ. In cluster 4, the expression levels of most genes increased in the early drought period and decreased in the later period, suggesting that these TFs may regulate the expression of early drought-responsive genes in Z. bungeanum. In addition, the average FPKM values suggested that the expression levels of ZbNAC057 (99.52) and ZbNAC090 (44.15) were significantly higher than those of other genes, which averaged 4.30. The total expression levels of the 109 ZbNAC genes were displayed in violin diagrams ( Figure 7C). Most genes had a low expression level under control conditions and were significantly induced under drought stress. Taken together, the 87 ZbNAC genes in clusters 2, 3, and 4 may positively regulate the drought tolerance of Z. bungeanum and were therefore selected for further study.

Comprehensive Analysis of Drought-Related ZbNAC Genes
Correlation analysis was performed for 87 genes in clusters 1, 2, and 3 with FPKM values greater than 1 to explore their possible relationships ( Figure 8A). The results showed that 742 pairs of genes showed positive correlations (p < 0.05), and only 12 pairs showed negative correlations (p < 0.05). GO analysis assigned GO terms to these genes from three categories: biological process (BP), cellular component (CC), and molecular function (MF) ( Figure 8B). Among BP terms, metabolic process, cellular process, and biological regulation were most common; among CC terms, cell, cell part, and organelle were most common; only one MF term, binding, was assigned to the selected genes. A coexpression network was then constructed based on Pearson correlation coefficients. Ten hub genes (ZbNAC007, ZbNAC017, ZbNAC018, ZbNAC043, ZbNAC047, ZbNAC072, ZbNAC079, ZbNAC087, ZbNAC100, and ZbNAC108) were identified using the Cytohubba program in Cytoscape 3.9.0 ( Figure 8C). Three hub genes (ZbNAC087, ZbNAC072, and ZbNAC100) were the most dominant in the network because of their high connectivity with other genes.  The AtNAC homologs of the ten hub ZbNAC proteins were predicted from the phylogenetic tree (Figure 1) and are listed in Table S2. To further examine the functions of the ten hub genes, a protein-protein interaction network of ten hub ZbNAC proteins was predicted based on their AtNAC homologs using STRING (Figure 9). NAC094 (homolog of ZbNAC018) was directly connected to MC5, with a score of 0.909. RD26 (homolog of ZbNAC047 and ZbNAC079) was directly connected to ZFHD1 (0.860) and HAI1 (0.804). NAP (homolog of ZbNAC007 and ZbNAC072) was directly connected to SAG12 (senescence-associated gene 12 encoding protein), SAG13 (senescence-associated gene 13 encoding protein), and HAI1 (highly ABA-induced 1). Notably, HAI1 was involved in a cluster network that also contained genes encoding sucrose non-fermenting-related kinase 2 (SnRK2), PYR, ABA receptors (PYLs), regulatory components of ABA receptor (RCAR), and open stomata 1 (OST). GO analysis annotated the genes in the predicted network with 32 GO terms, including abscisic acid binding, signaling receptor activity, abscisic acidactivated signaling pathway, response to abscisic acid, and signal transduction (Table S3). KEGG analysis showed that 14 genes were annotated with the pathways mitogen-activated protein kinase (MAPK) signaling pathway-plant and plant hormone signal transduction (Table S4). Moreover, qRT-PCR showed that the relative expression levels of the ten hub genes increased under drought stress (Figure 10), and the expression data from qRT-PCR were consistent with the RNA-seq data.

Discussion
NAC TFs are one of the largest plant-specific TF families and play significant roles in plant growth, development, and adaptation to environmental stress [25]. Genome-wide identification of NAC family members has been performed in many plant species but has not yet been available in Z. bungeanum owing to the absence of a reference genome. Recently, a chromosome-scale assembly of Z. bungeanum genome with high quality was published in our research [24], which provided the valuable opportunities for the further

Discussion
NAC TFs are one of the largest plant-specific TF families and play significant roles in plant growth, development, and adaptation to environmental stress [25]. Genome-wide identification of NAC family members has been performed in many plant species but has not yet been available in Z. bungeanum owing to the absence of a reference genome. Recently, a chromosome-scale assembly of Z. bungeanum genome with high quality was published in our research [24], which provided the valuable opportunities for the further study of Z. bungeanum. Here, we firstly performed the comprehensive genome-wide identification of NAC family members and characterized their gene expression under drought stress in Z. bungeanum.
Based on the presence of the NAC domain structure, 109 non-redundant ZbNAC TFs were identified, close to the 108 NACs in Dactylis glomerata [26]; more than the 82 in Cucumis melo [27], 80 in Fagopyrum tataricum [28], and 96 in Manihot esculenta [29]; and fewer than the 151 in O. sativa [30], 152 in Z. mays [31], and 152 in Glycine max [32]. These results suggest that the evolutionary expansion of the NAC family differed among species. There are a large number of chromosomes (n = 68) in the Z. bungeanum genome, owing to frequent fusion/fission events after whole-genome duplication [24]. The 109 NAC genes are distributed unevenly on most of the chromosomes (46/68), and this may be closely related to their functional diversity.
Gene duplication is a potential driving force for species evolution as well as gene functional diversity [33]. NAC TFs in 160 plant species are reported to have evolved from common ancestors and to have undergone numerous duplication events [34]. Segmental and tandem duplications are considered the main types of duplication events in eukaryotes [35]. In general, more segmental duplications than tandem duplications have been found in the NAC family, such as the 84 segmentally duplicated pairs and 5 tandemly duplicated pairs in Panicum miliaceum [36], the 17 segmentally duplicated pairs and 2 tandemly duplicated pairs in Musa acuminata [37], and the 116 segmentally duplicated pairs and 1 tandemly duplicated pair in peanut [23]. In Z. bungeanum, 42 pairs of segmental duplicates and 5 pairs of tandem duplicates were identified, indicating that segmental duplication was more common during evolution of the ZbNACs. Interestingly, five tandem duplicates were found together on chromosome 2. We speculate that these five genes may have a common ancestor and may therefore share similarities in structure and function. Ka/Ks = 1, <1, and >1 are indicative of neutral selection, purifying selection, and accelerated evolution under positive selection, respectively [38]. In this research, most duplicates showed Ka/Ks <1, indicating that most ZbNACs evolved primarily under purifying selection. C. sinensis, Z. bungeanum, Z. jujube, and Arabidopsis are dicots, whereas maize, rice, and wheat are monocots. Here, the dicots exhibited greater collinearity with Z. bungeanum than the monocots, indicating that NAC TFs experienced extensive evolution and duplication after the divergence of monocots and dicots. Feng et al. [24] found that Z. bungeanum shared a close phylogenetic relationship with C. sinensis and that they diverged from their most recent common ancestor approximately 35.3 million years ago. Consistent with this finding, the greatest number of collinear NAC genes were found between these two species.
Exons and introns are important features of gene evolution, and exons are an important part of protein-coding genes [39]. The insertion or deletion of introns can result in a loss of function or gene mutation. In addition, genes with fewer introns may respond more quickly, and genes with large-scale introns possibly have specific functions [32]. The number of introns in the ZbNAC genes varied from 1 to 8, similar to numbers reported in chickpea (1-6) [40] and sunflower (2-7) [41]. Family members with close phylogenetic relationships generally have similar gene structures [42], and this phenomenon was also found in our research. For instance, ZbNAC068, ZbNAC069, and ZbNAC070, which were most closely related in the same group, harbored five exons of similar location and length. The two ZbNAC genes with only one intron (ZbNAC073 and ZbNAC047) may have a faster response speed and may be located upstream in the response pathway. ZbNAC080 contained eight introns, and the length of each exon was relatively short. Thus, we suspect that the insertion of introns led to mutation in gene structure and that ZbNAC080 may have specific functions. The motifs in TF proteins are often associated with DNA binding, protein interaction, and transcriptional activity [28,43]. A total of 10 conserved motifs were identified in the 109 ZbNAC TFs. Motifs 1-5 were present in almost all NAC members in groups 1-3, indicating that they probably have significant biological functions. Motif 5 was the most common, suggesting that it may be located in a conserved domain of the ZbNAC proteins. All the conserved motifs were located in the N termini of the protein sequences, indicating that their existence is essential to NAC function [2].
The structure and regulatory mechanisms of gene promoters are closely related to specific plant traits [44]. The physiological functions of promoter cis-elements indicate possible responses of these genes to internal and external factors [44]. Twenty-four types of promoter elements were identified in ZbNACs and were divided into three categories: abiotic stress response elements, hormone response elements, and developmental regulation elements. Numerous elements were involved in abiotic stress and hormone response, consistent with NAC genes in jujube [20], implying that most ZbNACs may be associated with plant hormone signaling and abiotic stress responses. LTR (response to low temperature) and CCTG (response to drought stress) accounted for a large proportion of abiotic stress elements, suggesting a vital role for ZbNACs in low temperature and drought stress. Consistent with this notion, substantial evidence suggests that the NAC gene family plays a critical role in plant responses to abiotic stress. For example, overexpression of ONAC022 improved drought and salt tolerance in rice [45], and overexpression of MbNAC25 from Malus baccata increased the resistance of transgenic Arabidopsis to cold and salinity [46]. There were large numbers of ABA elements and MeJA elements in the ZbNAC promoters. Likewise, ABA and MeJA elements were also abundant in NAC promoters of Ipomoea trifida [12]. ABA and MeJA are well known as important phytohormones for plant response and adaptation to various abiotic stresses, including salinity, drought, and cold [47]. Previous research has shown that ABA and MeJA are easily induced by drought stress [48], and exogenous application of ABA and MeJA can enhance plant drought tolerance [49][50][51]. The results of this study suggest that most ZbNACs may participate in drought resistance via regulation of ABA and MeJA.
Z. bungeanum is a commercial crop that is widely grown in arid and semi-arid regions, and the above analyses revealed the possible function of ZbNACs in drought response. Therefore, exploring the expression of these ZbNACs under drought is necessary for a comprehensive analysis of the ZbNAC family. In the present research, RNA-seq was performed with two Z. bungeanum cultivars of different drought tolerance under progressive drought stress. The 109 ZbNAC genes were divided into four groups based on their expression levels, and three clusters (79 genes, 70%) were upregulated by drought, suggesting that most ZbNAC genes are involved in drought tolerance of Z. bungeanum. Likewise, drought also increased the expression levels of most NAC genes in peanut [52] and pepper [53]. Correlation analysis showed that the majority of upregulated genes were significantly positively correlated, suggesting that they may play a synergistic role in drought tolerance. GO analysis showed that their encoded proteins were involved in binding, indicating that NAC genes may regulate the drought response by binding to target genes. Co-expression analysis has been widely used to screen key genes from large gene sets [54][55][56]. Here, ten hub genes with a high degree of connectivity to other genes were identified in the constructed co-expression network. We speculate that these ten ZbNACs are pivotal genes that dominate the interactions of transcription factors, and they may have important regulatory functions in drought response. The qRT-PCR results showed that the expression levels of these ten ZbNACs increased with drought, and most of their expression levels were higher in HJ than in FJ. This result not only verified our speculation but also confirmed the reliability of our RNA-seq data.
Plant drought response mechanisms are complex and include ROS scavenging, stomatal closure, root elongation, production of antioxidant enzymes, and accumulation of antioxidant substances [57]. MAPK and ABA signaling pathways play an important role in plant responses to abiotic stress, including drought stress and salt stress [58]. Under drought stress, plants perceive the drought signal by the ABA signaling pathway, which consists of PYLs, Clade A protein phosphatase 2Cs (PP2Cs), and SnRK2 [59]. PYL-PP2C signaling can activate MAPK signaling cascades to respond to drought stress [60]. Our protein-protein interaction network predicted that NAP (homolog of ZbNAC007 and ZbNAC072) and RD26 (homolog of ZbNAC047 and ZbNAC079) interacted with HAI1, a member of the PP2C family. HAI1 could interact with SNRK2.2, SNRK 2.3, PYL4, RCAR3, and OST1. In addition, KEGG enrichment analysis showed that the target genes of these NAC TFs were enriched in the MAPK and ABA signaling transduction pathways. These results indicate that NAP and RD26 may regulate drought tolerance by affecting the MAPK and ABA signal transduction pathways. Consistent with this notion, previous reports indicate that NAP and RD26 mediate drought stress-responsive signaling in Arabidopsis [61,62].
The protein-protein interaction network also showed that NAC094 (ZbNAC018 homolog) could interact with MC5, which is associated with programmed cell death (PCD) under abiotic stress [63]. PCD has been reported to promote the development of lateral roots with enhanced drought tolerance, which are important for post-stress recovery [64]. Thus, the NAC094 protein may also be important for drought tolerance. Because of the conservation of NACs, NAC proteins encoded by homologous genes usually share common biological functions [11,21]. Consequently, the corresponding homologous genes ZbNAC007, ZbNAC047, and ZbNAC072, and ZbNAC079 may participate in the regulation of MAPK and ABA signaling transduction pathways, and ZbNAC018 may be involved in PCD in Z. bungeanum under drought. With regard to subcellular localization, these five ZbNAC proteins were all located in the nucleus, further emphasizing their importance in regulating the expression of target genes under drought stress. Taken together, ZbNAC007, ZbNAC018, ZbNAC047, ZbNAC072, and ZbNAC079 are potential candidate genes for improving the drought resistance of Z. bungeanum, but further verification research needs to be performed in the future.

Identification of NAC Genes in Z. bungeanum
The Z. bungeanum genome resource was published in our previous research (BioProject ID PRJNA524242) [24]. The hidden Markov model (HMM) file of the NAM domain (PF02365) was obtained from the Pfam database (http://pfam.sanger.ac.uk/, accessed on 12 February 2022). HMMER 3.0 was used to identify the NAC proteins in the Z. bungeanum genome with a cut-off value of 0.01 [65]. Redundant NAC proteins were deleted using an online tool (https://web.expasy.org/decrease_redundancy, accessed on 12 February 2022). Subsequently, the putative NAC proteins were further confirmed by BLASTP and the Conserved Domain Database (CDD, http://www.ncbi.nlm.nih.gov/cdd/, accessed on 12 February 2022). Protein physicochemical parameters, including number of amino acids, molecular weight, and pI value, were predicted using ExPasy website tools (http://web. expasy.org/protparam/, accessed on 12 February 2022). The subcellular localization of the NAC proteins was predicted using WoLF PSORT (https://wolfpsort.hgc.jp/, accessed on 12 February 2022) [66].

Chromosomal Distribution and Gene Duplication
The chromosomal positions of the ZbNAC genes were extracted from the GFF file of the Z. bungeanum genome. The chromosomal map with gene positions was produced using TBtools (v. 1.098696) [67]. The ZbNAC genes were named according to their chromosomal positions, ZbNAC001-ZbNAC109. MCScanX software with default parameters was used to analyze gene duplication events among 109 ZbNAC genes [68]. The genome file for rice was downloaded from Phytozome (http://phytozome.jgi.doe.gov/pz/portal.html#, accessed on 15 February 2022) [69], and the Arabidopsis genome file was obtained from the Arabidopsis Information Resource (https://www.arabidopsis.org/, accessed on 15 February 2022). The genome sequences and annotation files of wheat, maize, and jujube were acquired from NCBI (https://www.ncbi.nlm.nih.gov/, accessed on 16 February 2022). The synteny analysis between Z. bungeanum and other species was performed using TBtools (v. 1.098696, JAVA, China), and with parameters CPU for BLATP: 4, evalue: 1× 10 -10 , Num of BlastHits: 5.

Sequence Alignment and Phylogenetic Analysis of ZbNAC TFs
The NAC protein sequences from Z. bungeanum and Arabidopsis were aligned using ClustalW in MEGA 7.0, and an unrooted phylogenetic tree was constructed using the neighbor-joining (NJ) method in MEGA 7.0 with pairwise deletion and 1000 bootstrap replicates. The resulting phylogenetic tree was further processed with the online tool iTOL (https://itol.embl.de/, accessed on 18 February 2022).

Gene Structure, Conserved Motif, and Cis-element Analyses of the ZbNAC TFs
The online Gene Structure Display Server 2.0 (GSDS, http://gsds.gao-lab.org/, accessed on 18 February 2022) was used to determine exon-intron structures by comparing CDSs and genomic sequences of the ZbNACs. Conserved motifs were identified in the ZbNAC proteins using the online program MEME (http://meme-suite.org/tools/meme, accessed on 18 February 2022) with a maximum of 10 motifs and a range of motif widths from 6 to 200. The integrative analysis of gene structure and conserved motifs was visualized using TBtools [67].
The 2000-bp sequences upstream of the ZbNAC start codons were extracted and submitted to PlantCARE (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/, accessed on 22 February, 2022) for cis-element analysis. The PlantCARE output file was used to illustrate the cis-element distribution in the ZbNAC promoters using TBtools.

Plant Materials and RNA-seq
Mature seeds of FJ (Z. bungeanum cv. 'Fengjiao') and HJ (Z. bungeanum cv. 'Hanjiao') were collected from the Z. bungeanum Experimental Station of Northwest Agriculture and Forestry University in Fengxian, Shaanxi province, China (33 • 59 6.55 N, 106 • 39 29.38 E). The seeds were sown in a mixture of perlite, vermiculite, and chernozem, and the resulting seedlings were cultivated at 25 ± 2 • C and 60-70% humidity in an experimental greenhouse at Northwest Agriculture and Forestry University in Yangling, Shaanxi province, China. Three months after germination, healthy seedlings of similar size were selected for drought treatment. After stopping the water supply, leaves were sampled at 0 d (D1), 6 d (D2),

Conclusions
In this research, 109 ZbNAC TFs were identified in Z. bungeanum and were classified into four subgroups with NAC proteins from Arabidopsis. These NAC TF genes were unevenly distributed on 46 chromosomes. Conserved motifs and gene structures showed similarity in NAC TFs clustering together in phylogeny tree. The inter-species synteny analysis of NAC genes revealed the close phylogeny relationship between Z. bungeanum and C. sinensis. Meanwhile, promoter cis-element analysis implied the multiple functions of NAC in abiotic stresses. Moreover, the expression levels were significantly stimulated in most of ZbNACs in both Z. bungeanum cultivars under progressive drought stress, therein five ZbNAC genes (ZbNAC007, ZbNAC018, ZbNAC047, ZbNAC072, and ZbNAC079) were potential candidate genes for drought tolerance in Z. bungeanum. In general, the genomewide analysis of NAC family was first performed in Z. bungeanum and the expression pattern of ZbNAC genes was carefully explored under drought in this study. The results improved our understanding of NAC family in plants. Additionally, these candidate ZbNAC genes provide pivotal clues for molecular breeding of Z. bungeanum in the future.