Cotton is a major economic crop that serves as a principal source of natural fiber and a raw material of oil. Cotton is also an ideal model plant for research on polyploidization. The Gossypium genus experienced two major events. In cotton ancestry, the A-genome diploids diverged from the D-genome diploids ~5–10 million years ago (MYA). Afterward, allopolyploid Gossypium species, including G. hirsutum, formed through the interspecific hybridization between the A-genome ancestor resembling G. arboreum and the D-genome ancestor resembling G. raimondii around 1–2 MYA (Li et al., 2015; Zhang et al., 2015). Cotton allopolyploidization produces thousands of duplicated genes with different expression levels (Hu et al., 2015; Wang et al., 2016). The allopolyploid cotton species G. hirsutum differs greatly from the putative donor species G. arboreum and G. raimondii in plant morphology and economic traits (Paterson et al., 2012; Li et al., 2014). Transcription factors play important regulatory roles in the aforementioned networks. Many transcription factors, such as TCP, WRKY, and MYB, regulate numerous critical biological processes in cotton (Hao et al., 2012; Yu et al., 2012; Lu et al., 2017).

The NAC (NAM, ATAF, and CUC) family is one of the largest families of plant-specific transcription factors (Ooka et al., 2003; Olsen et al., 2005). The first reported NAC gene (NAM) is related to the formation of the shoot apical meristem and primordium in Petunia (Souer et al., 1996). In general, the NAC family has a highly conserved N-terminal region (NAC domain) and a relatively divergent C-terminal transcriptional activation region (TAR) (Puranik et al., 2012). The NAC family regulates plant growth and development processes, including lateral root formation (Guo et al., 2005), leaf senescence (Guo and Gan, 2006), flower morphogenesis (Sablowski and Meyerowitz, 1998), cellular metabolism (Kim et al., 2009), seed development (Kim et al., 2008), fruit ripening (Shan et al., 2012), secondary wall synthesis (Mitsuda et al., 2005), and hormonal signaling (He et al., 2005). Moreover, NAC genes respond to many biotic and abiotic stresses, such as pathogen disease (Wang et al., 2009), drought (Mao et al., 2015), salt (Hu et al., 2006), and temperature (Fang et al., 2015). NAC transcription factors are also associated with crop yield and quality (Uauy et al., 2006; Liang et al., 2014; Zhao et al., 2015).

NAC genes have been studied in various plant species, such as Arabidopsis thaliana (Jensen et al., 2010), Oryza sativa (Nuruzzaman et al., 2010), Zea mays (Fan et al., 2014), Glycine max (Pinheiro et al., 2009), Solanum tuberosum (Singh et al., 2013), Musa acuminata (Shan et al., 2012), Eucalyptus grandis (Hussey et al., 2015), and Vitis vinifera (Wang et al., 2013). Several NAC genes have been isolated in cotton. GhNAP influences cotton yield and its fiber quality by regulating leaf senescence (Fan et al., 2015a), and GhATAF1 responds to salt stress and fungal infection by coordinating phytohormone signaling networks (He et al., 2016). However, genome-wide analysis of the NAC family in cotton is lacking.

The completion of genome sequencing in G. arboreum, G. raimondii, and G. hirsutum has opened an opportunity to investigate gene families in cotton (Wang et al., 2012; Li et al., 2014, 2015). A comprehensive analysis of the NAC family during allopolyploidization is expected to accelerate molecular breeding in cotton. Previous studies conducted whole-genome annotation of the NAC family in G. arboreum and G. raimondii (Shang et al., 2013, 2016) but only identified the NAC family in G. hirsutum by using EST scanning (Meng et al., 2009; Huang et al., 2013; Shah et al., 2013, 2014). Thus, a systematic research of the molecular evolution of the NAC family in cotton is needed. Genome-wide and comparative genomic analyses of the NAC genes in G. arboreum, G. raimondii, and G. hirsutum revealed asymmetric evolution and expansion during cotton polyploidization, as evidenced by their biased subfamily distribution, selective gene loss, unequal gene localization, and biased orthologous gene expression. This study unraveled the evolution of NAC genes in polyploid cotton and elucidated how NAC genes from the different progenitor genomes interact with each other.

With cotton being an important crop, sequencing of its genome facilitates our understanding of eudicot and polyploid plant evolution. The present research analyzed the NAC expansion and evolutionary history from the diploid to the allotetraploid cotton.

In this study, we identified 142 GaNACs, 142 GrNACs, and 211 GhNACs in cotton species (Table S1). The number of NAC genes was larger in each Gossypium than in three other eudicots (Table S2). This finding may be due to the extra whole genome duplication in the Gossypium lineage (Paterson et al., 2012; Li et al., 2015; Zhang et al., 2015). Although the A-related genome is larger than the D-related genome, the number of NAC members is similar in the two genomes. The reason may be related to more transposable elements in the A-related genome (Li et al., 2014, 2015). Through phylogenetic analysis, 15 NAC subfamilies were clustered in cotton (Figure 1, Figures S1–S3). The OsNAC7 and ONAC022 subfamilies had the most genes, whereas the AtNAC3 and ANAC063 subfamilies had the least NACs (Figure S4). Similar distribution of the NAC family existed in maize (Fan et al., 2014) and banana (Shan et al., 2012). The distribution of OG, gene structure, and conserved motif confirmed the similar classification (Figure 3, Tables S3, S4). Meanwhile, G. hirsutum contains almost twice as many genes in most NAC subfamilies as its diploid ancestors, with the exception of the AtNAC3, ANAC011, TIP, ONA003, and ANAC063 subfamilies. This indicates that G. hirsutum might have experienced gene loss in those latter subfamilies (Figure 2). This phenomenon may be related to the diploidization process after polyploidization (Otto, 2007; Soltis and Soltis, 2009). After polyploidization, the ployploid genome may reshuffle extensively to regain the diploid heredity (Wang et al., 2016). Moreover, the NAC family in cotton has a typical NAC structure (Figure 3A, Figure S5). The lengths of the first and second exons in the gene structure were relatively conserved partially because of the conserved encoding of the NAC domain (Figure 3B). However, exon length and intron insertions were different in the third exon, especially for the ANAC011, ONAC003, and TIP subfamilies. The gain and loss of introns might have resulted in the different structures of NAC genes in cotton. Besides, the putative NLS existed in subdomain D. Most of the NAC genes in cotton were predicted to be nuclear proteins through PSORT and Cello analyses (Table S5). This has been confirmed for some GhNACs through subcellular localization (Fan et al., 2015a; He et al., 2016). In addition, tissue-specific expression patterns were relatively conserved in the NAC subfamily. Different NAC subfamilies may have different biological functions, mainly due to the distinct structures, especially in the TAR region (Figure 4, Figure S5).

The expansion of NAC genes was investigated in cotton. We found 34 duplicated NAC gene pairs in cotton (Table 1). The duplication events had strong expansion preference for some NAC subfamilies, including ONAC022, NAM, and NAP. Thus, the duplication event is not random across NAC subfamilies during the cotton evolution. The retention may be related to the genomic fractionation in cotton (Wang et al., 2016). Meanwhile, orthologous genes of the duplicated NACs were isolated in five other eudicots (Table S8). Due to their similar duplication ages, some duplicated GhNACs in the SENU5, ONAC022, and NAP subfamilies might have originated via duplications in their diploid ancestors. However, we could not find their corresponding duplicated GhNACs in other duplicated GaNACs or GrNACs. Thus, the NAC subfamily might have various retention rates during whole genome reshuffling after interspecific hybridization. On the basis of the chromosomal location, we found that the main duplication events is different during the cotton evolutionary history (Figures S6, S7, Table 1). The segmental duplication of NAC genes dominated the expansion in the paleohexaploidization event and might be related to the paleopolyploidy (Li et al., 2015; Zhang et al., 2015). Most of the duplicated genes after the cotton-specific decaploidy resulted from tandem duplication events, and might result from the chromosomal breakages and rearrangements after polyploidization (Buggs et al., 2012). However, the segmental and tandem duplication events were highly prevalent after interspecific hybridization, and might be different from previously reported polyploidy-related duplication events. Moreover, the Ka/Ks ratios were <1 in all duplicated gene pairs, indicating that the NAC genes in cotton have mainly experienced purifying selection (Table 1). The expression profiles of some duplicated genes revealed similar results (Figures 5C–H). These results indicate that most of the duplicated genes of the NAC family might have retained some essential functions during sequent evolution.

Cotton experienced two major events as evidenced by the Ks distribution (Figure 6D). Two Ks peaks were between 0.045 and 0.042 and between 0.068 and 0.083, and their corresponding times may correspond to the interspecific hybridization and the divergence between their diploid ancestors. Our Ks values are a little higher than the previous reports (Li et al., 2015; Zhang et al., 2015). The main reason may be fewer orthologous gene sets in our study. Then, cotton has undergone two major whole genome duplication events: the paleohexaploidization event and the cotton-specific decaploidy (Paterson et al., 2012; Li et al., 2015; Zhang et al., 2015). In our study, 67% of the duplication events might have occurred in the aforementioned periods, expecially in G. raimondii and G. arboretum (Figure 5A, Figure S8). The whole duplication events resulted in the decaploid ancestor in cotton (Wang et al., 2016). Half of the duplication events in G. hirsutum might have occurred after interspecific hybridization. The results suggest that G. hirsutum might have experienced a complex and distinct evolutionary history from its diploid ancestors. Phylogenetic analysis also indicated that the duplication events of GaNACs and GrNACs might have mainly occurred in the paleohexaploidization event and the cotton-specific decaploid and that most of the duplicated GhNACs might have occurred after interspecific hybridization (Figure 5B).

G. hirsutum originated from the interspecific hybridization between AA-genome and DD-genome species (Li et al., 2015). In the present study, orthologous genes of GhNACs were identified in the diploid ancestors (Figure 6C, Table S7). Compared with NAC loss in the different genomes, the loss rate of NAC genes was higher in G. hirsutum than in diploid species. During the polyploidization, many chromosomal breakages and rearrangements led to gene loss and gene retentions (Buggs et al., 2012). The high loss rate may exist in other allopolyploid Gossypium species (Liu et al., 2015). In addition, more NACs were lost in the At subgenome than in the Dt subgenome during the formation of G. hirsutum, which led to more orthologous gene pairs and synetic blocks of the NAC family in the Dt subgenome and the DD genome (Figure 6C, Table S9). In the Dt subgenome, multilocus interactions are largely preponderant, which highlights the complex of NAC genes in cotton (Waghmare et al., 2016). Moreover, due to the different functions of the At and Dt subgenomes, the NAC family in G. hirsutum may be primarily associated with stress tolerance (Zhang et al., 2015). In previous studies, some GhNAC genes exhibited upregulated expression under stress (Meng et al., 2009; Shah et al., 2013, 2014; He et al., 2016). Thus, NAC genes in cotton have different loss rates in different genomes.

Orthologous exchange was identified in the At and Dt subgenomes (Figure 6B). A large number of GhNACs were transferred from Dt to At, which may increase the expression of NACs in the current At subgenome relative to the ancestral AA genome. Similar results were found in the Adh locus (Small and Wendel, 2002). The orthologous exchange of NAC genes might distinguish G. hirsutum from its ancestor species in important phenotypic changes, including plant morphology and economic traits. This exchange is also confirmed by the genetic linkage map (Li et al., 2015). Moreover, gene loss and orthologous exchange showed NAC subfamily preference, and different NAC subfamilies exhibited different rates of gene loss and orthologous exchange. Furthermore, Ka and Ks analyses of the NAC genes showed that the Dt subgenome might have evolved slower than the DD genome while the At subgenome and the AA genome might have had similar evolutionary rates (Figures 6E,F). Analysis of Ka/Ks ratio revealed that NAC genes might have experienced less positive selection in the At subgenome than in the Dt subgenome (Figure 6G). This result may partly explain why most At or Dt subgenome shows the highest orthologous with the corresponding diploid AA homologous chromosomes (Li et al., 2015) Thus, the At and Dt subgenomes had undergone asymmetric evolution after interspecific hybridization, which might be caused by different levels of selection pressure and introgressed chromatin (Waghmare et al., 2016). In addition, GhNAC redundancy created by allotetraploidy might have allowed relatively relaxed purifying selection in both At and Dt subgenomes. In the allotetraploid cotton, different genome-derived NACs possibly interacted with each other. Furthermore, most of the orthologous gene pairs between GhNAC and its orthologous GaNACs or GrNACs have the similar expression model in different tissues, which indicated limited functional divergence after interspecific hybridization (Figure 7). However, most of the NAC genes exhibited biased D-ortholog expression, which might have led to subfunctionalization of the NAC genes. This result suggests that the NAC genes in the DD genome or the Dt subgenome play important roles in plant development. The asymmetric expression was also observed in the LTP and miRNA families (Xie and Zhang, 2015; Li, F. et al., 2016; Li, X. et al., 2016).

In this report, we isolated 142 GaNACs, 142 GrNACs, and 211 GhNACs in cotton. The main duplication event of NAC genes was different during the cotton evolutionary history. The duplication events mainly occurred in the paleohexaploidization event (35%) and the cotton-specific decaploidy (32%). Moreover, some duplicated GhNACs might have originated from their diploid ancestor, and another might have occurred after interspecific hybridization. Meanwhile, 15 NAC subfamilies in cotton were clustered with NAC subfamily preference. In addition, NAC genes in the At and Dt subgenomes have different gene loss rates, orthologous exchange, evolutionary rates, and expression levels. Taken together, the findings of the present study could broaden our understanding on the molecular evolution and expansion history of the NAC family in cotton during polyploidization.

KF and WenL designed the research. KF, FL, JC, and ZL performed the experiments. KF, FL, JC, WeiL, and SC analyzed the data. KF, JL, and WenL wrote the paper with contributions from all the authors.