Genome-Wide Identiﬁcation and Expression Analysis of the NAC Gene Family in Kandelia obovata , a Typical Mangrove Plant

: The NAC ( NAM , ATAF1/2 , and CUC2 ) gene family, one of the largest transcription factor families in plants, acts as positive or negative regulators in plant response and adaption to various environmental stresses, including cold stress. Multiple reports on the functional characterization of NAC genes in Arabidopsis thaliana and other plants are available. However, the function of the NAC genes in the typical woody mangrove ( Kandelia obovata ) remains poorly understood. Here, a comprehensive analysis of NAC genes in K. obovata was performed with a pluri-disciplinary approach including bioinformatic and molecular analyses. We retrieved a contracted NAC family with 68 genes from the K. obovata genome, which were unevenly distributed in the chromosomes and classiﬁed into ten classes. These Ko NAC genes were differentially and preferentially expressed in different organs, among which, twelve up-regulated and one down-regulated Ko NAC genes were identiﬁed. Several stress-related cis-regulatory elements, such as LTR (low-temperature response), STRE (stress response element), ABRE (abscisic acid response element), and WUN (wound-responsive element), were identiﬁed in the promoter regions of these 13 Ko NAC genes. The expression patterns of ﬁve selected Ko NAC genes ( Ko NAC6, Ko NAC15, Ko NAC20, Ko NAC38, and Ko NAC51) were conﬁrmed by qRT-PCR under cold treatment. These results strongly implied the putative important roles of Ko NAC genes in response to chilling and other stresses. Collectively, our ﬁndings provide valuable information for further investigations on the function of Ko NAC genes.


Introduction
Transcription factors (TFs) are of immense importance due to their crucial impact on controlling the transcription rate by binding to the cis-regulatory elements, resulting in activation or inhibition of the transcription level of target genes [1]. There are numerous types of TF families in plants, among which the NAC (NAM, ATAF1/2, and CUC2) family serves as one of the largest plant-specific TF families and is named after the Petunia hybrida E. Vilm. NO APICAL MERISTEM (NAM) [2] and Arabidopsis thaliana (L.) Heynh. genes ATAF1/2 and CUP-SHAPED COTYLEDON 2 (CUC2) [3]. A typical NAC protein contains an N-terminal conserved NAC domain for DNA binding and nuclear localization and a variable C-terminal region with transcriptional regulatory activity [4].

Expression Analysis of KoNAC Genes Based on Public RNA-Seq Data
Two previously released RNA-seq data sets of K. obovata were introduced here to analyze the expression profiles of KoNAC genes. The expression patterns of KoNAC genes in eight organs (root, stem, leaf, flower, pistil, stamen, sepal, and fruit) were obtained according to the previously published transcriptomic data under the NCBI BioProject accession number PRJNA416402 (https://www.ncbi.nlm.nih.gov/bioproject) (accessed on 31 March 2022) [31]. The expression levels of KoNAC genes in response to cold stress were determined based on the publicly released data from the NCBI BioProject under accession number PRJNA678025. These two RNA-seq data were remapped back to the K. obovata genome used here [35]. All expression data were normalized as fragments per kilobase of transcript per million fragments mapped (FPKM) values [48]. The differentially expressed genes (DEGs) related to chilling stress were defined under the criteria of fold change (FC) ≥ 1.5. The expression profiles of KoNAC genes were visualized as heatmaps using TBtools [47].

Plant Materials and Treatment
The healthy mature propagules of the typical viviparous mangrove plant K. obovata were sampled from Guangxi Maoweihai Mangrove Nature Reserve, Qinzhou, China (21 • 37 23 N, 108 • 44 13 E) and cultured in the Mangrove Germplasm Resources Center (MGRC) of Guangxi Forestry Research Institute (GFRI) ( Figure S1). The seedlings were grown in plastic pots containing sand and cultivated in a growth chamber at 28 • C and 75% humidity with a photoperiod of 14 h light/10 h darkness, and watered weekly with half-strength Hoagland's nutrient solution [49]. At the eight-leaf stage, the seedlings were treated under low temperature (4 • C) for 0 h, 6 h, 12 h, and 24 h, respectively. All treatments were performed with three replicates. The leaves were harvested, immediately frozen in liquid nitrogen, and stored at −80 • C for RNA extraction.

Cis-Regulatory Element Analysis of the KoNAC Genes
The upstream 1500 bp promoter sequences from the ATG start codon of the KoNAC genes were retrieved from the K. obovata genome, and the cis-regulatory elements in the promoter regions were predicted using Plant CARE (https://bioinformatics.psb.ugent.be/ webtools/plantcare/html/) (accessed on 31 March 2022) [50] and displayed by TBtools [47].

Genome-Wide Identification of the K. obovata NAC Genes
Two independent strategies for retrieval of KoNAC genes from the K. obovata genome, HMM search and BLASTP, were used here. Taken together, 68 putative KoNAC genes were identified and confirmed by PfamScan and NCBI-CDD. Based on their chromosome location, these KoNAC genes were named KoNAC1 to KoNAC68 and unevenly distributed on 17 chromosomes (Chrs), with no KoNAC gene present on Chr18 (Table 1, Figure 1). Detail-wise, nine KoNAC genes were located on Chr08, six genes were located on both Chr03 and Chr12, and five KoNAC genes each were located on Chr02, Chr04, Chr06, Chr09, and Chr17, while only one KoNAC gene each was found on Chr07, Chr14, and Chr16.    Moreover, every KoNAC gene contained one or more introns with an average length of 371 bp, while the proteins encoded by KoNAC genes ranged from 162 amino acid (aa) residues (KoNAC61) to 638 aa (KoNAC60) in length, with an average length of 344 aa. The pI values varied from 4.20 (KoNAC45) to 10.07 (KoNAC11), over half of the members (39/68) exhibiting pI > 7 (Table 1).

Phylogenetic Analysis and Classification of KoNAC Proteins
To illustrate the phylogenetic relationship among K. obovata and Arbidopsis thaliana NAC proteins, a neighbor-joining phylogenetic tree was constructed with 68 KoNAC proteins and 105 AtNAC proteins (Table S1). The result showed that the 173 NAC proteins could be classified into ten classes, namely, Class I to Class X ( Figure 2). Obviously, Class VII, with 21 KoNACs and 24 ANACs was the largest class, followed by Class X with 12 KoNACs and 14 ANACs. Other classes contained no more than 10 KoNACs each. Specially, no KoNAC belonged to Class I (Figure 3).  To better understand the phylogenetic relationship and classification of KoNAC genes, the gene structure and motif organization of the 68 KoNAC genes were analyzed. Each KoNAC gene had one or more introns and contained no more than six exons, while over half of the KoNAC genes (41/68) contained three exons (Figure 4c). Additionally, a total of 10 conserved motifs were queried within all K. obovata NAC proteins. Most motifs were located within the N-terminal region (Figure 4b), and motif 1, motif 2, motif 4, and motif 5 were the common elements in KoNAC genes. Clearly, these results showed that the KoNAC genes in the same phylogenetic cluster harbored similar gene structures and motif compositions (Figure 4), which further supported the evolutionary relationship of KoNAC genes demonstrated above.

Collinearity Analysis of KoNAC Genes
It is well-known that P. trichocarpa is a typical model plant for functional genomics and molecular studies in woody species. Moreover, K. obovata (Rhizophoraceae) and P.

Collinearity Analysis of KoNAC Genes
It is well-known that P. trichocarpa is a typical model plant for functional genomics and molecular studies in woody species. Moreover, K. obovata (Rhizophoraceae) and P. trichocarpa (Salicaceae) belong to the same order, Malpighiales (https://www.ncbi.nlm. nih.gov/Taxonomy/Browser/wwwtax.cgi) (accessed on 18 March 2022). Therefore, to better investigate the evolutionary relationship of NAC genes, the collinearity analysis was performed based on the genomes of A. thaliana and P. trichocarpa ( Figure 5). There were 16,355 collinear gene pairs between K. obovata and A. thaliana identified, among which 52 orthologous gene pairs between KoNACs and AtNACs were obtained (Figure 5a, blue lines). Meanwhile, a total of 26,594 collinear gene pairs between K. obovata and P. trichocarpa were available, among which 54 orthologous gene pairs between KoNACs and PtNACs were determined (Figure 5b, purple lines). Taken together, 49 common KoNAC genes shared homologous relationships with both A. thaliana and P. trichocarpa NAC genes (Table S3), implying these genes might function in a similar manner. In the meantime, there were 11 non-orthologous KoNAC genes (KoNAC13, KoNAC16, KoNAC17, KoNAC21, KoNAC40, KoNAC50, KoNAC53, KoNAC63, KoNAC65, KoNAC67, and KoNAC68) compared to NAC genes of A. thaliana and P. trichocarpa. These genes displayed different structure (Figure 4), among which four genes (KoNAC16, KoNAC17, KoNAC21, KoNAC40) clustered as the subgroup of Class II (Figure 2).

Expression patterns of KoNAC genes in different organs
To gain an insight into the function of NAC genes in K. obovata, the expressio of all KoNAC genes in various organs, including root, stem, leaf, flower, pistil, sepal, and fruit were determined based on previously published RNA-seq data o vata. Noticeably, the expression patterns of KoNAC genes were not in a constitutiv whereas they were differentially and preferentially expressed in different organs 6, Table S4). For example, 28 out of 68 KoNAC genes were highly expressed in roo 18 KoNAC genes were preferentially expressed in leaves, and 13 KoNAC gen mainly expressed in fruits.

Expression Patterns of KoNAC Genes in Different Organs
To gain an insight into the function of NAC genes in K. obovata, the expression levels of all KoNAC genes in various organs, including root, stem, leaf, flower, pistil, stamen, sepal, and fruit were determined based on previously published RNA-seq data of K. obovata. Noticeably, the expression patterns of KoNAC genes were not in a constitutive mode, whereas they were differentially and preferentially expressed in different organs ( Figure 6, Table S4). For example, 28 out of 68 KoNAC genes were highly expressed in roots, while 18 KoNAC genes were preferentially expressed in leaves, and 13 KoNAC genes were mainly expressed in fruits. sepal, and fruit were determined based on previously published RNA-seq data of K. obovata. Noticeably, the expression patterns of KoNAC genes were not in a constitutive mode, whereas they were differentially and preferentially expressed in different organs ( Figure  6, Table S4). For example, 28 out of 68 KoNAC genes were highly expressed in roots, while 18 KoNAC genes were preferentially expressed in leaves, and 13 KoNAC genes were mainly expressed in fruits.

Expression Analysis of KoNAC Genes under Cold Treatment
To gain more insight into the function of KoNAC genes, the expression profiles of these genes under cold treatment were detected based on the public transcriptomic data of K. obovata. There were 13 KoNAC genes differentially expressed in response to chilling stress, among which one down-regulated KoNAC gene (KoNAC51) and 12 up-regulated KoNACs (KoNAC6, KoNAC11, KoNAC15, KoNAC20, KoNAC24, KoNAC26, KoNAC32, KoNAC35, KoNAC38, KoNAC41, KoNAC62, and KoNAC68) were available (Figure 7a, Table S5). Specifically, four different genes, KoNAC6, KoNAC15, KoNAC20, and KoNAC38, were largely upregulated with higher and more significant values after treatment. The expression levels of these four up-regulated and one down-regulated KoNACs were confirmed by qRT-PCR (Figure 7c,d), implying that these KoNAC genes might act as positive or negative regulators in response to chilling stress.

Stress-Related Cis-Regulatory Elements Identified in KoNAC Genes
To obtain more evidence for the differentially expressed KoNAC genes on stress responses, the cis-regulatory elements in the promoter regions of these 13 KoNAC genes were predicated. Consequently, 8 well-known stress-related elements were available (Figure 7b, Table S6). LTR (low-temperature response; CCGAAA), a core cis-acting element involved in cold stress response, was present in the majority of the detected KoNAC genes. STRE (stress response element; AGGGG) and ARE (antioxidant response element; AAACCA) were two types of regulatory elements in rapid response to anaerobic stress and environmental stimuli. ABRE (abscisic acid response element; ACGTG), ERE (ethylene response element; ATTTTAAA), and TGACG-motifs were responsible for stress induction by three major stress-related hormones, ABA, ethylene, and methyl jasmonate (MeJA), respectively. Additionally, two biotic stress-responsive elements, WRE3 (wound-response element 3; CCACCT) and WUN-motif (wound-responsive element; AAATTACT) were found in several detected promoters, as well. ble S5). Specifically, four different genes, KoNAC6, KoNAC15, KoNAC20, and KoNAC38, were largely upregulated with higher and more significant values after treatment. The expression levels of these four up-regulated and one down-regulated KoNACs were confirmed by qRT-PCR (Figure 7c,d), implying that these KoNAC genes might act as positive or negative regulators in response to chilling stress.

Discussion
The NAC gene family, one of the largest TF families in plants, was used as positive or negative regulators in response to environmental stimuli including cold stress [19]. Multiple investigations on the functional characterization of NAC genes were reported for A. thaliana [5,53,54], P. trichocarpa [14], and other plants [19]. However, the function of the NAC genes in the typical woody mangrove K. obovata responding to abiotic stresses remains largely unknown. Here, we identified a contracted NAC gene family with 68 members from the K. obovata genome. These KoNAC genes were differentially and preferentially expressed in various organs, and 13 KoNAC genes were differentially expressed under cold treatment based on the publicly available RNA-seq data.
KoNAC proteins, in company with A. thaliana NAC proteins were categorized into 10 classes according to phylogenetic analysis. Obviously, the K. obovata NAC family exhibited a significant contraction in number compared to the NAC families in A. thaliana [5], P. trichocarpa [14], and other plants [15,16], and the decreased KoNAC genes in class V, class VIII, and class IX mainly contributed to the contraction (Figure 3). These results are consistent with the previous findings [35], and the contraction might relate to the evolu-tionary adaption to the intertidal zones. Moreover, 49 genes from the contracted KoNAC family shared orthologous relationships with the NAC genes of A. thaliana and P. trichocarpa, implying these genes might have similar functions [55]. Additionally, compared to AtNAC and KoNAC genes, there existed 11 non-orthologous KoNAC genes, among which, four genes (KoNAC16, KoNAC17, KoNAC21, and KoNAC40) clustered in class II (Figure 2). Another non-orthologous gene, KoNAC68, was induced under cold treatment, implying it might potentially function in response to cold stress in K. obovata (Figure 7a). More attention should be paid to these non-orthologous genes, and functional investigations of these genes will provide valuable knowledge about mangrove species.
To explore the function of NAC genes in K. obovata, the expression patterns of all KoNACs in various organs were determined. In contrast to the constitutive expression patterns of other gene families [56][57][58], KoNAC genes were differentially and preferentially expressed in different organs. For instance, there were 28 KoNAC genes expressed highly in roots, 18 KoNAC genes expressed preferentially in leaves, and 13 KoNAC genes expressed mainly in fruits. Referentially, the root-expressed gene OsNAC2 modulated root development in rice by involving the crosstalk of auxin and cytokinin pathways [59]. The A. thaliana rosette-expressed gene, ANAC087, positively regulated rosette development and leaf senescence [60]. FaRIF, a strawberry NAC gene, was reported as one key regulator controlling fruit ripening [61]. This evidence implied that these organ-specific KoNAC genes might function as key regulators in organ development. Moreover, the organ expression patterns of NAC genes in K. obovata were not similar to that in A. thaliana and other plants. For instance, KoNAC46 and KoNAC54 were primarily expressed in roots ( Figure 6), however, ANAC048 and ANAC074, the closest orthologues of these two KoNAC genes (Figure 2), respectively, were expressed in different organs. ANAC048 was involved in vascular development [62], and ANAC074 positively regulated programmed cell death of stigmatic tissue in A. thaliana [63]. Therefore, functional characterization of the KoNAC genes primarily expressed in roots distinct from other plants should be deeply covered in the future.
To better understand the roles of KoNAC genes, the expression analysis under chilling stress was performed based on the public transcriptomic data. In total, 13 out of 68 KoNAC genes were differentially expressed under cold treatment. Among them, KoNAC51 was the only downregulated gene, whereas its closest homologue KoNAC35 was up-regulated after treatment, implying these two class V genes might function oppositely in response to cold stress. Half of the up-regulated genes (KoNAC6, KoNAC15, KoNAC20, KoNAC26, KoNAC32, and KoNAC38) belonged to the class VII subgroup (Figures 2 and 7a). Particularly, KoNAC6, KoNAC15, and KoNAC26 clustered together and shared high sequence similarity to their closest orthologs ANAC002, ANAC081, and ANAC102 in A. thaliana. ANAC002 (ATAF1) was reported to serve as dual regulators responsive to abiotic and biotic stresses [64][65][66]. ANAC081 (ATAF2) was rapidly induced by pathogen attack and involved in plant defense [67,68], while ANAC102 was responsive to low-oxygen and high-light stresses [69,70]. Overexpression of MlNAC5, another closest ortholog of KoNAC26 and ANAC002 from Miscanthus lutarioriparius L.Liu ex S.L.Chen & Renvoize, led to enhanced tolerance to cold and drought stresses in A. thaliana [71]. Additionally, three closest orthologs of KoNAC32, ANAC019, ANAC055 and ANAC072, were required for drought tolerance in A. thaliana [72], among which ANAC019 and ANAC055 displayed a dual function in regulating ABA response and jasmonate response [73,74]. A. thaliana ANAC042, the closest ortholog of KoNAC38, conferred stress tolerance through regulating phytohormone metabolism and signaling [75][76][77]. Moreover, various stress-related cis-regulatory elements were identified from the promoters of these KoNAC genes (Figure 7b). The LTR element is an indispensable cis-acting element in plant response to low temperature [78,79]. Deletion of the LTR element will result in complete loss of promoter activity under cold stress [79]. STRE is a common cis-regulatory element in eukaryotes, and involved in response to multiple environmental stimuli [80]. ARE is an antioxidant response element in rapid response to anaerobic stress [81]. Meanwhile, ABRE, ERE, and TGACG motifs are three major types of elements related to plant hormones (ABA, ethylene, and MeJA) [82][83][84]. Among them, ABRE and TGACG motifs are enriched in the majority of the detected KoNAC genes, implying these KoNAC genes might respond to stresses via hormone-mediated pathways [85]. Additionally, both WRE3 and WUN motifs are biotic stress-responsive elements and present in several detected promoters as well, implying the KoNAC genes might function in response to biotic stresses [86,87]. Taken together, these findings suggest that these KoNAC genes may be involved responses to other abiotic or biotic stresses in addition to the cold response, providing auxiliary evidence for these KoNAC genes in response to abiotic and biotic stresses. To know more about the function of KoNAC genes, further investigation and more proof are required.

Conclusions
In the present study, a pluri-disciplinary work concerning comprehensive analysis of the KoNAC gene family was performed. We identified a contracted NAC TF family containing 68 genes from the genome of the typical mangrove plant K. obovata based on bioinformatic analysis. These KoNAC genes were unevenly distributed in 17 chromosomes of K. obovata. The NAC genes of K. obovata and A. thaliana were classified into ten classes, while no KoNAC gene belonged to class I. Obviously, the decreased members of class V, class VIII, and class IX mainly contributed to the contraction of the KoNAC family. KoNAC genes were differentially and preferentially expressed in different organs. Among them, 13 KoNAC genes were rapidly induced by chilling stress. The expression patterns of five selected KoNAC genes (KoNAC6, KoNAC15, KoNAC20, KoNAC38, and KoNAC51) were confirmed by qRT-PCR. Additionally, several stress-related cis-acting elements were detected in the promoter regions of these KoNAC genes, implying KoNAC genes might participate in multiple stress responses. Summarily, our findings will provide positive references for further investigations on functional characterization of KoNAC genes in stress responses.
Supplementary Materials: The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/cimb44110381/s1, Figure S1: Morphological features of K. obovata; Table S1: The basic information of A. thaliana NAC genes; Table S2: The primers of KoNAC genes in this study; Table S3: Common orthologous gene pairs of K. obovata between A. thaliana and P. trichocarpa; Table S4: The FPKM values of KoNAC genes expressed in different organs; Table S5: The FPKM values of differentially expressed KoNAC genes under cold treatment; Table S6: Number of stress-related cis-regulatory elements in the promoter regions of the differentially expressed KoNAC genes.