The genome data and annotation files of chestnut (N11-1) were downloaded from the Castanea Genome Database (http://castaneadb.net/), and the published NAC gene family members in A. thaliana were downloaded from TAIR (https://www.arabidopsis.org/) (Jensen et al., 2010). The Hidden Markov Model (HMM) file of the NAC domain (PF02365) was obtained from the Pfam database (http://pfam-legacy.xfam.org/). Firstly, the protein sequences of the NAC genes in A. thaliana were used as reference sequences to perform a Basic Local Alignment Search Tool (BLAST) program on all protein sequences of chestnut with the following parameters: evalue 1e-5 -outfmt 6. Subsequently, the HMM file of the NAM domain was used to search against the proteins of chestnut using HMMER 3.0 with an E-value of 1e-5. Then, we used Batch-CDD to confirm that the NAM domains existed in the candidate genes. Finally, 115 NAC genes were identified in the chestnut genome. ExPasy website (http://web.expasy.org/protparam/) was used to calculate the physicochemical properties of NAC proteins in chestnut.

A total of 209 NAC protein sequences from chestnut (115) and A. thaliana (94) (Jensen et al., 2010) were merged for phylogenetic analysis. MEGA 7.0 (Kumar et al., 2016) was used to construct the phylogenetic tree of NAC genes with the maximum likelihood. The “Find Best DNA/Protein Models (ML)” function in MEGA7.0 was used to find the best amino acid substitution model (partial deletion 95%), and the final parameters were as follows: Jones–Taylor–Thornton (JTT) model; Gamma Distributed (G); Partial deletion 95%; 1,000 bootstrap replications. Additionally, the neighbor-joining method was also used to construct a phylogenetic tree with the following parameters: Poisson model, pairwise deletion, and 1,000 bootstrap replications. MEME websites (http://meme.nbcr.net/meme/intro.html) were used to predict the conserved motifs of CmNAC proteins, with the following parameters: the number of repetitions, zero or one; and the maximum number of motifs, 10. The CmNAC proteins were submitted to NCBI-CDD Search to predict their conserved domains with default parameters. TBtools (Chen et al., 2020) was used to visualize the gene structure, conserved domains, and motifs of CmNAC members. We also constructed phylogenetic trees containing only members of the CmNAC gene family with methods the maximum likelihood method and neighbor-joining method, using the same strategy and method.

The position information and gene density file of the target sequence were extracted using TBtools software, and visualized using the Gene Location program. The Multiple Collinearity Scan toolkit (MCScanX) software was used to conduct collinear analysis of the genomes of chestnut and A. thaliana, rice, oak (Quercus robur), and grape (Vitis vinifera) to explore the evolution of CmNAC genes, with default parameters (Wang et al., 2012). Notably, we identified the CmNAC genes formed by whole-genome-duplication (WGD) events, as we did before (Yu et al., 2023b). Specifically, we drew the homologous collinear gene dot-plot within the chestnut genome. The “add_ka-and_ks_to collinearity” in MCScanX was used to obtain the non-synonymous (Ka) and synonymous substitution sites (Ks) values of homologous gene pairs, and the median Ks values of collinear blocks were calculated by writing the script (Yu et al., 2023b). The collinear blocks in the homologous gene dot-plot were colored differently based on different median Ks values, which will help to locate the blocks formed by different duplication events. Finally, based on the distribution of Ks corresponding to the WGD event that occurred in the chestnut genome before (Yu et al., 2023b), combined with the complementarity of the collinear blocks, we identified the CmNAC genes formed by the WGD event. In addition, we also obtained CmNAC members from other duplication models, such as segmental, proximal, dispersed, and tandem duplication, from the result file of collinear analysis of the chestnut genome. The 2000 bp nucleic acid sequences upstream of the transcription initiation of 115 CmNAC genes were submitted to the PlantCARE website to predict the composition of cis-acting elements, and the results were classified into functional categories.

The publicly available RNA-seq data from different tissues of chestnut were used to study the expression of CmNAC genes and explore their potential functions. Specifically, RNA-seq data of chestnut buds of 20, 25, and 30 days post-anthesis, nuts of 70, 82, and 94 days post-anthesis, and fertile/abortive ovules on 15-July, 20-July, and 25-July were analyzed. The accession numbers of the above RNA-seq data were in Supplementary Table S1, and all of these were three biological replicates. We used the Kallisto software to quantify RNA-seq data into Transcripts Per Kilobase of exon model per Million mapped reads (TPM) (Bray et al., 2016). The “Normalized” function in TBtools was used to normalize gene expression, and we obtained the heat-maps of gene expression based using TBtools.

We collected the nuts of “Yanlong” chestnuts 70, 82, and 94 days after flowering for real-time quantitative PCR experiments, which were stored at −80°C. Specifically, we used the RNAprep pure Plant Kit to extract and isolate RNA, and used the PrimeScript RT Master Mix (Takara, Beijing) to reverse transcribe RNA into single-stranded cDNA. ABI 7500 Real-Time PCR system (Applied Biosystems Inc., Foster City, CA, United States) with TB Green Premix Ex Taq (Takara) was used to conduct the RT-PCR experiments. The specific primer information was shown in Supplementary Table S2, in which the 18S gene of chestnut was used as the reference gene. The relative gene expression values were calculated using the comparative 2^-△△CT method (Livak and Schmittgen, 2001).

After strict screening by BLASTP and HMMER 3.0 programs, and using Batch-CDD to ensure the existence of the NAM domain, we finally identified 115 NAC gene family members in the chestnut genome. In addition, 115 CmNAC genes were renamed CmNAC1 to CmNAC115 based on their relative position. CmNAC1∼CmNAC113 were unevenly distributed on Chr1∼Chr12, while CmNAC114 and CmNAC115 were distributed on two scaffolds not attached to chromosomes (Figure 1; Supplementary Table S3). CmNACs had the largest number of genes distributed on Chr1 and Chr10, with 20 and 16 members, respectively. CmNACs had the least number of genes distributed on Chr11 and Chr12, with only three members. The amino acids number, molecular weight, instability index, and pI of CmNAC gene family members varied greatly (Supplementary Table S3). Specifically, the number of amino acids in CmNACs ranged from 93 (CmNAC64) to 934 aa (CmNAC60), and the molecular weight was between 10,992.77 Da (CmNAC64) and 106,329.77 Da (CmNAC60). The pI values ranged from 4.26 (CmNAC89) to 10.08 (CmNAC64), of which 83 CmNAC proteins were acidic, 31 were basic proteins, and CmNAC45 was a neutral protein. The instability index values of 25 CmNAC proteins were less than 40, and the protein properties were stable, accounting for approximately 21.74% of CmNAC proteins. CmNAC proteins were all hydrophilic proteins, due to the grand average of hydropathicity (GRAVY) values being less than zero. In summary, it is inferred that most of the CmNAC family proteins were unstable alkaline hydrophilic proteins. The subcellular localization predicted that CmNAC proteins were distributed in the nucleus, chloroplast, cytosol, Golgi apparatus, peroxisome, extracellular, mitochondrion, and endoplasmic reticulum, with 69.57% of the proteins located in the nucleus and a few CmNAC proteins located in other organelles (Supplementary Table S3).

To clarify the genetic relationships among members of the CmNAC genes, a phylogenetic tree with the maximum likelihood method was constructed using a total of 209 NAC protein sequences from chestnut and A. thaliana (Figure 2). The results showed that 115 CmNAC proteins were divided into 16 subgroups, namely, ONAC003, ANAC001, SENU5, NAP, AtNAC3, ATAF, TERN, ONAC022, NAC1, NAM, OsNAC7, ANAC011, NAC2, OsNAC8, TIP, and 29 unclassified CmNAC proteins, based on their homology with NAC proteins in A. thaliana (Ooka et al., 2004). Interestingly, no CmNAC protein was classified into the ANAC063 subgroup. Except for the UN subgroup, the TIP and NAP subgroups contained the largest number of CmNAC genes, with 17 and 15 members, respectively. The OsNAC8 and AtNAC3 subgroups contained the least number of CmNAC genes, both with only one member. Furthermore, we also constructed phylogenetic trees of NAC genes in chestnut and A. thaliana with the neighbor-joining method (Supplementary Figure S1). These results indicated that CmNAC proteins exhibit a certain degree of diversity, which was similar to the phylogenetic results of NAC proteins reported in many species (Jin et al., 2020; Jia et al., 2021; Du et al., 2022).

We analyzed the gene structure, motif distribution, and promoter cis-elements to further understand CmNAC genes. The phylogenetic tree constructed using only 115 CmNAC protein sequences was consistent with the phylogenetic tree constructed using 209 NAC protein sequences from chestnut and A. thaliana (Figure 2; Supplementary Table S1–3). Specifically, Gene Structure Display Server 2.0 (GSDS) and MEME (http://meme.nbcr.net/meme/intro.html) were used to analyze the distribution of exons and introns in the CmNAC genes, and conserved motifs (Figure 3). The structure analysis of the CmNAC genes showed a high degree of difference in the number of introns and exons, with approximately 1–10 introns present, and the gene structure of CmNAC transcription factors in the same subfamily was relatively conserved. In addition, the introns of most CmNAC genes were phase zero, and this further indicated that they were conserved in gene structure. Diversified motif compositions were detected in the protein sequence of the CmNAC genes. Interestingly, some motifs can be detected in almost all CmNACs, which may be related to the shared domains of the family. In addition, CmNAC members within the same subgroup had a relatively consistent motif distribution. For example, members in the NAC1 subgroup had a completely consistent motif distribution, while motif 9 mainly existed in the TIP subgroup. Some subgroups had individual members missing motifs that other members all contained. For example, the CmNAC115 in the OsNAC7 subgroup lacked motif 3, which was preserved in other members of the OsNAC7 subgroup (Figure 3B). Indeed, the specific functions of these motifs required further in-depth research (Yang et al., 2023).

The cis-acting elements analysis was necessary to understand the functions of CmNAC genes. A total of 2,887 cis-acting elements were identified by predicting the upstream 2000 bp sequences of 115 CmNAC gene promoters, and these cis-acting elements involve various types of functions (Figure 4; Supplementary Table S4). For example, TGA-element or TGA-box related to auxin response was identified in the promoter regions of 31 CmNAC genes. Gibberellin-responsive elements were identified in the promoter regions of 56 CmNAC genes, and cis-regulatory elements involved in endosperm expression were identified in the promoter regions of 20 CmNAC genes. These results suggested that members of the CmNAC gene family were closely related to the growth and development of chestnut. Additionally, a large number of cis-acting elements related to plant resistance to stress were identified, such as wound-responsive elements, cis-acting elements involved in defense and stress responsiveness (Figure 4; Supplementary Table S4) (Sun et al., 2022).

Gene duplication is one of the reasons for the formation of gene families and the diversification of gene functions (Liu and Widmer, 2014; Quan et al., 2019; Yu et al., 2023a). We conducted collinear analysis on the chestnut genome using MCScanX to explore duplication events among members of the CmNAC gene family (Figure 5). Fifty-three dispersed duplicate genes were found among 115 CmNAC genes (Supplementary Table S5). In addition, five and four CmNAC genes were identified to be proximal and tandem duplications, respectively. Interestingly, 32 of the CmNAC genes were considered WGD or segmental duplication. A homologous dot-plot of the chestnut genome and color distinguished them based on the Ks value of homologous genes was drawn to preliminarily determine homologous genes from different duplication events (Figure 6). Furthermore, referring to our previous distribution of Ks values of WGD of the chestnut genome and the complementary relationship between homologous gene fragments (Yu et al., 2023b), we ultimately confirmed 9 pairs of 18 CmNAC genes from the chestnut WGD event. The final results indicated that dispersed duplication and WGD were the main drivers of CmNAC gene expansion (Supplementary Table S5).

The collinear relationships between chestnut and four representative species, namely, A. thaliana, rice, grape, and oak, were analyzed to further understand the evolution of the CmNAC genes (Figure 7). The results showed that 65 collinear blocks were found between the genomes of chestnut and grape, with an average block length of 49.23 (Supplementary Table S6). The number of collinear blocks found between the genomes of chestnut and A. thaliana, rice, and oak was 66, 28, and 38, respectively, with average block lengths of 21.08, 10.61, and 26.29, respectively (Supplementary Table S7–9). Furthermore, we have discovered 71 orthologous gene pairs containing a total of 50 CmNAC genes between the genomes of chestnut and grape (SSupplementary Table S10). However, only 42 (73 pairs), 19 (29 pairs), and 36 (39 pairs) CmNAC genes were found between chestnut and A. thaliana, rice, and oak genomes, respectively (Supplementary Table S11–13). These results indicated that better collinearity was preserved between the genomes of chestnut and grape. In addition, we calculated the Ka/Ks values of gene pairs formed by different duplication models (such as tandem, WGD, and segmental duplication), which suggested that they experienced purification selection pressure during their evolution process (Ka/Ks < 1) (Supplementary Table S14) (Anisimova et al., 2001; Roth and Liberles, 2006).

To further understand the possible functions of NAC genes in different tissue development processes of chestnuts, we fully utilized RNA-seq from the NCBI database of chestnut ovules (fertile/abortive), buds, and nuts (Figure 8). Some CmNAC genes exhibited high expression levels throughout the entire developmental stage of fertile ovules, while they were almost not expressed throughout the entire developmental process of abortive ovules (Figure 8A). For example, the TPM values of CmNAC113 in the three developmental stages of fertile ovules were 20.71, 13.17, and 44.92, respectively, while the TPM values were 5.82, 1.70, and 2.15 in the same stage of abortive ovules, respectively. On the contrary, some CmNAC genes were almost not expressed throughout the entire developmental stage of fertile ovules, but exhibited higher expression levels in abortive ovules. For example, the TPM values of CmNAC109 in the three developmental stages of fertile ovules were 1.149, 0.67, and 0.19, respectively, while the TPM values were 26.74, 7.12, and 25.76 in the same stage of abortive ovules, respectively. Notably, some cis-regulatory elements related to endosperm expression were found in the promoter regions of these CmNAC genes (Figure 4). These results indicated that these NAC genes in chestnut were related to the fertility of chestnut endosperm.

Furthermore, we analyzed the expression levels of CmNAC genes in chestnut buds and nuts at different developmental stages (Figure 8B). Interestingly, we found that some CmNAC genes maintained high expression levels, and their expression levels continued to increase, during chestnut bud development. For example, the expression level of CmNAC14 in chestnut buds 20 days after flowering was 66.12 TPM, while the expression levels in chestnut buds 25 and 30 days after flowering were 126.49 and 148.57 TPM, respectively. The expression level of CmNAC38 in chestnut buds 20 days after flowering was 90.05 TPM, while the expression levels in chestnut buds 25 and 30 days after flowering were 427.51 and 871.72 TPM, respectively. Meanwhile, gibberellin-responsive element and cis-acting regulatory element related to meristem expression were identified in the promoter regions of CmNAC14 and CmNAC38 (Figure 4). Notably, gibberellin is a very important plant hormone that participates in many biological processes, such as the development and differentiation of plant buds (Richards et al., 2001). In addition, during the development of chestnut nuts, some CmNAC genes remained highly expressed, such as the TPM values of CmNAC49 at 70, 82, and 94 days after flowering of chestnut nuts, which were 122.94, 166.93, and 123.73, respectively (Figure 8C). Some CmNAC genes exhibited sharp changes in expression with the development of chestnut nuts. For example, CmNAC54 was almost not expressed in nuts at 70 and 82 days after flowering (TPM values of 0.17 and 0.19, respectively), while the expression level increased sharply at 94 days after flowering (chestnut maturity), reaching a TPM value of 72.39. Interestingly, cis-acting elements involved in abscisic acid and MeJA responsiveness were identified in CmNAC49 and CmNAC54 (Figure 4). The expression patterns of 12 CmNAC genes that we were interested in during the development of chestnut nuts were validated using qRT-PCR experiments (Figure 9A). These results showed that the qRT-PCR experiment results and RNA-seq analysis showed consistent expression patterns of CmNAC genes. For example, the expression level of CmNAC38 continuously increased until it reached its highest level at maturity (94 days after flowering). CmNAC47 has the lowest expression level on 70 days after flowering, the highest expression level on 82 days after flowering, and a significant decrease in expression level at maturity (94 days after flowering). These results suggested that some CmNAC genes may play important roles in the development of chestnut buds and nuts. Additionally, to determine the subcellular localization of CmNAC genes through experiments, we randomly selected three CmNAC members, namely, CmNAC49, CmNAC74, and CmNAC105, and instantaneously expressed GFP-CmNAC49, GFP-CmNAC74, and GFP-CmNAC105 fusion proteins in A. thaliana protoplasts for subcellular localization (Figure 9B). Confocal microscopy analysis showed that the GFP-fused CmNAC49 and CmNAC74 displayed fluorescence signals distributed in the nucleus, and GFP-fused CmNAC105 was localized to the nucleus and peroxisome. These results were consistent with the subcellular localization prediction of CmNAC through bioinformatics analysis.