Genome-Wide Identiﬁcation and Analysis of the NF-Y Transcription Factor Family Reveal Its Potential Roles in Salt Stress in Alfalfa ( Medicago sativa L.)

: Nuclear factor Y (NF-Y) is a heterotrimeric transcription factor that plays an important role in various biological processes in plants, such as ﬂowering regulation, drought resistance, and salt stress. However, few in-depth studies investigated the alfalfa NF-Y gene family. In this study, in total, 60 MsNF-Y genes, including 9 MsNF-YAs , 26 MsNF-YBs , and 25 MsNF-YCs, were identiﬁed in the alfalfa genome. The genomic locations, gene structures, protein molecular weights, conserved domains, phylogenetic relationships, and gene expression patterns in different tissues and under different stresses (cold stress, drought stress, and salt stress) of these NF-Y genes were analyzed. The illustration of the conserved domains and speciﬁc domains of the different subfamilies of the MsNF-Y genes implicates the conservation and diversity of their functions in alfalfa growth, development, and stress resistance. The gene expression analysis showed that 48 MsNF-Y genes (7 MsNF-YAs , 22 MsNF-YBs , and 19 MsNF-YCs ) were expressed in all tissues at different expression levels, indicating that these genes have tissue expression speciﬁcity and different biological functions. In total, seven, seven, six, and eight MsNF-Y genes responded to cold stress, the ABA treatment, drought stress, and salt stress in alfalfa, respectively. According to the WGCNA, molecular regulatory networks related to salt stress were constructed for MsNF-YB2 , MsNF-YB5 , MsNF-YB7 , MsNF-YB15 , MsNF-YC5 , and MsNF-YC6 . This study could provide valuable information for further elucidating the biological functions of MsNF-Ys and improving salt tolerance and other abiotic stress resistance in alfalfa.


Introduction
Alfalfa (Medicago sativa L.) is one of the top-quality forages for livestock, such as dairy cattle, because of its high nutrient content and crude protein content and has been hailed as the 'king of forages' [1,2]. However, alfalfa is frequently affected by various abiotic stresses, such as drought, cold, and salt, during alfalfa growth and development, which seriously affect the yield and quality of alfalfa [3]. Alfalfa has formed complex gene regulatory networks in response to increasingly deteriorating growth environments over a long period of evolution [3]. In addition, with the reporting of the alfalfa reference genome, several transcription factor families, such as WRKY, CBF, SPL, and MYB, have been identified to be involved in abiotic stress processes in alfalfa [3]. However, our knowledge of the complex transcriptional regulatory network is limited.
The nuclear factor Y (NF-Y) transcription factor, also known as the CCAAT-binding factor (CBF) or the heme activator protein (HAP), is widespread in eukaryotes. NF-Y consists of three distinct subunits, NF-YA, NF-YB, and NF-YC [4]. The subfamily (NF-YA, NF-YB, and NF-YC) members can be identified by the sequence length and conserved motifs or domains. In general, NF-YA subunits have a longer sequence than the NF-YB and NF-YC sequences [5]. The conserved region in the NF-YA subunits usually consists The protein lengths of the 60 MsNF-Y genes showed a wide distribution, ranging from 136 AA (MsNF-YB5) to 400 AA (MsNF-YC14) ( Table 1 and Table S1). The protein length of the MsNF-YA members ranged from 312 AA to 334 AA, while the protein length of the MsNF-YB members ranged from 136 AA to 237 AA, and the MsNF-YC members were distributed from 160 AA to 400 AA. Among the three MsNF-Y subfamilies, the MsNF-YA subfamily member proteins had the longest average length of 328 AA, whereas MsNF-YB had the shortest protein with an average length of 182. 3

Analysis of the Motifs and Conserved Domains of MsNF-Y Family Members
To analyze the motif composition of these MsNF-Y proteins, MEME software was used to identify the motifs. The results showed that different subfamilies had different motif compositions (Figure 1). In the MsNF-YA group, motif 1 and motif 2 were found in all members ( Figure 1A). MsNF-YA8 and MsNF-YA9 did not have motif 3. In the MsNF-YB subfamily, all members, except for MsNF-YB5, contained the three motifs identified ( Figure 1B). In the MsNF-YC subfamily, the composition of the motifs in these proteins was relatively diverse. There were only 10 MsNF-YC proteins with one motif and five MsNF-YC proteins with two motifs ( Figure 1C).  To analyze the conserved domains in the different NF-Y subfamilies, multiple alignments were conducted by using DNAMAN software. The results showed that there were a few conserved domains in the MsNF-Y proteins, indicating that these proteins belong to the same subgroups ( Figure 2). For example, the MsNF-YA proteins have conserved domains, including one domain for the protein-protein interaction of NF-YB/C and another for the DNA-binding domain (Figure 2A). The MsNF-YB proteins have one domain for DNA binding and another domain for interacting with NF-YA and NF-YC proteins ( Figure 2B). Among the MsNF-YC proteins, the results are the same (one domain for DNA binding and another domain for the interaction of NF-YA/B) but with lower conservation compared with that in the MsNF-YA and MsNF-YB subfamilies, implying a probably more abundant functional role in alfalfa ( Figure 2C). To analyze the conserved domains in the different NF-Y subfamilies, multiple alignments were conducted by using DNAMAN software. The results showed that there were a few conserved domains in the MsNF-Y proteins, indicating that these proteins belong to the same subgroups ( Figure 2). For example, the MsNF-YA proteins have conserved domains, including one domain for the protein-protein interaction of NF-YB/C and another for the DNA-binding domain (Figure 2A). The MsNF-YB proteins have one domain for DNA binding and another domain for interacting with NF-YA and NF-YC proteins ( Figure 2B). Among the MsNF-YC proteins, the results are the same (one domain for DNA binding and another domain for the interaction of NF-YA/B) but with lower conservation compared with that in the MsNF-YA and MsNF-YB subfamilies, implying a probably more abundant functional role in alfalfa ( Figure 2C).

Phylogenetic Relationships and Gene Structure of the MsNF-Y Family Members
To understand the phylogenetic relationships among the MsNF-Ys, phylogenetic trees were constructed by using MEGA 5.0 software based on the protein sequences. The

Genome Distribution and Gene Duplication Events of the MsNF-Ys in the Alfalfa Genome
To show the genome distribution of the 60 MsNF-Y genes in the alfalfa genome, the TBtools program was used. The results showed that the 60 MsNF-Y genes were distributed unevenly on the 32 alfalfa chromosomes, except for chr4.1, chr4.2, chr4.3, chr4.4, chr5.4, chr6.1, chr6.2, chr6.3, chr6.4, and chr7.4 ( Figure 4A). Six NF-Y genes were on chr7.3, while on the other chromosomes there were no more than five NF-Y genes. tained eight, seven, six, and four exons, respectively ( Figure 3B). Based on the phylogenetic tree and gene structure, we found that the same tree groups had the same gene structure ( Figure 3).

Genome Distribution and Gene Duplication Events of the MsNF-Ys in the Alfalfa Genome
To show the genome distribution of the 60 MsNF-Y genes in the alfalfa genome, the TBtools program was used. The results showed that the 60 MsNF-Y genes were distributed unevenly on the 32 alfalfa chromosomes, except for chr4.1, chr4.2, chr4.3, chr4.4, chr5.4, chr6.1, chr6.2, chr6.3, chr6.4, and chr7.4 ( Figure 4A). Six NF-Y genes were on chr7.3, while on the other chromosomes there were no more than five NF-Y genes. In total, 60 duplication events were found among these MsNF-Y genes in the alfalfa genome ( Figure 4). Two duplication events (MsNF-YC9 and MsNF-YC10 on chr3.1 and MsNF-YC19 and MsNF-YC20 on chr7.3) were formed in the NF-YC genes on the same chromosome, and 58 duplication events were formed by different NF-Y genes on different chromosomes ( Figure 4B and Table S2). Among the 58 duplication events, 9 duplication In total, 60 duplication events were found among these MsNF-Y genes in the alfalfa genome ( Figure 4). Two duplication events (MsNF-YC9 and MsNF-YC10 on chr3.1 and MsNF-YC19 and MsNF-YC20 on chr7.3) were formed in the NF-YC genes on the same chromosome, and 58 duplication events were formed by different NF-Y genes on different chromosomes ( Figure 4B and Table S2). Among the 58 duplication events, 9 duplication events were formed by the 9 MsNF-YA genes, 25 duplication events were formed by the 24 MsNF-YB genes, and 24 duplication events were formed by the 21 MsNF-YC genes. To understand the evolution type among these duplication events on different chromosomes, the Ka, Ks, and Ka/Ks values were calculated in this study. The results showed that among the 58 duplication events, 41 duplication events were estimated as purifying selection, 11 events were neutral evolution, and 6 events were positive selection (Table S2).

cis-Element Analysis in the Promoter Regions of MsNF-Y Genes
To clarify the transcriptional regulation of the MsNF-Y genes, the cis elements in the promoter regions were identified by using the PlantCare program. In total, 1679 cis elements belonging to 21 different categories were found in the promoter regions of the 60 MsNF-Y genes ( Figure 5A). Light-responsive elements existed in the promoter regions of all MsNF-Ys, and there were 94, 390, and 291 light-responsive elements in MsNF-YA, MsNF-YB, and MsNF-YC, respectively ( Figure 5B). In addition to light-responsive elements, MeJAresponsive elements are the second largest category in the promoter regions of MsNF-Ys. In total, 38 MsNF-Y genes (5 MsNF-YAs, 22 MsNF-YBs, and 11 MsNF-YCs) contained MeJA-responsive elements in the promoter regions ( Figure 5B). In total, 160 abscisic acid response elements (ABREs) were found in the promoter regions of 51 MsNF-Y genes (2 MsNF-YAs, 25 MsNF-YBs, and 24 MsNF-YCs). There were 55 auxin-responsive elements and 39 gibberellin-responsive elements in the promoter regions of 29 and 27 MsNF-Y genes, respectively ( Figure 5A). The promoter regions of 24, 11 and 21 MsNF-Y genes harbored MYB transcription-factor-binding sites related to drought inducibility, flavonoid biosynthetic gene regulation, and light-responsive elements, respectively. Many stress response elements were also found in the promoter regions of the MsNF-Ys ( Figure 5B). The promoter regions of 19, 22, and 10 MsNF-Y genes contained 29 low temperature response elements (LTRs), 35 defense and stress responsiveness elements (TC-rich repeats), and 11 dehydration, low temperature, and salt stress elements (DREs), respectively. These results indicate that MsNF-Ys not only participate in the growth and development of alfalfa but also play a key role in various stresses.

Expression Analysis of MsNF-Y Genes in Different Tissues from Alfalfa
To investigate the expression pattern of these MsNF-Y genes, the global gene expression in different tissues (root, flower, nodule, leaf, ES1, and PES) of alfalfa was collected and analyzed. The results showed that the 60 MsNF-Y genes had significantly different expression patterns and indicated that these genes played different roles in alfalfa development ( Figure 7). Among the nine MsNF-YA genes, two MsNF-YAs (MsNF-YA5 and MsNF-YA7) were not expressed in the six tissues in this study. MsNF-YA1 and MsNF-YA3 exhibited high expression levels in the PES, and MsNF-YA4 had a higher expression level in the leaf. MsNF-YA2/9/6/8 had a similar expression pattern and were mainly expressed in the roots and nodules ( Figure 7A). Among the twenty-six MsNF-YB genes, four genes (MsNF-YB19/22/24/25) were not expressed in these tissues ( Figure 7B). MsNF-YB4 was only expressed in the roots, and MsNF-YB6/8/26 were only expressed in the nodules. Based on the expression patterns, these MsNF-YB genes were clustered into different groups. MsNF-YB1/MsNF-YB13 were mainly expressed in ES1 and PES, and MsNF-YB14/15 were mainly expressed in ES1, PES, and leaves. MsNF-YB9/12/10 were mainly expressed in the leaf, while MsNF-YB17/3/18/16 were mainly expressed in the flower, and MsNF-YB6/8/26 were mainly expressed in the nodule ( Figure 7B). Among the twenty-five MsNF-YC genes, six genes (MsNF-YC3, MsNF-YC9, MsNF-YC17, MsNF-YC19, MsNF-YC20, MsNF-YC22, and MsNF-YC24) were not expressed in these tissues ( Figure 7C). MsNF-YC7 was only expressed in the root, while MsNF-YC10 is only expressed in the leaf. The other MsNF-YC genes were expressed in at least two different tissues ( Figure 7C). These results indicate that these MsNF-Y genes perform different functions in alfalfa growth and development.

Expression Analysis of MsNF-Y Genes in Alfalfa under Different Abiotic Stresses
To clarify the role of the MsNF-Y genes under different abiotic stresses, we collected transcriptome data under cold, ABA, drought, and salt stress. The results showed that multiple MsNF-Y genes are involved in different abiotic stress processes (Figure 8). The expression levels of seven MsNF-Y genes (five MsNF-YB genes and two MsNF-YC genes) significantly changed under cold stress ( Figure 8A). Among the seven genes, the expression levels of six genes decreased as the cold stress treatment time increased, and only one gene (MsNF-YC1) showed an increase as the cold stress treatment time increased ( Figure 8A). Under the ABA treatment, the expression levels of seven genes (one MsNF-YA gene, two MsNF-YB genes, and four MsNF-YC genes) were significantly upregulated ( Figure 8B).
The expression levels of six MsNF-Y genes (three MsNF-YB genes and three MsNF-YC genes) were significantly upregulated under drought stress ( Figure 8C). MsNF-YB21 and MsNF-YB14 expression abundance was significantly upregulated in the early stage of drought stress. The gene expression levels of the other four genes (MsNF-YC5, MsNF-YB7, MsNF-YC14, and MsNF-YC7) were significantly upregulated in the later stage of drought stress. Under salt stress, the expression levels of eight MsNF-Y genes (five MsNF-YB genes and three MsNF-YC genes) were significantly upregulated ( Figure 8D). MsNF-YB2 gene expression was significantly upregulated in the early stage of salt stress. MsNF-YC5 and MsNF-YB21 gene expression first increased and then decreased as the salt treatment time increased. The expression levels of the other five genes were significantly upregulated in the later stage of salt stress ( Figure 8D).

Expression Analysis of MsNF-Y Genes in Different Tissues from Alfalfa
To investigate the expression pattern of these MsNF-Y genes, the global gene expression in different tissues (root, flower, nodule, leaf, ES1, and PES) of alfalfa was collected and analyzed. The results showed that the 60 MsNF-Y genes had significantly different expression patterns and indicated that these genes played different roles in alfalfa development (Figure 7). Among the nine MsNF-YA genes, two MsNF-YAs (MsNF-YA5 and MsNF-YA7) were not expressed in the six tissues in this study. MsNF-YA1 and MsNF-YA3 MsNF-YB6/8/26 were mainly expressed in the nodule ( Figure 7B). Among the twenty-five MsNF-YC genes, six genes (MsNF-YC3, MsNF-YC9, MsNF-YC17, MsNF-YC19, MsNF-YC20, MsNF-YC22, and MsNF-YC24) were not expressed in these tissues ( Figure 7C). MsNF-YC7 was only expressed in the root, while MsNF-YC10 is only expressed in the leaf. The other MsNF-YC genes were expressed in at least two different tissues ( Figure 7C). These results indicate that these MsNF-Y genes perform different functions in alfalfa growth and development.

Expression Analysis of MsNF-Y Genes in Alfalfa under Different Abiotic Stresses
To clarify the role of the MsNF-Y genes under different abiotic stresses, we collected transcriptome data under cold, ABA, drought, and salt stress. The results showed that multiple MsNF-Y genes are involved in different abiotic stress processes (Figure 8). The expression levels of seven MsNF-Y genes (five MsNF-YB genes and two MsNF-YC genes) significantly changed under cold stress ( Figure 8A). Among the seven genes, the expression levels of six genes decreased as the cold stress treatment time increased, and only one gene (MsNF-YC1) showed an increase as the cold stress treatment time increased ( Figure  8A). Under the ABA treatment, the expression levels of seven genes (one MsNF-YA gene, two MsNF-YB genes, and four MsNF-YC genes) were significantly upregulated ( Figure  8B). The expression levels of six MsNF-Y genes (three MsNF-YB genes and three MsNF-YC genes) were significantly upregulated under drought stress ( Figure 8C). MsNF-YB21 and MsNF-YB14 expression abundance was significantly upregulated in the early stage of drought stress. The gene expression levels of the other four genes (MsNF-YC5, MsNF-YB7, MsNF-YC14, and MsNF-YC7) were significantly upregulated in the later stage of drought  Combining these different abiotic stress results, we also found multiple MsNF-Y genes involved in two or more stresses. MsNF-YB14 responds to both cold and drought stresses, MsNF-YB15 responds to both cold and salt stresses, and MsNF-YB7, MsNF-YC5, and MsNF-YC14 respond to both drought and salt stresses. Combining these different abiotic stress results, we also found multiple MsNF-Y genes involved in two or more stresses. MsNF-YB14 responds to both cold and drought stresses, MsNF-YB15 responds to both cold and salt stresses, and MsNF-YB7, MsNF-YC5, and MsNF-YC14 respond to both drought and salt stresses.

Important Candidate MsNF-Y Genes in Salt Stress
To further identify candidate genes among the eight differentially expressed MsNF-Y genes in salt stress, a weighted gene coexpression network analysis (WGCNA) was conducted in this study. The results showed that these genes were divided into 14 coexpression modules ( Figure 9A). The correlation analysis between the modules and the salt stress treatment time revealed that a few modules were significantly correlated with salt stress ( Figure 9B). The magenta and black modules were significantly correlated with the early stage (Tables S1 and S2) of salt stress (R = 0.66, p = 1.2 × 10 −3 ; R = 0.58, p = 5.7 × 10 −3 ). The pink, purple, and tan modules were significantly correlated with the middle stage (Tables S3 and S4)  Among these significant modules, six MsNF-Y genes were included. MsNF-YB2, which responds to salt stress at the early stage, was found in the magenta module; MsNF-YC5, which responds to salt stress at the middle stage, was found in the purple module; and the other four MsNF-Y genes (MsNF-YB5, MsNF-YB7, MsNF-YB15, and MsNF-YC6) were located in the turquoise module.
To further understand the coexpression network of the six MsNF-Y genes, the genes that were synergistically expressed under salt stress were identified. The results showed that MsNF-YB2 was mainly coexpressed with the DEAD genes under salt stress at the Among these significant modules, six MsNF-Y genes were included. MsNF-YB2, which responds to salt stress at the early stage, was found in the magenta module; MsNF-YC5, which responds to salt stress at the middle stage, was found in the purple module; and the other four MsNF-Y genes (MsNF-YB5, MsNF-YB7, MsNF-YB15, and MsNF-YC6) were located in the turquoise module.
To further understand the coexpression network of the six MsNF-Y genes, the genes that were synergistically expressed under salt stress were identified. The results showed that MsNF-YB2 was mainly coexpressed with the DEAD genes under salt stress at the early stage ( Figure 10A and Table S3). Regarding MsNF-YB5, there were 43 coexpressed genes, including 33 protein kinase genes, 6 MYB transcription factors and 4 DEAD genes ( Figure 10B and Table S4). In total, 22 coexpressed genes containing 12 AP2 domain genes, 7 bHLH transcription factors, and 3 PP2C genes with MsNF-YB7 were found ( Figure 10C and Table S5). Regarding MsNF-YB15, there were 30 coexpressed genes, which were mainly MYB transcription factors ( Figure 10D and Table S6). In total, 38 genes were found to be coexpressed with MsNF-YC5 ( Figure 10E and Table S7). Regarding MsNF-YC6, 17 genes were coexpressed. Among them, 10 genes encode CBL-interacting serine/threonine-protein kinases, indicating that MsNF-YC6 may be involved in the Ca 2+ signaling process under salt stress ( Figure 10F and Table S8). These results provide information for a regulatory network analysis of these MsNF-Y genes. 7 bHLH transcription factors, and 3 PP2C genes with MsNF-YB7 were found ( Figure 10C and Table S5). Regarding MsNF-YB15, there were 30 coexpressed genes, which were mainly MYB transcription factors ( Figure 10D and Table S6). In total, 38 genes were found to be coexpressed with MsNF-YC5 ( Figure 10E and Table S7). Regarding MsNF-YC6, 17 genes were coexpressed. Among them, 10 genes encode CBL-interacting serine/threonineprotein kinases, indicating that MsNF-YC6 may be involved in the Ca 2+ signaling process under salt stress ( Figure 10F and Table S8). These results provide information for a regulatory network analysis of these MsNF-Y genes.

RT-qPCR Analysis of MsNF-Y Genes under Salt Conditions
To confirm the salt-induced expression profiles of six MsNF-Y genes, a RT-qPCR was further applied. As shown in Figure 11, the expression profiles of the candidate genes were consistent with RNA-Seq results. The expression of MsNF-YB2 was significantly upregulated at S1 stage and then decreased gradually. The expressions of MsNF-YB5, MsNF-YB7, MsNF-YB15, and MsNF-YC6 were increased gradually after salt treatment and reached the highest at S6 stage. The expression of MsNF-YC5 was comparatively high at S4 and S5 stage. The RT-qPCR results verified the results of RNA-Seq and indicated that these six MsNF-Y genes identified by our analyses were important candidate genes and may play roles in response to salt stress in alfalfa.

RT-qPCR Analysis of MsNF-Y Genes under Salt Conditions
To confirm the salt-induced expression profiles of six MsNF-Y genes, a RT-qPCR was further applied. As shown in Figure 11, the expression profiles of the candidate genes were consistent with RNA-Seq results. The expression of MsNF-YB2 was significantly up-regulated at S1 stage and then decreased gradually. The expressions of MsNF-YB5, MsNF-YB7, MsNF-YB15, and MsNF-YC6 were increased gradually after salt treatment and reached the highest at S6 stage. The expression of MsNF-YC5 was comparatively high at S4 and S5 stage. The RT-qPCR results verified the results of RNA-Seq and indicated that these six MsNF-Y genes identified by our analyses were important candidate genes and may play roles in response to salt stress in alfalfa.

Discussion
As a gene family important for plant growth and development and abiotic stress resistance, NF-Y transcription factors have been identified and functionally analyzed genome-wide in many plant species. For example, 36 NF-Y gene members (10 NF-YAs, 13 NF-YBs, and 13 NF-YCs) were identified in Arabidopsis [25], 78 NF-Y gene members were identified in soybean (21 NF-Yas, 32 NF-Ybs, and 25 NF-YCs) [19], 34 NF-Y gene members were identified in rice (11 NF-YAs, 11 NF-YBs, and 12 NF-YCs) [26], and 50 NF-Y gene members were identified in maize (14 NF-Yas, 18 NF-Ybs, and 18 NF-YCs) [27]. However, because alfalfa is autotetraploid, genome assembly is difficult, and it was not until 2020 that the alfalfa reference genome was assembled successfully [1]. In this study, in total, 60 MsNF-Y members (9 NF-YAs, 26 NF-YBs, and 25 NF-YCs) were identified in the alfalfa reference genome, and their gene structures, sequence features, conserved domains, and expression patterns were systematically analyzed, particularly the expression changes of MsNF-Ys under different abiotic stresses. These findings provide valuable information for a subsequent functional analysis and molecular regulatory network construction of a single MsNF-Y gene.
Many studies provide evidence suggesting that NF-Y transcription factors are widely involved in plant flowering time regulation and seed development. The overexpression of AtNF-YB1 in Arabidopsis resulted in delayed flowering in plants [28], while a phylogenetic analysis found that MsNF-YB1/2/3/4 had high homology with AtNF-YB1 ( Figure 6B). Furthermore, through an expression pattern analysis, MsNF-YB3 was found to be specifically expressed in flowers ( Figure 7B); thus, we inferred that MsNF-YB3 might be an im-

Discussion
As a gene family important for plant growth and development and abiotic stress resistance, NF-Y transcription factors have been identified and functionally analyzed genome-wide in many plant species. For example, 36 NF-Y gene members (10 NF-YAs, 13 NF-YBs, and 13 NF-YCs) were identified in Arabidopsis [25], 78 NF-Y gene members were identified in soybean (21 NF-Yas, 32 NF-Ybs, and 25 NF-YCs) [19], 34 NF-Y gene members were identified in rice (11 NF-YAs, 11 NF-YBs, and 12 NF-YCs) [26], and 50 NF-Y gene members were identified in maize (14 NF-Yas, 18 NF-Ybs, and 18 NF-YCs) [27]. However, because alfalfa is autotetraploid, genome assembly is difficult, and it was not until 2020 that the alfalfa reference genome was assembled successfully [1]. In this study, in total, 60 MsNF-Y members (9 NF-YAs, 26 NF-YBs, and 25 NF-YCs) were identified in the alfalfa reference genome, and their gene structures, sequence features, conserved domains, and expression patterns were systematically analyzed, particularly the expression changes of MsNF-Ys under different abiotic stresses. These findings provide valuable information for a subsequent functional analysis and molecular regulatory network construction of a single MsNF-Y gene.
Many studies provide evidence suggesting that NF-Y transcription factors are widely involved in plant flowering time regulation and seed development. The overexpression of AtNF-YB1 in Arabidopsis resulted in delayed flowering in plants [28], while a phylogenetic analysis found that MsNF-YB1/2/3/4 had high homology with AtNF-YB1 ( Figure 6B). Furthermore, through an expression pattern analysis, MsNF-YB3 was found to be specifically expressed in flowers ( Figure 7B); thus, we inferred that MsNF-YB3 might be an important candidate gene for flowering time regulation in alfalfa. AtNF-YB9 regulates seed development by integrating light signals with hormonal signals, and the homologous gene AtNF-YB6 is involved in the morphogenesis of seed embryos by mediating the ABA signaling pathway [28]. In the present study, we found that MsNF-YB20/22/23/25 and AtNF-YB6/9 were clustered in one evolutionary branch ( Figure 6B); therefore, the four MsNF-YB genes might also be involved in the seed development process in alfalfa.
The NF-Y transcription factor is not only involved in plant growth and development but also plays an important role in plant abiotic stress [4]. The overexpression of AtNF-YB1 and ZmNF-YB2 in Arabidopsis and maize, respectively, significantly improved drought resistance in the transgenic lines [10,28]. Similarly, the overexpression of soybean GmNF-YA3 or poplar PdNF-YB7 in Arabidopsis can reduce leaf water loss and improve plant drought resistance [29,30]. The overexpression of AtNF-YA1 can significantly improve the resistance of Arabidopsis to salt stress [28]. GmNF-YA16, GmNF-YB2, GmNF-YC13, and GmNF-YC14 in soybean can significantly enhance the resistance of plants to drought and salt stress [19,22,29]. In this study, MsNF-YB14 and MsNF-YB15 were homologous genes of GmNF-YB2, while MsNF-YC5 and MsNF-YC6 were homologous genes of GmNF-YC13 ( Figure 6). Their expression was induced by drought and salt stress ( Figure 8). Therefore, MsNF-YB14/15 and MsNF-YC5/6 may also perform the function of regulating alfalfa resistance to drought and salt stress.
NF-Y transcription factors regulate drought and salt stress in plants by mediating multiple signaling and plant hormone pathways. The regulatory pathways mainly include two types. One type mediates the ABA signaling pathway to regulate expression changes in downstream genes. For example, PdNF-YB21 can interact with PdFUS3 to jointly activate the expression of downstream PdNCED3 genes, thereby promoting the synthesis of ABA and ultimately improving drought tolerance in poplars [31]. In soybean, GmNF-YC14 can interact with GmNF-YB2 and GmNF-YA16 to form a complete, active NF-Y transcriptional complex, and this complex can regulate the ABA signaling pathway mediated by the ABA receptor PYR/PYL to enhance drought resistance and salt tolerance in soybeans [22]. WGCNA, a method commonly used to identify coexpression regulatory networks, has been shown to be highly effective in many studies [32]. In this study, we used WGCNA and found that MsNF-YB2 was coexpressed with four DEAD genes in the ABA signaling pathway under salt stress, while MsNF-YB7 was coexpressed with three PP2C genes in the ABA pathway under salt stress ( Figure 10). Another regulatory pathway of NF-Y genes alters the stress resistance ability of plants by mediating independent ABA pathway genes. After overexpressing AtNF-YB1 in Arabidopsis, although drought resistance was improved in the transgenic lines, there was no significant change in drought resistance genes related to the ABA pathway, suggesting that NF-Y family members exist in a drought resistance pathway independent of the ABA signaling pathway [10,28]. In this study, many genes independent of the ABA pathway were found to be coexpressed with the MsNF-Y genes under drought and salt stresses. For example, 6 MYB transcription factors were coexpressed with MsNF-YB5, 16 MYB transcription factors were coexpressed with MsNF-YB15, and 10 CBL genes were coexpressed with the MsNF-YC6 under salt stress ( Figure 10). These results could greatly enhance our understanding of the molecular regulatory network of MsNF-Y genes in the future, especially under abiotic stress.
In summary, we provide genome-wide results of the alfalfa NF-Y genes and expression patterns and coexpression regulatory networks in different tissues under different abiotic stresses. The information described here could contribute to further studies investigating alfalfa NF-Y gene families, especially in the context of abiotic stress.

Identification of NF-Y Genes in the Alfalfa Genome
The reference genome of Xinjiangdaye was obtained from https://figshare.com/ (accessed on 14 November 2021). The HMMER3.0 program was used to identify the NF-Y gene family members in alfalfa based on e < 10 −10 [33]. The protein sequences and the conserved domains of the NF-Y subfamilies were collected as described in Siefers et al. [25]. The candidate NF-Y gene family members in alfalfa were screened via CD-search (https://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi, accessed on 22 December 2021) by default to analyze the presence of the conserved domains. The molecular weight (MW) and isoelectric point (pI) were analyzed by ExPASy (https://web.expasy.org, accessed on 29 December 2021).

Motif Analysis and Sequence Alignments
A motif analysis was conducted with the online program MEME (https://meme-suite. org/, accessed on 10 January 2022) by using the protein sequences of these NF-Y gene family members [34]. The gene structure was obtained by using TBtools software [35]. DNAMAN software was used to conduct the multiple sequence alignments.

Phylogenetic Tree Construction, Genome Localization and Gene Duplication
The NF-Y gene family sequences from Arabidopsis, rice, soybean, and Tartary buckwheat were collected from previous studies [21,[24][25][26]. The neighbor-joining method with 1000 replicated bootstrap values was used to conduct a phylogenetic analysis by using MEGA 7.0 software based on the protein sequences. A duplication event analysis, including the parameters Ka and Ks, and a genome localization analysis were conducted by TBtools [35].

cis-Element Analysis
The 2 kb sequences upstream of ATG, which is the translation start codon of all NF-Y gene family members, were downloaded from the Xinjiangdaye reference genome by using TBtools [35]. Then, the cis elements of these sequences were analyzed using the online program PlantCARE, http://bioinformatics.psb.ugent.be/webtools/plantcare/ html/ (accessed on 8 February 2022) [36].

Weighted Gene Coexpression Network Construction
We analyzed and constructed the weighted gene coexpression network (WGCNA) of alfalfa under salt stress using TBtools with the default parameters [32,35]. First, using transcriptome data (salt stress) collected at seven different time points at 0 h (CK), 0.5 h (S1), 1 h (S2), 3 h (S3), 6 h (S4), 12 h (S5), and 24 h (S6), all expressed genes were filtered based on FPKM values (more than 1) for the subsequent analysis. Then, the correlation between different modules and the salt stress treatment time was obtained using a cluster analysis and correlation analysis. Finally, coexpression network drawing was performed using Cytoscape software [41].

Expression Analysis by qRT-PCR
To further verify the expression of candidate genes under salt stress, we collected the root tips of Zhongmu No.1 under NaCl treatment at 0 h (CK), 0.5 h (S1), 1 h (S2), 3 h (S3), 6 h (S4), 12 h (S5), and 24 h (S6) according to a previous study [38]. All the samples were three biological repetitions. Total RNAs were extracted using an RNAprep Pure Plant Kit (Tiangen Biotech) and cDNAs were prepared using Tran-Script-Uni One-Step gDNA Removal and cDNA Synthesis SuperMix (TRAN). Then, we performed qPCR using a SYBR Green RT-PCR Kit (Takara). The primers were designed by Primer Premier 5 software, and the primer sequences were presented in Table S9. The qRT-PCR program was: 95 • C for 30 s; 40 cycles of 95 • C for 5 s; 58 • C for 34 s; and 95 • C for 15 s. MsActin was used to as the internal reference gene for data normalization analysis. For each analysis, three technical replicates from three biological replicates were conducted. Ultimately, we calculated the relative expression of MsNF-Ys by the 2 −∆∆Ct method.

Conclusions
In our study, in total, 60 MsNF-Y genes, including 9 MsNF-YAs, 26 MsNF-YBs, and 25 MsNF-YCs, were identified in the alfalfa reference genome, and conserved domain, gene structure, genomic location, cis element, and expression pattern analyses were conducted. The expression analysis showed that 48 MsNF-Y genes were expressed in all tissues at different levels, showing that these MsNF-Y genes play an important role in alfalfa growth and development. In addition, we analyzed the expression of MsNF-Y genes under different abiotic stresses (cold, drought, and salt). The results indicated that a few MsNF-Y genes were involved in these abiotic stresses. Finally, the co-expression networks under salt stress of MsNF-YB2, MsNF-YB5, MsNF-YB7, MsNF-YB15, MsNF-YC5, and MsNF-YC6 were constructed by using a WGCNA. In future studies, these candidate genes should be overexpressed or gene knocked out in alfalfa to verify their biological functions and the yeast two hybrid system or Chip-seq should be used to identify the interaction network for these candidate genes. Overall, our results could provide valuable information for further elucidating the biological functions of MsNF-Ys and improving salt tolerance and other abiotic stress resistance in alfalfa.