Genome-wide analysis of the nuclear factor Y (NF-Y) gene family in bermudagrass (Cynodon dactylon L.) identified an NF-Y heterotrimer comprising CdNF-YA-1A1/CdNF-YB-7B2/CdNF-YC-1A2

Nuclear factor Y (NF-Y), also known as heme activator protein (HAP) or CCAAT-binding factor (CBF), is a highly conserved transcription factor (TF) in all eukaryotes^[1]. Functional NF-Ys are comprised of three distinct subunits: NF-YA (HAP2/CBF-B), NF-YB (HAP3/CBF-A), and NF-YC (HAP5/CBF-C). The core region of NF-YA protein is composed of two conserved α-helix domains, whereas three conserved α-helix collectively form a histone-fold domain (HFD) in both NF-YB and NF-YC proteins^[2]. To perform their regulatory functions, NF-YBs and NF-YCs firstly form heterodimers in the cytoplasm and enter the nucleus to further form heterotrimers with NF-YAs^[3]. The resulting heterotrimeric NF-Ys can bind to the CCAAT-box cis-element in the promoter of target genes to regulate their expression. In some cases, the heterodimers of NF-YBs and NF-YCs could also interact with other classes of TFs rather than NF-YAs to function as transcription regulators^[3].

In yeast and animals, each subunit of NF-Y is encoded by a single gene. In contrast, the three subunits of NF-Y are encoded by multiple genes in plants^[4]. For example, the model plant Arabidopsis (Arabidopsis thaliana) has ten AtNF-YA, ten AtNF-YB, and ten AtNF-YC genes, whereas rice (Oryza sativa) has ten OsNF-YA, 11 OsNF-YB, and 12 OsNF-YC genes^[5]. Different heterotrimeric NF-Ys can be formed from these subunits with different tissue and stress-induced expression specificities, thus playing essential regulatory roles in the growth, development, and stress response of plants^[6]. In Arabidopsis, two NF-YA proteins (AtNF-YA2 and AtNF-YA6), associated with AtNF-YB2 and AtNF-YC3, could bind to the CCAAT-box in the promoter sequence of Flowering Locus T gene to regulate its expression, which in turn controls the flowering time^[7]. The heterotrimeric NF-Y composed of AtNF-YA2, AtNF-YB3, and AtNF-YC10, could enhance the expression of many heat stress-inducible genes and is required for the proper heat stress response^[8]. In rice, aleurone layer-specific OsNF-YB1 could interact with OsNF-YC11 and OsNF-YC12 to regulate the expression of downstream genes, thus contributing to the grain filling and endosperm development^[9].

As an important warm-season turfgrass species that originated from Africa, bermudagrass (Cynodon dactylon L.) has great economic value^[10,11]. Elite bermudagrass cultivars are widely planted in sports fields, public parks, and home lawns to produce uniform and high-quality turfs in tropical and subtropical areas around the world^[12,13]. The usage of bermudagrass as a turfgrass species is mainly derived from its fast propagation rate, perennial growth mode, and strong tolerance to multiple stresses^[14]. Elucidation of the mechanisms regulating the growth, development, and stress response of bermudagrass is vital for breeding new bermudagrass cultivars. Two previous studies reported the isolation of an OsNF-YC4 homologous gene from hybrid bermudagrass (C. dactylon × C. transvaalensis) and found that the CdtNF-YC1 gene could increase the drought and salinity tolerance of both transgenic rice and seashore paspalum (Paspalum vaginatum)^[15,16]. However, systematic analysis of the NF-YC gene family is still deficient in bermudagrass, which impedes the application of NF-YC genes in the molecular breeding studies of bermudagrass.

In this study, we performed a comprehensive genome-wide analysis of the NF-Y gene family in bermudagrass. The phylogenetic relationships, gene structure characteristics, protein domains, and expression patterns of the identified CdNF-Y genes were all analyzed to understand their functional divergence. Furthermore, the NF-Y heterotrimer consisting of CdNF-YA-1A1, CdNF-YB-7B2, and CdNF-YC-1A2 was also identified using the combination of yeast two-hybrid (Y2H) assay, yeast three-hybrid (Y3H) assay, and bimolecular fluorescence complementation (BiFC) assay. Results of this study established a solid foundation to further characterize the function of NF-Y genes in bermudagrass.

C. dactylon cultivar Yangjiang and Nicotiana benthamiana accession LAB were used in this study. The C. dactylon plants were ordinarily grown at turfgrass plots of the Yangzhou University (longitude and latitude：32°35' N, 119°40' E; soil type: 80% river sand, 20% peat soil) under standard turf management conditions. For heat stress treatment, the C. dactylon plants were grown in a growth chamber with a 16/8 h light/dark cycle at 25 °C for one week, then transferred into a growth chamber of 42 °C for 12 h^[17]. The N. benthamiana plants were grown in a growth chamber with a 16/8 h light/dark cycle at 22 °C.

Database search and sequence retrieval

The translated amino acid sequences of NF-Y genes from A. thaliana, Brachypodium distachyon, and O. sativa were downloaded from the NCBI database (www.ncbi.nlm.nih.gov)^[18−20] and were BLAST searched against the C. dactylon protein database (76,879 sequences, 33,188,769 residues) annotated from the genome^[21] to identify putative bermudagrass NF-Y proteins with the following parameters: score value ≥ 100 and E-value < 10⁻⁵. The identified homologous proteins were further BLAST searched against the protein sequence repositories of NCBI to validate the correct identities. After manually removing the redundant sequences, the finally obtained CdNF-Y proteins and the corresponding encoding genes were named according to their chromosomal locations.

Phylogenetic and collinearity analysis

The amino acid sequences of NF-Y proteins were aligned using the ClustalW algorithm of MEGA X^[22]. The aligned amino acid sequences were used to construct the maximum-likelihood (ML) phylogenetic tree with the following parameters: complete deletion and 2,000 replicates of bootstrap^[23]. The collinearity of A. thaliana, B. distachyon, O. sativa, and C. dactylon NF-Y genes was analyzed and plotted using the MCScan toolkit of TBtools software with default parameters^[24]. DnaSP 6 software was used to calculate non-synchronous (Ka) and synchronous (Ks) values of the identified gene pairs^[25].

Gene structure, protein domain, and conserved motif analysis

The mRNA and genomic DNA sequences of CdNF-Y genes were analyzed using the gene structure display server (https://gsds.gao-lab.org) to obtain the gene structures^[26]. The amino acid sequences of CdNF-Y proteins were searched against the InterPro database (www.ebi.ac.uk/interpro/search/sequence) and MEME (https://meme-suite.org/meme/tools/meme) with default parameters to acquire the sequence composition and location of conserved protein domains and motifs, respectively.

Gene expression and promoter cis-regulatory element analysis

Raw sequencing data of transcriptomic analysis of bermudagrass plants under different stress treatments were downloaded from the NCBI SRA database (www.ncbi.nlm.nih.gov/sra)^[27−30] and mapped to bermudagrass genes using the HISAT2 software (http://daehwankimlab.github.io/hisat2) with default parameters^[31]. The numbers of mapped reads for each gene were converted to FPKM (fragments per kilobase of transcript per million mapped fragments) values. CdNF-Y gene expression data were extracted from the gene expression profiles and visualized as heatmaps using the MultiExperiment Viewer (MeV) software. The 2,000 bp nucleotide sequences upstream of CdNF-Y genes were extracted from the bermudagrass genomic DNA sequence as putative promoter sequences, which were searched against the PlantCARE database (https://bioinformatics.psb.ugent.be/webtools/plantcare/html) to identify the putative cis-regulatory elements^[32].

RT-qPCR analysis

Total RNA was extracted from different tissues of normally growing bermudagrass plants and leaves of bermudagrass plants under heat stress treatment using the Plant RNA Midi Kit (Omega, Norcross, USA). cDNA was synthesized using the SuperScript II reverse transcriptase kit (Takara, Beijing, China). Quantitative PCR was performed on a Mini Opticon Real-Time PCR System (Bio-Rad, Hercules, USA) using the Hieff® qPCR SYBR Green Master Mix (Yeasen Biotechnology, Shanghai, China), and the CdSUD1 gene was used as a standard control (Supplementary Table S1). The relative gene expression level was calculated using the 2^−ΔΔCᴛ method^[33]. Results of three biological replicates were statistically analyzed using Tukey's multiple comparison test with the SPSS 16.0 software, which was declared at the 0.05 level of probability.

Gene cloning and sequence analysis

Total RNA was extracted from the leaves of normally growing bermudagrass plants and reverse transcribed to cDNA as described above. The full-length coding sequences (CDSs) of CdNF-YA-1A1, CdNF-YB-7B2, and CdNF-YC-1A2 were PCR amplified from the synthesized cDNA using the primers listed in Supplementary Table S1 and sequenced to confirm the clone composition. The translated amino acid sequences of cloned CdNF-YA-1A1, CdNF-YB-7B2, and CdNF-YC-1A2 genes and homolog genes from other plant species were aligned using the DNAMAN software with default parameters.

Subcellular localization analysis

The full-length CDSs of CdNF-YA-1A1, CdNF-YB-7B2, and CdNF-YC-1A2 without the stop codon were amplified by PCR, digested with XbaI and KpnI, and cloned into the pBI121-GFP vector. The resulting construct was transformed into Agrobacterium tumefaciens GV3101, which was injected into leaves of three-week-old N. benthamiana plants along with the control GV3101 strain carrying the pBI121-mCherry-NLS plasmid. GFP and mCherry fluorescence signals were observed after 3 d under an LSM 880NLO confocal microscope (Carl Zeiss, Jena, Germany).

Transactivation activity assay

The full-length CDSs of CdNF-YA-1A1, CdNF-YB-7B2, and CdNF-YC-1A2 were amplified by PCR, digested with EcoRI and BamHI, and cloned into the pGBKT7 vector. The resulting constructs, along with the empty vector (negative control) and pGBKT7-CdTB1 plasmid (positive control)^[34], were separately transformed into yeast strain AH109 using the LiAc/PEG method and grown on SD-Trp medium^[35]. The positive transformants were cultured on SD-Trp-His-Ade selective medium supplemented with 2.5 mM 3-AT to test the transcriptional activation activity.

Yeast two-hybrid assay

The full-length CDSs of CdNF-YB-7B2 and CdNF-YC-1A2 were amplified by PCR, digested with EcoRI and BamHI, and cloned into the pGADT7 vector. The resulting pGADT7-CdNF-YB-7B2 and pGADT7-CdNF-YC-1A2 constructs were individually combined with the above-constructed pGBKT7-CdNF-YA-1A1 and pGBKT7-CdNF-YB-7B2 plasmids and cotransformed into the yeast strain AH109 using the LiAc/PEG method and grown on SD-Leu-Trp medium^[35]. The positive transformants were then cultured on SD-Leu-Trp-His-Ade selective medium supplemented with 2.5 mM 3-AT to determine the protein-protein interactions.

Yeast three-hybrid assay

The full-length CDSs of CdNF-YB-7B2 and CdNF-YC-1A2 were amplified by PCR, digested with EcoRI/BamHI and NotI/BglII, respectively, and jointly cloned into the pBridge vector. The full-length CDS of CdNF-YA-1A1 was amplified by PCR, digested with EcoRI and BamHI, and cloned into the pGADT7 vector. The resulting pBridge-CdNF-YB-7B2-CdNF-YC-1A2 and pGADT7-CdNF-YA-1A1 constructs were combined and cotransformed into the yeast strain AH109 using the LiAc/PEG method and grown on SD-Leu-Trp-Met medium^[35]. The positive transformants were then cultured on SD-Leu-Met-Trp-His-Ade selective medium supplemented with 2.5 mM 3-AT to determine the protein-protein interactions.

Bimolecular fluorescence complementation assay

The full-length CDSs of CdNF-YA-1A1, CdNF-YB-7B2, and CdNF-YC-1A2 without the stop codon were amplified by PCR and cloned into the pSm35s-nYFP and pSm35s-cYFP vectors via the Gibson Assembly method^[36]. The full-length CDS of CdNF-YA-1A1 without the stop codon was amplified by PCR, digested with XbaI and KpnI, and cloned into the pBI121-mCherry vector. The resulting constructs were individually transformed into A. tumefaciens GV3101. The A. tumefaciens strains carrying pSm35s-nYFP, pSm35s-cYFP, pSm35s-nYFP-CdNF-YA-1A1, pSm35s-nYFP-CdNF-YB-7B2, pSm35s-cYFP-CdNF-YA-1A1, pSm35s-cYFP-CdNF-YC-1A2, and pBI121-mCherry-CdNF-YA-1A1 plasmids were mixed and injected into the leaves of 3-week-old N. benthamiana plants. YFP and mCherry fluorescence signal was observed after 3 d under the LSM 880NLO confocal microscope.

Using NF-Y protein sequences from A. thaliana, B. distachyon, and O. sativa as references, a total of 88 NF-Y genes, including 32 NF-YA, 31 NF-YB, and 25 NF-YC genes, were successfully identified in the C. dactylon genome through BLAST searches and manual validations (Supplementary Table S2). Because the C. dactylon genome is comprised of four haplotypes of subgenomes (A1, A2, B1, and B2) with high sequence similarity^[21], the 88 CdNF-Y genes could be classified into 28 homoeologous gene groups with one to four allelic genes, and were named as CdNF-YA-1A1 to CdNF-YA-9B2, CdNF-YB-1A1 to CdNF-YB-11A2, and CdNF-YC-1A1 to CdNF-YC-8B2, respectively, according to their chromosomal locations and sequence similarities (Supplementary Fig. S1, Supplementary Table S2).

Phylogenetic analysis indicated that 88 CdNF-Y genes, along with 30 AtNF-Y, 33 OsNF-Y and 31 BdNF-Y genes, were distinctively divided into three subfamilies, which was consistent with the gene identification and classification result (Fig. 1a). In line with the species divergence time, the clade distances between CdNF-Y genes and homologous OsNF-Y/BdNF-Y genes were shorter than those of CdNF-Y genes and homologous AtNF-Y genes (Fig. 1a). Accordingly, collinearity analysis indicated that CdNF-Y, OsNF-Y and BdNF-Y genes formed 23 homologous gene pairs, whereas only five homologous gene pairs could be detected among the three grass species and A. thaliana (Fig. 1b). Inter-species collinearity analysis also identified 108 homologous gene pairs from the 88 CdNF-Y genes (Fig. 1b). Interestingly, Ka/Ks ratios of 88.8% (96/108) homologous gene pairs were less than one (Supplementary Table S3), suggesting these CdNF-Y genes might have suffered robust purifying selective pressure during evolution.

Figure 1. 

Genomic characterization of NF-Y genes in bermudagrass. (a) Phylogenetic relationship. (b) Collinearity of NF-Y genes in C. dactylon, B. distachyon, O. sativa, and A. thaliana.

Sequence divergence of NF-Y genes in bermudagrass

The three subfamilies of CdNF-Y genes had different gene sequence characteristics. Notably, the average ORF length of CdNF-YA genes (831.8 bp) was longer than that of CdNF-YB and CdNF-YC genes (637.2 and 646.44 bp, respectively) (Supplementary Table S2). Accordingly, gene structure analysis revealed that CdNF-YA genes had complicated exon-intron organizations with at least four introns, whereas most CdNF-YB and CdNF-YC genes were intronless or only contained a single intron (Fig. 2a, b). Exceptionally, three homoeologous gene groups of CdNF-YB genes (CdNF-YB3, CdNF-YB5, and CdNF-YB6) clustering as one clade in the phylogenetic analysis also had two to three introns (Fig. 2a, b).

Figure 2. 

Sequence characteristics of bermudagrass NF-Y genes. (a) Phylogenetic relationship. (b) Gene structure of CdNF-Y genes. (c) Conserved protein domain of CdNF-Y proteins.

In agreement with the phylogenetic analysis result, protein domain analysis indicated that the three subfamilies of CdNF-Y proteins harbored the conserved CBFB_NFYA, HFD_NFYB, and HFD_NFYC domains, respectively (Fig. 2c). However, the location of the conserved protein domain in the protein sequence exhibited variation among the CdNF-Y proteins of different homoeologous gene groups. Notably, CdNF-YB-6A1, CdNF-YB-6B1, and CdNF-YB-6B2 proteins further harbored a VQ domain after the HFD_NFYB domain (Fig. 2c). Accordingly, motif analysis indicated that motif 1, representing conserved domains of CBFB_NFYA and HFD_NFYB/YC, was identified in nearly all (84/88) CdNF-Y proteins, whereas motifs 3, 7, and 14, representing unique domains of CBFB_NFYA, HFD_NFYB, and HFD_NFYC, respectively, were only identified in corresponding CdNF-YA, CdNF-YB, and CdNF-YC proteins (Supplementary Fig. S2, Supplementary Table S4). Except for these conserved motifs, many motifs (9, 10, 12, 13, 15, 16, 17, 19, and 20) were only identified in a few CdNF-Y proteins belonging to some specific homoeologous gene groups (Supplementary Fig. S2, Supplementary Table S4), suggesting these CdNF-Y proteins might have specialized sequence characteristics.

Expression profiling of NF-Y genes in bermudagrass

To gain insights into the putative biological functions, the expression patterns of 88 CdNF-Y genes were examined using the transcriptome sequencing data. The result indicated that some CdNF-Y genes, including CdNF-YA-3A1, CdNF-YA-6A1, CdNF-YA-6B2, CdNF-YB-3B1, CdNF-YB-3B2, CdNF-YC-2B1, CdNF-YC-6A1, and CdNF-YC-6B2, were highly expressed in all organs and maintained the high expression abundance under stress treatments (Fig. 3, Supplementary Table S5). By contrast, some other CdNF-Y genes were preferentially expressed in specific organs or exhibited expression variance under specific stress treatment. For example, CdNF-YA-1B1, CdNF-YB-1A1, and CdNF-YB-5A2 were preferentially expressed in the inflorescence, rhizome, and root, respectively (Fig. 3, Supplementary Table S5). On the other hand, the expression of CdNF-YA-4A1 was decreased under heat stress treatment and increased under mild cold (20 °C) stress treatment, whereas the expression of CdNF-YC-2A1 was induced by salt stress treatment (Fig. 3, Supplementary Table S5). In agreement with the transcriptome sequencing results, RT-qPCR analyses of six randomly selected CdNF-Y genes indicated that CdNF-YA-1A1, CdNF-YB-7B2 and CdNF-YC-1A2 were preferentially expressed in the leaf, CdNF-YA-4B2 was highly expressed in the stolon and rhizome, CdNF-YB-9B2 was highly expressed in the root, whereas CdNF-YC-7A1 was highly expressed in the inflorescence (Fig. 4).

Figure 3. 

Expression of bermudagrass NF-Y genes in six organs and under different stress treatments. Left panel, phylogenetic tree of CdNF-Y genes. Right panel, heatmap showing the average expression abundance of the transcriptome dataset.

Figure 4. 

Expression analysis of the six bermudagrass NF-Y genes using RT-qPCR. The expression of (a) CdNF-YA-1A1, (b) CdNF-YA-4B2, (c) CdNF-YB-7B2, (d) CdNF-YB-9B2, (e) CdNF-YC-1A2, and (f) CdNF-YC-7A1 was quantified using the CdSUD1 gene as a reference. Error bars represent SE. Different letters indicate significant differences determined by Tukey's multiple comparison test.

To better understand the different expression pattern of CdNF-Y genes, the 2,000 bp upstream promoter sequences were further analyzed. The result indicated that multiple phytohormone-responsive and stress-responsive cis-regulatory elements were widely distributed in the promoter sequences of 88 CdNF-Y genes (Supplementary Fig. S3). Specifically, two tandem low-temperature responsive cis-elements were identified in the CdNF-YA-4A1 gene promoter sequence, which is in accordance with the gene expression variation under heat and cold stress treatments (Fig. 3). Interestingly, except CdNF-YA-2A2, CdNF-YA-2B2, CdNF-YA-6A2 and CdNF-YB-1A1, all other 84 CdNF-Y genes harbored at least one CCAAT-box cis-element in the promoter sequences, suggesting their expression were feedback regulated by themselves or other CdNF-Y genes (Supplementary Fig. S3).

Molecular characterization of the CdNF-YA-1A1/CdNF-YB-7B2/CdNF-YC-1A2 heterotrimer

RT-qPCR analysis revealed that CdNF-YA-1A1, CdNF-YB-7B2 and CdNF-YC-1A2 genes exhibited similar tissue expression patterns (Fig. 4). Moreover, the expression abundance of the three genes also decreased similarly under heat stress treatments (Fig. 3, Supplementary Fig. S5, Supplementary Table S5). These results suggested that the encoding proteins of the three genes might function together as an NF-Y heterotrimer. To verify this suspicion, the coding sequences of three genes were amplified through RT-PCR (Fig. 5a). Sequence analysis indicated that the encoded CdNF-YA-1A1, CdNF-YB-7B2, and CdNF-YC-1A2 proteins all had conserved DNA-binding domains (Fig. 5b, Supplementary Fig. S4), suggesting the three proteins might function as transcription factors. Accordingly, subcellular localization analysis indicated that CdNF-YA-1A1 protein was located in the nucleus, whereas both CdNF-YB-7B2 and CdNF-YC-1A2 proteins were located in the nucleus and cytoplasm simultaneously (Fig. 6a). However, yeast transactivation activity assay indicated that the three proteins alone had no transactivation activity (Fig. 6b).

Figure 5. 

Cloning of the CdNF-YA-1A1, CdNF-YB-7B2, and CdNF-YC-1A2 genes. (a) Gel electrophoresis showing the PCR amplification of CdNF-YA-1A1, CdNF-YB-7B2, and CdNF-YC-1A2 genes. The nucleic acid sequence of (b) CdNF-YA-1A1, (c) CdNF-YB-7B2, and (d) CdNF-YC-1A2 genes, and the deduced amino acid sequence of the encoding proteins.

Figure 6. 

Subcellular localization and transactivation activity analysis of the CdNF-YA-1A1, CdNF-YB-7B2, and CdNF-YC-1A2 proteins. (a) Subcellular location of CdNF-YA-1A1, CdNF-YB-7B2, and CdNF-YC-1A2 proteins in tobacco leaf epidermis cells. (b) Transactivation activity assay of CdNF-YA-1A1, CdNF-YB-7B2, and CdNF-YC-1A2 proteins in yeast cells.

Protein sequence analysis also indicated that CdNF-YA-1A1 protein had a NF-YB/C interaction domain, whereas CdNF-YB-7B2 and CdNF-YC-1A2 had a NF-YC/A interaction domain and a NF-YB/A interaction domain, respectively (Supplementary Fig. S4). In line with these observations, Y2H and BiFC assays collectively indicated that CdNF-YB-7B2 could interact with CdNF-YC-1A2 to form a heterodimer (Fig. 7). Y2H and BiFC assays also revealed that either CdNF-YB-7B2 or CdNF-YC-1A2 protein alone could not interact with CdNF-YA-1A1 (Fig. 7), however, Y3H assay and BiFC assays indicated that the CdNF-YB-7B2/CdNF-YC-1A2 heterodimer could further interact with CdNF-YA-1A1 to form a heterotrimer (Fig. 8).

Figure 7. 

Interaction analysis of the CdNF-YA-1A1, CdNF-YB-7B2, and CdNF-YC-1A2. (a) Yeast two-hybrid assay. (b) Bimolecular fluorescence complementation assay showing the positive interaction between CdNF-YB-7B2 and CdNF-YC-1A2 proteins.

Figure 8. 

Interaction analysis of the CdNF-YA-1A1 protein and the CdNF-YB-7B2/CdNF-YC-1A2 heterodimer. (a) Yeast three-hybrid assay. (b) Bimolecular fluorescence complementation assay showing the positive interaction between CdNF-YA-1A1 and the heterodimer of CdNF-YB-7B2/CdNF-YC-1A2.

Since its origination from Africa, bermudagrass has become a worldwide distributed perennial turfgrass species of great economic value. Unlike cereal grass species such as rice, wheat, and maize, bermudagrass not only has the erect-growing shoot but also has two types of prostrate growing stems, stolon and rhizome, which lay the foundation for its successful application in the turf industry since the two specialized stems provide bermudagrass with fast regeneration ability under frequent clipping^[37]. Moreover, some bermudagrass ecotypes are highly tolerant to abiotic stresses, especially salt and heat stresses^[38,39]. Elucidation of the growth and development regulatory mechanism, as well as the stress-responsive mechanism of bermudagrass plants, could not only promote the breeding of new turfgrass cultivars with required phenotypes but also assist the agronomic trait improvement of cereal grass species^[40]. NF-Ys are highly conserved transcription factors in regulating plant growth, development, and stress responses. In the past decades, several studies reported the characterization of NF-Y gene families in many important plant species, including A. thaliana^[5], N. tabacum^[41], O. sativa^[20], B. distachyon^[19], Triticum turgidum^[42], Zea mays^[43], Solanum lycopersicum^[44], S. tuberosum^[45], Vitis vinifera^[46], Medicago sativa^[47], Populus trichocarpa^[48], Camellia sinensis^[49], and Panicum virgatum^[50]. In the current study, we performed the genome-wide analysis of the NF-Y gene family in bermudagrass for the first time. Results of this study provided a solid basis for the in-depth functional study and breeding application of NF-Y genes in bermudagrass.

A total of 88 CdNF-Y genes, including 32 CdNF-YA, 31 CdNF-YB, and 25 CdNF-YC genes, were identified in the C. dactylon genome (Fig. 1, Supplementary Table S2). At first glance, the number of CdNF-Y genes is much larger than that of other plants, especially O. sativa (33 OsNF-Ys), B. distachyon (31 OsNF-Ys), Z. mays (52 ZmNF-Ys), and P. virgatum (47 PvNF-Ys), which also belong to the grass family^{[19,20,43,50]}. However, since bermudagrass is an allotetraploid plant with four sets of chromosomes, the 88 CdNF-Y genes are allelic copies of 28 unique NF-Y genes, which are actually fewer than those of the above-mentioned grass species. Accordingly, phylogenetic analysis revealed that orthologs of OsNF-YC8, OsNF-YC9, OsNF-YC10, OsNF-YC11, OsNF-YC12, BdNF-YC9, BdNF-YB5, and BdNF-YB10 were missing in bermudagrass (Fig. 1). Two whole-genome duplication events occurring at about 5.38 and 0.77 million years ago were suspected to successively form the A/B and A1/A2 (B1/B2) subgenomes of bermudagrass^[21]. Chromosome doubling and then rapid sequence arrangement of the genome would lead to massive gene loss in the polyploidization process^[51], which might explain the observed loss of NF-YB and NF-YC genes in bermudagrass. On the other hand, gene structure analysis revealed that all CdNF-YA genes had multiple exons and introns, whereas most of the CdNF-YB genes and all the CdNF-YC genes had zero or one intron (Fig. 2a), which is a sign of rapid evolution^[52]. Protein domain analysis further indicated that one homoeologous gene groups, CdNF-YB6, were mutated to harbor a surplus VQ domain in the encoded proteins (Fig. 2b). These results collectively implied that evolutionary selection pressure might be stricter on CdNF-YA genes than that of CdNF-YB and CdNF-YC genes.

Gene expression information could provide important clues for understanding the function of specific genes^[53]. Through analyzing the transcriptome sequencing data, the organ-specific and stress-responsive expression profiles of 88 CdNF-Y genes were successfully acquired (Fig. 3). Previous functional analysis of CdtNF-YC1, the homologous gene of OsNF-YC4 in hybrid bermudagrass, indicated that it was highly expressed in the root and rhizome, and could be induced by short-term chilling and salt stresses^[15]. In agreement with this result, the four allelic genes of CdNF-YC5, which showed the highest homology to OsNF-YC4, were all highly expressed in root and rhizome (Fig. 3, Supplementary Table S5). Furthermore, the expression of one allelic gene, CdNF-YC-5B2, was also induced by mild cold and short-term salt stresses (Fig. 3, Supplementary Table S5). In addition to supporting the previous studies, the expression profiling results further revealed that eight CdNF-Y genes (CdNF-YA-3A1, CdNF-YA-6A1, CdNF-YA-6B2, CdNF-YB-3B1, CdNF-YB-3B2, CdNF-YC-2B1, CdNF-YC-6A1, and CdNF-YC-6B2) were highly expressed in all organs of bermudagrass, whereas other CdNF-Y genes were only preferentially expressed in one or a few organs. Moreover, the constitutively expressed CdNF-Y genes were also stably expressed under different stress treatments (Fig. 3). These observations strongly suggested that these CdNF-Y genes might be housekeeping genes to sustain the growth of bermudagrass plants, whereas other CdNF-Y genes might play more specific roles in the growth and development of different organs as well as delicate responses to diverse stress conditions. Similar differentiation of housekeeping and specific NF-Y genes was also observed in O. sativa^[20], B. distachyon^[19], C. sinensis^[49], and P. virgatum^[50], suggesting the functional differentiation of NF-Y genes might be a universal phenomenon in plants.

Functional NF-Ys usually form NF-YA/NF-YB/NF-YC heterotrimers to perform transcriptional regulatory functions^[54]. In Arabidopsis, the AtNF-YA5/AtNF-YB9/AtNF-YC9 heterotrimer regulates the expression of chlorophyll A/B binding protein genes in etiolated seedlings^[55], while the AtNF-YA2/AtNF-YB3/AtNF-YC10 heterotrimer participates in the expression regulation of heat stress-inducible genes^[8]. Similarly, the OsNF-YA8/OsNF-YB9/OsNF-YC10 heterotrimer composes a natural gene on-off system to modulate thermotolerance in rice^[56]. In Z. mays, two heterotrimers, ZmNF-YA3/ZmNF-YB9/ZmNF-YC13 and ZmNF-YA3/ZmNF-YB10/ZmNF-YC13, synergistically regulate the plant architecture^[57]. In S. lycopersicum, the SlNF-YC3 protein could interact with SlNF-YB3a, SlNF-YB3b, and SlNF-YB5c, and further form NF-Y heterotrimers with SlNF-YA3a to regulate arbuscular mycorrhizal symbiosis, whereas SlNF-YC1a, SlNF-YC1b, SlNF-YC1d, and SlNF-YC9 proteins interact with SlNF-YB8a, SlNF-YB8b, and SlNF-YB8c proteins, and further recruit SlNF-YA1b and SlNF-YA9 to form functional NF-Y complexes and regulate flavonoid biosynthesis^[58,59]. Here, a new NF-Y heterotrimer comprising CdNF-YA-1A1/CdNF-YB-7B2/CdNF-YC-1A2 was successfully identified in bermudagrass by Y2H, Y3H, and BiFC assays (Figs. 7, 8). Notably, the CdNF-YA-1A1 and CdNF-YB-7B2 genes were homologous with AtNF-YA2 and AtNF-YB3, respectively (Fig. 1). It is possible that the CdNF-YA-1A1/CdNF-YB-7B2/CdNF-YC-1A2 heterotrimer might function as AtNF-YA2/AtNF-YB3/AtNF-YC10 to regulate heat stress response in bermudagrass. However, the functional mechanism of the two heterotrimers might be different since the expression of the three CdNF-Y genes was all down-regulated under heat stress treatment (Supplementary Fig. S5), whereas the expression of the AtNF-YB3 gene was induced by heat stress^[8]. In the future, transgenic bermudagrass plants with altered expression abundance of the three genes could be analyzed to characterize their detailed biological functions.

In this study, we successfully identified 88 NF-Y genes in bermudagrass. The sequence characteristics and expression profiles of 88 CdNF-Y genes were systematically analyzed. As observed in other plants, 88 CdNF-Y genes underwent different evolutionary selection pressures and showed various expression biases in bermudagrass organs and stress conditions. Moreover, an NF-Y heterotrimer consisting of CdNF-YA-1A1, CdNF-YB-7B2, and CdNF-YC-1A2 was further identified using Y2H, Y3H, and BiFC assays. Although the detailed function of the NF-Y heterotrimer remains unclear, gene expression data suggested it might participate in the heat stress response. In summary, the results of this study established a solid foundation for the in-depth functional study and breeding application of NF-Y genes in bermudagrass.