Genome-Wide Identification and Expression Profile Analysis of the NF-Y Transcription Factor Gene Family in Petunia hybrida

Nuclear Factor Ys (NF-Ys) are a class of heterotrimeric transcription factors that play key roles in many biological processes, such as abiotic stress responses, flowering time, and root development. The petunia (Petunia hybrida) is a model ornamental plant, and its draft genome has been published. However, no details regarding the NF-Y gene family in petunias are available. Here, 27 NF-Y members from the petunia genome were identified, including 10 PhNF-YAs, 13 PhNF-YBs, and 4 PhNF-YCs. Multiple alignments showed that all PhNF-Y proteins had clear conserved core regions flanked by non-conserved sequences. Phylogenetic analyses identified five pairs of orthologues NF-YB proteins from Petunia and Arabidopsis, and six pairs of paralogues NF-Y proteins in Petunia. Analysis of the gene structure and conserved motifs further confirmed the closer relationship in each subfamily. Bioinformatics analysis revealed that 16 PhNF-Ys could be targeted by 18 miRNA families. RNA-seq results showed that expression patterns of PhNF-Ys among four major organs (leaf, stem, flower, and root) were clustered into six major groups. The stress response pattern of PhNF-Ys was identified under cold, heat, drought, and salinity treatments. Based on the RNA-seq data, we found that 3 genes responded to drought, 4 genes responded to salt, 10 genes responded to cold, and 9 genes responded to hot. In conclusion, this study provides useful information for further studying the functions of NF-Ys in stress response.


Introduction
Abiotic stress, such as temperature, salinity, and drought, can cause great damage to plant growth and development. To survive, plants have evolved strategies to tolerate stress, including developmental, physiological, morphological, ecological, biochemical, and molecular mechanisms [1]. Members of many transcription factor (TF) families, such as AP2/ERF, MYB/MYC, zinc-finger protein, and NAC, play important roles at molecular levels by increasing or decreasing the expression of stress responsive genes [2]. Compared to those TF families, Nuclear Factor Y (NF-Y) is a relatively late reported TF family that participates in abiotic stress response.
NF-Y is a heterotrimeric complex with three distinct subunits, NF-YA, NF-YB, and NF-YC, that binds to the CCAAT element in the promoters of genes [3]. All three subunits contain the conserved DNA binding domains and subunit interaction domains [4]. To form heterotrimeric complexes, NF-YB and NF-YC, who possess H2B and H2A motifs, respectively, form a tight histone dimer in the cytoplasm, PhNF-YA1-PhNF-YA10, PhNF-YB1-PhNF-YB13, and PhNF-YC1-PhNF-YC4. The characteristics of the PhNF-Y sequences are listed in Table 1. The lengths of gene sequences ranged from 327 bp to 1065 bp. The lengths, the molecular weights (MWs), and the isoelectric point (pI) values of these proteins ranged from 108 aa to 327 aa, from 11.92 kDa to 39.76 kDa, and from 4.25 to 10.39, respectively.

Multiple Alignment and Conserved Motifs of NF-Y Proteins
Multiple alignments showed that all PhNF-Y proteins had clear conserved core regions flanked by non-conserved sequences. The core regions of most petunia NF-YA proteins contained two conserved domains, one for the NF-YB/C interaction (20AA) and the other for binding CCAAT domain in the promoter (21AA) (Figure 1a). The two domains were separated by a conserved linker (11AA). However, the PhNF-YA6 and PhNF-YA4 lacked the conserved DNA binding domain of the NF-YA subfamily, which was necessary for the recognition of the CCAAT domain. In addition, a conserved 12 AA sequence with an unknown function was identified in front of the NF-YB/C interaction domain. We also noticed that most amino acid residues required in mammals and yeast [3,41] were remained in most petunia NF-YA proteins. For these required amino acid residues, the arginine (R24) was often replaced by alanine or glycine, and the alanine (A31) was mutative in plant NF-YA proteins [31]. In addition, eight of the petunia NF-YA proteins contained three predicted nuclear localization signals [5,42] (black boxes in Figure 1a). and the reported functions of AtNF-Ys, the functions of PhNF-Y members could be concluded and further tested. Notably, five pairs of NF-YB orthologues (AtNF-YB7 and PhNF-YB7; AtNF-YB4 and PhNF-YB4; AtNF-YB10 and PhNF-YB10; AtNF-YB6 and PhNF-YB6; AtNF-YB3 and PhNF-YB3) were found in the NF-YB group, suggesting a significant possibility of the similar biological functions between the NF-YB orthologues. Moreover, six pairs of paralogues were identified: PhNF-YB1 and PhNF-YB8; PhNF-YB11 and PhNF-YB12; PhNF-YA1 and PhNF-YA9; PhNF-YA5 and PhNF-YA6; PhNF-YA2 and PhNF-YA10; PhNF-YA7 and PhNF-YA8.  The conserved region of Petunia NF-YB proteins had structural similarities with H2B histone. This region possessed characteristic domains for DNA-binding and NF-YA/C interaction (Figure 1b). The required amino acid residues were highly conserved in PhNF-YB proteins, particularly the functionally important protein residues. The characteristic intra-chain arginine-aspartate (R58-D65) bidentate pairs, which were necessary in stabilizing the NF-YB/C heterodimer, were absolutely conserved [43]. Notably, similarly with Arabidopsis and soybean, the disulfide bond between C33 and C37 was not presented, which played a crucial role in the translocation of NF-YB/C into nuclei in humans and Aspergillus nidulans [43]. This evidence further confirmed the conserved histone structure and changed the model of nucleus translocation of NF-YB proteins in plants.
The PhNF-YC proteins contained highly conserved NF-YA/B interaction domains and had structural similarities with H2A histone [43]. The required amino acids were conserved with few alterations by similar functional residues [44] (Figure 1c). Similarly with NF-YB subunits, the arginine (R54) and aspartate (D61) pairs were present in all petunia NF-YCs [43]. Another specific feature-tryptophan (W46) at the end of α2, which made the specific interaction of NF-YC and NF-YB, was also absolutely conserved in all PhNF-YCs.

Gene Structures and Conserved Motifs of PhNF-Y Members
To further understand the evolutionary relationships in the petunia NF-Y family, the exonintron structures were analyzed using the GSDS website ( Figure  3; http://gsds.cbi.pku.edu.cn/index.php). Most of the PhNF-YAs had five or six exons with similar distribution, an exception being PhNF-YA6. Half of the PhNF-YBs had no intron, and the other half were variable, but the members involved in the same clade had similar numbers of exons and introns. For the PhNF-YC subfamily, PhNF-YC1 and PhNF-YC3 had only one exon, and PhNF-YC2 and PhNF-YC4 had two exons. The gene structures of the NF-Y members were positively associated with their phylogenetic relationship and were consistent with the previous reports [8,43].
The putative conserved motifs of 27 PhNF-Ys were identified using Multiple Em for Motif Elicitation (Figure 4; MEME 5.1.0). Ten conserved motifs were recognized in each subfamily, and the lengths of the motifs ranged from 6 to 50 amino acids. For the PhNF-YA subfamily members, Motif 1, 2, and 3 were presented in 8 of the 10 NF-YAs. PhNF-YA4 lacked Motif 1 and 3, and NF-YA8 lacked Motif 3. Furthermore, the members of the same clade contained similar motif numbers and distribution patterns, indicating that they might have similar functions. All PhNF-YB proteins contained Motif 11, 12, and 13, and the variation of motifs were consistent with the phylogenetic relationships. In the case of the PhNF-YC subfamily, Motif 21, 22, and 23 were observed in all four members, and Motif 25, 26, 27, and 28 were found in PhNF-YC2 and PhNF-YC4.

Gene Structures and Conserved Motifs of PhNF-Y Members
To further understand the evolutionary relationships in the petunia NF-Y family, the exon-intron structures were analyzed using the GSDS website ( Figure 3; http://gsds.cbi.pku.edu.cn/index.php). Most of the PhNF-YAs had five or six exons with similar distribution, an exception being PhNF-YA6. Half of the PhNF-YBs had no intron, and the other half were variable, but the members involved in the same clade had similar numbers of exons and introns. For the PhNF-YC subfamily, PhNF-YC1 and PhNF-YC3 had only one exon, and PhNF-YC2 and PhNF-YC4 had two exons. The gene structures of the NF-Y members were positively associated with their phylogenetic relationship and were consistent with the previous reports [8,43].
The putative conserved motifs of 27 PhNF-Ys were identified using Multiple Em for Motif Elicitation (Figure 4; MEME 5.1.0). Ten conserved motifs were recognized in each subfamily, and the lengths of the motifs ranged from 6 to 50 amino acids. For the PhNF-YA subfamily members, Motif 1, 2, and 3 were presented in 8 of the 10 NF-YAs. PhNF-YA4 lacked Motif 1 and 3, and NF-YA8 lacked Motif 3. Furthermore, the members of the same clade contained similar motif numbers and distribution patterns, indicating that they might have similar functions. All PhNF-YB proteins contained Motif 11, 12, and 13, and the variation of motifs were consistent with the phylogenetic relationships. In the case of the PhNF-YC subfamily, Motif 21, 22, and 23 were observed in all four members, and Motif 25, 26, 27, and 28 were found in PhNF-YC2 and PhNF-YC4.

miRNA Target Site Prediction and Validation
It is well known that microRNAs can regulate target genes though cleavage of mRNA or translation inhibition [45]. In the petunias, 58 mature miRNA sequences that belonged to 39 miRNA families have been reported. Twenty-five were conserved at a wider taxonomic range, 19 were specific for the Solanaceae, and 13 variants were unique to the petunia [46].
All these sequences were used to predict target transcript candidates of PhNF-Ys. The results showed that 16 PhNF-Ys could be regulated by 18 miRNA families (Table S1). Eight PhNF-YAs, 6 PhNF-YBs, and 2 PhNF-YCs had target sites of 14, 5, and 3 miRNA families, respectively. In most cases, one miRNA family had one or two target genes in the NF-Y family. However, eight NF-YA members had target sites of miRNA169 ( Figure 5), which was consistent with previous reports [47,48]. These results indicated that the miR169-NF-YA module existed in Petunia and that PhNF-YAs had more uncovered relationships with other miRNA families.

miRNA Target Site Prediction and Validation
It is well known that microRNAs can regulate target genes though cleavage of mRNA or translation inhibition [45]. In the petunias, 58 mature miRNA sequences that belonged to 39 miRNA families have been reported. Twenty-five were conserved at a wider taxonomic range, 19 were specific for the Solanaceae, and 13 variants were unique to the petunia [46].
All these sequences were used to predict target transcript candidates of PhNF-Ys. The results showed that 16 PhNF-Ys could be regulated by 18 miRNA families (Table S1). Eight PhNF-YAs, 6 PhNF-YBs, and 2 PhNF-YCs had target sites of 14, 5, and 3 miRNA families, respectively. In most cases, one miRNA family had one or two target genes in the NF-Y family. However, eight NF-YA members had target sites of miRNA169 ( Figure 5), which was consistent with previous reports [47,48]. These results indicated that the miR169-NF-YA module existed in Petunia and that PhNF-YAs had more uncovered relationships with other miRNA families.

Expression Profiles of the PhNF-Ys in Different Tissues
To investigate the biological functions of PhNF-Y genes in the petunias, the expression patterns of PhNF-Y genes were analyzed in four major tissues (leaf, stem, root, and flower) of plants under normal conditions. Based on the RNA-seq data, the corresponding expression data of 23 PhNF-Y genes were collected, and the other four NF-Y genes were not detected and were thus omitted from the analysis. As shown in Figure 6

Expression Profiles of the PhNF-Ys in Different Tissues
To investigate the biological functions of PhNF-Y genes in the petunias, the expression patterns of PhNF-Y genes were analyzed in four major tissues (leaf, stem, root, and flower) of plants under normal conditions. Based on the RNA-seq data, the corresponding expression data of 23 PhNF-Y genes were collected, and the other four NF-Y genes were not detected and were thus omitted from the analysis. As shown in Figure 6

Expression Profiles of the PhNF-Ys under Abiotic Stress
Previous reports have verified that NF-Y members could respond to various abiotic stresses [49]. We performed RNA-seq to identify the expression levels of PhNF-Ys under several stresses ( Figure  7). Based on the RNA-seq data, 9 NF-Y genes not detected were omitted from the analysis. The results showed that PhNF-YA6 and -10 and PhNF-YB13 were regulated by drought stress (Figure 7a). PhNF-YA6 and PhNF-YA10 reached a peak at 6 h and subsequently decreased to a normal level after 12 h of treatment, while PhNF-YB13 decreased sharply after 12 h of treatment. With the salt treatment (Figure 7b), the expression levels of four PhNF-Ys (PhNF-YA5, -6, and-−10 and PhNF-YB3) were regulated by salt. PhNF-YA5, -6, and -10 were induced during the early stages, reached the peak at 12 h, and then decreased. The expression level of PhNF-YB3 reached the top at 1 h and then quickly decreased to a normal level. Under cold stress (Figure 7c), 10 PhNF-Ys were upregulated or downregulated at different timepoints. For NF-YA subfamily members, PhNF-YA1 and PhNF-YA7 were obviously downregulated and PhNF-YA4, -5, -6, and -10 were markedly upregulated at one or more timepoints, suggesting these genes might have obtained divergent functions under cold stress during evolution. In addition, three NF-YB (PhNF-YB1, -3, and 2) members were significantly inhibited after 12 h of treatment. PhNF-YB13 showed a gradual increase from 3 to 12 h. With hot stress (Figure 7d), four genes were induced, and five genes decreased. PhNF-YA1, -2, and -10 were increased within 1 h and PhNF-YC1 exhibited a gradual increase. PhNF-YA4, -5, and -6 and PhNF-YB13 decreased immediately in 1 h. The expression of PhNF-YA7 was the lowest at 6 h. These results suggested that the PhNF-Ys had a faster response to hot stress. Notably, we found that more NF-Y

Expression Profiles of the PhNF-Ys under Abiotic Stress
Previous reports have verified that NF-Y members could respond to various abiotic stresses [49]. We performed RNA-seq to identify the expression levels of PhNF-Ys under several stresses (Figure 7). Based on the RNA-seq data, 9 NF-Y genes not detected were omitted from the analysis. The results showed that PhNF-YA6 and -10 and PhNF-YB13 were regulated by drought stress (Figure 7a). PhNF-YA6 and PhNF-YA10 reached a peak at 6 h and subsequently decreased to a normal level after 12 h of treatment, while PhNF-YB13 decreased sharply after 12 h of treatment. With the salt treatment (Figure 7b), the expression levels of four PhNF-Ys (PhNF-YA5, -6, and -10 and PhNF-YB3) were regulated by salt. PhNF-YA5, -6, and -10 were induced during the early stages, reached the peak at 12 h, and then decreased. The expression level of PhNF-YB3 reached the top at 1 h and then quickly decreased to a normal level. Under cold stress (Figure 7c), 10 PhNF-Ys were upregulated or downregulated at different timepoints. For NF-YA subfamily members, PhNF-YA1 and PhNF-YA7 were obviously downregulated and PhNF-YA4, -5, -6, and -10 were markedly upregulated at one or more timepoints, suggesting these genes might have obtained divergent functions under cold stress during evolution. In addition, three NF-YB (PhNF-YB1, -3, and 2) members were significantly inhibited after 12 h of treatment. PhNF-YB13 showed a gradual increase from 3 to 12 h. With hot stress (Figure 7d), four genes were induced, and five genes decreased. PhNF-YA1, -2, and -10 were increased within 1 h and PhNF-YC1 exhibited a gradual increase. PhNF-YA4, -5, and -6 and PhNF-YB13 decreased immediately in 1 h. The expression of PhNF-YA7 was the lowest at 6 h. These results suggested that the PhNF-Ys had a faster response to hot stress. Notably, we found that more NF-Y genes participated in temperature stress response, suggesting that these genes might be key candidates for improving tolerance to temperature stress.
Plants 2020, 9, x FOR PEER REVIEW 9 of 16 genes participated in temperature stress response, suggesting that these genes might be key candidates for improving tolerance to temperature stress.

Discussion
Based on the petunia draft genome and the sequences of Arabidopsis NF-Y family members, we initially identified 10 NF-YAs, 13 NF-YBs, 4 NF-YCs, 3 NC2βs, and 1 NC2α, Since NC2α/β could not pair with NF-YB/C [50], these four members were removed from further analyses. For the PhNF-YC subfamily, only 4 members were found, less than 10 members of the Arabidopsis NF-YC gene family. The reduced numbers of PhNF-YC genes could result in a reduction of heterotrimers by two-thirds in theory. Smaller amounts of NF-YC members could also decrease the high redundancy of the NF-Y family in plants.
Multiple alignments showed that most PhNF-Y proteins contained the conserved regions responsible for subunit interaction and DNA binding (Figure 1), which were also found in other plants, mammals, and yeasts [32,38]. Previous studies have shown that the DNA binding domain in the C-terminus of NF-YA could bind to the CCAAT sites [41]. In this study, PhNF-YA4 and PhNF-YA6, which lacked the conserved DNA binding domain, might fail to recognize the CCAAT sites. The functions of these two members were worth further investigation. Additionally, although most required AA were conserved, PhNF-YC4 still exhibited some alterations in the NF-YA interaction domain. In consideration of the absence of homologues of the PhNF-YC4 gene in Arabidopsis, the ability to interact with the NF-YA subunit and the biological function of PhNF-YC4 required further testing.
To explore the functions of these PhNF-Ys, exon-intron structures and evolutionary analysis were performed (Figures 2 and 3). The number and distribution of introns and exons were similar to previous results [4,8,37], suggesting that the functions of PhNF-Ys might be similar to homologues genes in other species. In phylogenetic tree analysis, genes showing a close evolutionary relationship usually have a similar function, so we can predict the functions of PhNF-Ys based on the known functions of AtNF-Ys. For example, AtNF-YB3 was associated with the development of root tissue and flowering time [51], suggesting that its homologues gene, PhNF-YB3, might have similar functions. The higher expression level of PhNF-YB3 in leaf tissue further supported this possibility. AtNF-YA1 and AtNF-YA9 were reported to be critical regulators of embryo development [21,52], suggesting that the PhNF-YA genes clustered in the same clade might have similar functions. PhNF-YC2 and PhNF-YC4 were distant to NF-YCs in Arabidopsis, suggesting that these genes might have a unique role in determining specific traits in petunias. Moreover, the corresponding heterotrimeric complexes might be quite different from that in Arabidopsis.
Previous studies have shown that miRNA169 could direct the cleavage of the mRNA of NF-YA to regulate abiotic stress responses in different plant species [26,47,53]. In our study, we identified eight PhNF-YA genes possessing complementary sequences of different miRNA169 members, further supporting the existence of the miRNA169/NF-YA module in petunias. It was noted that some NF-YA genes contained the binding sites of other 13 miRNA family, including miR156, 159, 168,171, 172, 395, 482, 6164, 8016, 6149, 171(V), and 172(v). The functions of some miRNA were well studied in plants. For example, miR156 was a key component of the aging pathway through directly regulating its targets SQUAMOSA PROMOTER BINDING PROTEIN-LIKE (SPL) at a posttranscriptional level [45]. miRNA164 is involved in plant growth and abiotic stress responses through guiding the cleavage of the mRNAs of NAC genes [54]. The miR172 regulated flowering time and vegetative phase change through clearing the mRNA of AP2 and AP2-like genes [45]. In our study, PhNF-YA1, -2, -5, -7, and -10 showed an obvious stress response. PhNF-YA3 had a higher expression level in flower tissue. The functions of these PhNF-YA genes and the underlying relationships of PhNF-YA and miRNA were worth further investigation. In addition, PhNF-YB2, -8, -10, and -12 contained binding sites of miRNA156. In our previous study, CmNF-YB8 was shown to regulate flowering time via directly binding to the promoter of cmo-MIR156 [55]. In this study, the expression level of PhNF-YB8 was relatively constant. The function of PhNF-YB8 and the relationship of PhNF-YB8 and miR156 in petunias might differ from the situation in chrysanthemum. For the NF-YC subfamily, we first revealed that PhNF-YC2 contained target sites of miRNA164 and miRNA171(v), and PhNF-YC3 contained target sites of miRNA394. Both miR164 and miR171 were involved in stress response [45], and miR394 played roles in flower development [45]. The deeper relationships of NF-YC and miRNAs need to be explored.
Many NF-Y genes have been reported to participate in drought and salt stress responses. In previous studies, NF-Y genes regulated drought tolerance via different mechanisms, including the ABA-dependent pathway, the ABA-independent pathway, and even novel mechanisms. In our study, only PhNF-YA6 and -10 and PhNF-YB13 responded to drought. In previous studies, AtNF-YA6 involved ABA signaling [34], and GmNF-YA21, the homologues of AtNF-YA10, played roles in drought response [34]. Our results indicated that PhNF-YA6 and PhNF-YA10 were highly expressed in root tissue, which had a close relationship with stress tolerance, supporting the idea that PhNF-YA6 and PhNF-YA10 might be associated with drought response. PhNF-YB13 showed distant relationships with other NF-YB members, and had a high expression in leaves, roots, and stems, so its function requires further study. Salt stress showed a similar pattern of drought stress, and PhNF-YA5, -6, and -10 were quickly induced and were highly expressed in roots. In Arabidopsis, the induction of AtNF-YA5 under drought occurred at both the transcriptional and posttranscriptional levels [15]. The tissue-specificity of PhNF-YA5 expression indicated that it had higher expression in roots and stems. In our study, PhNF-YA5 was weakly induced after 3 h of drought stress. We speculated that PhNF-YA5 might regulate salt tolerances through regulating root or stem development. It has been well known that the mechanisms of drought response and salt response overlap to some degree. It is easier to understand that PhNF-YA5, -6, and -10 responded to both drought and salt. PhNF-YB3 was quickly induced by salt. AtNF-YB3, the homologues of PhNF-YB3, was identified to regulate flowering time [34]. In our study, PhNF-YB3 showed higher expression in leaf and flower tissue. We concluded that PhNF-YB3 should play roles in flowering time regulation and salt response, which requires further evidence.
There is little evidence that NF-Y genes are involved in cold and hot responses [30]. In our study, 10 PhNF-Y genes were obviously induced or suppressed by cold responses, and 9 PhNF-Y genes quickly responded to hot responses. Among these genes, seven (PhNF-YA1, -4, -5, -6, -7, and -10 and PhNF-YB13) were involved in both cold and hot responses. The expression of PhNF-YA4, -5, and -6 increased under cold responses and decreased under hot responses. The performances of PhNF-YA7 and -10 and PhNF-YB13 were similar under hot and cold responses. Moreover, except for PhNF-YA4, six other genes showed high expression in roots, indicating the underlying mechanisms that might be overlapped to a certain extent. Notably, we found that PhNF-YA6 and PhNF-YA10 responded to all stresses, including drought, salt, and cold and hot stress, indicating their multifunctional roles under different abiotic stresses. In our studies, more PhNF-Y genes responded to hot and cold stresses than to drought and salt. We concluded that PhNF-Y genes might play unknown, important roles in temperature stress. However, we could not rule out the possibility that different treatment concentrations and times result in dramatic expression changes under drought or salt treatment. In conclusion, the results of the stress response experiments, together with the tissue-specific expression, suggested that many PhNF-Y genes might respond to stresses through regulating root or leaf development. These PhNF-Y genes deserve further research and might be utilized for improving the abiotic stresses tolerance of petunias.

Identification, Alignments, and Phylogenetic Analysis of NF-Y Members
The sequences of 36 Arabidopsis NF-Y proteins were obtained from the Arabidopsis Information Resource (TAIR, https://www.arabidopsis.org/). The full length of all Arabidopsis NF-Y proteins were used as queries to blast petunia NF-Y protein sequences from Sol Genomics Network (SGN, https://www. solgenomics.net/). Alignments of full-length amino acid sequences of identified PhNF-Ys and AtNF-Ys were performed using ClustalX (http://www.clustal.org/) and BioEdit. Based on the sequence alignments generated by ClustalX, all potentially redundant NF-Y sequences were discarded. Phylogenetic trees were constructed using MEGA6 [56] and the neighbor-joining method with 1000 bootstrap replicates. The dendrograms were drawn by the online tool of iTOL (https://itol.embl.de/). Molecular weight (MW) and isoelectric point (pI) of the deduced amino acid sequences were predicted with the Sequences Manipulation Suite (https://www.genscript.com/sms2/index.html).

Conserved Motifs and Gene Structure Analyses
PhNF-Y conserved motifs were identified using the Multiple Em for Motif Elicitation (MEME, http://www.meme.sdsc.edu/meme/meme.html). The parameter settings were as follows: the maximum number of motifs was set to 10, the width range was 6-55, and other factors were at default selections [8].
Information about Genomic DNA Sequences and Coding Sequences (CDS) were obtained from the SGN database. The exon-intron structures of the PhNF-Ys were constructed using Gene Structure Display Server 2.0 (http://gsds.cbi.pku.edu.cn/).

miRNA Target Site Prediction
Genomic DNA sequences of 27 AtNF-Ys were used to predict the possible targets of all mature miRNA members using the online tool PsRNA server (http://plantgrn.noble.org/psRNATarget/). Parameters are set as follows: The maximum expectation value is 3.5, the length for complementarity scoring (hsp size) is 21 bp, target accessibility (UPE) is 25, flanking sequence length is 17 bp in upstream and 13 bp in downstream, and the translational inhibition range is 9-11NT.

Plant Material, Tissue, and Stress Treatment
Petunia "Ultra" seeds were grown in 7 cm diameter pots containing sterile nutrient soil and cultured in a room at 22 ± 3 • C under an LD cycle (14 h light/10 h dark). The 10-leaf stage seedlings were exposed to varieties of abiotic stresses in chambers, including cold stress (4 • C), hot stress (40 • C), drought stress (20% polyethylene glycol-6000 (PEG-6000)), and salt stress (500 mM NaCl). For the temperature treatment, seedlings were exposed at either 4 • C or 40 • C in a chamber with an LD cycle (14 h light/10 h dark), after which the second fully expanding leaves from the top were sampled. For the NaCl and PEG 6000 assays, seedlings were transferred to 22 • C water containing the stress agent in chambers at 22 • C under an LD cycle (14 h light/10 h dark), and the second fully expanding leaves from the top were sampled [57]. Samples were collected at 0 h, 1 h, 3 h, 6 h, and 12 h, immediately frozen in liquid nitrogen, and stored at −80 • C. Each treatment was replicated three times.
In the tissue expression pattern experiment, samples of roots, stems, and leaves were collected from seedlings without bloom. The first flowers of three individual plants were collected as flower samples. The samples were frozen in liquid nitrogen and stored at −80 • C.

Transcriptome Data Analysis
Total RNA samples were extracted from the above samples using TRIZOL (Invitrogn, Carlsbad, CA, USA) following the manufacturer's protocol [58]. The integrity of the RNA samples was confirmed by electrophoresis on 1% agarose gels, and the amount and quality of RNA were detected by NanoDrop 2000 Spectrophotometer (Thermo Scientific, Waltham, MA, USA) [37]. RNA-Seq libraries were prepared and sequenced using a HiSeq 2000 (Illumina) at the Shanghai Origingene Biotechnology Co. Ltd. (http://www.origin-gene.com/). RNA-seq data were processed and assembled as previously described [55]. After removing low-quality reads, adapters, and barcode sequences, high-quality clean reads were assembled de novo into contigs. The resulting contigs were blasted against Petunia axillaris draft genome sequence databases (https://www.solgenomics.net/organism/ Petunia_axillaris/genome) with a cutoff E value of 1 × 10 −5 . The Reads Per Kilobase of exon model per Million mapped reads (FPKM) values of all genes were determined with a treatment of log 2 (X), and subsequently used for generating heatmaps. Finally, a series of heat maps was generated using TBtools [59].
Supplementary Materials: The following are available online at http://www.mdpi.com/2223-7747/9/3/336/s1. Table S1: Identification of mature miRNAs and targeted NF-Y genes in Petunia. Table S2: Spectrophotometric analyses of RNA samples from plants under different stress treatments. Table S3: Spectrophotometric analyses of RNA samples from different tissues. Figure S1: Electropherograms of RNA samples.