To identify HcNAC genes, we downloaded the files of H. citrina including the genome sequences, coding sequences (CDS), and protein sequences from NCBI (https://www.ncbi.nlm.nih.gov/assembly/GCA_017893485.1) (Qing et al., 2021). Hidden Markov model (HMM) files were constructed based on the NAM domain (PF02365) retrieved from the Pfam website (https://pfam.xfam.org/). The conserved NAM domain was utilized to search for HcNAC protein sequences by a program of HMMER3.0 (E-value ≤ 10⁻⁵). Then, the physicochemical property prediction of HcNAC protein sequences was detected using TBtools software (Chen et al., 2020). The secondary structure prediction was performed by the online website (https://npsa-prabi.ibcp.fr/NPSA/npsa_sopma.html). WoLF PSORT software online (https://wolfpsort.hgc.jp/) was used to predict the subcellular localization of HcNAC proteins.

We obtained the genome annotation GFF3 file of the NAC genes in night lily from the NCBI database, which was used for mapping the HcNAC gene chromosomal positions using one small program (Gene Location Visualize from GTF/GFF) of TBtools. A. thaliana, O. sativa, and C. lanatus NAC proteins were obtained from the Arabidopsis Information Resource (http://www.Arabidopsis.org/), rice genome annotation (http://rice.plantbiology.msu.edu/), and Cucurbit Genomics Database (http://cucurbitgenomics.org/), respectively (Ooka et al., 2003; Nuruzzaman et al., 2010; Lv et al., 2016). Combined with identified HcNAC proteins, a phylogenetic tree was constructed using the MEGA11.0.13 integrated tool by the Neighbor-Joining method (Kumar et al., 2016). The tree nodes were evaluated by 5,000 bootstrap replicates. HcNAC gene structure was analyzed using the Visualize Gene Structure of TBtools applets.

Gene-conserved motif prediction was performed via the MEME tool (http://meme-suite.org/index.html) with default settings. NCBI conserved domain database was used to predict the conserved domains of the HcNACs. Genome synteny analysis was made as described previously (Sun et al., 2017). We used Advanced Circos and dual synteny plot of TBtools software to show homologous gene pairs. Nonsynonymous (Ka) and synonymous (Ks) rates among protein sequences were used to assess the DNA sequence evolution. To appraise the divergence of duplicated night lily NAC genes, the selective strength was estimated by calculating the Ka/Ks ratio between paralogous gene pairs using the Simple Ka/Ks Calculator Tool (NG) in TBtools. Ka/Ks larger than 1 indicates positive selection, Ka/Ks less than 1 indicates purifying selection, and Ka/Ks equal to 1 indicates neutral mutation (Zhang and Yu, 2006).

Gene functional enrichment analysis was performed to reveal the biological processes, cellular components, and molecular functions of the HcNAC genes using STRING (https://cn.string-db.org/). The results were visualized using microscopic letter website (http://www.bioinformatics.com.cn/). For the cis-acting element analysis, about 2,000 base pairs of promoter regions upstream from the initiation codon of HcNAC genes were extracted and then analyzed by the PlantCARE database (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/) to hunt for promoter cis-elements (Lescot et al., 2002).

In this study, H. citrina was used as treated plant species and cultivated in the experimental base. We chose F1 hybrid population 116 as the treated material based on its moderate resistance to abiotic stresses, obtained by Dongzhuang Huanghua as the female parent and Chonglihua as the male parent. Scapes were taken as explants for tissue cultivation (Zuo et al., 2024). These obtained seedlings were cultured in a 28°C growing box with the 16 h light/8 h dark condition and treated with Hoagland solution respectively containing 20% PEG6000 and 250 mM NaCl after 4 weeks old (Cao et al., 2024). Roots were collected at different time points with the same interval and frozen at a lower temperature than −70°C immediately for later use.

To investigate the transcriptional dynamics of the HcNAC genes in different tissues, bud, tender leaf, mature leaf, tender scape, mature scape, tender root, and mature root were collected and then sent to BMK Biotechnology Company for transcriptome sequencing analysis (NCBI accession: PRJNA1154842). We also chose collected roots of five groups (0, 24, 48, 72, and 108 h) under drought stress (NCBI accession: PRJNA1154840) and five groups (0, 24, 48, 72, and 96 h) under salinity stress for RNA-Seq analysis (NCBI accession: PRJNA1154841). All samples were designed for three replicates. Raw data were performed quality control checks through FastQC and filtered by Trimmomatic 0.36 using the quality control results (Bolger et al., 2014). HISAT2 2.2.1 was used to map the paired-end clean reads to the H. citrina reference genome (Kim et al., 2015). DESeq2 1.30.1 were used to determine differentially expressed genes with the analysis of variance (ANOVA) method (p.adjust < 0.05, |Log₂FC| ≥ 1.5) (Love et al., 2014). The fragments per kilobase per million mapped fragments (FPKM) values of HcNAC genes were summed and used for measuring the transcript abundance of HcNACs. The heat map was generated with log₂FPKM values to visualize different expression levels.

RT-qPCR analysis was carried out according to previously reported with little modification (Hu et al., 2016). The cDNA was synthesized by reverse transcription of high-quality RNA using the TIANScript II RT Kit by the corresponding instructions (TianGen, Beijing, China). RT-qPCR–specific primers, designed using Primer Premier 5, were shown in Supplementary Table S11 . The gene expression levels were detected using a LightCycler480 II (Roche, Basel, Switzerland). FastKing One-Step SYBR Green Kit (TianGen, Beijing, China) was used for RT-qPCR reactions, and the calculation of the gene expression levels was analyzed through the 2^−ΔΔCt Method. Significant difference in gene expression levels between the two treatment groups was analyzed by a one-way ANOVA analysis of variance method with Duncan test at a p-value < 0.05. Gene expression was analyzed using three independent biological repeats. To normalize expression levels of the selected HcNAC genes, HcACTIN gene was used as an internal control (Cao et al., 2024).

HcNAC35 and HcNAC71 fragments with KpnI and XhoI restriction sites were obtained by conventional PCR. The plasmid pSuper1300 with Green Fluorescent Protein (GFP) was digested with these two enzymes to get the linear vector fragment, which was subsequently ligated with the target fragments using pEASY-Basic Seamless Cloning and Assembly Kit (TransGen, Beijing, China). Constructed fusion expression plasmids including pSuper1300-HcNAC35-GFP and pSuper1300-HcNAC71-GFP were transformed into tobacco leaf cells by the Agrobacterium-mediated transient transformation. After 48–60 h of dark culture, GFP fluorescence signals were captured using a laser fluorescence microscope (Zeiss LSM800). Relative primers were listed in Supplementary Table S11 .

The full-length CDS with the BamHI and KpnI restriction sites of HcNAC35 were amplified for constructing an overexpression vector and then cloned into a 1305.4-with Green Fluorescent Protein (GFP) vector to carry out sequencing Tongchuan verification. The resulting construct 1305.4-HcNAC35-GFP was transformed into Agrobacterium tumefaciens strain EHA105 to perform genetic transformation using an optimized transformation system in watermelon Tongchuan (TC) (Cao et al., 2022). Transgenic plants were screened by GFP observation and PCR analysis. Wild-type (WT) TC and transgenic plants were transplanted and cultivated in the same growing conditions. Primers referred to in the experiment were given in Supplementary Table S11 .

The levels of MDA (malondialdehyde) were assessed using the thiobarbituric acid method as described in previous study (Hu et al., 2019). Approximately 0.2 g leaf samples of both salinity-stressed and unstressed plants were collected and ground with 2 ml of ice-cold 0.5% TCA (trichloroacetic acid). The homogenates were centrifuged at 7,000 rpm for 15 min at 4°C. A 1.5-ml volume of the supernatant and 1.5 ml of 8% TCA containing 0.5% thiobarbituric acid were mixed, and boiled for 10 min. The homogenate cooled to ambient temperature was centrifuged at 7,000 rpm for 15 min. Absorbance was measured at 450, 532, and 600 nm.

The level of O₂ ⁻ was also determined as previously described (Hu et al., 2019). Leaf samples (0.3 g per sample) were ground with 2 ml of ice and precooled 50 mM phosphate buffer solution (PBS, pH 7.8). The homogenates were then centrifuged at 10,000 rpm for 20 min at 4°C. A 1 ml of the supernatant, 250 μL of hydrochloride hydroxylamine, and 750 μL of 65 mM PBS (pH 7.8) were mixed and incubated for 1 h at 25°C. Then, a 1 ml of the reaction mixture, 1 ml of 7 nM α-naphthylamine, and 1 ml of 17 nM paminobenzenesulfonic acid were mixed and incubated for 30 min at 25°C. Finally, the absorbance was measured at 530 nm. The O₂ ⁻ content was calculated using a standard curve relating O₂ ⁻ concentration to absorbance. H₂O₂ was assessed using Micro Hydrogen Peroxide (H₂O₂) Assay Kit (Solarbio, Beijing, China). Absorbance was recorded at 450, 532, and 600 nm using an Infinite M200 microplate reader (Tecan, Männedorf, Switzerland).

Multiple sequence alignment analysis was performed, followed by a total of 113 complete non-redundant NAC genes identified based on the PlantTFDB database and the reported H. citrina genome (Qing et al., 2021). The chromosomal distributions of genes showed that they were unequally dispersed across LG1-LG11 with gene counts ranging from 4 to 19. These genes were designated HcNAC1-HcNAC113 according to their chromosomal positions ( Figure 1 ). The full length of encoded HcNAC proteins varied significantly from 94 to 1,285 amino acid residues. Physicochemical properties were further analyzed based on their protein sequences, with molecular weights ranging from 10.7 to 141.9 kDa, and theoretical pI values from 4.42 to 10.4 ( Supplementary Table S1 ). The secondary structure including alpha helix, beta-turn, random coil, and extended strand was listed in Supplementary Table S2 . The subcellular localization prediction indicated that most HcNAC proteins were localized in the nucleus region, but the others were mainly in either chloroplast or cytoplasm.

To explore the evolutionary relationships among H. citrina NAC TFs, a phylogenetic tree was established using the NAC protein sequences from H. citrina, C lanatus, O. sativa, and A. thaliana ( Figure 2 ). The protein sequences (ClaNACs, AtNACs, OsNACs) involved in tree construction are listed in Supplementary Table S3 . The phylogenetic analysis indicated that 113 HcNAC proteins tightly clustered with ClaNAC, AtNAC, and OsNAC proteins possessed 16 subgroups except the NAC-B and NAC-G subgroup, which suggested evolutionary conservation of the HcNAC gene family. The NAC-L subgroup comprised 17 HcNACs, 7 ClaNACs, 7 OsNACs, and 14 AtNACs, belonging to the divided largest subfamily. The NAC-D subgroup consisted of six OsNAC genes, but only two HcNACs, two ClaNACs, and two AtNAC genes. In the NAC-B subgroup, 27 OsNACs, 1 ClaNAC, and 1 AtNAC gene were found while no NAC gene was in H. citrina. The NAC-G subgroup contained six AtNACs and two OsNAC genes, but no ClaNACs and HcNAC genes. The types and number of genes within each subfamily varied greatly. HcNAC TF was absent in the NAC-B and NAC-G subgroups, which implies that these groups might be lost in night lily during evolution.

To further investigate the evolutionary conservation of HcNAC genes, multiple sequence alignment was conducted to explore the homologous domain features and frequency of amino acids, exhibiting high-sequence homology in the same subgroup and showing a bit of difference in the conservatisms of some amino acids among HcNAC proteins of different subgroups ( Figure 3A ). According to the above results of sequence alignment, we carried out conserved domain prediction ( Figure 3B ). NAM domain for DNA binding was detected in 112 NAC proteins. In addition, HcNAC16 had an incomplete subdomain named NAM superfamily. The results indicated that strong sequence conservation existed in their evolutionary process. The NAM domain was the core region for investigating the biological functions of NAC proteins in H. citrina.

An analysis of gene motifs may help to better understand the diversity of the HcNAC proteins. The 113 HcNAC TFs were divided into 16 subgroups in the NJ phylogenetic tree. Among them, the NAC-L subgroup was the highest in numbers and NAC-D was the lowest with only two members ( Figure 4A ), consistent with the clustered results in Figure 2 . As shown in Figure 4B , 10 different conserved motifs were obtained in HcNACs using MEME online software, and 53 of 113 HcNAC proteins contained 8 common motifs except for motifs 3 and 9, suggesting their important biological functions to be determined. However, motif 9 was exclusive to the NAC-F subgroup (HcNAC2, HcNAC27, HcNAC29, HcNAC60, HcNAC78, HcNAC80, and HcNAC81), indicating the specific functions of different subgroups might be owing to specific motifs. Intron/exon compositions were analyzed to gain the evolution information of HcNAC members ( Figure 4C ). Of the 113 HcNACs, more than 50% contained three code sequences. The number of intron regions was mainly 2 and varied from 1 to 12, proving that significant variation existed in the gene structure of HcNAC genes. Additionally, genes within each subgroup were similar in intron–exon structures. These results suggested that HcNACs clustered in the same subgroup shared similar conserved motifs and exon–intron organizations.

As effective methods for inferring the evolutionary history of species genome, intra- and inter-species collinear analysis are indispensable to address the study of NAC family expansions (Qiao et al., 2019). Through the annotation and intragenomic synteny analysis of HcNAC genes, 78 syntenic pairs were identified, among which two pairs of tandemly duplicated genes were detected ( Figure 5A ). It can be seen that tandem duplication and segment duplication might be major driving forces that form the expansion of the HcNAC gene family. Comparative syntenic maps of night lily with one monocot (O. sativa) and one dicot (A. thaliana) were constructed ( Figure 5B ). As a result, a comparison between night lily and O. sativa showed that 85 orthologous gene pairs were presented, whereas 21 HcNAC genes formed collinearity pairs with NAC genes from A. thaliana. The number of collinear NAC genes was higher in O. sativa than in A. thaliana, reflecting the closer relationship between night lily and monocots. The Ka, Ks, and Ka/Ks ratio, measurements of the protein conservation, were used to determine whether selective pressures occurred on genes encoding the HcNAC proteins ( Figure 5C and Supplementary Table S4 ). The Ka/Ks value of only one HcNAC gene pair was greater than 1, which signified strong positive selection. The remaining genes were subjected to purification selection considering Ka/Ks ratios less than 0.5. The results revealed that the HcNAC gene family was mainly affected by purifying selection in evolutionary selection.

Gene Ontology (GO) analysis of 113 HcNACs were subsequently performed to exposit biological processes, molecular functions, and cellular components ( Figure 6 ). The enrichment results showed that the cellular component was mainly the nucleus, accounting for 108 of 113. HcNAC34, HcNAC42, HcNAC87, HcNAC88, and HcNAC94 were not enriched. For the category of biological process, all other NAC proteins except HcNAC33, HcNAC34, and HcNAC42 of night lily were enriched in the regulation of transcription. Furthermore, some HcNAC proteins were enriched in the regulation of nucleic acid–templated transcription, regulation of RNA biosynthetic process, and response to abiotic stresses. The main enriched function was functional DNA binding, accounting for 109 of 113 total functions, and either HcNAC34, HcNAC42, HcNAC80, or HcNAC106 were enriched. The next most enriched function was DNA-binding TF activity (total of 11 of 113). Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis showed that most HcNACs except HcNAC106 were not involved in any pathway ( Supplementary Table S5 ). The HcNAC protein GO enrichment information was provided in Supplementary Table S5 .

Seven different tissues, including bud, tender leaf, mature leaf, tender root, mature root, tender scape, and mature scape were used to check the tissue-specific expressions of the identified 113 HcNACs. RNA-Seq data were shown in Supplementary Table S6 . Based on the analysis, the expression of HcNAC genes was detected in selected tissues, whereas their expression levels varied considerably ( Supplementary Figure S1 ). Except for eight of the HcNACs not having expression data, the expressions of remaining genes in roots were relatively higher than the other tissues. Furthermore, 21 of 105 genes showed high transcript abundances in all tissues. To validate the RNA-Seq result, 16 HcNACs were selected for RT-qPCR analysis ( Figure 7A ). The expression levels were consistent with those from RNA-Seq data ( Figure 7B ). The results described herein represented that HcNACs might perform distinct functions through different tissues and be worth further studies.

To elucidate the potential transcriptional regulation of HcNAC genes, cis-acting elements of the upstream promoter regions (2,000 base pairs) were classified and analyzed ( Figure 8 and Supplementary Table S8 ). A total of 115 abundant regulatory elements were detected via the online PlantCARE database, HcNAC111 hosted 295 cis-elements whose number was the largest ( Supplementary Table S7 ). We selected 21 typical responsive cis-elements to further unveil transcriptional regulation of HcNAC genes ( Figure 7 ; Supplementary Table S8 ). All HcNAC genes had regulatory elements linked to stress response, involving defense and stress, wound, drought, and low-temperature responsive elements. In addition, elements of CAT-box, circadian, GCN4_motif, and O2-site were identified as key regulars of plant development. Light-response covered TCT-motif, TGA-elements, G-box, and SpI elements. Hormone-related elements were associated with auxin, abscisic acid, gibberellin, salicylic acid, and MeJA, respectively, CGTCA-motif, ABRE, AuxRR-core, GARE-motif, and P-box. Our analysis revealed the presence of different types of cis-acting elements in the HcNAC gene family, which potentially regulate plant growth and development, and in response to hormone and abiotic stress.

To better illustrate the possible roles of HcNACs in abiotic stress, we collected night lily roots to carry out RNA-Seq sequencing to identify the expression patterns of HcNAC genes after drought (PEG) treatment for 0–108 h and salt (NaCl) treatment for 0–96 h. The expressions of night lily HcNACs under drought stress were shown in Supplementary Tables S9 , S10 under salinity stress. In 16 subgroups, a total of 12 HcNACs of night lily were detected with no expressions when suffering from drought while 11 were under salinity stress ( Figures 9A, B ). As the treatment time increased from 0 to 108 h under drought stress and 96 h under salinity stress, 10 HcNAC genes were upregulated and 5 genes were downregulated, among which HcNAC57 and HcNAC86 belonging to the NAC-A subgroup exhibited the opposite trends under two stresses. The results from RNA-Seq data revealed that the HcNAC genes were likely to function distinguishably in response to different abiotic stresses.

According to the analysis of cis-elements of upstream 2,000 bp sequence of night lily NAC genes ( Figure 8 and Supplementary Table S8 ), most of HcNACs in 16 subgroups contained six general stress response elements: MYC, stress response element (STRE), WUN-motif, MYB-binding sites (MBS), TC-rich repeat (cis-acting factor involved in defense and stress response), and low-temperature response motif (LTR). They accounted for 40%, 34%, 1%, 12%, 5%, and 8% of the total number of cis-elements identified in abiotic and biotic stress categories, respectively. In 16 genes which showed high expression in different tissues and abiotic stresses, the number of stress response elements ranged from 2 to 14. Among them, seven HcNAC genes (HcNAC34, HcNAC35, HcNAC45, HcNAC57, HcNAC65, HcNAC75, and HcNAC89) contain more stress response elements than other genes. HcNAC63 in the NAC-E subgroup did not contain MYC response element involved in abiotic stress response. Three genes (HcNAC71, HcNAC86, and HcNAC102) did not contain STRE response element. It was speculated that HcNACs might be involved in diverse abiotic stress responses and regulatory pathways.

Several NAC genes have been reported to regulate the growth and development of plants against abiotic stress (Huang et al., 2013). To gain insights into the putative functions of HcNAC genes in growth and development, the transcription levels of 16 genes, which showed significant expression disparities under PEG treatment, were scrutinized by RT-qPCR ( Figure 10A ). HcNAC65, HcNAC86, and HcNAC102 were downregulated. On the contrary, seven HcNAC genes (HcNAC34, HcNAC35, HcNAC47, HcNAC63, HcNAC75, HcNAC89, and HcNAC112) showed an increasing trend with time extension after treatment. Furthermore, HcNAC35 showed the most significant upregulation by drought stress. Therefore, HcNAC35 was considered a candidate gene for abiotic stress responses, and further studies of the gene expression during PEG treatment were needed. The results hinted at the consistency of RT-qPCR and transcriptome results under drought stress ( Figures 10A, B ), which provided a crucial reference for functional gene selection in the NAC gene family of H. citrina.

For salinity stress, we also performed RT-qPCR experiments to confirm the data accuracy of the selected 16 genes. As observed in Figure 11A , the expression of HcNAC57 displayed a decrease of twofold to threefold obviously while HcNAC67 was downregulated slightly, being suppressed at all processing times. In contract, eight genes (HcNAC35, HcNAC47, HcNAC65, HcNAC71, HcNAC75, HcNAC89, HcNAC102, and HcNAC112) exhibited upregulation after NaCl treatment. Notably, the expression of HcNAC35 and HcNAC71 were significantly increased by threefold to sevenfold at 24–96 h of treatment in response to salinity stress. Therefore, the consistency was equally manifested between RNA-Seq and qPCR results under salinity stress ( Figures 11A, B ), indicating that the sequencing data were highly reproducible. Overall, the expression patterns of HcNAC genes under different stress conditions proved that down- or up-expressed genes might engage in various stress responses.

Subcellular localization is known to be used to predict protein function. As one of the gene families classified into the TFs, most HcNAC genes in addition to HcNAC35 and HcNAC71 were predicted to locate in the nucleus ( Supplementary Table S2 ). To test the accuracy of predicted results, we selected HcNAC35 and HcNAC71, which exhibited distinct patterns in expression under drought and salt treatment, to construct HcNACs-GFP fusion proteins. As shown in Figure 12 , the results of the empty vector showed that strong fluorescence signals could be detected in the nucleus, cytoplasm, and cell membrane when pSuper1300-GFP was transformed. The signal location of HcNAC35-GFP and HcNAC71-GFP was only within the nucleus, which was the same as the predictions. It is speculated that HcNAC35 and HcNAC71 may play a role in the nuclear regulation as typical TFs, and the deep specific functions remain to be uncovered in future studies.

Previous experiments have implied that HcNAC35 exhibited significant upregulation under drought and salt treatments. To investigate the biological function of HcNAC35 in watermelon abiotic stress responses, we generated HcNAC35-OE plants in the TC watermelon genetic background. Based on expression levels of HcNAC35, two independent transgenic lines were chosen for further analysis. The RT-qPCR results showed that the transcript abundance of HcNAC35 in the OE-1 and OE-2 plants was about 4.82-fold and 6.25-fold higher than WT plants, respectively ( Figure 13A ).

To further assess whether alteration of HcNAC35 expression affects drought and salinity stress tolerance, we treated 4-week-old HcNAC35-OE plants and WT seedlings with 25% PEG and 250 mM NaCl, respectively ( Figure 13B ). The obvious wilting was observed in the shoot tips (STs) and leaves of WT seedlings after 12 h under 250 mM NaCl treatment, while STs and leaves of OE seedlings showed little damage. Additionally, visible damage was also caused in WT seedlings when 0–24 h under 25% PEG treatment. The phenotype observation indicated that HcNAC35 might play a critical role in response to abiotic stresses, especially salinity stress.

As indicators of membrane peroxidation, contents of MDA were measured to further evaluate the damage degree under salinity stress. MDA levels increased in all treated plants, but 30% and 38% lower in the OE plants than those of WT plants when 0–12 h after salt treatment ( Figure 13C ), indicating that the overexpression of HcNAC35 contributed to decreased damage under salinity stress. The levels of superoxide radicals (O₂ ⁻) were determined to explore the change of sensitivity caused by altered redox status in OE plants under salinity stress. HcNAC35-overexpressing plants showed lower O₂ ^- and H₂O₂ contents during salinity stress than the TC plants ( Figures 13D, E ). These results indicated that overexpression of HcNAC35 increased the tolerances of watermelon to salt and drought stresses to different degrees.