Identification and Comprehensive Analysis of the Nuclear Factor-Y Family Genes Reveal Their Multiple Roles in Response to Nutrient Deficiencies in Brassica napus

Nuclear Factor-Y (NF-Y) transcription factors play vital roles in plant abiotic stress response. Here, the NF-Y family in Brassica napus, which is hyper-sensitive to nitrogen (N) deprivation, was comprehensively identified and systematically characterized. A total of 108 NF-Y family members were identified in B. napus and categorized into three subfamilies (38 NF-YA, 46 NF-YB and 24 NF-YC; part of the Arabidopsis NF-YC homologous genes had been lost during B. napus evolution). In addition, the expansion of the NF-Y family in B. napus was driven by whole-genome duplication and segmental duplication. Differed expression patterns of BnaNF-Ys were observed in response to multiple nutrient starvations. Thirty-four genes were regulated only in one nutrient deficient condition. Moreover, more BnaNF-YA genes were differentially expressed under nutrient limited environments compared to the BnaNF-YB and BnaNF-YC subfamilies. Sixteen hub genes responded diversely to N deprivation in five rapeseed tissues. In summary, our results laid a theoretical foundation for the follow-up functional study of the key NF-Y genes in B. napus in regulating nutrient homeostasis, especially N.


Introduction
Nuclear Factor-Y (NF-Y) transcription factors were discovered in yeast, mammals, plants, and other eukaryotes [1][2][3]. The NF-Y family members can specifically recognize and bind to the CCAAT sequence in the target gene's promoters [4]. The NF-Y transcription factor complex is composed of protein subunits from three unique subfamilies: NF-YA, NF-YB and NF-YC. A comparison of the protein sequences of three NF-Y subunits in different species, reveals that each subunit contains a conserved structural region. The NF-YA subunit has a carbon terminal region containing 56 amino acids, the NF-YB subunit a middle region containing 90 amino acids, and the NF-YC subunit a nitrogen (N) terminal region encompassing 84 amino acids [3]. The NF-Y family has been functionally annotated in a number of species. In yeast and animals, all NF-Y subunits contain only one or two members, but in plants, the number of genes encoding NF-Y subunits is expanded significantly [2]. In Arabidopsis thaliana, there are 36 NF-Y genes including 10, 13 and 13 members in the NF-YA, NF-YB, and NF-YC subunits, respectively. In theory, these subunits can be integrated into 1690 distinct complexes [5]. Likewise, in rice (Oryza sativa),

Physiological Response of B. napus to N, P and K Deficiencies
In order to explore the response of B. napus to diverse nutrient deficiencies, the B. napus cultivar "ZS11" was cultivated for 11 days under low N (LN, 100 µM), low P (LP, 5 µM) and low K (LK, 5 µM) conditions ( Figure 1). Compared with normal nutrient supply (CK), the growth of B. napus under nutrient limitation was significantly inhibited, especially under N deficient conditions (Figure 1a-c). As expected, shoot dry weight under LN, LP and LK conditions was reduced to 70.11%, 56.40% and 29.01% of that under CK condition, respectively. However, no difference was observed in roots (Figure 1b,c). Moreover, the taproot growth was the most sensitive under LN condition, which was about 59.91% longer than that under CK. The SPAD value of leaves decreased significantly under LN condition, while the opposite trend was observed under LP condition (Figure 1i). We further determined the photosynthetic characters of B. napus under the same conditions. Compared with CK, the transpiration rate, net photosynthesis rate, internal CO 2 and stomatal conductance were notably reduced under LN and LP conditions, indicating that the lack of N and P greatly inhibited B. napus photosynthesis (Figure 1d-g). For K deficiency, it led to 33.72% and 19.67% reduction in net photosynthesis rate and internal CO 2 while the stomatal conductance was significantly increased, by 70.04% compared to CK. However, little effect was observed on the transpiration rate (Figure 1d-g). In general, B. napus has diverse responses to different nutrient deficiencies, and is especially sensitive to N deficiency and less affected by LK.

Identification and Molecular Characterization of the NF-Y Family Genes
Previous reports showed that the NF-Y family genes could improve nutrient uptake and assimilation [37]. Thus, in order to explore the role of the NF-Y family genes in regulating nutritional deficiency tolerance in B. napus, we comprehensively identified the family members. Through multiple BLAST search against the homologue of 36 Arabidopsis NF-Y family members, a total of 108 NF-Y genes were detected in the whole genome of B. napus (Table S1). The NF-Y genes of B. napus included 38 genes of NF-YA subfamily, 46 genes of NF-YB subfamily and 24 genes of NF-YC subfamily (Tables S2 and S3). Meanwhile, 54 members of the NF-Y family (18 NF-YA, 24 NF-YB and 12 NF-YC) in B. rapa were identified. In B. oleracea 62 (22 NF-YA, 27 NF-YB and 13 NF-YC) genes were identified (Tables S2 and S3). Surprisingly, we found that several copies of the NF-YC subfamily members were lost in the whole genomes of three Brassica species, including NF-YC1/3/5/6/7/8/12 (Table S1). We then summarized the NF-Y family genes in 14 species. Besides the three Brassica species, Glycine max has the highest number of NF-Y genes, with 68 (21 NF-YA, 32 NF-YB and 15 NF-YC) members, while Prunus persica has the least NF-Y genes, with 24 (6 NF-YA, 12 NF-YB and 6 NF-YC) members (Table S2). Moreover, multiple sequence alignment was conducted based on the specific conserved protein sequences of three NF-Y subfamilies [3]. The results showed that the NF-YA and NF-YB subfamilies of the three species were relatively conserved, while the NF-YC subfamily varied largely, which may be related to the loss of some NF-YC subfamily genes during evolution ( Figure 2).
The gene ID, gene length (coding sequence), gene location, intron number, protein length and predicted subcellular localization of these genes are listed in Table S3. In general, the gene length of BnaNF-Ys varied from 324 bp to 2674 bp, and the genes encoding polypeptides varied from 107 to 642 amino acids in length. The intron number of BnaNF-Ys ranged from 0 to 7. The subcellular localization predicted by WoLF PSORT indicated that 86 BnaNF-Y proteins were localized in nucleus, 17 in chloroplast, two in mitochondria, two in cytoplasm, and one in plasma (Table S3)  The gene ID, gene length (coding sequence), gene location, intron number, protein length and predicted subcellular localization of these genes are listed in Table S3. In general, the gene length of BnaNF-Ys varied from 324 bp to 2674 bp, and the genes encoding polypeptides varied from 107 to 642 amino acids in length. The intron number of BnaNF-Ys ranged from 0 to 7. The subcellular localization predicted by WoLF PSORT indicated that 86 BnaNF-Y proteins were localized in nucleus, 17 in chloroplast, two in mitochondria, two in cytoplasm, and one in plasma (Table S3). The grand average of hydropathy (g) Stomatal conductance; (h) Taproot length; (i) SPAD value. CK: normal nutrient supply; LN: N deficiency (100 µM N); LP: P deficiency (5 µM P); LK: K deficiency (5 µM K). Error line represents variance (n = 4), different letters (a-d) indicate significant differences per Student's t test, * p < 0.05, ** p < 0.01.
(GRAVY) value of the BnaNF-Y proteins varied from −1.358 to 0.027, most of them negative except BnaC06.NF-YA3b and BnaC04.NF-YA3, indicating that most of the BnaNF-Ys are hydrophilic (Figure 3a). The molecular weight (MW) of the BnaNF-Y proteins varied from 12.211 kDa to 73.177 kDa (Figure 3b), with isoelectronic points (pI) values which ranged from 4.34 to 10.55. The pI of most NF-YA subfamily genes was more than seven, but the pI of most NF-YB and NF-YC subfamily genes were less than seven (Figure 3c).

Phylogenesis of the NF-Y Genes in B. napus
In order to explore the genetic relationships of the NF-Y genes between Arabidopsis and B. napus, the protein sequences of 36 AtNF-Ys and 108 BnaNF-Ys were used to construct a phylogenetic tree in MEGA7.0 ( Figure S1). The results showed that all of the NF-Y family members can be divided into three subgroups. Ten genes in Arabidopsis and 38 genes in B. napus were included in the NF-YA subfamily, 13 genes in Arabidopsis and 46 genes in B. napus were included in the NF-YB subfamily, and 13 genes in Arabidopsis and 24 genes in B. napus were included in the NF-YC subfamily. The NF-Y genes in the same subfamily were closely related to each other. Furthermore, some genes of the NF-Y group

Phylogenesis of the NF-Y Genes in B. napus
In order to explore the genetic relationships of the NF-Y genes between Arabidopsis and B. napus, the protein sequences of 36 AtNF-Ys and 108 BnaNF-Ys were used to construct a phylogenetic tree in MEGA7.0 ( Figure S1). The results showed that all of the NF-Y family members can be divided into three subgroups. Ten genes in Arabidopsis and 38 genes in B. napus were included in the NF-YA subfamily, 13 genes in Arabidopsis and 46 genes in B. napus were included in the NF-YB subfamily, and 13 genes in Arabidopsis and 24 genes in B. napus were included in the NF-YC subfamily. The NF-Y genes in the same subfamily were closely related to each other. Furthermore, some genes of the NF-Y group in Arabidopsis did not have homologous genes in B. napus, while some other genes of the Arabidopsis NF-Y groups had two or more homologous genes in B. napus. This may be due to the loss of certain genes and the homologous replication of chromosomes during evolution. To distinguish each subfamily of the NF-Y genes, all the BnaNF-Y members were renamed according to the Arabidopsis homologous genes following the international nomenclature from BnaA02.NF-YA1 to BnaC07.NF-YC13 ( Figure S1, Table S3).

Structural Analysis and Conserved Motif Analysis of the NF-Y Genes in B. napus and Arabidopsis
The gene structure of the 144 NF-Y family members in B. napus and Arabidopsis were further analyzed. The results showed that the gene structure varied largely among different subfamilies, while it was relatively conserved in the same family ( Figure 4, Table S3). The number of introns varied from 0 to 7 with an average number of 2.6. Twenty-two NF-Y genes had no intron, while BnaA03.NF-YB8 had eight introns, the largest number among all the members. Generally, the BnaNF-Y genes had similar numbers of introns as their counterparts in Arabidopsis. In addition, we analyzed the motif structures of the NF-Y family genes using their amino acid sequences in MEME. A total of 20 conserved motifs were identified, and the number of amino acids in each motif ranged from 6 to 50 (Figures 4 and S2). Similar to gene structure, nearly all the NF-Y genes in the same subfamily had relatively conserved motif composition and arrangement; however, it varied to a substantial degree in different subfamilies, especially among the NF-YA, NF-YB and NF-YC subfamilies. In addition, each subfamily contains specific motifs. In detail, the NF-YA subfamily specifically contains motifs 2, 4, 9, 14, 16 and 20, NF-YB subfamily specifically contains motif 1, 3, 7, 8, 13 and 18, and NF-YC subfamily specifically contains motif 6, 12, 15, 17 and 19, and they all have motif 10 ( Figure 4 and Figure S2). Taken together, the results of the structural analysis and conserved motif analysis strongly demonstrate the reliability of the genetic identification results.
To unravel the evolutionary processes of the 108 NF-Y genes in oilseed rape, the gene expansion patterns including tandem and segmental duplication were further analyzed using the coding sequences. A total of 144 segmental duplication events were identified using MCScanX methods ( Figure 5). The NF-Y genes in the A subgenome of B. napus had high homology with the NF-Y genes in the C subgenome. For example, two genes (BnaA03.NF-YB7 and BnaA03.NF-YC11) located on chromosome A03 had strong homology with two genes (BnaC03.NF-YB7 and BnaC03.NF-YC11a) located on the chromosome C03. Two genes (BnaA09.NF-YB2 and BnaA09.NF-YB6) located on chromosome A09 are highly homologous to two genes (BnaC09.NF-YB2 and BnaC09.NF-YB6) on chromosome C09. However, only two tandem duplication events were observed. These results indicate that the major driving force for the expansion of the NF-Y family genes in B. napus might be segmental duplication events.
Furthermore, a comparative syntenic map between the NF-Y genes in B. napus and Arabidopsis was constructed to further deduce the expansion mechanisms of the B. napus NF-Y family members. The results showed that there were strong orthologs among B. napus, B. rapa, B. oleracea and A. thaliana ( Figure 6). In detail, 110 and 187 pairs of collinear relationships were observed in the A subgenome of B. napus with A. thaliana and B. rapa, respectively. Additionally, the C subgenome of B. napus contained 114 and 219 pairs of syntenic relationships with A. thaliana and B. oleracea, respectively. In order to better understand the evolutionary constraints acting on the NF-Y genes between B. napus and A. thaliana, the nonsynonymous (Ka)/synonymous (Ks) ratios were estimated. The results showed that the Ka/Ks ratios of all the NF-Y gene pairs between B. napus and A. thaliana were less than one (Figure 3d-f, Table S4), suggesting that the NF-Y genes in rapeseed might have suffered robust purifying selective pressure during evolution.

Expression Patterns of the BnaNF-Y Genes in Response to Different Nutrient Deficiencies
To check the functions of the BnaNF-Y genes in response to various nutrient limitations, the expression profiles of the 108 BnaNF-Y members were analyzed under N, P or K deficiencies as well as sufficient nutrient supply conditions using RNA-seq data. The results demonstrated that the expression levels of the BnaNF-Y genes changed largely in response to different nutrient deficiencies and varied among different subfamilies (Table S5,  Figures 7 and S4a). In leaf, the expression of 61, 23 and 16 genes were significantly changed under N, P and K starvations, respectively ( Figure S3b). In root, 32, 25 and 20 genes were affected dramatically by N, P and K deficiencies, respectively ( Figure S3c). Some of the differentially expressed genes (DEGs) responded specifically to a certain nutrient depriva-tion, while others changed in response to two or three nutrient deprivations ( Figure S3b,c). Among all the DEGs, 86.89% (53/61), 69.57% (16/23) and 56.25% (9/16) of the BnaNF-Y genes were upregulated by N, P and K deficiencies in leaf, respectively ( Figure S3d). In root, 81.25% (26/32), 68.00% (17/25) and 40.00% (8/20) of the DEGs were enhanced by N, P and K deficiencies, respectively ( Figure S3e). The rest of the DEGs in leaf and root were downregulated by different nutrient limitations ( Figure S3d,e). Further analysis of the DEGs showed that 34 genes were regulated only in one nutrient deficiency condition in both leaf and root, including 26 under N deficiency, three under P deficiency, and five under K deficiency (Figure 7a and Table S5). A total of 26 DEGs responded specifically to a certain environment in one organ (Figure 7a). Fifteen DEGs were specifically upregulated by N deficiency in leaf, and six DEGs were specifically induced by N starvation in both leaf and root (Figure 7a), indicating their key roles in N deficiency response. Twelve BnaNF-Y genes were upregulated in one nutrient deficiency treatment but downregulated in another ( Figure S6

Expression Patterns of the BnaNF-Y Genes in Response to Different Nutrient Deficiencies
To check the functions of the BnaNF-Y genes in response to various nutrient limitations, the expression profiles of the 108 BnaNF-Y members were analyzed under N, P or K deficiencies as well as sufficient nutrient supply conditions using RNA-seq data. The results demonstrated that the expression levels of the BnaNF-Y genes changed largely in response to different nutrient deficiencies and varied among different subfamilies (Table  S5, Figures 7 and S4a). In leaf, the expression of 61, 23 and 16 genes were significantly changed under N, P and K starvations, respectively ( Figure S3b). In root, 32, 25 and 20 genes were affected dramatically by N, P and K deficiencies, respectively ( Figure S3c). Some of the differentially expressed genes (DEGs) responded specifically to a certain nutrient deprivation, while others changed in response to two or three nutrient deprivations ( Figure S3b and S3c). Among all the DEGs, 86.89% (53/61), 69.57% (16/23) and 56.25% (9/16) of the BnaNF-Y genes were upregulated by N, P and K deficiencies in leaf, respectively ( Figure S3d). In root, 81.25% (26/32), 68.00% (17/25) and 40.00% (8/20) of the DEGs were enhanced by N, P and K deficiencies, respectively ( Figure S3e). The rest of the DEGs in leaf and root were downregulated by different nutrient limitations ( Figure S3d and S3e). Further analysis of the DEGs showed that 34 genes were regulated only in one nutrient deficiency condition in both leaf and root, including 26 under N deficiency, three under P deficiency, and five under K deficiency (Figures 7a and S5). A total of 26 DEGs responded specifically to a certain environment in one organ (Figure 7a). Fifteen DEGs were specifically upregulated by N deficiency in leaf, and six DEGs were specifically induced by N starvation in both leaf and root (Figure 7a), indicating their key roles in N deficiency response. Twelve BnaNF-Y genes were upregulated in one nutrient deficiency Since the NF-Y family is a heterotrimeric complex with three distinct subfamilies, we further analyzed the expression pattern of each subfamily in response to N, P and K limitations. The results showed that 68.42% (26/38), 36.96% (17/46) and 33.33% (8/24) of the NF-YA, NF-YB and NF-YC genes, respectively, were differentially expressed under nutrient limitation conditions (Figure 7b-d). To further characterize the response of the three BnaNF-Y subfamilies to N, P and K starvations in two organs, we calculated the ratios of DEGs under different contexts and visualized the data using Sankey diagrams ( Figure 8). In terms of the three NF-Y subfamilies, the differentially expressed NF-YA, NF-YB and NF-YC genes accounted for 41% (73/177), 31% (55/177) and 28% (49/177), respectively (Figure 8a). For the two organs, 44% (77/177) DEGs were detected in root, and 56% (100/177) in leaf (Figure 8b). Moreover, 53% (93/177), 27% (48/177) and 20% (36/177) of the total DEGs responded to N, P and K deficiencies, respectively (Figure 8c). These results suggest that the NF-YA subfamily may play a more important role in response to nutrient deprivation as compared to the NF-YB and NF-YC subfamilies, especially in N deficient environments.   NF-YC genes accounted for 41% (73/177), 31% (55/177) and 28% (49/177), respectively (Figure 8a). For the two organs, 44% (77/177) DEGs were detected in root, and 56% (100/177) in leaf (Figure 8b). Moreover, 53% (93/177), 27% (48/177) and 20% (36/177) of the total DEGs responded to N, P and K deficiencies, respectively (Figure 8c). These results suggest that the NF-YA subfamily may play a more important role in response to nutrient deprivation as compared to the NF-YB and NF-YC subfamilies, especially in N deficient environments.

Discovery of the NF-Y Hub Genes in B. napus and Expression Analysis of these Genes among Different Tissues under N Limitation
To reveal the coexpression relationships among the BnaNF-Y genes, we calculated the interaction weight of the target genes based on the fragments per kilobase million (FPKM) values of the RNA-seq data under N, P and K starvation conditions ( Figure S7). Generally, we identified 1091 pairs of coexpression relationships among the 108 BnaNF-Y family genes. Then, the BnaNF-Y genes with the ten strongest interaction relationships with other BnaNF-Y members were considered as the hub genes, as displayed in Figure S7. As shown in Figure 8, more than half (53%) of the DEGs were detected under N deficiency. Hence, to understand the detailed expression profiles of the BnaNF-Y family genes in response to N limitation, we further investigated the expression levels of the 16 core genes in root, hypocotyl, basal node, petiole, fully expanded leaf and new leaf of B. napus using quantitative real-time PCR (qRT-PCR). The results illustrated that the mRNA levels of these genes in oilseed rape were basically consistent with the results from the RNA-Seq data ( Figure 9, Table S5). Some genes, such as BnaC05.NF-YA2b, BnaCnn.NF-YA3 and BnaC06.NF-YA3a, were specifically induced by N starvation in root and some had high expression levels in fully expanded leaf, such as BnaC08.NF-YC4 and BnaA09.NF-YC9. Additionally, the core NF-Y genes were expressed differently in tissues other than root and leaf in response to N deficiency (Figure 9). For example, the expression levels of BnaA01.NF-YA6, BnaC08.NF-YC9, BnaC05.NF-YC9, BnaC09.NF-YB12, BnaA01.NF-YC11 and BnaC01.NF-YC11 were significantly induced in petiole by N deficiency. RNA was extracted separately from five tissues including root, hypocotyl, petiole, fully expanded leaf and new leaf under N sufficient and deficient conditions. +N: N sufficient condition; -N: N free (0 μM N) condition. Each sample included three independent biological replicates. * and ** indicate significant differences at p < 0.05 and p < 0.01 per student's t test, respectively.

Subcellular Localization of Four NF-Y Hub Genes in Rapeseed
We further selected four NF-Y hub genes in rapeseed (BnaA03.NF-YA2, BnaC01.NF-YA6, BnaC03.NF-YC11b and BnaC05.NF-YA2b) to explore their subcellular localization. The full length cDNA of four BnaNF-Y genes was isolated and integrated with PMDC83 Figure 9. The expression levels of the 16 selected BnaNF-Y genes in five tissues of B. napus in response to nitrogen (N) deficiency by qRT-PCR. Fourteen-day-old rapeseed seedlings were cultivated in a nutrient solution free of N for six days. RNA was extracted separately from five tissues including root, hypocotyl, petiole, fully expanded leaf and new leaf under N sufficient and deficient conditions. +N: N sufficient condition; −N: N free (0 µM N) condition. Each sample included three independent biological replicates. * and ** indicate significant differences at p < 0.05 and p < 0.01 per student's t test, respectively.

Subcellular Localization of Four NF-Y Hub Genes in Rapeseed
We further selected four NF-Y hub genes in rapeseed (BnaA03.NF-YA2, BnaC01.NF-YA6, BnaC03.NF-YC11b and BnaC05.NF-YA2b) to explore their subcellular localization. The full length cDNA of four BnaNF-Y genes was isolated and integrated with PMDC83 vector to form 35S::BnaNF-Ys-GFP protein structure. Then, these constructed vectors were injected into tobacco leaves. Meanwhile, we used nuclear dye 4 ,6-diamidino-2phenylindole (DAPI) for auxiliary labeling of nuclei. The results suggested that all of the GFP fusion protein signals were completely merged with the DAPI signals ( Figure 10). Consistent with their characteristics as transcription factors, these BnaNF-Y genes were located in the nucleus. Combined with our above analysis, these results jointly emphasized the important roles of NF-Y genes in regulating plant responses to nutrient deficiencies, especially N deprivation.

Discussion
The NF-Y genes are trimeric transcription factors in plants, containing the NF-YA, NF-YB and NF-YC subfamilies which include multiple genes [2]. Recently, the NF-Y families have been identified and characterized in several plant species, including Arabidopsis [5], rice [6], soybean [7], sorghum [10], barley [13], apple [15] and more. Previously, the NF-Y genes were identified and characterized in B. napus based on an expressed sequence tag database and the completed B. napus genome sequence in the CoGe database [42,43]. However, there had been no comprehensive genome-wide investigation or systematic characterization in Brassica species. In the present study, 62, 54 and 108 NF-Y genes were identified in B. rapa, B. oleracea and B. napus, respectively (Table S1). Especially in B. napus, the number of NF-Y family members is far greater than that found in other species (Table  S2) or in the two previous reports on B. napus [42,43]. Therefore, substantial changes via evolution and selection must have occurred during the evolution processes of Brassica species. Moreover, we found that such purifying selective pressure may have been a strong driving force in the evolution of Brassica species (Figure 3e, Table S4). Interestingly,

Discussion
The NF-Y genes are trimeric transcription factors in plants, containing the NF-YA, NF-YB and NF-YC subfamilies which include multiple genes [2]. Recently, the NF-Y families have been identified and characterized in several plant species, including Arabidopsis [5], rice [6], soybean [7], sorghum [10], barley [13], apple [15] and more. Previously, the NF-Y genes were identified and characterized in B. napus based on an expressed sequence tag database and the completed B. napus genome sequence in the CoGe database [42,43]. However, there had been no comprehensive genome-wide investigation or systematic characterization in Brassica species. In the present study, 62, 54 and 108 NF-Y genes were identified in B. rapa, B. oleracea and B. napus, respectively (Table S1). Especially in B. napus, the number of NF-Y family members is far greater than that found in other species (Table S2) or in the two previous reports on B. napus [42,43]. Therefore, substantial changes via evolution and selection must have occurred during the evolution processes of Brassica species. Moreover, we found that such purifying selective pressure may have been a strong driving force in the evolution of Brassica species (Figure 3e, Table S4). Interestingly, and similar to other species, the Ka/Ks ratios of all the NF-Y gene pairs were less than one in B. napus, suggesting that the NF-Y family genes may be conserved in function [10]. The gene replication pattern analysis showed that 144 segmental duplication events occurred in the B. napus genome; however, only two tandem duplication pairs were observed ( Figure 5). In rice, the NF-Y family expansion was mainly driven by segmental duplication events, which indicates that the NF-Y family might have similar expansion patterns in dicotyledons and monocots [6,32].
Moreover, the conserved domain analysis of the NF-Y subfamilies in three Brassica species and Arabidopsis indicates that the NF-YC subfamily excepted, the NF-YA and NF-YB subfamilies are relatively conservative (Figure 2). Both the NF-YA and NF-YC subfamilies in B. oleracea showed large differentiation compared with others, which further provides strong evidence for an obvious selection during the evolution of Brassica species. Similarly, the gene structure and motifs differed across different subfamilies in B. napus, which is consistent with previously reports in sorghum [10], citrus [12] and tea [44]. However, the NF-Y genes in the same subfamily, such as NF-YB11s, had a similar composition and arrangement of structure and motifs, suggesting that the function of genes in each subfamily may be relatively conserved in B. napus (Figure 4, Table S3). Generally, these results put forward new insights into the phylogenetic mechanism of the NF-Y family in Brassica species.
Over the past decade, the NF-Y genes as the key regulators of drought resistance in plants have been widely investigated [34]. Previous studies showed that over-expression of the NF-Y genes could enhance drought tolerance in Arabidopsis [45,46] and Populus [47]. However, fewer studies focused on the relationships between the NF-Y genes and plant stress resistance under nutrient limitation conditions. Here, we investigated the expression of the NF-Y family in B. napus in response to multiple nutrient limitations. The results suggest that the mRNA levels of the BnaNF-Y family members display various expression patterns under distinct nutrient deficiencies, which may be related to their complex structure and motif compositions (Figures 7 and S4a, Table S5). In soybean, NF-YA12 was upregulated by abscisic acid, NaCl and cold stresses [7]. Likewise, several NF-Y genes were induced under multiple stress conditions in sorghum [11]. Interestingly, except for the BnaNF-Y genes (whose expression levels were regulated by different nutrient deficiencies) we also found that the expression changes of 34 genes were stress-specific (Figures 7 and S5). Among them, the expression levels of seventeen, one and three BnaNF-Y genes were specifically regulated in leaf under N, P and K deprivations, respectively, while only five genes which specifically responded to a certain nutrient deficiency were observed in roots (Figure 7). Thus, that these NF-Y genes responded specifically to a certain tissue under a specific nutrient deficiency may indicate a unique role in improving stress resistance in B. napus. In addition, reports show that ABA-dependent or independent signaling pathways are involved in the response of many NF-Y genes to abiotic stresses [31,48]. In the current research, the expression levels of 12 BnaNF-Y genes were inversely regulated under multiple nutrient deprivations, indicating that they might participate in completely different regulatory pathways in order to regulate plant growth. However, whether the regulatory pathways associated with the BnaNF-Y genes are ABA-dependent or independent still needs to be further verified.
In sorghum, the expression profiles of the three NF-Y subfamily members differed largely in response to various abiotic stresses, suggesting that different NF-Y subfamily members might have specific subfunctionalization [11]. Similar results were observed here, in that the BnaNF-YA family genes had the most obvious response to multiple nutrient starvations (Figures 7b-d and 8). Consistent with previous studies, the NF-YA family members were the genes most related to regulation of plant nutrient deficiency responses. For example, three Arabidopsis NF-YA genes (AtNF-YA2, AtNF-YA7 and AtNF-YA10) might be post-transcriptionally regulated by miR169 to regulate plant tolerance to P deficiency [35]. In addition, the absorption of N and P as well as grain yield were significantly enhanced in the TaNFYA-B1 over-expression lines under N and P deficient levels in wheat [37]. Therefore, compared with other subfamily members the NF-YA subfamily might act as a more critical player in response to nutrient deficiencies in B. napus. Among these, the responses of the BnaNF-YA family to nutritional starvations were slightly changed compare to that in Arabidopsis. For instance, the expressions of AtNF-YA2, AtNF-YA3 and AtNF-YA5 were be induced by P deficiency [49], while their homologous genes in B. napus were nearly unaffected by P limitation (Table S5). However, the three NF-YA genes were upregulated by N deprivation in both Arabidopsis and B. napus [36], indicating that they still retained some important functions in B. napus (Table S5). N deficiency significantly inhibited the growth of B. napus compared with other nutrient deficiencies, and strongly regulated the expression of most NF-Y genes (Figures 1 and 8). Thus, the results fully suggest that the BnaNF-Y family genes may have important roles in response to N limitation in oilseed rape, especially the NF-YA subfamily (Figure 8).
The homologous genes of other species contained more than one copy compared with Arabidopsis, and might have experienced functional differentiation. In B. napus, the complex evolutionary processes of the NF-Y family genes give rise to its functional diversification under multiple nutrient deficiencies. To investigate the expression changes of the BnaNF-Y genes in response to different nutrient limitations in more details, we identified 10 hub genes based on the coexpression network analysis ( Figure S7). Besides the NF-YA and NF-YC subfamilies that have been reported to be closely related to nutrient deficiencies, another four BnaNF-YB genes also displayed strong relationships ( Figure S7). Reports showed that AtNF-YA2, AtNF-YA3, AtNF-YA5 and AtNF-YA8 were post-transcriptionally regulated by miR169a to enhance plant adaptability to low N availability in soil, while AtNF-YC4 was related to N and carbon distribution in plants [36,49]. Furthermore, the sensitivity of the NF-Y family genes to N starvation was higher than that to P and K starvations in B. napus (Figure 8). Thus, we further investigated the mRNA levels of the BnaNF-Y hub genes and other BnaNF-Y genes specifically responding to N deficiency in different tissues. In line with other species, most NF-Y genes revealed intense tissue specificity, and seven BnaNF-YA subfamily genes were mainly induced by N deficiency in root (Figure 9). Such tissue-specific expression patterns have been observed in a few plants [11,50]. Interestingly, BnaA10.NF-YC11 was significantly induced in hypocotyl by N deficiency, but was inhibited in root ( Figure 9), suggesting that this gene might play different roles in specific tissues under N deficiency. In addition, unlike other genes, BnaC08.NF-YC4 was significantly inhibited in petiole and fully expanded leaf ( Figure 9). Meanwhile, BnaNF-Ys as a transcription factor was located in the nucleus ( Figure 10). This evidence indicates that the BnaNF-Y genes might participate in modulating nutritional stress tolerance through positive or negative mechanisms, which have been testified to in Arabidopsis [35].

Identification of the NF-Y Family Genes in B. napus
The NF-Y family members of B. napus were obtained using BLASTP in the CNS-Genoscope genomic database (http://www.genoscope.cns.fr/brassicanapus/, accessed on 1 September 2019) based on the protein sequences of 36 Arabidopsis NF-Y genes from TAIR database (https://www.arabidopsis.org/, accessed on 1 September 2019). The homologous genes in three Brassica species were obtained based on E-value < 1 × 10 −5 . The redundant gene sequences of the identified Brassica species were removed manually. In order to ensure the reliability of identified genes, the Hidden Markov Model (HMM) in the Pfam database (http://pfam.xfam.org/, accessed on 1 September 2019) was employed with E-value < 0.001 to eliminate out-of-standard genes, which also did not contain the complete CBFD-NFYs domain (PF02045 and PF00808). The full genomic sequences, coding sequences and protein sequences of B. napus NF-Y members were derived from the B. napus genome database. For B. rapa and B. oleracea, the same methods described above were used to obtain the sequence information of NF-Y family members from the Brassica Database (BRAD, http://brassicadb.org/brad, accessed on 1 September 2019).

Chromosomal Location and Gene Duplication Analysis of the NF-Y Genes in B. napus
The chromosomal locations of the NF-Y genes were determined from the B. napus genome database and were displayed by TBtools [51]. Gene duplication events and collinearity relationships of the NF-Y genes in Arabidopsis and three Brassica species were analyzed using Multiple Collinearity Scan toolkit (MCScanX) [52]. The standard for identifying potential gene duplication was defined by the following two points: (a) the comparison region covers more than 75% of the gene; and (b) the similarity of the aligned regions is more than 75%. Likewise, the gene duplication events and collinearity relationships of the NF-Y genes between B. oleracea and B. napus or between B. rapa and B. napus were analyzed using the same methods described above. MCScanX was employed to construct the map which was used to demonstrate the syntenic relationships between the B. napus NF-Y family genes and its ancestors in Arabidopsis, B. rapa and B. oleracea.

Phylogenetic Relationships and Evolutionary Pressure Analysis of the B. napus NF-Y Genes
MEGA7 software was employed to conduct the protein sequence alignments with default parameters, and the evolutionary tree analysis tool was used to build phylogenetic trees with the neighbor joining (NJ) method based on 1000 bootstraps [53]. In order to analyze the selection pressure of BnaNF-Y family evolution, values of the Ka, Ks and Ka/Ks ratio between Arabidopsis and B. napus were calculated based on the sequence alignments. Pairwise alignments of the gene CDS without stop codon were performed to calculate Ka/Ks ratio by the Ka/Ks calculator in TBtools software [51]. Generally, a Ka/Ks ratio greater than one indicates positive selection, while a Ka/Ks ratio less than one suggests a functional constraint; equal to one indicates neutral selection [54].

Molecular Characterization Analysis of the B. napus NF-Y Family Genes
The full genomic sequences, coding sequences and protein sequences of the BnaNF-Y genes downloaded from the CNS-Genoscope genomic database were used to analyze gene structure by multiple alignments. The gene structure was then visualized in TBtools [51]. The potential conserved motifs of the NF-Y family genes in Arabidopsis and B. napus were detected using Multiple Expectation Maximization for Motif Elicitation program (MEME, http://meme-suite.org/tools/meme, accessed on 3 June 2019) with the following four parameters: (a) the maximum input data was 6000; (b) the analysis mode was anr; (c) the optimum motif width ranged from 6 to 50; and (d) the maximum motif number was 20 [55]. Moreover, the MW and pI of BnaNF-Y proteins were calculated using ExPASy tool (http://www.expasy.org/tools/, accessed on 3 June 2019), and the GRAVY values of the proteins were analyzed using the PROTPARAM tool (http://web.expasy.org/protparam, accessed on 3 June 2019). WoLF PSORT server (https://wolfpsort.hgc.jp/, accessed on 3 June 2019) was employed to predict the subcellular localizations of the BnaNF-Y proteins.

Heat Map and Coexpression Networks of the BnaNF-Y Family Genes Using RNA-seq Data
The heat map and gene coexpression network analysis were performed based on the RNA-seq data. In the RNA-seq experiment, 14-day-old rapeseed seedlings were treated with N-free, P-free and K-free nutrient solution for 10 days. Then, the fully expanded leaf and root were harvested separately for RNA extraction. After harvest, all samples were frozen immediately in liquid N and then stored at −80 • C. Each treatment included three independent biological replicates. Illumina HiSeq 2500 platform (Illumina, San Diego, CA, USA) was applied to perform transcriptome sequencing. The FPKM values indicating the transcript abundance of all genes under control and N/P/K deficient conditions were obtained based on gene length and the reads mapped to the gene. Moreover, DEGs were identified based on p value < 0.05 between control and treatments. The interactions and weight values of the target genes were analyzed using DeGNServer online tools (http://plantgrn.noble.org/DeGNServer/Analysis.jsp, accessed on 9 October 2019) based on the FPKM values. Cytoscape software was employed to visualize the coexpression networks of the BnaNF-Y genes [56].

Plant Materials and Growth Conditions
A B. napus cultivar named "ZS11" was planted in an illuminated growth room (22 • C; 16 h light/8 h dark) with a light intensity 300-320 mmol protons m −2 s −1 . The seeds were first disinfected and soaked in deionized water in the dark for two days, then transferred to 0.5 mM CaCl 2 solution for three days to promote root development. Finally, modified Hoagland's solution with pH 5.8 was used for seedling growth, containing 5. After nine days' growth, plants were treated with different nutrient concentrations for six days, including N, P and K deficiencies. CaSO 4 and K 2 SO 4 were used to substitute Ca(NO 3 ) 2 and KNO 3 in N starvation treatment, respectively. K 2 SO 4 was replaced by KH 2 PO 4 in P starvation treatment. In K starvation treatment, NaH 2 PO 4 and NaNO 3 were used to replace KH 2 PO 4 and KNO 3 , respectively. Tissue samples were collected after the treatments and stored at −80 • C until RNA extraction. Three biological replicates were prepared for each treatment. For phenotypic analysis, the B. napus cultivar "ZS11" was cultivated for 11 days under low N (LN, 100 µM), low P (LP, 5 µM) and low K (LK, 5 µM) conditions ( Figure 1). Then, the taproot root length, shoot and root dry weight were measured. SPAD value was measured by SPAD meter (SPAD 502 plus, Minolta, Japan). Photosynthetic parameters were detected by Photosynthesis System (CIRAS-3, America).

Expression Analysis of the BnaNF-Y Genes under Multiple Nutrient Deficiencies
A RNeasy Plant Mini Kit (Qiagen) was used to extract the total RNA of each sample, and a First Strand cDNA synthesis kit (Toyobo) was employed to synthesize cDNA. The qRT-PCR was carried out in a 10 µL volume containing 2.5 µL cDNA, 0.2 µL primers, 5 µL Hieff qPCR SYBR Green Master Mix (Low Rox). The thermal cycle was as follow: 95 • C 5 min; 40 cycles of 95 • C for 10 s, 60 • C for 20 s, 72 • C for 20 s. qRT-PCR was performed using the QuantStudio 6 Flex instrument (Life Technologies, USA). The relative gene expression levels were calculated according to the 2 −∆∆Ct method [57]. The housekeeping gene EF1-α (Accession number DQ312264) was used as an internal standard [58]. The primers used for qRT-PCR in this study are listed in Table S6.

Subcellular Localization Analysis of the BnaNF-Y Genes
For subcellular localization, the coding regions of four BnaNF-Y genes without stop codons amplified from "ZS11" were ligated separately into the vector (PMDC83) driven by the CaMV35S promoter. DAPI was used to stain nuclei; they were cultured on LB agar supplemented with selection antibiotics and incubated at 28 • C for two days. The confluent Agrobacterium with the target vector were resuspended in infiltration buffer (10 mM MgCl 2 , 0.5M MES (pH 5.6), 100 mM acetosyringone) to an optical density (600 nm) of 1.0, and incubated at room temperature without shaking for 3 h before infiltration. Approximately 2 mL of the Agrobacterium mixture was then infiltrated into the back of 3-4 young leaves of tobacco. The subcellular localization assay was performed 2 d after inoculation [59]. A fluorescence microscope was used to record the images.

Statistical Analysis
Student's t test was used to perform the statistics in SPSS 22 (IBM, Chicago, IL, USA). Significance of differences was defined as * p < 0.05, ** p < 0.01.

Supplementary Materials:
The following are available online at https://www.mdpi.com/article/10 .3390/ijms221910354/s1, Figure S1: Phylogenetic tree of the NF-Y family in B. napus and Arabidopsis. The analysis includes 144 protein sequences, including 108 from B. napus (white star) and 36 from Arabidopsis (light green circle). The different colors mean different subfamilies. Figure S2: The sequence information of the NF-Y family protein for each motif in Arabidopsis and B. napus. Figure S3: Chromosomal distribution of the NF-Y family members in B. napus. The localization of each BnaNF-Y gene on the chromosome was represented via colored fonts. Bar = 5 Mb. Figure S4: The expression map of all B. napus NF-Y genes: (a) The heat map of all B. napus NF-Y genes. The color scale is shown above the expression map. The data in the heat map was normalized (z-scores) using FPKM values and then visualized by Tbtools. The higher the value of the color scale, the higher the expression, and vice versa; (b) Venn diagram of the differentially expressed genes (DEGs) among three nutrient deficiencies in leaf; (c) Venn diagram of the DEGs among three nutrient deprivations in root; (d) The number of total, up-and downregulated DEGs in response to three nutrient deficiencies in leaf; (e) The number of total, up-and downregulated DEGs in response to three nutrient limitations in root. CK: sufficient nutrient supply; -N: nitrogen deficiency condition; -P: phosphorus deficiency condition; -K: potassium deficiency condition. Figure S5: Heat map of the B. napus NF-Y family genes differentially expressed in only one treatment based on the RNA-seq data. The color scale is shown above each subfamily. The data in the heat map was normalized (z-scores) using FPKM values and then visualized by TBtools. The higher the value of the color scale, the higher the expression, and vice versa. CK: sufficient nutrient supply; -N: nitrogen deficiency condition; -P: phosphorus deficiency condition; -K: potassium deficiency condition. * and ** indicate significant differences at p < 0.05 and p < 0.01 per student's t test, respectively. Figure S6: Heat map of the B. napus NF-Y family genes differentially expressed in two or three treatments based on the RNA-seq data. The color scale is shown above each subfamily. The data in the heat map was normalized (z-scores) using FPKM values and then visualized by Tbtools. The higher the value of the color scale, the higher the expression, and vice versa. CK: sufficient nutrient supply; -N: nitrogen deficiency condition; -P: phosphorus deficiency condition; -K: potassium deficiency condition. * and ** indicate significant differences at p < 0.05 and p < 0.01 per Student's t test, respectively. Figure S7: Co-expression networks of the BnaNF-Y genes. Cycle nodes represent genes, and the size and color of the nodes represents the power of the interrelation among the nodes by degree value. The width and the color of the lines between two nodes represent interactions between two genes. The hub NF-Y genes are located in the center of the network, while the 10 most coexpressed genes are displayed in the network. Table  S1: Statistics of the NF-Y genes in each subfamily in Arabidopsis and three Brassica species. Table  S2: Copy number variations of the NF-Y family in higher plants. Table S3: Characteristics of the NF-Y family members in B. napus and subcellular localization prediction. Table S4: Non-synonymous (Ka) and synonymous (Ks) nucleotide substitution rates for the NF-Y coding loci in Arabidopsis and B. napus. Table S5. FPKM values of the NF-Y family genes in leaf and root of B. napus under N, P and K deficiencies based on the RNA-seq data. Table S6: Primers used for quantitative real-time PCR in this study.

Data Availability Statement:
The data presented in this study are available on request from the corresponding authors.

Conflicts of Interest:
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.