Identiﬁcation and Expression Analysis of the NAC Gene Family in Co ﬀ ea canephora

: The NAC gene family is one of the largest families of transcriptional regulators in plants, and it plays important roles in the regulation of growth and development as well as in stress responses. Genome-wide analyses have been performed in diverse plant species, but there is still no systematic analysis of the NAC genes of Co ﬀ ea canephora Pierre ex A. Froehner. In this study, we identiﬁed 63 NAC genes from the genome of C. canephora . The basic features and comparison analysis indicated that the NAC gene members increased via duplication events during the evolution of the plant. Phylogenetic analysis divided the NAC proteins from C. canephora , Arabidopsis and rice into 16 subgroups. Analysis of the expression patterns of CocNACs under cold stress and co ﬀ ee bean development indicated that 38 CocNACs were di ﬀ erentially expressed under cold stress; six genes may play important roles in the process of cold acclimation, and four genes among 54 CocNACs showing a variety of expression patterns during di ﬀ erent developmental stages of co ﬀ ee beans may be positively related to the bean development. This study can expand our understanding of the functions of the CocNAC gene family in cold responses and bean development, thereby potentially intensifying the molecular breeding programs of Co ﬀ ea spp. plants.


Introduction
The NAC transcription factor, named after the proteins of NAM (Petunia No Apical Meristem), ATAF1/ATAF2 (Arabidopsis transcription activation factors), and CUC2 (Cup-shaped cotyledon), is one of the largest plant-specific transcription factors (TFs) [1][2][3]. The NAC domain (around 160 amino acids) is located in the N-terminal region of the NAC family proteins and can be divided into five subdomains (A to E) that provide the DNA binding ability of NAC TFs [4,5]. In addition to the conserved NAC-domain, a variable transcription regulatory region was found in the C-terminal of NAC family proteins; this region is suggested to produce regulation diversity of transcriptional activation activities [4].
The quality of the coffee depends on the biochemical composition of the coffee bean, which is controlled by genetic and environmental factors influencing the accumulation of key components of the bean [32][33][34]. To investigate genetic control of these processes, studying the changes in transcripts during bean development as well as in response to stress are required, but only ESTs (Expressed Sequence Tags) are available [35]. However, the possibility that the NAC family is associated with bean development arises from studies of other plants. Nacregulated seed morphology 1 and -2 (NARS1 and NARS2, also known as ANAC056 and ANAC018, respectively) regulate embryogenesis by regulating the development and degeneration of ovule integuments in Arabidopsis [36]. At least three soybean NACs (GmNAC1-GmNAC3) are expressed at particular stages of seed development [37]. Two maize NACs (ZmNAC128 and ZmNAC130) regulate carbohydrate and protein levels in seeds [38]. Recently, it was found that six of 112 strawberry NAC genes are positively correlated with fruit ripening: FaNAC006, FaNAC021, FaNAC022, FaNAC035, FaNAC042, and FaNAC092 [39]. In the Poaceae family of monocotyledonous plants, the six-barley d-9 NACs in a subclade are specifically expressed in the developing grain [40], and wheat [41], rice [42] and maize homologs [43] fall into the same clade and are expressed specifically in developing grains. Therefore, the study of CocNACs in relation to coffee bean development will provide insights into the regulation of coffee bean development.
Coffee beans are the seeds of the coffee plant and the source for coffee, one of the major commodities in the word. Coffea canephora (diploid, 2n = 2x = 22 chromosomes, 710 Mb) (Rubiaceae), known as robusta coffee, provides approximately 30% of the coffee beans for the worldwide coffee market [44,45]. C. canephora is considered as one of the parents of Coffea arabica (tetraploid), which provides about 70% of the coffee bean production in the world [46]. C. canephora can provide intense flavor and bitterness as well as greater abiotic and biotic resistance compared to C. arabica [47,48]. To produce a good quality beverage, coffee cultivation has spread towards the highlands in both Brazil and China, where the temperature occasionally reaches 0 • C [45,49]. In addition to cold tolerance, another important trait, coffee bean development, is poorly understood with few published studies [34,50,51]. In order to examine both cold tolerance and coffee bean development, we carried out RNA_seq during various developmental stages of coffee beans, and we performed genome-wide analysis of C. canephora NAC TFs (CocNACs). Our results may provide insights regarding the molecular significance of CocNACs under cold stress and in the development of coffee beans.

Plant Material Collection and Cold Treatment
Following the previous research [46,52], fruits of Coffea canephora at different stages of development (SG (small green fruit), LG (large green fruit), Y (yellow fruit), and R (red fruit)) were harvested from coffee trees cultivated in the open field of the Germplasm Repository of Coffee (Coffea spp.), RuiLi City, Ministry of Agriculture (RuiLi, Yunnan, China). Fruits were collected from three harvests (2016-2018), and 50 fruits were sampled from five independent trees for each of the differential developmental stages. As soon as sampling was finishing, the coffee beans were isolated and frozen immediately using liquid nitrogen, then stored at −80 • C until use.
For cold treatment, one-year-old seedlings of C. canephora were used. The seedlings were provided by the Germplasm Repository of Coffee (Coffea spp.), RuiLi City, Ministry of Agriculture (RuiLi, Yunnan, China). Briefly, the coffee plants were grown under a long-day cycle (16 h light/8 h cark) in a greenhouse at 24 • C before cold treatment. Cold treatments were carried out as described previously description [45] with some modifications. The plants were submitted to 24 • C/20 • C (day/night) for seven days in a growth chamber, followed by seven days at 13 • C/8 • C (day/night) to express cold acclimation ability, and then treated with 4 • C/4 • C (day/night) for three days. During the cold treatment, the following parameters were fixed: humidity (60%), luminosity (600-650 µmol m −2 s −1 ), and photoperiod (16 h light/8 h dark). Five plants were taken out at the end of each treatment time point, and the two top pairs of recent mature leaves were sampled. After sampling, the leaves were frozen immediately in liquid nitrogen and then stored at −80 • C until use. The cold treatments were carried out with three independent biological replicates.

RNA Extraction from Coffee Beans, Library Construction, Sequencing and Assembly
Total RNA was isolated from each coffee bean sample by using the RNAprep Pure Plant Plus Kit (TIANGEN BIOTECH CO., LTD, Beijing, China) according to the manufacturer's protocols. Total RNA (20 µg) was sent to BGI (Beijing Genomics Institute, Shenzhen) where the libraries were constructed and sequenced using a BGISEQ-500. To enrich mRNA, the total RNA was fractionated using oligo-dT attached magnetic beads. Then, short mRNA fragments (approximately 200 bp) were generated by adding the fragmentation buffer. The first strand cDNA was synthesis primed by random N6 primers, and the second strand was consequently synthesized by adding the dNTPs, buffer, DNA polymerase I and RNase H. The short double stranded cDNA fragments were then prepared for ligation to the adapters and PCR amplification by adding an extra single nucleotide adenine to the 3 end. Next, double stranded PCR products were denatured by heating and circularized by the splint oligo sequence. The ssCir DNA (single strand circle DNA) was formatted as the final sequencing library, and the library was amplified with phi29 to make DNA nanoball (DNB). Then, the libraries were sequenced using the BGISEQ-500 sequencer. All RNA-seq raw data from this study have been submitted to the NCBI Sequence Read Archive (SRA) under accession No. PRJNA561881.
The sequencing reads were mapped to the reference genome of Coffea canephora (http://coffeegenome.org/) using HISAT2 [53]. Before mapping, the quality of RNA-seq reads was analyzed by SOAPnuke [18]; the adapter sequence and low-quality reads were removed by Trimmomatic [54]. For gene expression, the clean reads were mapped to genes by Bowtie2 [53], and the FPKM (fragments per kilo bases of exons for per million mapped reads) value for transcripts were calculated by RSEM [55].

Data Collection and Identification of NAC Genes
The Coffea canephora genome sequences were downloaded from Coffee Genome Hub (http://coffeegenome.org/) [56]. For comparative analysis, all the predicted protein sequences of Arabidopsis thaliana, Carica papaya, Populus trichocarpa, Vitis vinifera, Oryza sativa, Amborella trichopoda, and Physcomitrella patens were retrieved from Phytozome 13 (https://phytozome-next.jgi.doe.gov/). To identify the NAC genes, all the retrieved proteins were used as queries to search against the profile Hidden Markov Model (HMM) by HMMER (version 3.1b2, http://hmmer.org/) [57] with the Pfam HMM library Pfam 32.0 [58], and the proteins with NAC domains (PF02365) with E-values below 0.1 were isolated. For Arabidopsis and rice, the isolated NAC genes were compared with previously published sequences to obtain the nomenclature [4,6].

Phylogenetic Analysis
The multiple sequence alignment of the full-length proteins sequences of NAC genes was carried out using the MUSCLE program with default parameters [59]. The unrooted polygenetic tree was constructed by using the NJ (Neighbor-joining) approach via MEGA 6 with parameters set as follows: Jones-Taylor-Thornton (JTT) model, bootstrap values of 1000 replicates and pairwise deletion [60].

Chromosome Location, Gene Structure, and Conserved Motifs in the CocNAC Family
The information concerning the detailed chromosome locations of CocNAC genes was obtained from the Coffee Genome Hub with the Tool 'Advanced Search'. These data were then integrated and modified by using Inkscape software (version 0.92, Inkscape's contributors, https://inkscape.org/). Gene structures of CocNACs were drawn by using the GSDS (Gene Structure Display Server, version 2.0, http://gsds.cbi.pku.edu.cn/) [61]. The conserved motifs of CocNAC proteins were analyzed by using the MEME (Multiple Em for Motif Elicitation) tool (version 5.0.4, http://meme-suite.org/) [62] with the following parameters: maximum motif width: 50; minimum motif width: 6; maximum number of motifs, 10; distribution of motif occurrences, zero or one per sequence.

qRT-PCR
RNA was isolated using the RNAprep Pure Plant Plus Kit (TIANGEN BIOTECH CO., LTD, Beijing, China), according to the manufacturer's protocols. The 1st-strand cDNA was synthesized using the PrimeScript™ RT reagent Kit with gDNA Eraser (TAKARA BIO INC., Shiga, Japan). The cDNA was diluted to 12.5 ng/µL for PCR. A QuantStudio TM 7 Flex Real Time PCR System (Applied Biosystems ® , Foster City, CA, USA) was used for PCR analysis. The reactions were prepared in a total volume of 20 µL containing 2 µL cDNA, 1.0 µL of each primer at 10 µM, 10 µL of SYBR ® Premix Ex Taq™ II (Tli RNaseH Plus) (TAKARA BIO INC., Shiga, Japan), and 6 µL distilled water. PCR was performed with the following program: 95 • C for 30 s, followed by 40 cycles of 95 • C for 5 s, and 56 • C for 30 s. All of the primer sequences that were used in this study are listed in Table S1. The 2 −∆∆C T method was used for PCR data analysis [63], and the data were re-represented by a heatmap, as described previously [64,65] using MeV software [66] using the average mean values of three independent biological replicates. All qRT-PCR analyses was carried out using a sample from the three biological repeats with triplicate technique repeats for each sample.

Identification of Putative NAC Genes in Coffea canephora
To identify the members of the CocNAC gene family, the entire Coffea canephora (V1.0) genome sequence and 25,574 protein sequences were downloaded from the Coffee Genome Hub (http: //coffee-genome.org/) [56]. After HMM (Hidden Markov Models) search, 63 CocNAC gene models were isolated. These gene models were named from CocNAC1 to CocNAC63 according to their position on the C. Canephora chromosomes 1-11 and from top to bottom ( Table 1). The protein sizes of CocNACs varied from 87 amino acids (a.a.) (CocNAC10, Cc01_g19370) to 664 a.a. (CocNAC39, Cc07_g12760), with an average length of 350 a.a. The relative molecular weight (Mw) varied from 9.91 kDa (CocNAC10, Cc01_g19370) to 75.49 kDa (CocNAC39, Cc07_g12760). The locus ID, chromosome location, protein length, introns, molecular weight (Mw), and PI (isoelectric point) of the resulting genes are shown in Table 1. CocNACs were widely distributed on the C. Canephora chromosomes with a non-homogeneous distribution: except for eight CocNACs (CocNAC56, CocNAC57, CocNAC58, CocNAC59, CocNAC60, CocNAC61, CocNAC62, and CocNAC63), all other CocNACs could be mapped on chromosomes Ch1-Ch11 (Table 1 and Figure 1). No CocNAC was found on chromosome Chr3, and more than 61.9% of all CocNACs were located on five chromosomes: chromosome Chr1 (ten CocNACs), chromosome Chr2 (ten CocNACs), chromosome Chr5 (seven CocNACs), chromosome Chr7 (six CocNACs), and chromosome Chr10 (six CocNACs) (Figure 1). For comparative genomic analysis, the NAC genes of Arabidopsis thaliana, Carica papaya, Populus trichocarpa, Vitis vinifera, Oryza sativa, Amborella trichopoda, and Physcomitrella patens were obtained from previous publications [4,6,10,11] or identified in this study ( Figure 2). Finally, 115, 80, 169, 74, 145, 46, and 33 NACs were found from Arabidopsis thaliana, Carica papaya, Populus trichocarpa, Vitis vinifera, Oryza sativa, Amborella trichopoda, and Physcomitrella patens, respectively. The evolutionary relationship of each species and the number of NAC genes are shown in Figure 2. The results showed that the number of NAC genes in Coffea canephora is larger than in Amborella trichopoda and Physcomitrella patens and less than in Arabidopsis thaliana, Populus trichocarpa, and Oryza sativa. This observation implies that the number of NAC genes was not dependent the genome size, and that various forms of regulation might be present in NAC function.

Phylogenetic Analysis of the CocNACs
To obtain insight into the evolutionary relationships of the CocNAC genes, a phylogenetic analysis of CocNACs, ANACs (Arabidopsis thaliana NAC proteins) and OsNACs (Oryza sativa L. NAC proteins) was carried out using MEGA6 [60]. The ANACs and OsNACs were also isolated from the Arabidopsis and rice genomes (Araport11) as described above for C. canephora. After HMM search, 120 ANAC and 145 OsNAC proteins were obtained (Table S2). For rice, the "LOC_" prefix was removed from the RGAP (Rice Genome Annotation Project) locus IDs as the description of previous research [6]. After comparing with the phylogenetic analysis of rice and Arabidopsis, the CocNACs, ANACs and OsNACs were divided into 16 subgroups, NAC-a to NAC-p ( Figure 3). The CocNAC genes were found to be unequally distributed in all the subgroups except the NAC-l subgroup, in which all the NAC members were from Arabidopsis. No Arabidopsis NAC gene was present in subgroup NAC-n. Multiple sequence alignments of NAC proteins were carried out using the MUSCLE program with default parameters. The un-rooted phylogenetic tree was constructed by MEGA 6.

Gene Structure and Protein Motif Analysis of C. canephorae NAC Gene Family
In order to identify the structural diversity of the CocNACs, the exon-intron content in the coding sequences of the individual CocNACs was analyzed ( Figure 4B). Additionally, we built a separate unrooted phylogenetic tree using the protein sequences of all the CocNACs, which divided the CocNACs into 13 subfamilies, subgroup I to XIII ( Figure 4A). All CocNAC genes contained intron(s), implying possible alternative splicing. With few exceptions, CocNAC genes in the same group contained similar numbers of exons and introns. CocNAC genes in subgroups I, II, IV, V, X, and XI contained two or three introns. More variable gene structures were found in subgroups VI, IX, XII and XIII. Several genes, such as CocNAC25, CocNAC44, CocNAC22 and CocNAC30, included relatively large introns.
The MEME program was used to investigate the structure of CocNAC proteins. In total, 10 motifs, motifs 1-10, were detected from the CocNAC proteins. As expected, similar motif compositions were found among the closely related CocNAC members ( Figure 4C; Figure S1). Motif 3, motif 4, motif 2/5, motif 1/7, and motif 6 comprise the subdomains A, B, C, D, and E of NAC, respectively ( Figure S1). All of the CocNACs in subgroups III, IV and IX lacked motif 4, but IX contained motif 8. Most of subgroup XIII contained motif 9, and two out of three CocNACs in subgroup IV contained motif 10. This pattern suggests that the specific motifs may be related to specific functions of different subgroups.

qPCR Analysis of CocNAC Genes under Low Temperature Treatment
To investigate the response of CocNACs to low temperature treatments, one-year-old C. canephora seedlings were subject to 13 • C (acclimation) and 4 • C (cold temperature), and the top two pairs of recent mature leaves were sampled for qRT-PCR. After removing the genes with Ct values above 35 cycles, 38 CocNAC genes were found to be differentially responsive to low temperature in C. canephora ( Figure 5). Four genes including CocNAC18 were up-regulated by 13 • C treatment and sustained the up-regulation at 4 • C ( Figure 5A). Six genes including CocNAC36 were found to be down-regulated by 13 • C treatment when it was considered a cold acclimation process ( Figure 5B). Seven genes including CocNAC47 showed preferential up-regulation by 4 • C treatment ( Figure 5C). Four genes were downregulated by 13 • C/8 • C, but slightly upregulated by 4 • C/4 • C cold treatments ( Figure 5D). Seventeen genes including CocNAC23 were continuously down-regulated during the treatment ( Figure 5E). Grouping by the expression pattern was not consistent with classification by sequence motifs (Figure 4C), indicating that genes showing similar evolutionary relationships might be differentially responding to temperatures.

CocNAC Genes Involved in Coffee Bean Development
Since several NAC genes are involved in seed size and quality as mentioned in the Introduction, identification of CocNACs expressed specifically during coffee bean development is a promising topic. To explore the patterns of CocNACs expression and highlight the roles of CocNACs in developing coffee beans, four differential developmental stages of coffee beans were collected, and RNA_seq data were analyzed regarding the CocNAC expression patterns. Based on the FPKM (fragments per kilobases of exons for per million mapped reads) values, 54 out of 63 CocNAC genes showed a variety of expression patterns during different developmental stages of coffee beans, in which transcription patterns were distributed into five clades ( Figure 6A). Clade I contained four genes, including CocNAC04, that showed preferential expression in Y (coffee bean from yellow fruit) and slightly upregulated in R (coffee bean from red fruit). Clade II Included nine genes highly expressed in R. Clade III Included 15 genes, displaying the highest expression in Y and downregulated in R. Clade IV included 16 genes whose relative transcription levels were high at the LG stage (coffee bean from large green fruit), and six of these also were highly expressed at the SG stage (coffee bean from small green fruit). Clade V was composed of ten genes displaying the highest expression in the SG stage, and five were highly accumulated until the LG stage. To confirm the expression pattern of detected CocNACs in developing coffee beans, ten genes were randomly selected for qRT-PCR analysis, and each showed similar expression patterns with the RNA sequencing results ( Figure 6B). Based on the expression patterns, these genes showed a coffee bean development-associated expression pattern with either development onset or the dynamic process. Therefore, some genes might be highlighted as positive or negative regulators associated with coffee bean development: three genes from clade II (CocNAC30, CocNAC44, and CocNAC31) were continuously upregulated and were considered as the positive regulators, while six genes from clade IV (CocNAC45, CocNAC38, CocNAC25, CocNAC29, CocNAC33, and CocNAC23) and five genes (CocNAC62, CocNAC03, CocNAC26, CocNAC39, and CocNAC12) from clade V were mostly downregulated and negatively correlated with coffee bean development.

CocNAC Gene Identification and Evolutionary Analysis in C. canephora
Plant-specific NAC transcription factors (TFs) are one of the largest families of transcriptional regulators in plants, and they play important roles in the regulation of growth and development as well as in plant stress responses [67][68][69]. Due to the importance of NAC genes in plants, genome-wide identification and classification of NAC genes have been carried out in many plant species, but little is known about this gene family in the tropical tree C. canephora. In this study, we identified 63 CocNAC genes based on the C. canephora genome sequence [56] and named these genes as CocNAC01 to CocNAC63 based on their chromosome locations (Table 1 and Figure 1). To obtain insights into the evolution of CocNACs, comparative genomic analysis was carried out among C. canephora, A. thaliana, C. papaya, P. trichocarpa, V. vinifera, O. sativa, A. trichopoda, and P. patens (Figure 2). Clearly, the size of the NAC gene family in C. canephora is similar to those of C. papaya and V. vinifera, larger than in A. trichopoda, and P. patens, but smaller than in A. thaliana, P. trichocarpa, and O. sativa. The size difference might due to whole-genome polyploidization events during the evolution and differentiation of each species; for example, the C. canephora genome experienced three genome duplication events, and the Arabidopsis genome went through five such events [56,70]. Furthermore, the number of NAC genes did not match with the genome size of each species, indicating the stable features of NAC genes in different species during evolution.
Based the analysis of basic characteristics of CocNAC genes, the gene sizes and isoelectric points (PI) were very different, but relatively conservative gene structure was observed (Table 1 and Figure 4B). Of the 63 CocNACs, 34 had two introns, while the number of introns varied from one to six. Similar results were observed in F. tataricum [9], pepper (Capsicum annuum L.) [71], tobacco [72], and P. trichocarpa [10]. The presence of introns in CocNACs implied the presence of alternative splicing (AS) depending upon developmental and environmental changes. AS is a novel means of regulating the environmental fitness of plants [73], especially for controlling abiotic stress responses such as cold stress [74]. More introns produce more alternatively spliced mRNA isoforms, resulting in production of more diverse proteins [75].
The unrooted tree was constructed using NAC protein sequences from C. canephora, Arabidopsis and rice to explain the phylogenetic relationships existing among the three species, and all the NAC protein sequences were grouped into 16 subgroups (Figure 3). Most of the subgroups contained NAC members similar to those of Arabidopsis and rice previously published [4,6]. Consistent with previous studies, the bootstrap values were somewhat low due to the large number of NAC proteins [10,76]. Remarkably, the NAC-l subgroup did not include any CocNAC members, while no Arabidopsis NAC members were found in subgroup NAC-n, indicating that the NACs were differentially acquired and retained in each species during evolution (Figure 3).
The MEME program was used to examine the structure of CocNAC proteins. Ten conserved motifs were found, and the conserved motifs were present in the N-terminal of the NAC proteins ( Figure 4C). The five subdomains (A to E) of NAC domains were represented by motif 3, 4, 2/5, 1/7, and 6, respectively ( Figure S1). In subgroup IX, the motifs 2 (subdomain C), 3 (subdomain A), and 4 (subdomain B) were replaced by motif 8 (Figure 4C), indicating that NACs in this subgroup might have different functions. This pattern was similar to the results for F. tataricum [9].

Identification of CocNACs Responding to Cold Stress in C. canephora
Low temperature can cause irreversible damage and can impact plant growth and development. Many plants can improve their cold tolerance by cold acclimation [77]. Under cold stress, the COR (cold-regulated) genes, governed by multiple TFs, can be induced to change the features of plant physiology and biochemistry [77,78]. Among these TFs, CBFs (C-repeat binding factors) are well known [79]. However, it has been reported that only 10%-25% of COR genes are mediated by CBFs, indicating that other TFs may be involved in cold tolerance [80].
Low temperature negatively impacts coffee plant growth and yield [49] and also influences biochemical composition of green Arabica coffee beans [34]. In particular, cold stress affects photosynthesis, antioxidative systems, membrane lipid composition, and expression of genes involved in these aspects of plant physiology [44,81]. This situation led to the identification of TFs responding or controlling cold stress in C. canephora. As mentioned in the Introduction, NAC TFs play both positive and negative roles in cold stress. Apple NAC 29 (MdNAC029) is a negative regulator that acts by repressing CBF genes under cold stress [82], while banana MaNAC1 is a positive regulator of cold stress that acts by interacting with MaCBF1 [28].
A large number of NAC TFs confer cold tolerance in various plant species. In our study, 38 of 63 CocNACs differentially responded to cold stress in terms of transcript levels ( Figure 5). Six genes were downregulated by 13 • C treatment, indicating that they might be involved in the cold acclimation process ( Figure 5B). Seven genes were preferentially up-regulated by the 4 • C treatment ( Figure 5C), suggesting that they might play roles in cold acclimation. Seventeen genes were continuously downregulated by the cold treatment ( Figure 5E), and they may act as negative regulators during cold treatments. Notably, the functions of some CocNACs (CocNAC18, CocNAC23 and CocNAC47) can be predicted based their orthologs in the model plant Arabidopsis [83][84][85]. Overexpression of ANAC035, the ortholog of CocNAC18, confers cold tolerance by induction of COR genes [83], and the CocNAC47 ortholog ANAC072 functions as a transcription activator under abiotic stress [84]. CocNAC23 seems to be function in cold stress as well as in coffee bean development (Figures 4 and 5). Its Arabidopsis ortholog, ANAC083/VNI2, plays an integrator role in the signals for leaf senesce and COR genes [85]. Taken together, our analysis extends the knowledge regarding NACs functioning for cold responses in plants. After elucidation of these genes' functions, some could be employed for production of cold-tolerant coffee plants.

CocNAC Genes May Play an Important Role in Coffee Bean Development
Coffee beans, the seeds of the coffee plant, contain all the precursors essential for the generation of the coffee aroma during roasting [34,51]. Establishing the relationships between the beverage quality and chemical composition of green coffee beans is the main target of research into coffee quality. Characterization of gene profiles during coffee bean development and identification of key quality-related genes can assist in the coffee breeding industry [51]. In Arabidopsis, seed development is a complex process, and many of TFs are involved [86]. Regarding NAC TFs and seed development, promising results have been obtained from the study of the Poaceae family of monocotyledonous plants [40]. Notably, the six barley NACs in subclade NAC-d-9 are specifically expressed in the developing grain [40]. Wheat homologs of these genes show strong and exclusive expression patterns in wheat grain tissues [41]. In addition, two rice NAC TFs, ONAC020 and ONAC026 that fall into the same clade, are highly expressed specifically in grains during maturation, determining grain size and weight [42]. Furthermore, two maize NACs that also fall into the same clade I are specifically expressed in grains [43]. However, this is not the case for other plant systems such as Arabidopsis, implying the possible presence of other CocNACs for bean development.
Similar to previous reports, several CocNACs appeared to be involved in coffee bean development. Fifty-four out of 63 CocNAC genes were highlighted by RNA_seq data as being involved in coffee bean development, and these could be divided into five clades ( Figure 6A), indicating the variety of functions of CocNACs during coffee bean development. Ten of those genes were confirmed by qRT-PCR ( Figure 6B). CocNAC29 showed continuous down-regulation during coffee bean development and was considered as one of the negative regulators. CocNAC29 belongs to subgroup NAC-f, and its Arabidopsis homolog is ANAC056 (also called NARS1) (Figures 3 and 6A). ANAC056 (NARS1) and ANAC018 (NARS2) redundantly regulated embryogenesis in Arabidopsis [36] and both were grouped into NAC-f. In subgroup NAC-f, six CocNACs (CocNAC16, 07, 08, 32, 29, and 53) were found, and five (CocNAC07, 08, 32, 29, and 53) of them were expressed in coffee beans (Figures 3 and 6A), indicating that all the NAC-f members have conserved functions during seed development. Os01g01470 (ONAC020) and Os01g29840 (ONAC026) were reportedly related to seed size and weight in rice [42], and both were grouped into subgroup NAC-k (Figures 3 and 6A). Furthermore, we found that CocNAC57, closest to Os01g01470 and Os01g29840, was highly expressed in the LG and Y stages. As is well known, the coffee bean size is determined at the stages of LG and Y, indicating that CocNAC57 might also have a similar function. Notably, the heat map III in Figure 6 may include important CocNACs (CocNAC10, CocNAC32, CocNAC57 and CocNAC60) for coffee bean development, but none of them except CocNAC57 have been characterized. According to the expression patterns of CocNAC genes, it can be hypothesized that these genes may have conserved functions during coffee bean development; the specific functions of CocNACs need to be further verified.

Conclusions
In this study, we performed a comprehensive analysis of gene structure, chromosomal location, and phylogenetic relationships of the NAC gene family in C. canephora and examined the expression patterns of CocNAC genes in developing coffee beans by RNA_seq followed by qRT-PCR, and under cold stress using qRT-PCR. Analysis of 63 CocNAC genes provided evolutionary relationships, putative function with distinct motifs and diversity of regulation of gene expression. Several genes appeared to be positive or negative regulators under cold stress. At least four CocNACs may be positive regulators of coffee bean development. Our systematic analysis of CocNAC genes lays the foundation for functional characterization of CocNAC genes and development of cold tolerance in coffee plants.
Supplementary Materials: The following are available online at http://www.mdpi.com/2073-4395/9/11/670/s1, Figure S1. Sequence logos for conserved motifs identified in CocNAC proteins by MEME tools, Table S1. List of primers used in this study, Table S2. NAC genes from Arabidopsis thaliana and Oryza sativa L.