Genome-Wide Identification and Characterization of NAC Family in Hibiscus hamabo Sieb. et Zucc. under Various Abiotic Stresses

NAC transcription factor is one of the largest plant gene families, participating in the regulation of plant biological and abiotic stresses. In this study, 182 NAC proteins (HhNACs) were identified based on genomic datasets of Hibiscus hamabo Sieb. et Zucc (H. hamabo). These proteins were divided into 19 subfamilies based on their phylogenetic relationship, motif pattern, and gene structure analysis. Expression analysis with RNA-seq revealed that most HhNACs were expressed in response to drought and salt stress. Research of quantitative real-time PCR analysis of nine selected HhNACs supported the transcriptome data’s dependability and suggested that HhNAC54 was significantly upregulated under multiple abiotic stresses. Overexpression of HhNAC54 in Arabidopsis thaliana (A. thaliana) significantly increased its tolerance to salt. This study provides a basis for a comprehensive analysis of NAC transcription factor and insight into the abiotic stress response mechanism in H. hamabo.


Introduction
NAC transcription factor is one of the largest plant gene families. The most notable structural feature of NAC transcription factor is the presence of a highly conserved NAC domain (about 150-160 amino acids) at the N-terminus of the protein, while the C-terminus is a modified transcriptional regulatory region (TAR) [1]. It participates in the regulation of plant stress responses and life processes such as plant leaf senescence, flower formation, fruit maturation and coloring, seed development, and root development. For example, tomato SlNAC6 is involved in drought stress response and the reproductive process [2]. Pepper CaNAC035 is a positive regulator of abiotic stress tolerance [3], and CaNAC064 modulates plant tolerance to cold stress [4]. Rice ONAC127 and ONAC129 are involved in the grain filling process [5], and OsNAC2 participates in ABA-induced leaf senescence [6]. Wheat TaNAC30 negatively regulates resistance to stripe rust [7], and TaNAC2-5A controls the nitrate response and increases wheat yield [8]. Because of the development of highthroughput sequencing technology, an increasing number of plants have completed wholegenome sequencing. Therefore, the identification research of NAC has entered an explosive period. Research on the identification and expression of related gene families has been carried out in Triticum aestivum [9], Fagopyrum tataricum [10], Betula pendula [11], Prunus mume [12], Panicum virgatum [13], and other species.
Studying the responses of plants to abiotic stress is helpful in better protecting and utilizing plants. Hibiscus hamabo Sieb. et Zucc. (H. hamabo) is a semi-mangrove plant belonging to Hibiscus Maliace. It is extremely resistant to salt and alkalis, as well as drought and barren environments, and it has developed roots. Moreover, it is a pioneer tree species with excellent wind resistance and embankment protection in coastal areas. It is a woody plant that integrates medicinal, edible, health care, ornamental, greening, water and soil protection, and fibrous raw material values [14]. Our team has completed the whole genome sequencing of H. hamabo (not yet published). On this basis, we have carried out a genome-wide identification of the NAC transcription factor gene family, and we have verified the function of the HhNAC54 gene. This study has certain significance for the analysis of the molecular response mechanism of H. hamabo under abiotic stress, and it can also provide an excellent gene reserve for the molecular selection of salt-tolerant plants.

Identification and Analysis of HhNAC Genes in H. hamabo
Based on the genomic information of H. hamabo, 185 putative NAC genes were obtained with HMM. Three genes without the NAC and NAM domain were eliminated according to the results from the SMART and Pfam databases, leaving 182 HhNAC genes (Table S1). The lengths of protein sequences, chromosome locations, start and stop sites, protein molecular weights, and isoelectric points (pI) were examined (Table S2). The results showed that 182 genes were randomly distributed on 45 chromosomes, except chromosome 35, and they were named HhNAC1-HhNAC182 according to their positions on the chromosomes. The protein lengths were between 203 (HhNAC79) and 799 (HhNAC170) amino acids, and the protein molecular weights were between 23,021.5 (HhNAC79) and 88,094.4 Da (HhNAC170).
NAC domains were found to cover a range of 40-50 amino acids with a high similarity compared to HhNAC protein domains. To further study the evolutionary relationships, we compared 182 HhNAC and AtNAC protein sequences containing an NAC domain ( Figure 1). The 182 HhNACs could be classified into 19 subfamilies (A-S), of which the C subfamily contained the most members (22), whereas the O, N and F subfamilies contained only two members each. In addition, HhNAC subfamilies were expanded or contracted to different degrees compared with A. thaliana. For example, 0 members belonged to the Q subfamily in A. thaliana, but the number increased to 11 in H. hamabo. This situation also appeared in B, R, C, H, I, J, E, and P subfamilies and others, whereas in the O and N subfamilies, the number of HhNACs decreased significantly.
The motif compositions of 182 HhNAC proteins were analyzed using an online MEME tool to better understand their conserved structures and evolutionary relationships and ultimately revealed 20 distinct motifs ( Figure 2). Corresponding with the result of the phylogenetic analysis, the components of the same subfamily were formed by similar motifs. Except in the widely distributed Motifs 2, 3, 5, and 13, HhNAC members in the same subfamily have similar motifs, such as subfamilies A, B, C, and S, suggesting that the function of HhNAC proteins in the same subfamily may be similar. In addition, some motifs only appear in a specific subfamily, such as Motif 7 and Motif 10, which only appear in subfamily S, and Motif 17, which only appears in subfamily P. These motifs may be important in specific subfamilies, which are worthy of further study. A few HhNAC proteins with different motifs were classified in the same branch, indicating that these HhNACs may have undergone functional differentiation in the evolution process. The structural study of HhNAC is further performed in this research (  The motif compositions of 182 HhNAC proteins were analyzed using an online MEME tool to better understand their conserved structures and evolutionary relationships and ultimately revealed 20 distinct motifs ( Figure 2). Corresponding with the result of the phylogenetic analysis, the components of the same subfamily were formed by similar motifs. Except in the widely distributed Motifs 2, 3, 5, and 13, HhNAC members in the same subfamily have similar motifs, such as subfamilies A, B, C, and S, suggesting that the function of HhNAC proteins in the same subfamily may be similar. In addition, some motifs only appear in a specific subfamily, such as Motif 7 and Motif 10, which only appear in subfamily S, and Motif 17, which only appears in subfamily P. These motifs may be important in specific subfamilies, which are worthy of further study. A few HhNAC proteins with different motifs were classified in the same branch, indicating that these HhNACs may have undergone functional differentiation in the evolution process. The structural study of HhNAC is further performed in this research ( Figure 2). A typical structure containing two exons existed in most HhNAC members (104 HhNACs), and the intron number of the remaining HhNACs ranged from 0 to 17, indicating the functional differentiation of HhNAC family members. Similarly conserved motifs and gene structures in HhNACs of the same subfamilies further confirmed the reliability of evolutionary relationships. According to the annotation information, 182 HhNAC genes were annotated on 45 chromosomes ( Figure 3). HhNAC genes are distributed on 45 chromosomes, except chromosome 35, among which chromosome 1 contained the highest number of genes (9 HhNACs), whereas some chromosomes only contained one HhNAC gene, such as Chromosomes 37, 44, and 46. In addition, a positive correlation was not found between the number of genes and the length of the chromosome. In order to explore the evolutionary process of HhNAC genes, intragenomic synteny analysis was conducted, and a total of 18 pairs of HhNAC genes were found to have a collinear relationship ( Figure 3). The substitution rate (Ka/Ks) of the 18 pairs of gene synonymous (Ks) and non-synonymous (Ka) mutations was calculated, which ranged from 0.16 to 0.47, with an average of 0.28, suggesting that purifying the selection was key to the evolution of the HhNAC family.
In order to understand the phylogenetic mechanism of HhNACs, we performed synteny analysis of HhNACs with NAC genes in A. thaliana and Populus trichocarpa. There were 24 and 26 HhNAC genes in a syntenic relationship with A. thaliana and P. trichocarpa, respectively ( Figure 4). These genes may be important in evolution of the NAC gene family. In addition, some syntenic gene pairs were anchored to highly conserved syntenic blocks, spanning more than 100 genes, indicating that there may be an evolutionary relationship between H. hamabo and two other plants.  tionary process of HhNAC genes, intragenomic synteny analysis was conducted, and a total of 18 pairs of HhNAC genes were found to have a collinear relationship ( Figure 3). The substitution rate (Ka/Ks) of the 18 pairs of gene synonymous (Ks) and non-synonymous (Ka) mutations was calculated, which ranged from 0.16 to 0.47, with an average of 0.28, suggesting that purifying the selection was key to the evolution of the HhNAC family.

Expression Patterns of HhNAC Genes in Different Abiotic Stresses
To further explore the biological function of HhNACs, the expression level of 182 HhNAC genes were analyzed under drought and salt stress, which was based on the transcriptomic data ( Figure 5). Only a few HhNACs were not expressed in all samples (FPKM = 0), and the expressed genes (FPKM > 0) showed significant differences in expression pattern during different treatments and time points. Some genes exhibited the same pattern of expression under drought and salt stress, as HhNAC70 and HhNAC168 of the E subfamily, HhNAC146 of the G subfamily, and HhNAC41, HhNAC73 and HhNAC118 of the D subfamily were all upregulated under the two stresses, whereas HhNAC17 of the J subfamily, HhNAC44 of the P subfamily, and HhNAC100 of the R subfamily were downregulated under both stress conditions. Several genes showed opposing expression trends under the two types of stress. HhNAC103 of subfamily B was induced under PEG treatment and reduced under NaCl treatment. HhNAC50 and HhNAC158 of subfamily E were induced under NaCl treatment and reduced under PEG treatment. The results suggest that they may participate in the stress response process in different ways.
In order to understand the phylogenetic mechanism of HhNACs, we performed synteny analysis of HhNACs with NAC genes in A. thaliana and Populus trichocarpa. Ther were 24 and 26 HhNAC genes in a syntenic relationship with A. thaliana and P. trichocarpa respectively (Figure 4). These genes may be important in evolution of the NAC gen family. In addition, some syntenic gene pairs were anchored to highly conserved syntenic blocks, spanning more than 100 genes, indicating that there may be an evolu tionary relationship between H. hamabo and two other plants.

Expression Patterns of HhNAC Genes in Different Abiotic Stresses
To further explore the biological function of HhNACs, the expression level of 18 HhNAC genes were analyzed under drought and salt stress, which was based on th transcriptomic data ( Figure 5). Only a few HhNACs were not expressed in all sample (FPKM = 0), and the expressed genes (FPKM > 0) showed significant differences in ex pression pattern during different treatments and time points. Some genes exhibited th same pattern of expression under drought and salt stress, as HhNAC70 and HhNAC168 o the E subfamily, HhNAC146 of the G subfamily, and HhNAC41, HhNAC73 and HhNAC118 of the D subfamily were all upregulated under the two stresses, wherea HhNAC17 of the J subfamily, HhNAC44 of the P subfamily, and HhNAC100 of the R subfamily were downregulated under both stress conditions. Several genes showed op posing expression trends under the two types of stress. HhNAC103 of subfamily B wa induced under PEG treatment and reduced under NaCl treatment. HhNAC50 and HhNAC158 of subfamily E were induced under NaCl treatment and reduced under PEG treatment. The results suggest that they may participate in the stress response process in different ways. Plant responses to drought and salt stress are also influenced by ABA and MeJA signals. Nine genes with the largest fold changes under the PEG and NaCl treatments were selected for qRT-PCR analysis, so that H. hamabo responses to salt, drought stresses, and ABA and MeJA treatments could be studied at six time points (Figure 6). The expression trends of these nine genes at the same time point in all samples showed consistency with the transcriptome, supporting the transcriptome result's dependability. Plant responses to drought and salt stress are also influenced by ABA and MeJA signals. Nine genes with the largest fold changes under the PEG and NaCl treatments were selected for qRT-PCR analysis, so that H. hamabo responses to salt, drought stresses, and ABA and MeJA treatments could be studied at six time points (Figure 6). The expression trends of these nine genes at the same time point in all samples showed consistency with the transcriptome, supporting the transcriptome result's dependability. Meanwhile, some genes showed different expression trends in the early stage of the stress response. For example, HhNAC45 and HhNAC61 had higher expression levels at 6 h under salt stress, which began to decrease after 6 h, thus indicating that they may play different roles in the early and late stages of the stress response. Some HhNAC genes could be significantly induced or suppressed by the four treatments, such as HhNAC54 and HhNAC61. Furthermore, several genes exhibited opposing expression trends when responding to distinct environments. HhNAC127, for example, was upregulated under PEG and NaCl stress, but downregulated under ABA and MeJA treatments. These findings suggested that the HhNAC family participated in multiple stress responses and had diverse regulatory mechanisms. 022, 23, 3055 8 of 13 Figure 6. HhNAC genes were analyzed by qRT-PCR. The significance analysis was carried out using Student's t-test (* p < 0.05, ** p < 0.01).

HhNAC54 Overexpression Increases Transgenic A. thaliana Salt Tolerance
Because of the expression pattern under the four treatments, the HhNAC54-overexpressing plants were reinforced with A. thaliana to study whether the expression of HhNAC54 could improve the salt tolerance (Figure 7). Five T3 transgenic plants were obtained through resistance screening. The PCR was set for detecting and determining the expression of HhNAC54 in T3 plants. Three independent transgenic lines were selected for subsequent studies (OE1-3). Transgenic plants and wild types (WT) were grown on a 1/2 MS medium with progressively elevated salt concentrations. There was no obvious difference in the root length of transgenic plants and WT planted on normal 1/2 MS medium after 10 d. Root growths of transgenic plants were less inhibited by high salinity than those of wild type, showing that resistance to salt stress in the transgenic A. thaliana was increased, and HhNAC54 had a function in improving salt tol- Figure 6. HhNAC genes were analyzed by qRT-PCR. The significance analysis was carried out using Student's t-test (* p < 0.05, ** p < 0.01).

HhNAC54 Overexpression Increases Transgenic A. thaliana Salt Tolerance
Because of the expression pattern under the four treatments, the HhNAC54-overexpressing plants were reinforced with A. thaliana to study whether the expression of

Discussion
The NAC transcription factor mainly exists in terrestrial plants [15]. NAC family was identified in H. hamabo for the first time in this study, and the member quantity (182) was more than that in A. thaliana (117), rice (151) [16], and soybean (152) [17]. The phenomenon of multiple members may be caused by the extensive duplication and diversification of the genome during evolution [18]. The gene family can be augmented in a variety of ways, including genome-wide replication, tandem replication, fragment replication, and others [18].
HhNAC members can be divided into 19 subfamilies, and almost all subfamilies contained both HhNAC and AtNAC proteins. Differences in the number of gene subfamily members between H. hamabo and A. thaliana may be due to genes that are copied or lost. Genes clustered in the same subfamily were closely related, which indicates that they might have similar functions. The evolutionary relationship between HhNAC members can also be analyzed by comparing conserved motifs. In this study, the HhNACs in the same subfamily shared common motifs, which indicates that they were highly conserved, further proving the reliability of HhNAC classification.
In recent years, many transcription factors have been identified as enhancing plant adaption to stress [19,20]. A large amount of evidence indicated that the NAC family played a key role in determining plant responses to abiotic stress. For example, the overexpression of NAC enhanced rice tolerance to salt and drought stress [21]. The wheat TaSNAC4 homeologous gene TaSNAC4-3A was induced under a series of abiotic stresses, and the overexpression of TaSNAC4-3A in A. thaliana can endow it with a tolerance to drought stress [22]. The expression trends of HhNAC genes were analyzed based on transcriptome data from roots exposed to drought and salt stress. The expression trends of most HhNAC genes after drought and salt stress exhibited significantly different levels of expression. Some genes exhibited the same expression trend under the two types of stresses, indicating that they might play similar roles in responding to the two stresses. Some genes showed opposite expression trends, indicating that they might play completely different roles. qRT-PCR analyses in this study were performed to further verify gene expression levels under stresses. Moreover, the additional three time points were set in order to compare the expression patterns in the early and late stages of the response to stresses.

Discussion
The NAC transcription factor mainly exists in terrestrial plants [15]. NAC family was identified in H. hamabo for the first time in this study, and the member quantity (182) was more than that in A. thaliana (117), rice (151) [16], and soybean (152) [17]. The phenomenon of multiple members may be caused by the extensive duplication and diversification of the genome during evolution [18]. The gene family can be augmented in a variety of ways, including genome-wide replication, tandem replication, fragment replication, and others [18].
HhNAC members can be divided into 19 subfamilies, and almost all subfamilies contained both HhNAC and AtNAC proteins. Differences in the number of gene subfamily members between H. hamabo and A. thaliana may be due to genes that are copied or lost. Genes clustered in the same subfamily were closely related, which indicates that they might have similar functions. The evolutionary relationship between HhNAC members can also be analyzed by comparing conserved motifs. In this study, the HhNACs in the same subfamily shared common motifs, which indicates that they were highly conserved, further proving the reliability of HhNAC classification.
In recent years, many transcription factors have been identified as enhancing plant adaption to stress [19,20]. A large amount of evidence indicated that the NAC family played a key role in determining plant responses to abiotic stress. For example, the overexpression of NAC enhanced rice tolerance to salt and drought stress [21]. The wheat TaSNAC4 homeologous gene TaSNAC4-3A was induced under a series of abiotic stresses, and the overexpression of TaSNAC4-3A in A. thaliana can endow it with a tolerance to drought stress [22]. The expression trends of HhNAC genes were analyzed based on transcriptome data from roots exposed to drought and salt stress. The expression trends of most HhNAC genes after drought and salt stress exhibited significantly different levels of expression. Some genes exhibited the same expression trend under the two types of stresses, indicating that they might play similar roles in responding to the two stresses. Some genes showed opposite expression trends, indicating that they might play completely different roles. qRT-PCR analyses in this study were performed to further verify gene expression levels under stresses. Moreover, the additional three time points were set in order to compare the expression patterns in the early and late stages of the response to stresses.
Based on the results of the transcriptome data and qRT-PCR analysis, HhNAC54 was selected for further functional analysis using transgenic A. thaliana. The changes in the root length of the transgenic A. thaliana were significantly lower than those of the WT, which indicates that HhNAC54 had a function in improving salt tolerance, and that the regulation mechanism of HhNAC54 is worth further studying.

Plants Materials and Treatments
H. hamabo was taken from the Institute of Botany, Chinese Academy of Sciences, Jiangsu Province. After 20 days of vernalization, the plump and healthy H. hamabo seeds were soaked in concentrated sulfuric acid for 15 min and then washed with tap water repeatedly. The treated seeds then were sown onto plugs, the substrate was composed of peat and vermiculite (1:1) in the pH range of 5-6, and they were cultured in a light incubator (Volume: 1008 L) with a photoperiod of 16/8 h (day (12,000 lux intensity white light LEDs)/night), and a temperature of 24/20 • C under a relative humidity of 65%. The H. hamabo seedlings that grew 8~10 true leaves were treated. The seedlings were taken out of the container, cleaned, transferred to a hydroponic culture in 1/2 MS solution for 7 d, and then subjected to abiotic stress treatment in the container with same light intensity, photoperiod, and temperature. Treatment 1: 400 mM NaCl. Treatment 2: 15% polyethylene glycol (PEG 6000). Treatment 3: 50 µM abscisic acid (ABA). Treatment 4: 1 mM methyl jasmonate (MeJA). Treatments 3 and 4 were spraying seedlings evenly. After the treatment, the leaves were collected at 1 h, 2 h, 6 h, 12 h, and 24 h, frozen quickly in liquid nitrogen, and then stored at −80 • C for later use. Each treatment has three biological replicates.

Identification and Analysis of NAC Genes in H. hamabo
Based on the whole genome data of H. hamabo (not yet published), NAM domain (PF02365) and NAC (PF01849) Hidden Markov Model (HMM) profiles from the Pfam database (http://pfam.xfam.org, accessed on 16 March 2021) were used as probes to screen NAC genes in H. hamabo by HMMER (v. 3.1) software, and the E value was set to E < 10 −3 [23,24]. ClustalW software (v2.1) was used for performing a multiple sequence alignment with the obtained NAC candidate sequences to remove redundant sequences [25,26]. MEGA7.0 software was used to compare the NAC protein sequences of H. hamabo and A. thaliana [27,28]. The neighbor-joining method was selected for constructing the phylogenetic tree. The parameter of random sampling repetition times was set to 1000 bootstrap, and the system default values were used for other parameters.
The GSDS website (http://gsds.cbi.pku.edu.cn, accessed on 16 March 2021)) was used to draw a schematic diagram of the NAC gene structure [29]. MEME (http://memesuite.org/, accessed on 16 March 2021) was made to find the conservative motif and draw its protein secondary structure [30]. The parameter is set to 10 motifs. TBtools software (v. 1.089, JAVA, China) was used to draw the gene structure map and for chromosome mapping [31]. MCScanX software was used to analyze collinearity, tandem repeats, and the homologous evolution among different species, and Circos software was used for drawing [32,33]. KaKs_Calculator (v.2.0) software was used to calculate the Ks and Ka/Ks values between gene pairs [34].

Quantitative Real-Time PCR
Quantitative real-time PCR was determined using the mature method [35]. The PCR instrument and reagent were a StepOnePlus real-time PCR system (Applied Biosystems) and a SYBR Green Master Mix (Bimake, Houston, TX, USA), respectively. Total RNA was extracted with Plant Rneasy Mini Kit (Qiagen, Hilden, Germany). The first cDNA strand was synthesized by A PrimeScript ® RT kit (TaKaRa, Dalian, China). Genscript online design software (https://www.genscript.com/tools/pcr-primers-designer, accessed on 16 March 2021) was used to design the primer pairs (Table S3). A clustered heat map was drawn using Tbtools software. The reference gene was Actin (HhACT) [14]. Three independent experiments were conducted. Relative transcript abundances were analyzed using the 2 −∆∆CT method [36].

Identification of Transgenic A. thaliana
HhNAC transgenic lines had obtained as described previously [37]. In brief, the HhNAC ORF was inserted into the pCAMBIA1305 vector using the ClonExpress II One Step Cloning Kit (Vazyme, Nanjing, China). Then, pCAMBIA1305-HhNAC was transformed into Agrobacterium tumefaciens GV3101 using the freeze-thaw method [38]. Col-0 of A. thaliana was transformed via the floral dip method [38]. Hygromycin screened (150 mg/L) homozygous T 3 transgenic A. thaliana lines were used for functional identification.

Conclusions
This is the first study to identify the NAC family in H. hamabo at the genome level. A total of 182 HhNACs were obtained, which were then classified into 19 groups based on phylogenetic relationship, motif pattern, and gene structure analysis. The expression patterns of HhNACs under drought and salt stress were analyzed using RNA-seq and qRT-PCR, and the stress response gene HhNAC54 was obtained for further functional analysis. Overexpression of HhNAC54 in A. thaliana significantly increased the tolerance to salt. This study provided a basis for a comprehensive analysis of the NAC gene family and abiotic stress response mechanism in H. hamabo.