Identification and Characterization of Abiotic Stress–Responsive NF-YB Family Genes in Medicago

Nuclear factor YB (NF-YB) are plant-specific transcription factors that play a critical regulatory role in plant growth and development as well as in plant resistance against various stresses. In this study, a total of 49 NF-YB genes were identified from the genomes of Medicago truncatula and Medicago sativa. Multiple sequence alignment analysis showed that all of these NF-YB members contain DNA binding domain, NF-YA interaction domain and NF-YC interaction domain. Phylogenetic analysis suggested that these NF-YB proteins could be classified into five distinct clusters. We also analyzed the exon–intron organizations and conserved motifs of these NF-YB genes and their deduced proteins. We also found many stress-related cis-acting elements in their promoter region. In addition, analyses on genechip for M. truncatula and transcriptome data for M. sativa indicated that these NF-YB genes exhibited a distinct expression pattern in various tissues; many of these could be induced by drought and/or salt treatments. In particular, RT-qPCR analysis revealed that the expression levels of gene pairs MsNF-YB27/MtNF-YB15 and MsNF-YB28/MtNF-YB16 were significantly up-regulated under NaCl and mannitol treatments, indicating that they are most likely involved in salt and drought stress response. Taken together, our study on NF-YB family genes in Medicago is valuable for their functional characterization, as well as for the application of NF-YB genes in genetic breeding for high-yield and high-resistance alfalfa.


Introduction
Nuclear factor Y (NF-Y), also called heme-activated protein (HAP) or CCAAT binding factor (CBF), can be found in almost all eukaryotes [1,2]. The NF-Y Transcription Factors binds to cis-elements with the conserved core sequence CCAAT to activate or inhibit the expression of related functional genes in metabolism [3]. The NF-Y family protein consists of three different subunits: NF-YA, NF-YB, and NF-YC [4]. These three subunits are distinguished by conserved structure and protein length, with NF-YAs being longer than NF-YB and NF-YC. In general, NF-YA contains two structural domains, A1 and A2 [5], while the protein structures of NF-YB and NF-YC are similar to H2B and H2A histones, respectively [5]. Subunits NF-YB and NF-YC form dimers in the cytoplasm and then bind to the NF-YA protein and form a trimer in the nucleus [6]. Interestingly, in yeast and mammals, each NF-Y subunit is encoded by a single gene; but in plants, it is encoded by multiple genes, and the number of genes encoding individual subunits are also different from species to species [1].
Recent studies have shown that individual NF-Y subunits in plants are involved in many important growth processes, especially in embryogenesis [7], seed maturation [8], chloroplast synthesis [9], and tissue division [10]. It is worth noting that a large body of

Identification of NF-YB Family Genes in Medicago
The conserved mode of NF-YB and the NF-Y proteins from A. thaliana and O. sativa were used as queries to identify NF-YB genes in M. truncatula and M. sativa. In total, 21 MtNF-YB and 28 MsNF-YB genes were identified and designated according to their location on chromosomes, as shown in Table 1. Previously, 24 MtNF-YB members were deposited in PlantTFDB with V4 genome sequencing results, and three of them were found to be duplicate sequences with the V5 genome sequencing data in the present study. Therefore, these 21 MtNF-YB were named according to their positions on the chromosome and were listed along with the corresponding gene locus deposited in PlantTFDB in Table 1. MtNF-YB and MsNF-YB genes encode proteins ranging from 91 to 257 and 78 to 289 amino acids in length, respectively. The detailed information of NF-YB family genes in Medicago, including sequence ID, chromosome location, amino acid length (aa), protein isoelectric point (PI) value and protein molecular weight (MW) (kDa) are listed in Table 1. Moreover, the corresponding NF-YB homologous genes of M. truncatula or M. sativa were identified based on sequence alignment. Subcellular location analysis showed that NF-YB proteins from M. truncatula and M. sativa were predicated to be located in nuclear (20 out of 21 MtNF-YB; 22 out of 28 MsNF-YB) or extracellular areas. Recent studies demonstrated that miRNAs located in the nucleus could act on CREs (e.g., promoters and enhancers) to regulate gene expression (activation or suppression) [27]. Those localized outside the cell may be closely related to cellular signal reception and transduction [28].

Multiple Sequence Alignment of NF-YB Genes in Medicago
To further investigate the conserved regions of NF-YB subfamilies in Medicago, multiple protein sequence alignments were analyzed using MEGA-X, and displayed via jalview ( Figure 1). The results illustrated that the core regions for MtNF-YB and MsNF-YB were 92 and 88 amino acids in length, respectively, which were close to the average length reported previously [29]. In addition, most NF-YBs from M. truncatula and M. sativa also contained DNA binding α1, NF-YA interaction α2, and NF-YC interaction α3/αC, which have a similar structure as the histone-fold motif (HFM) of the core histone H2B [30]. Among them, MtNF-YB3 lacks the αC domain, MsNF-YB21 and 23 lack α3 and αC domains, MsNF-YB25 lacks the α1 domain, and MsNF-YB24 lacks α1 and α3 domains, which may be due to the genome sequences.

Phylogenetic Analysis of NF-YB Proteins
To better understand the evolutionary relationships of the NF-YB gene, a neighborjoining (NJ) phylogenetic tree was constructed with the NF-YB proteins from M. truncatula, M. sativa, A. thaliana, P. vulgaris, P. sativum, and T. pratense. The details are shown in Table S2. The NF-YB proteins were divided into six clades, designated as A, B, C, D, E, and F. The F clade was the largest group, containing 69 NF-YB proteins, whereas the B and C clade were the smallest, consisting of only one member (only MsNF-YB25 in group B and PcNF-YB2 in group C), indicating that NF-YB proteins were distributed unevenly in different clades ( Figure 2). Among them, three AtNF-YB members were found in each of clusters A and E. Notably, only two members from M. sativa (MsNF-YB2, 13) and M. truncatula (MtNF-YB2, 13) were present in cluster D, demonstrating that these genes may play specific roles in the evolutionary process ( Figure 2).

Analysis of Gene Structure and Conserved Motifs of NF-YB Genes in Medicago
To comprehensively study the function of NF-YB genes, we performed analysis on gene structure and conserved motifs. We performed another phylogenetic analyses for NF-YB from M. sativa and M. truncatula and found they were regrouped into five clusters. As shown in Figure 3a, all NF-YB proteins were classified into five clades that were consistent with the phylogenetic relationships, as illustrated in Figure 2. The MEME analysis tool was used to predict conserved motifs in NF-YB genes (Figure 3b). A total of 20 motifs were identified; we found that all 21 MtNF-YB members contained motif 1. In total, 23 out of 28 MsNF-YB members contained motif 1, except for MsNF-YB24, 23, 21, 5, and 15. In addition, most NF-YB contained similar motifs, for example, motifs 2 and 3 were widely distributed in most NF-YB. We also found that NF-YB proteins with close phylogenetic relationships exhibited similar motif arrangements, for example, motifs 5, 11, 16, 17, and 18 only present in cluster B. The specificity of these motifs may result in functional differences among NF-YB within each cluster. To elucidate the gene structure of NF-YB family genes, we compared the coding sequences with their corresponding genomic sequences to determine the positions of exons and introns. As shown in Figure 3c, the numbers of exons ranged from one to seven, where genes with one exon accounted for 53% of the total NF-YB genes, most of which were from B and C clades. These results indicate that the intron/exon distribution in NF-YB genes from Medicago is highly variable. In addition, 29 genes lacked one or two of the 5 -or 3 -untranslated region (UTR), particularly in cluster E, which may be due to the genome annotation.

Analyses of Chromosomal Distribution and Synteny of NF-YB Genes in Medicago
All NF-YB genes were unevenly distributed on chromosomes in M. truncatula and M. sativa (Figure 4a,b). The NF-YB genes of M. truncatula were widely distributed on chr1 with seven members, followed by chr4 with five members (Figure 4a). Notably, chr6 had no NF-YB (Figure 4a). Nine NF-YB genes of M. sativa were also distributed on chr1, and one to six genes were distributed on chr2-6 ( Figure 4b).

Analyses of Chromosomal Distribution and Synteny of NF-YB Genes in Medicago
All NF-YB genes were unevenly distributed on chromosomes in M. truncatula and M. sativa (Figure 4a,b). The NF-YB genes of M. truncatula were widely distributed on chr1 with seven members, followed by chr4 with five members (Figure 4a). Notably, chr6 had no NF-YB (Figure 4a). Nine NF-YB genes of M. sativa were also distributed on chr1, and one to six genes were distributed on chr2-6 ( Figure 4b).  Tandem duplication, segmental duplication, and whole-genome duplication are the main impetus for gene family expansion [31]. Two pairs of segmental duplication were found in M. truncatula (MtNF-YB8/MtNF-YB15 and MtNF-YB9/MtNF-YB14) and M. sativa Over the course of evolutionary history, duplicated genes have three potential evolutionary fates: non-functionalization, neo-functionalization, and sub-functionalization [32]. In comparing the non-synonymous (Ka) and synonymous substitution (Ks) rates of substitution (Ka/Ks), one could infer the magnitude of selective constraint and positive selection. Generally, Ka/Ks > 1, Ka/Ks = 1, and Ka/Ks < 1 indicate positive selection, neutral evolution, and purifying selection, respectively [33]. In the present study, the Ka, Ks, and Ka/Ks of NF-YB homologous gene pairs were estimated in Medicago (Supplementary file S1). This showed that the MsNF-YB1/MtNF-YB2 homologous gene pairs had Ka/Ks ratios of 1.821, indicating a high degree of positive selection. In sharp contrast to this gene pair, we found that the Ka/Ks ratios of other NF-YB gene homologous pairs were less than 0.5, and that the ratios of another five homologous pairs were even smaller than 0.1, suggesting that NF-YB genes have undergone purifying selection after segmental and whole-genome duplications.
Overall, the promoters of NF-YB genes contained various cis-acting elements with different numbers. In particular, almost all NF-YB genes contain a high number of ARE elements (p < 0.05 for 11 of them with Fisher's exact test), and they may play a crucial role in anaerobic induction response in roots of Medicago. Methyl jasmonate (MeJA), as a wounding-related phytohormone, is able to stimulate the expression of defense-like genes [34]. Interestingly, NF-YB genes with relatively more MeJA-responsive elements were grouped in cluster B (Figure 5b,c), with 12 in MsNF-YB18 (p < 0.01 with Fisher's exact test), indicating that NF-YB genes in this cluster play a specific role in resistance against wounding stress. It is generally known that three cis-elements, ABRE, MBS, and W-box, are related with responsiveness to drought-induced signaling and regulation of downstream gene expression [35]. Our results showed that many NF-YB genes contained more ABRE and W-box elements (p < 0.05 for nearly half of them with Fisher's exact test), while MBS elements were grouped in C and D, indicating that the NF-YB family gene plays a role in drought resistance in Medicago. Notably, the NF-Y family protein binds specifically to the CCAAT cis-acting elements, though it is not widely distributed in the NF-YB family of Medicago (Figure 5b,c).

Expression Patterns of NF-YB Genes in Different Tissues
We investigated the expression patterns of NF-YBs in various tissues of M. truncatula with the genechip dataset from the MtGEA web server, including roots, stems, leaves, flowers, pods, and seeds ( Figure 6a). Overall, half (7) of the 14 MtNF-YB genes were expressed at a relatively low level in these tissues, and the other half were expressed at a relatively high level (Figure 6a). Among the lowly expressed seven genes, MtNF-YB16 was highly expressed in seeds, MtNF-YB7 in roots, and MtNF-YB17 in leaves (Figure 6a). Among the highly expressed seven genes, MtNF-YB5 and 11 were lowly expressed in seeds, and MtNF-YB15 and 19 in leaves (Figure 6a).

Expression Patterns of NF-YB Genes in Different Tissues
We investigated the expression patterns of NF-YBs in various tissues of M. truncatula with the genechip dataset from the MtGEA web server, including roots, stems, leaves, flowers, pods, and seeds ( Figure 6a). Overall, half (7) of the 14 MtNF-YB genes were expressed at a relatively low level in these tissues, and the other half were expressed at a relatively high level (Figure 6a). Among the lowly expressed seven genes, MtNF-YB16 was highly expressed in seeds, MtNF-YB7 in roots, and MtNF-YB17 in leaves (Figure 6a). Among the highly expressed seven genes, MtNF-YB5 and 11 were lowly expressed in seeds, and MtNF-YB15 and 19 in leaves (Figure 6a).
For M. sativa, gene expression levels in six tissues were analyzed based on transcriptome data, including roots, elongated stems, pre-elongated stems, leaves, flowers, and nodules (Figure 6b). The transcriptome data for MsNF-YBs could be divided into two distinct categories: higher expression genes (from MsNF-YB11 to MsNF-YB6, Figure 6b from top to bottom) and lower expression genes (from MsNF-YB1 to MsNF-YB26, Figure 6b from top to bottom). Compared with other genes, MsNF-YB11 and MsNF-YB22 were highly expressed in all tissues (Figure 6b). Interestingly, the expression of four genes (MsNF-YB10, 8,19,26) in nodules was significantly higher than in other tissues or than the remaining MsNF-YB genes, and these genes may play a specific role in nodules. For M. sativa, gene expression levels in six tissues were analyzed based on transcriptome data, including roots, elongated stems, pre-elongated stems, leaves, flowers, and nodules (Figure 6b). The transcriptome data for MsNF-YBs could be divided into two distinct categories: higher expression genes (from MsNF-YB11 to MsNF-YB6, Figure 6b from top to bottom) and lower expression genes (from MsNF-YB1 to MsNF-YB26, Figure  6b from top to bottom). Compared with other genes, MsNF-YB11 and MsNF-YB22 were highly expressed in all tissues (Figure 6b). Interestingly, the expression of four genes (MsNF-YB10, 8,19,26) in nodules was significantly higher than in other tissues or than the remaining MsNF-YB genes, and these genes may play a specific role in nodules.

Expression Levels of NF-YB Genes in Medicago under Stress
It was reported that NF-YB genes are involved in a number of important processes, such as embryo and seed development, flowering, photosynthesis, and stress responses [36]. Therefore, we also analyzed their expression level with genechip data for M. truncatula, including roots under in vitro culture salinity and hydroponic salinity, as well as roots and shoots under drought treatment. One probe set was selected as representative for each MtNF-YB gene, and 14 out of 21 MtNF-YB genes had their corresponding probe set (Figure 7a and Supplementary File S3).

Expression Levels of NF-YB Genes in Medicago under Stress
It was reported that NF-YB genes are involved in a number of important processes, such as embryo and seed development, flowering, photosynthesis, and stress responses [36]. Therefore, we also analyzed their expression level with genechip data for M. truncatula, including roots under in vitro culture salinity and hydroponic salinity, as well as roots and shoots under drought treatment. One probe set was selected as representative for each MtNF-YB gene, and 14 out of 21 MtNF-YB genes had their corresponding probe set (Figure 7a and Supplementary file S3).
Meanwhile, transcriptome data of M. sativa for samples from roots under drought and NaCl treatment were also analyzed (Supplementary file S4). Due to extremely low expression of some genes, only 19 MsNF-YB genes were presented in Figure 7f In M. truncatula, 11 genes (MtNF-YB13, MtNF-YB6, MtNF-YB19, etc.) and 7 genes (MtNF-YB10, MtNF-YB15, MtNF-YB7, etc.) were highly induced in roots or in shoots under drought treatment, respectively (Figure 7a,b). In contrast, all genes were significantly up-regulated at different levels under NaCl treatment, except MtNF-YB5 and MtNF-YB12 genes, which were down-regulated under hydroponic NaCl treatment (Figure 7c,d). Taken together, as shown in the Venn diagram, ten genes were highly expressed under at least three stress treatments (Figure 7e  (i) Venn diagram of genes whose expression levels were increased by stress treatment; these genes were selected based on the homolog listed in Table 1.
Meanwhile, transcriptome data of M. sativa for samples from roots under drought and NaCl treatment were also analyzed (Additional File 4). Due to extremely low expression of some genes, only 19 MsNF-YB genes were presented in Figure 7f In M. truncatula, 11 genes (MtNF-YB13, MtNF-YB6, MtNF-YB19, etc.) and 7 genes (MtNF-YB10, MtNF-YB15, MtNF-YB7, etc.) were highly induced in roots or in shoots under drought treatment, respectively (Figure 7a,b). In contrast, all genes were significantly up-regulated at different levels under NaCl treatment, except MtNF-YB5 and MtNF-YB12 genes, which were down-regulated under hydroponic NaCl treatment (Figure 7c,d). Taken together, as shown in the Venn diagram, ten genes were highly expressed under at least three stress treatments (Figure 7e  (i) Venn diagram of genes whose expression levels were increased by stress treatment; these genes were selected based on the homolog listed in Table 1.

Validation of the Expression Profile of Stress-Responsive NF-YB Genes by RT-qPCR
In order to verify the expression profiles of NF-YB genes from genechip for M. truncatula and transcriptome data for M. sativa, RT-qPCR was carried out for the verification of six candidate gene pairs. The expression levels of these above-mentioned six genes in four tissues (roots, stems, leaves, and flowers) were further verified by RT-qPCR, with their lowest expression as a control, as shown in the model figure (Figure 8a and Supplementary file S5). Notably, the expression profile of almost all of these genes was consistent with genechip data (Figures 6a and 8a); for example, MtNF-YB7, MtNF-YB10, and MtNF-YB17 showed the highest expression in roots, leaves, and leaves, respectively. As an exception, MtNF-YB16 was expressed with a significantly higher level in flowers than in other tissues, whereas it was expressed at a similar level to the genechip data in the four tissues (Figures 6a and 8a).
Moreover, seedlings were treated with NaCl or mannitol at 1 h, 3 h, 6 h, 12 h, 24 h, 48 h and used for RT-qPCR analysis. RT-qPCR showed that MtNF-YB7 and MtNF-YB15 were highly induced by both NaCl and mannitol treatment from 3 h to 48 h (Figure 8b and Supplementary file S5). The expression level of MtNF-YB10 was induced by mannitol at 6 h and 12 h, and MtNF-YB13 at 1 h, 3 h, 6 h, 12 h, and 24 h, though the fold change was not more than 2.5-fold (Figure 8b). The expression level of MtNF-YB16 was induced by more than 15-fold and 50-fold at 1 h and 48 h by NaCl treatment (Figure 8b). However, the expression of MtNF-YB17 appeared to be repressed by both treatments, except for a slight increase by the mannitol treatment at 3 h (Figure 8b).
For M. sativa, the expression levels detected by RT-qPCR for MsNF-YB9, 11, and 28 were relatively consistent with the transcriptome data (Figures 6a and 9b): the expression level of MsNF-YB9 was significantly higher in roots than in the other three tissues; MsNF-YB11 was expressed at a relatively high level in all four tissues (Figure 9a and Supplementary file S5). Unlike transcriptome data, RT-qPCR data showed that MsNF-YB18 had a significantly higher expression level in leaves than in other tissues. Meanwhile, RT-qPCR data showed that MsNF-YB27 was significantly expressed in roots, while transcriptome data showed its expression was the highest in pre-elongated stems (Figure 9a). Since the treatment timepoints for RT-qPCRs were the same as for the transcriptome data, we compared the transcriptome data and RT-qPCR data with the correlation analysis. RT-qPCR data showed that MsNF-YB9, 27, and 28 showed a significant increase in gene expression under both NaCl and mannitol treatment (Figure 9b and Supplementary file S5), which was the same as for transcriptome data, and they were positively correlated (Figure 9c). Notably, MsNF-YB9/MtNF-YB7, MsNF-YB27/MtNF-YB15, and MsNF-YB28/MtNF-YB16 were homologous gene pairs, and their expression was consistent under the two stress treatments, indicating that these genes play the same role in M. truncatula and M. sativa. Moreover, the expression level of MsNF-YB11 was transiently elevated under NaCl treatment and then reduced greatly, which was inconsistent with its significant elevation as shown by the transcriptome data. The expression of MsNF-YB18 detected by RT-qPCR was consistent with the changes of the transcriptome data (Figure 9c).

Discussion
Previous studies have demonstrated that NF-YB proteins play an important role in plant resistance to abiotic stresses [13]. The function of NF-YB in abiotic resistance was previously analyzed in several plant species, including A. thaliana, rice, wheat, tung tree, soybean, canola, grape, and tomato [9]. However, a genome-wide identification and characterization of NF-YB genes in Medicago is still lacking. In the present study, we conducted an integrated investigation on the NF-YBs, and a total of 21 NF-YB members from M. truncatula and 28 NF-YB members from M. sativa were identified.
Multiple sequence alignment confirmed that almost all M. truncatula and M. sativa NF-YB members contained DNA binding as well as NF-YA interaction and NF-YC interaction domains (Figure 1), and they shared highly conserved amino acid length, which was consistent with model plants such as rice [37]. Phylogenetic analysis indicated that they are highly homologous to NF-YB proteins from Arabidopsis (Figure 2), but one particular sub-cluster E was only presented in Medicago, which may play a special role in evolution. The same sub-clusters are more closely related, and the type and position of their patterns are similar (Figures 2 and 3). These findings suggest that these motifs may be involved in the functional diversity of NF-YB proteins. In analyzing gene structure, we found that many NF-YB genes in Medicago had only one exon and no introns (Figure 3), which is consistent with the findings for Arabidopsis and Brassica napus. [30]. Previous studies have postulated that an intron-rich gene would lose multiple introns simultaneously by retrotransposition, thereby producing intron-less ancestral genes [38]. Thus, several NF-YB genes in Medicago may experience the loss of multiple introns during gene

Discussion
Previous studies have demonstrated that NF-YB proteins play an important role in plant resistance to abiotic stresses [13]. The function of NF-YB in abiotic resistance was previously analyzed in several plant species, including A. thaliana, rice, wheat, tung tree, soybean, canola, grape, and tomato [9]. However, a genome-wide identification and characterization of NF-YB genes in Medicago is still lacking. In the present study, we conducted an integrated investigation on the NF-YBs, and a total of 21 NF-YB members from M. truncatula and 28 NF-YB members from M. sativa were identified.
Multiple sequence alignment confirmed that almost all M. truncatula and M. sativa NF-YB members contained DNA binding as well as NF-YA interaction and NF-YC interaction domains (Figure 1), and they shared highly conserved amino acid length, which was consistent with model plants such as rice [37]. Phylogenetic analysis indicated that they are highly homologous to NF-YB proteins from Arabidopsis (Figure 2), but one particular sub-cluster E was only presented in Medicago, which may play a special role in evolution. The same sub-clusters are more closely related, and the type and position of their patterns are similar (Figures 2 and 3). These findings suggest that these motifs may be involved in the functional diversity of NF-YB proteins. In analyzing gene structure, we found that many NF-YB genes in Medicago had only one exon and no introns (Figure 3), which is consistent with the findings for Arabidopsis and Brassica napus. [30]. Previous studies have postulated that an intron-rich gene would lose multiple introns simultaneously by retrotransposition, thereby producing intron-less ancestral genes [38]. Thus, several NF-YB genes in Medicago may experience the loss of multiple introns during gene family diversification. Genomewide analyses have shown that the loss and gain of introns were extensive during the process of eukaryotic diversification.
Duplication and divergence play a critical role in the expansion and evolution of gene families [39]. We found two segmental duplications (MtNF-YB8/MtNF-YB15, MtNF-YB9/MtNF-YB14) among 21 MtNF-YB genes (Figure 4), while two segmental duplication events (MsNF-YB4/MsNF-YB11/MsNF-YB15) among 28 MsNF-YBs were found in M. sativa. The gain and loss of genes or the expansion or contraction of gene families is common following polyploidization [40]. Unlike M. truncatula, M. sativa has tandem duplicated gene pairs (MsNF-YB1/MsNF-YB2). Thus, the expansion of the MsNF-YB gene family could be an indication that MsNF-YB genes play roles in additional biological processes or have novel functions. Previous studies have demonstrated that NF-YB plays an important role in drought resistance [13]. In our studies, ARE (anaerobic induction), TGACG-motif (MeJAresponsive), and ABRE (abscisic acid-responsive) elements were found widely distributed in NF-YB genes, and all of these elements were involved in response to drought stress, suggesting that the NF-YBs of Medicago are increased in response to drought stress.
High-salinity or drought soil is the most serious abiotic stress [41]. It is urgent to improve the salinity and drought tolerance of alfalfa to increase yield. Genechip data and transcriptome data under NaCl and drought treatments, along with expression profiles in various tissues, suggested that the expressions of six homologous gene pairs were highly induced or drastically changed in Medicago (Figures 6 and 7). Previous studies have reported that NF-YB genes are also involved in plant developmental processes, including embryogenesis, flowering time, etc. [42]. In the present study, we identified the tissuespecific expression patterns of six NF-YB genes in Medicago in a variety of tissues, and the results showed that NF-YB genes were expressed ubiquitously, with the exception of a few genes that are expressed in specific tissues: MtNF-YB16 is specifically expressed in flowers (Figure 8a), and MsNF-YB9 in roots (Figure 9a). This observation was consistent with previous studies for NF-YB genes in rice [43], suggesting that NF-YB genes are multifunctional and are involved in a wide range of biological processes [44]. Subsequently, the expression pattern of all six genes was verified under NaCl and mannitol treatments in M. truncatula and M. sativa by RT-qPCR analyses (Figures 8b and 9b), and the results suggested that several genes are involved in NaCl and/or mannitol induction. Correspondingly, three gene pairs (MsNF-YB9/MtNF-YB7, MsNF-YB27/MtNF-YB15, and MsNF-YB28/MtNF-YB16) were significantly up-regulated under NaCl and mannitol treatments, and their expression patterns were the same in M. truncatula and M. sativa. This evidence indicated that these genes are likely key genes in response to abiotic stress.
In the phylogenetic analysis, NF-YB genes were divided into five branches in addition to the E branch, which includes Arabidopsis-specific NF-YB genes. Among them, NF-YB1, NF-YB2, NF-YB3, NF-YB6, and NF-YB9 have been extensively studied in Arabidopsis. Previous studies have shown that NF-YB1 not only regulates drought resistance but also interacts with CO (CONSTANS) to affect the transcript levels of two key integrators (FT and SOC1) in the flowering pathway, therefore adjusting flowering time [45]. Interestingly, MtNF-YB18 clustered with AtNF-YB1, indicating that it may have similar functions as AtNF-YB1. Moreover, MsNF-YB2, 27 and MtNF-YB11, 15 were observed to cluster with AtNF-YB2 (Figure 2), which has been reported to regulate the photoperiod-dependent flowering time [46]. However, in our study, the homologous gene pair MtGRF15/MsGRF27 showed significant high expression under both salt and mannitol treatments (Figures 8b and 9b), suggesting that this gene pair may also play a role in coping with abiotic stresses.
In Arabidopsis, LEAFY COTYLEDON1 (LEC1) is a central regulator that controls many different aspects of embryo development. Early in embryogenesis, LEC1 is required to maintain the fate of embryonic cells that constitute the suspensor and to specify the identity of cotyledons, the embryonic leaves [47]. While previous studies have shown that NF-YB proteins of Arabidopsis can be divided into two classes, LEC1-like (LEC1 or AtNF-YB9; LEC1-LIKE or AtNF-YB6) and non-LEC1-like. NF-YB9/LEC1 was the first NF-YB gene identified and studied in A. thaliana and has been shown to be required for the embryonic maintenance of cell fate, where the ectopic expression of NF-YB9 can induce somatic embryos from vegetative cells [16]. In addition, NF-YB9 has also been shown to play an essential role in embryogenesis and seed maturation [48]. LEC1 and LEC1-LIKE (NF-YB6) regulated embryo development by activating the expression of genes required for embryogenesis and cellular differentiation [16]. In the present study, MsNF-YB28 was grouped with AtNF-YB9, while MtNF-YB16 was grouped with AtNF-YB6 (Figure 2), suggesting they may share similar functions.
Furthermore, the homologous gene pairs (MtNF-YB16/MsNF-YB28) were all highly expressed under NaCl and mannitol treatment as evidenced by RT-qPCR (Figures 8b and 9b). Thus, the homologous pairs may be involved in regulating embryonic development, as well as in response to abiotic stress. Therefore, these two candidate homologous genes (MtGRF15/MsGRF27 and MtNF-YB16/MsNF-YB28) were considered as key candidates for in-depth study of NF-YB genes in Medicago.

Analyses of Sequence, Conserved Motif, and Structural Characterization
To exhibit the structural divergence of NF-YB proteins, the conserved motifs were performed with Multiple Em for Motif Elicitation (MEME) 5.0.2 online program (https: //meme-suite.org/meme/tools/meme) (accessed on 25 September 2020) [49]. The following parameters were employed: the maximum number of motifs was 20, minimum motif width was 10 (aa), and maximum motif width was 200 (aa). Subsequently, sequence alignment analysis of NF-YB protein sequences was carried out using jalview (https://issues.jalview.org/secure/Dashboard.jspa) (accessed on 25 September 2020). The visualization of exon-intron positions and conserved motifs was performed through TBtools software (Guangzhou, China) [50].

Analyses of Phylogenetic Relationship
NF-YB proteins from five plant species (M. truncatula, M. sativa, A. thaliana, and O. sativa) were used in a multiple alignment in ClustalW [51]. Phylogenetic trees were constructed by the neighbor-joining method using the program MEGA-X with a bootstrap of 1000 replicates. Meanwhile, subfamily clustering on the phylogenetic tree was determined based on Arabidopsis member clustering. Subsequently, EvolView (https://evolgenius.info/ evolview-v2/) (accessed on 29 September 2020) was used to view the phylogenetic tree.

Analysis of Chromosome Locations and Collinearity
The loci of NF-YB genes were obtained from the genome annotation data. TBtools was applied to map the chromosome locations for each gene. Next, these sequences were analyzed to identify collinearity blocks against the whole genome using MCSCAN [52]. Moreover, the intraspecific synteny relationship (M. truncatula and M. sativa) and interspecific synteny relationship (M. truncatula, M. sativa, A. thaliana, and O. sativa) were analyzed, and they were mapped to the chromosomes of M. truncatula and M. sativa using TBtools, respectively. Lastly, the synonymous (Ks) and nonsynonymous (Ka) substitution rates were estimated using TBtools [53].

Identification of Stress-Response Cis-Elements of the NF-YB Promoter
The promoter sequences (length, 2 kb) of NF-YBs were collected by TBtools, and the cis-elements in the promoters were predicted using PlantCARE (http://bioinformatics.psb. ugent.be/webtools/plantcare/html/) (accessed on 2 October 2020). The visualized models of cis-elements in the promoters were made with TBtools [54].

Analysis of Gene Expression Profile
Genechip data from roots and shoots and those under drought and salt stress conditions for MtNF-YB genes were downloaded from the M. truncatula Gene Expression Atlas. The expression levels of MtNF-YB genes in different tissues were also analyzed. Amazing HeatMap software was used to generate the heatmap [50]. The original transcriptome data from M. sativa under NaCl and mannitol treatments at 0 h, 1 h, 3 h, 6 h, 12 h, and 24 h (SRR7160314-15, 22-23, 25-49, 51-52, 56-57) were downloaded from the NCBI database [55]. Then, the data were converted to fastq files using the SRA-Toolkit v2.9 [56]. Raw reads were trimmed using Trimmomatic-0.39 [57]. Gene expression level was determined by mapping cleaned reads to the corresponding alfalfa reference genomes using the StringTie v2.1.3 package and the log2(FPKM) values [58].

Plant Materials and Treatments
M. truncatula (cv. Jemalong A17) and M. sativa (cv. Zhongmu No.1) plants used in this study were stored at the Institute of Animal Sciences of the Chinese Academy of Agricultural Sciences, Beijng, China. Roots, stems, leaves, and flowers of mature M. truncatula and M. sativa plants were collected separately for RNA extraction and RT-qPCR analysis. To investigate the expression pattern of NF-YB genes in response to NaCl and mannitol stress, seeds were germinated and transferred into the MS liquid medium (MS basal salts supplemented with 30 g/L sucrose), then kept in a growth chamber at 25 • C under a photoperiod using a 16/8 light/dark regime. When the third leaf was fully expanded, seedlings were transferred to fresh MS liquid medium supplied with 300 mM NaCl and 15% mannitol, respectively, and the whole plants were collected at 0 h, 1 h, 3 h, 6 h, 12 h, 24 h, and 48 h for each treatment, and plants without treatment were used as controls. The samples were frozen in liquid nitrogen and stored at −80 • C for subsequent analysis.

Gene Expression Analyses by RT-qPCR
Total RNAs were extracted using Eastep ® Super total RNA Extraction kit (Promega, Shanghai, China) according to the manufacturer's instructions. First-strand cDNA synthesis was performed using Trans ® Script One-Step gDNA Removal and cDNA Synthesis SuperMix (TransGen Biotech, Beijing, China) per the manufacturer's recommendations. RT-qPCRs were performed with a 2'RealStar Green Fast Mixture (GeneStar, Shanghai, China) on an ABI Q7 QuantStudioTM Real-Time PCR Detection System (Applied Biosystems, CA, USA). PCRs were performed with the following program: 94 • C for 30 s, followed by 40 cycles of 94 • C for 5 s and 60 • C for 34 s. Melting curve analysis was performed using ABI Q7 QuantStudio TM Real-Time PCR Software with the RT-qPCR data. The transcript levels of each gene were determined by relative quantification using the 2 -∆∆Ct method and normalized with the actin gene as a reference. Data are the average of three independent biological samples ± SE, and vertical bars indicate standard deviation. Student's t test was used to compare the difference between two treatments (n = 3, * p < 0.05, ** p < 0.01) [59]. The primer sequences used in this study are shown in Table S1.