Identification of NAC domain proteins in rose
The HMMER (Prakash et al. 2017) and the Pfam NAC domain (PF02365) were employed to identify the rose NAC family using the protein sequences obtained from the Arabidopsis NAC gene sequence. A total of 123 protein sequences were initially obtained by comparing and searching the rose genome database (https://www.rosaceae.org). A two-way comparison between TBtools (Chen et al. 2020) and the NCBI database was used to detect the NAC domain after removing the redundant sequences. A total of 116 RcNACs were identified and named from RcNAC001 to RcNAC116 according to their chromosome location. These 116 RcNACs were analyzed for their gene location, gene length, molecular weight (Mw), isoelectric point (pI), and subcellular localization (Table S2). The protein length of RcNACs ranged from 70 (RcNAC076) to 714 amino acids (aa, AA) (RcNAC102), with an average of 336 aa (Table S1). The pI ranged from 4.43 (RcNAC107) to 9.78 (RcNAC98), whereas the Mw ranged from 7.87 kDa (RcNAC076) to 79.99 kDa (RcNAC102). Subcellular localization analysis predicted nuclear localization for 49% of RcNACs and extracellular localization for 51%. These results suggest the crucial roles of RcNACs in the evolution of the rose genome.
Phylogenetic analysis of NAC proteins
A phylogenic tree was constructed using 736 NAC protein sequences from five species, including rose (116), kiwifruit (156), Arabidopsis (105), apple (228), and Oryza barthii (131). Based on the classification of NAC family proteins in Arabidopsis by Ooka et al. (2003), the NACs were divided into 16 subfamilies, including ANAC011, NAM, NAC1, OSNAC7, OSNAC8, TIP, NAC2, ONAC22, ANAC001, SEUN5, ONAC003, TERN, ANAC063, ATAF, NAP, and Unclassified. The ANAC001 subfamily was the largest with 23 members (~20% of the total RcNACs) in rose, followed by the NAM subfamily containing 11 RcNACs; OSNAC8 and ANAC063 subfamilies were the minor categories with only one member each. However, no plant species had all the subfamilies, and no OSNAC8 member was found in Arabidopsis (Table S3). Furthermore, the percentage of each subfamily was different among the plant species. For example, monocots like O. barthii (19%) showed a more significant percentage of the TIP subfamily than the dicots likeroses(4%), kiwifruit(6%), Arabidopsis(7%), andapple(4%) (Fig. 1). These findings indicate potential differences in levels of gene duplication or loss between monocots and dicots that might have occurred during evolution.
Gene structure and motif analysis of RcNACs
Further, the conserved motifs, exon-intron structure, and phylogeny were analyzed to understand better the structural diversity in rose NAC proteins (Fig. 2). A total of ten conserved motifs, mainly localized in the N-terminus, were identified for NAC proteins (Fig. S1). Motifs 1–7 were highly conserved and distributed in most rose NACs (~60%), while motifs 8–10 were distributed in a few subfamilies (ANAC001 and Unclassified). Most of the conserved motifs of rose NACs of the subfamilies, such as ANAC011, TIP, and NAC2, were similar, indicating a relatively close evolutionary relationship. However, the conserved motifs of RcNACs in subfamilies ANAC001, ONAC003, and ANAC063 were different, indicating a rather distant evolutionary relationship.
Meanwhile, the exon-intron organization showed that the intron number was diverse among RcNACs, ranging from one to nine (Fig. 2). All 116 RcNACs consisted of one or more exons. The maximum number of exons (seven) was found in RcNAC102 of the TIP subfamily, and the average number of exons in the 116 RcNACs was four. Furthermore, the genes in the same subfamily had similar gene structures, with similar intron and exon arrangement and numbers, indicating highly conserved biological functions of these genes.
Chromosome localization and collinearity analysis of RcNACs
Chromosome mapping revealed that the RcNACs were unevenly distributed within the rose genome composed of seven chromosomes. Among them, chromosome 1 contained the maximum RcNACs (31, 26.7%), followed by chromosome 5 (21, 18.1%); the least number of RcNACs was on chromosome 4 (6, 5.2%; Fig. 3A).
Furthermore, a collinearity analysis was performed using NACs from rose and four other Rosaceae plants, including apple, strawberry, mei, and almond. As shown in Fig. 3B, the collinearity between the NACs in rose and apple (142) was significantly higher than the other three species (mei, 79; strawberry, 66; and almond, 40). The RcNACs on chromosomes 5 and 7 had significantly more collinear gene pairs than the other chromosomes. Less than ten collinear gene pairs were identified between rose and almond, which may be related to the phylogenetic relationship between species (Fig. 3B). These results suggest that few RcNACs were possibly produced by gene and segmental duplication events, which probably played a major driving force for RcNAC gene evolution.
Relative transcript levels of stress-related RcNACs in different transcriptome datasets
NAC proteins play various roles in plant growth, development, and stress response. In this study, 44 stress-related RcNACs belonging to eight subfamilies (ATAF, NAC2, TIP, SNUE5, NAM, NAP, ONAC022, ONAC003) were selected for further analysis. Sequence homology of all these genes with already reported stress-related NACs indicated their close association with biotic or abiotic stresses (Shao et al. 2015). The expression levels of 44 RcNACs at the three stages of the rose flower transition (vegetative meristem, pre-floral meristem, floral meristem) and in the six organs (prickle, stamen, leaf, pistil and ovary, stem, and root) based on public RNA-seq data are shown in Fig. 4A. The transcript levels of RcNACs of the same or different subfamilies were highly variable among the organs or the transition stages, suggesting their multiple functions in growth and development. Interestingly, the expression levels of four NAC2s (RcNAC008, RcNAC086, RcNAC009, RcNAC088), three ATAFs (RcNAC034, RcNAC061, RcNAC091), two TIPs (RcNAC031, RcNAC102), and one SNEU5 (RcNAC056) in the various organs and flower transition processes were relatively high, demonstrating their roles in rose growth and development.
To further gain insight into the putative functions of these 44 RcNACs in response to abiotic stress, their expression patterns under drought stress in leaves (no drought, NDL; mild drought stress, MDL; severe drought stress, SDL) and roots (NDR, MDR, SDR) and salt stress in leaves (TL0, TL2, TL4, TL7, TL9) and roots (TR0, TR2, TR4) were evaluated (Fig. 4B). Overall, the expression levels of the RcNACs from eight subfamilies were different. These genes showed dramatic changes under various degrees of drought or salt stress in rose leaves and roots, further confirming their roles in response to abiotic stress. Meanwhile, 11 RcNACs of the same clade showed induced expression under drought or salt stress in rose leaves and roots. Specifically, three ATAFs (RcNAC034, RcNAC069, RcNAC091) showed higher expression levels under drought or salt stress. These expression profiles demonstrated functional convergence and differences of the RcNACs, which probably occurred during evolution.
Analysis of the cis-regulatory element (CRE) in stress-related RcNACs promoter
To analyze the regulation of stress-related RcNAC gene expression, the CREs of its promoter, a region of approximately 1.5 kb upstream of the initial transcript site, were tested. Li et al. (2021) classified the promoters into three categories based on functions, corresponding to plant growth and development (Box4, G-box, CAT-box, AE-box, GT1-motif, GATA-motif, I-box), plant hormone response (ABRE, CGTCA-motif, ERE, TGACG-motif, GARE-motif, TGA-element), and abiotic and biotic stress (TC-rich repeats, MBS, LTR, WUN-motif, W box, WRE3, STRE, MYB, MYC, ARE, Myb-binding site, TCT-motif).These CREs were randomly dispersed among the different RcNACs promoters. RcNAC034 had the most significant number of CREs (57), followed by RcNAC091 (52), while RcNAC081 had the least number of CREs (6). Among the different types of CREs, abiotic and biotic stress-related CREs, including MYB, MYC, Box4, G box, and ABRE, were significantly more than the other CREs (Fig. 5); 31 out of 44 RcNACs contained ABRE elements, suggesting RcNACs roles in regulating the response to abiotic stresses through ABA-dependent pathway.
Expression of selected RcNACs
According to the transcriptome expression and CREs analysis, 11 representative genes (RcNAC009, RcNAC031, RcNAC034, RcNAC056, RcNAC059, RcNAC062, RcNAC069, RcNAC086, RcNAC088, RcNAC091, RcNAC102) belonging to six subfamilies of the same clade (Fig. 4B) were selected for further analysis using RT-qPCR under ABA, PEG, and NaCl treatments (abiotic stress) in rose leaves and roots. As shown in Fig. 6, members of the different subfamilies showed differential expression patterns. RcNAC034 and RcNAC069 showed upregulated expression levels when treated with ABA in leaves and roots, while RcNAC059 and RcNAC031 did not. PEG upregulated RcNAC056 in leaves, and ABA upregulated in roots. The expression levels of both RcNAC031 and RcNAC102 showed no change under ABA, PEG, and NaCl treatments. However, RcNAC091 showed upregulated expression levels when treated with ABA, PEG, and NaCl. These expression profiles of RcNACs from different subfamilies indicated conserved or diverse roles in rose leaves and roots in response to external cues.
Silencing RcNAC091 decreases tolerance to dehydration stress
RcNAC091 showed an stongly induced expression pattern in response to abiotic stress (Fig. 6). Therefore, VIGS was used to silence the selected RcNAC091 for functional analysis. As shown in Fig.7, about 40% of RcNAC091 was silenced in rose (Fig. 7A). No obvious difference in phenotype was observed between TRV and TRV-RcNAC091 under dehydration 0 h. However, when dehydrated for 4 h and 8 h, discs of TRV-RcNAC091 showed more wilting symptoms than TRV controls (Fig. 7B). The discs areas of TRV-RcNAC091 under dehydration for 4 h and 8 h and rehydration for 3h were 510 mm², 395 mm², 556 mm², respectively, which were significantly lower than TRV controls at the indicated time points (574 mm², 459 mm², 615 mm²) (Fig. 7C). Further, TRV-RcNAC091 samples after rehydration for 3 h showed lower chlorophyll content (Fig. 7D), with a higher ion leakage than TRV controls (Fig. 7E). Furthermore, TB, DAB, and NBT staining used to investigate the cell death rate and ROS accumulated levels revealed deeper staining for TRV-RcNAC091 than TRV controls under rehydration for 3 h (Fig. 7F), indicating severely damaged cell membranes. These results suggest that RcNAC091 positively regulates dehydration stress response in rose plants.