Screening and identification of key genes regulating fall dormancy in alfalfa leaves.

Fall dormancy (FD) determines the adaptation of an alfalfa variety and affects alfalfa production and quality. However, the molecular mechanism underlying FD remains poorly understood. Here, 44 genes regulating FD were identified by comparison of the transcriptomes from leaves of Maverick (fall-dormant alfalfa) and CUF101(non-fall-dormant), during FD and non-FD and were classified them depending on their function. The transcription of IAA-amino acid hydrolase ILR1-like 1, abscisic acid receptor PYL8, and monogalactosyldiacylglycerol synthase-3 in Maverick leaves was regulated by daylength and temperature, and the transcription of the abscisic acid receptor PYL8 was mainly affected by daylength. The changes in the expression of these genes and the abundance of their messenger RNA (mRNA) in Maverick leaves differed from those in CUF101 leaves, as evidenced by the correlation analysis of their mRNA abundance profiles obtained from April to October. The present findings suggested that these genes are involved in regulating FD in alfalfa.

Introduction

Fall dormancy(FD) is defined as an adaptive growth characteristic of alfalfa in autumn in response to the reduction of daylength and drop in temperature. Alfalfa varieties are classified into three types, fall-dormant varieties (FD I-III), semi-dormant varieties(FD IV-VI), and non-dormant varieties(FD IIV-IX), according to plant regrowth after mowing in the autumn[1-3]. FD, which is one of the most important factors that influences plant adaptation, has a dramatic impact on the production of alfalfa[2].Therefore, FD is considered the primary evaluation index of alfalfa varieties in North America [4]. In this region, fall-dormant alfalfa is adapted to the cold climate, grows slowly after cutting in autumn, and has low yield, but exhibits strong cold-hardiness and overwintering ability, which non-dormant alfalfa does not possess[5]. The FD of fall-dormant alfalfa is stronger than that of semi-dormant and non-dormant alfalfa, and that of non-dormant alfalfa is the weakest [6-8]. Recently, it has been reported that environmental factors, such as photoperiod and temperature, regulate FD[9, 10]. Furthermore, photoperiod is considered a limiting FD-inducing factor, and complex interactions associated with daylength and temperature are known to occur[6].

Deeper insights into the molecular mechanisms underlying FD have been recently reported. Douglas identified genomic regions that control the FD of alfalfa, although the specific nucleotide sequences of these genomic regions remain unknown[11]. Wang et al. found that PHYA/PHYB mRNA content increases with the shortening of daylength and reduction in temperature[8]. Fan et al. identified 583 dormancy-related microRNAs via high-throughput sequencing, of which 28 microRNAs were found to play key roles in regulating dormancy[12]. Zhang et al. detected 2,064 differentially expressed genes in alfalfa leaves between dormant and non-dormant periods by transcriptome sequencing[13]. In addition, the reduction of indoleacetic acid (IAA), zeatin riboside (ZR), and gibberellins (GA3) in leaves of fall-dormant varieties is more significant than that of semi-dormant and non-dormant varieties. Further, increased levels of abscisic acid (ABA) were observed in semi-dormant and non-dormant varieties[14]. However, to the best of our knowledge, the genes, specifically those related to photoperiod and temperature, that regulate FD in alfalfa remain unknown.

In the present study, the candidate genes that regulate alfalfa FD were identified by comparing the leaf transcriptomes of Maverick (fall-dormant alfalfa) and CUF101 (non-fall-dormant alfalfa) at FD and non-dormant periods, the results of the bioinformatics analysis, and the data for gene function from previous studies. Finally, three genes regulating fall dormancy were identified by analyzing the changes in mRNA abundance of key candidate genes from April to October and by performing correlation analysis between their mRNA abundance and daylength and temperature. The present study establishes the foundation for highly targeted study of molecular mechanisms underlying FD.

Materials and methods

Plant materials and growth condition

Seeds of Maverick (FD I) and CUF101 (FD IX) of alfalfa standard varieties were introduced from the United States and planted in sandy loam soil at the Experimental Station of Henan Agricultural University (34°19 N, 113°35 E) by hand on October 1,2009, with 0.6m spacing between rows. The plants were irrigated regularly during drought, but they were not fertilized during growth. Weeds were controlled by hand or hoeing, and insects were controlled as required by hand. Leaf samples of the fall-dormant alfalfa variety collected in May and September are labeled DM and DS, respectively, and the samples of the non-dormant alfalfa variety collected in September are labeled NDS thereafter[13].

The plants were cut, then let grow for 14 days, the leaf samples were collected, which was done once a month. Three replicates of mature fresh leaves of both varieties were collected from the top of the plant between 08:00 and 09:00 on April 12, May 10, June 12, July 8, August 6, September 2, and October 18 in2011,and on April15, May11, June13, July9, August7, September5, and October16 in 2016.The samples were immediately frozen in liquid nitrogen, and stored at -80°C until use for RNA extraction and assessment of mRNA abundance of genes. Meanwhile, height of 30 randomly chosen plants was recorded to calculate the mean height for each plot. Leaf area was calculated based on the length and width of ten randomly chosen leaves of each plant [15]. Day length and temperature of each sampling day were also recorded.

Maverick and CUF101 grown under artificial growth conditions in a green house were exposed to different illumination times or different temperature. One group of plants including three plants of each variety was cultured at 24°C and a daily illumination (3000lux) of 8h, 12h, or 16h. Another group of three plants of each variety was grown under 3000 lux illumination for 48h and a temperature of 16°C, 24°C,or 32°C. Leaf samples of each variety from each treatment were collected and immediately frozen in liquid nitrogen and subsequently stored at -80°C until use for RNA extraction and mRNA abundance assessment.

Significant differences were analyzed by one-way ANOVA using SPSS 19.0 (IBM Corp., Armonk, NY, USA). Daylength and temperature of each sampling day were also recorded. Correlation between plant growth rate or mRNA abundance of genes of the two varieties and daylength or temperature was analyzed using double factor correlation analysis and one-sided t-test, respectively, both available in SPSS 19.0 (IBM Corp., USA). Correlation between mRNA abundance of genes in the two varieties and daylength or temperature was also analyzed. A linear chart of the obtained data was constructed using GraphPad Prism 5 software (GraphPad Software, Inc., San Diego, CA, USA).

Screening and bioinformatic analysis of candidate genes

After these sequencing reads were trimmed, paired-end reads were assembled by Trinity [16] as a single dataset (reference transcriptome) that was then annotated using BLASTX. The transcript was deposited in DDBJ/EMBL/GenBank under the accession number GAFF00000000. Differentially expressed genes (DEGs) were identified using Simbiot® platform, and the accuracy of the results was tested[13]. DEGs with absolute foldchange ≥4 and adjusted p-value ≤ 0.05 for DM vs. DS and DS vs. NDS were obtained using R software[17], and a custom script was implemented to identify common DEGs(co-DEGs) and unique DEGs for DM vs. DS and DS vs. NDS. Next, the total DEGs in DM vs. DS and DS vs. NDS were screened for unique genes of DS vs. NDS and unique genes of DM vs. DS, respectively. Some of the unique DEGs were also classified as co-DEGs based on the fold-change value and adjusted p-value. In addition, DEGs associated with drought resistance, insect resistance, and disease resistance, as well as DEGs with the same expression trend in DM vs. DS and DS vs. NDS, were excluded. The co-DEGs, unique DEGs of DM vs. DS, and unique DEGs of DS vs. NDS were confirmed and collectively designated ACDEGs.

ACDEGs were imported into Blast2GO for gene ontology (GO) annotation. The sequences of these genes were submitted to the Kyoto Encyclopedia of Genes and Genomes (KEGG) Automatic Annotation Server (KAAS; http://www.genome.jp/kaas-bin/kaas_main; Version 1.67x) to retrieve KEGG orthology (KO) assignments and KEGG pathways using single-directional best hit assignment method. The function of each ACDEG was determined by submitting the amino acid sequences of the ACDEGs to the UniProt database (http://www.uniprot.org/uniprot/) and consulting the literature. The ACDEGs were classified according to their common functions. Finally, key candidate genes involved in the regulation of alfalfa FD were identified based on foldchange in their expression, p-value, and their function in plant growth and development.

Assessment of mRNA abundance of key genes and correlation analysis between mRNA abundance of genes and daylength or temperature

The mRNA abundance of 44 candidate genes in Maverick leaves (three samples over a 7-month period, resulting in a total of 21 samples) was detected and analyzed from April to October in 2011. Key genes among the 44 candidate genes were selected based on whether their transcription was regulated by photoperiod or temperature and their role in plant growth and development. The mRNA abundance of six key genes in CUF101 leaves(3 samples per month over a 7-month period, resulting in a total of 21 samples) was detected and analyzed again in 2011. The mRNA levels of these genes in Maverick and CUF101 leaves (3 samples ×7months × 2 varieties = 42 samples) were detected and analyzed again in 2016.

The mRNA abundance of the six key genes in Maverick and CUF101 leaves grown under artificial culture conditions (3 samples × 3 treatments × 2 varieties × 2 factors = 36 samples) was detected and analyzed.

Total RNA in all samples was extracted in strict accordance with the TRIzol method (Invitrogen, Carlsbad, CA, USA). Reverse transcription was performed using a RT kit following the manufacturer’s instructions (Takara Bio, Inc.). RNA content and quality were analyzed in a Nano2000 Ultramicro spectrophotometer (ThermoFisher Scientific, Waltham, MA, USA), and each RNA sample was adjusted to the same concentration.

Quantitative reverse transcription polymerase chain reaction (qRT-PCR) primers for each gene were designed using Primer5.0 software (Premier Biosoft International, Palo Alto, CA, USA) following the established principles of qRT-PCR primer design (S1 Table). mRNA abundance of each gene in all samples was detected using the Roche SYBR Green fluorescent dye method and a Roche Cycler9.0 fluorescent PCR instrument (Roche Diagnostics GmbH, Manheim, Germany). The data were normalized against the expression of the housekeeping gene GAPDH, and the expression was calculated relative to a calibrator samples using the formula 2−ΔΔCt; the values were presented as mean±standard deviation[18]. Significant differences were assessed by one-way ANOVA, and the correlation of the mRNA abundance of each gene with daylength or temperature was analyzed using double factor correlation analysis and one-sided t-test, respectively, implemented in SPSS 19.0 (IBM Corp., USA). The obtained qRT-PCR data were used to construct a linear chart in GraphPad prism 5 (GraphPad Software, Inc., USA).

Results

Analysis of DEGs

The DEGs of DM vs. DS and DS vs. NDS were listed in S1 and S2 Files. A total of 538 genes were identified as DEGs in DM vs. DS, of which 122 were annotated and 416 were not annotated. Further, the expression of 337 DEGs was upregulated and201 DEGs were downregulated in DM vs. DS(S3 File). Similarly, of the 836 DEGs identified in DS vs. NDS,156 were annotated and 680 were not annotated; 545 DEGs were upregulated and 291 DEGs were downregulated (S4 File). 86 genes were co-DEGs of DM vs. DS and DS vs. NDS(S5 File); among these, 78 genes were un-annotated and 8 were annotated (S5 File). After a series of analyses and screening, in total, 489 co-DEGs of DM vs. DS and DS vs. NDS were obtained; of these, 405 genes were un-annotated and 84 genes were annotated. Further, 218genes were unique to DS vs. NDS, with163 genes being un-annotated and 55 genes annotated. Of the354genes unique to DM vs. DS, 290 genes were un-annotated and 64 genes were annotated. Overall, ACDEGs comprised 1,069 genes (S6 File).

GO enrichment analysis of ACDEGs

The GO enrichment analysis classified most ACDEGs into three significant, broad GO categories: “biological process”, “molecular function” and “cell component”. Within the category “biological process”, the DEGs were assigned to the terms: metabolic process (153 genes), cellular process (113 genes), response to stimulus (37 genes), single-organism process (55 genes), localization (36 genes), cellular component organization or biogenesis (26 genes), signaling (21 genes), developmental process (13 genes), and growth (6 genes). The DEGs assigned to “molecular function” category were enriched in cell (84 genes), organelle (65 genes), macromolecular complex (26 genes), membrane (8 genes), and extracellular region (6 genes). Finally, in the “cell component,” the genes were functionally assigned to catalytic activity (139 genes), binding (95 genes), transporter activity (15 genes), molecular structure (15 genes), and nucleic acid binding transcription factor activity (9 genes) (Fig 1).

Fig 1

Gene ontology functional classification of key candidate genes.

The functional assignments with biological processes, molecular functions, and cellular components are shown based on the number of proteins and the converted corresponding proportion.

Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis of ACDEGs

ACDEGs were significantly enriched in 33 pathways (p ≤0.05) according to the KEGG enrichment analysis. Ribosome pathway and metabolic pathways were the two main pathways, and DEGs were significantly enriched in carbohydrate metabolic pathways, such as the starch and sucrose metabolism pathway and the pentose and glucuronate interconversion pathway. The biosynthesis of secondary metabolites, such as cyanoamino acid metabolism, phenylpropanoid biosynthesis, amino sugar and nucleotide sugar metabolism, phenylalanine metabolism, and sesquiterpenoid biosynthesis, was also enriched. Genes involved in protein processing and endoplasmic reticulum pathway coded mainly for heat shock proteins. The RNA transport pathway, ubiquitin-mediated proteolysis pathway, tryptophan metabolism pathway, plant circadian rhythm, and plant hormone signal transduction were also significantly enriched (S7 File).

Based on the results of the KEGG and GO analysis of DEGs and the results of previous studies, 44 key candidate genes may be involved in the regulation of FD (Table 1).

Table 1

Key candidate genes involved in regulating alfalfa fall dormancy.

Classification	Number	Tags	Genes	Function [reference]	Expression trend of DM vs. DS	Expression trend of DS vs. NDS
Genes involved in carbohydrate metabolism and transport	1	comp1260_c0	MGDG synthase 3	Synthesis of photosynthetic membranes and chloroplast envelope [19]		Downregulated(-3.53)
	2	comp59094_c0	Beta-D-glucosidase	Catalyzes the hydrolysis of terminal non-reducing residues in beta-D-glucosides with a release of glucose[20]	Downregulated(-2.10)	Upregulated(3.05)
	3	comp25768_c0	UDP-sugar pyrophosphorylase	Catalyzes the conversion of various monosaccharide 1-phosphates to respective UDP-sugars in the salvage pathway[21]	Upregulated(2.14)	Downregulated(-3.19)
	4	comp53413_c0	Ribulose bisphosphate carboxylase/oxygenase (Rubisco) activase	Activates Rubisco	Downregulated(-2.87)	Upregulated(1.01)
Genes regulated by photoperiod and light	5	comp46970_c0	Cryptochrome-2(CRY-2)	Inhibits hypocotyl elongation [22] and stem and root growth[23]	Upregulated(3.50)
	6	comp1262167_c0	Protein FAR-RED IMPAIRED RESPONSE 1(FAR1)	Positively regulates chlorophyll biosynthesis via the activation of HEMB1 gene expression[24]; activates transcription; elongates hypocotyls and reduces expansion of cotyledons[25]		Upregulated(2.23)
	7	comp38709_c0	GATA transcription factor 12(GATA12)	Activates transcription, regulates some light-responsive genes and circadian rhythm[26–28]	Upregulated(2.40)	Downregulated(-2.41)
	8	comp12163_c0	Serine protease SPPA, chloroplastic	the light-dependent degradation of antenna and photosystem II in chloroplasts[29, 30];	Downregulated(-2.22)	Upregulated(2.23)
	9	comp46405_c0	Protein HEADING DATE 3A	Regulates dormancy [31]	Downregulated(-1.72)	Upregulated(4.38)
	10	comp56985_c0	Granule-bound starch synthase 1(GBSS1)	Regulated by photoperiod, may be accompanied by abolition of expression of starch synthesis genes [32]	Downregulated(-2.69)	Upregulated(5.74)
Phytohormone	11	comp36708_c0	Probable indole-3-pyruvate monooxygenase YUCCA3	Participates in indoleacetic acid (IAA) synthesis[33]	Upregulated(3.40)	Downregulated(-1.35)
	12	comp343499_c0	Methylesterase17	Hydrolyzes conjugates of IAA	Downregulated(-1.72)	Upregulated(3.38)
	13	comp522282_c0	IAA-amino acid hydrolase ILR1-like 1	Hydrolyzes conjugates of IAA[34]		Upregulated(2.36)
	14	comp15714_c0	Auxin response factor 6(ARF6)	Activates transcription, participates in transcriptional regulation of a variety of biological processes related to growth and development[35, 36]		Downregulated(-2.74)
	15	comp57456_c0	ABA receptor PYL8	Positively regulates abscisic acid (ABA) signaling pathway[37–40]	Upregulated(2.66)
	16	comp61323_c0	Protein early responsive to dehydration 15	Negatively regulates ABA responses and mediates stress-related ABA signaling[41, 42]	Upregulated(2.05)
Transportation	17	comp46282_c0	ATP-binding cassette (ABC) transporter C family member 7	Involved in the response to biotic stress		Upregulated(2.36)
	18	comp10285_c0	ABC transporter B family member 11	Essential regulator of plant growth	Upregulated(2.55)
	19	comp70285_c0	Monosaccharide-sensing protein 3	Involved in direct growth and differentiation or flexible response to changing environmental conditions (UniProt: Q9SD00)	Upregulated(2.13)
	20	comp40098_c0	Probable ion channel SYM8	A calcium channel[43]	Upregulated(3.15)
	21	comp343967_c0	Putative phospholipid-transporting ATPase 7	Transports phospholipids (UniProt: Q9LVK9)	Upregulated(3.15)	Downregulated(-3.32)
	22	comp928260_c0	Ras-related protein RABH1d	Involved in protein transport(UniProt: Q9SID8)	Downregulated(-0.83)	Upregulated(2.23)
	23	comp1120261_c0	Oligopeptide transporter 6	Translocates tetra- and pentapeptides[44, 45]	Upregulated(2.40)	Downregulated(-2.41)
Cell cycle, division, and growth	24	comp33056_c0	Cyclin-D5-2	Mediates the linking of extracellular and developmental signals to the cell cycle[46]		Downregulated(-3.60)
	25	comp67791_c1	Anaphase-promoting complex (APC)subunit 1	Marks target cell cycle proteins for degradation by 26s proteasome; involved in eukaryotic cell reproduction[47–49]	Downregulated(-1.27)	Upregulated(5.11)
	26	comp686013_c0	Protein S-acyltransferase 24 (PAT24)	Involved in cell growth regulation [50]; affects root hair formation and growth [51, 52]	Upregulated(2.20)
	27	comp32326_c0	Monocopper oxidase-like protein SKU5	Possibly involved in directional growth processes and cell wall expansion[53]	Upregulated(2.20)
	28	comp672512_c0	Prohibitin-3, mitochondrial	Loss of function of the homologs AtPHB3 causes mitochondrial welling, decreases meristematic cell production, increases cell division time, and reduces cell expansion rates, leading to severe growth retardation[54]	Downregulated(-0.83)	Upregulated(2.36)
	29	comp887354_c0	Kinesin-4	Involved in cell cycle[55] and cell division[56–61].	Downregulated(-2.92)	Upregulated(2.23)
	30	comp37666_c0	Structural maintenance of chromosomes protein 3(SMC3)	Essential protein for plant viability, required for cell division during embryogenesis, increased expression accelerates cell cycle [62, 63]		Upregulated(2.23)
Genes regulating transcription	31	comp424089_c0	Probable histone H2A.2	Plays a central role in transcription regulation, DNA repair, DNA replication, and chromosomal stability. Involved in gene regulation[64]		Upregulated(2.36)
	32	comp57595_c0	Histone H3.3		Upregulated(6.23)	Downregulated(-4.08)
	33	comp391402_c0	MADS-box transcription factor PHERES 1	Regulates dormancy[31, 65]	Downregulated(-3.20)	Upregulated(3.12)
	34	comp50413_c1	Squamosa promoter-binding protein 1(SBP1)	Activates transcription, and involved in leaf development, vegetative phase change, etc. [66, 67]	Upregulated(3.97)
Ubiquitination	35	comp53272_c0	Probable mediator of RNA polymerase II transcription subunit 37c	Coactivator regulates transcription of nearly all RNA polymerase II-dependent genes [UniProt: Q9LHA8—MD37C_ARATH]		Downregulated(-6.76)
	36	comp538610_c0	BTB/POZ domain-containing protein At1g67900	Mediates transcriptional repression [UniProt: Q9C9V6—Y1790_ARATH]	Upregulated(2.40)
	37	comp933338_c0	BTB/POZ domain-containing protein At3g19850	Mediates transcriptional repression [UniProt: Q9LT24—Y3985_ARATH]	Upregulated(2.40)
	38	comp29926_c0	F-box only protein 6	Substrate-recognition component of some SCF-type E3ubiquitin-ligase complexes, participates in regulation of auxin-mediated signaling pathway, leaf vascular tissue pattern formation, and simple leaf morphogenesis[68, 69]		Upregulated(5.50)
	39	comp655910_c0	UBX domain-containing protein 7	Acts in many cellular events such as ubiquitin-dependent degradation and membrane fusion[70]	Upregulated(2.87)	Downregulated(-3.00)
	40	comp909276_c0	Putative E3 ubiquitin-protein ligase LIN	Catalyzes polyubiquitination with ubiquitin-conjugating enzyme E2 UBC8 in vitro, involved in plant C/N response and early steps of the plant defense signaling pathway[71]	Upregulated(2.40)
	41	comp395328_c0	E3 ubiquitin-protein ligase UPL6	Mediates ubiquitination and subsequent proteasomal degradation of target proteins [UniProt:Q8RWB8—UPL6_ARATH]	Downregulated(-2.96)	Upregulated(1.68)
Receptor kinases	42	comp70176_c0	Wall-associated receptor kinase(WAK)5; Wall-associated receptor kinase-like (WAKL)2	Involved in cell expansion, elongation[72, 73]	Upregulated(2.36)
	43	comp41596_c0	WAK3			Downregulated(-3.00)
	44	comp403595_c0	Serine/threonine-protein kinase Nek1	Involved in sensing and repair of DNA strand breaks at the G1-S and G2-M transitions[74]	Upregulated(2.40)	Downregulated(-2.41)

DM and DS: Leaf samples of the fall-dormant alfalfa variety collected in May and September, respectively; NDS: leaf samples of the non-dormant alfalfa variety collected in September.

Daylength, temperature, leaf area, and plant height of the Maverick variety from April to October

Daylength and temperature first increased and then decreased from April to October: the longest day was in June, and the highest temperature was measured in July (Fig 2).

Fig 2

Daylength and temperature from April to October.

The plant height of Maverick did not differ significantly between April and August, but it decreased significantly in September and October as compared with that from April to August. The plant height of CUF101 did not differ significantly between July and October, and it was significantly higher than that of Maverick in the same period (Fig 3). The leaf area of Maverick was smaller than that of CUF101 from April to September, and the difference in leaf area between the two varieties reached its maximum in September. There was no significant difference in the leaf area of CUF101 from June to September, whereas the leaf area of Maverick decreased from June to September. In addition, the leaf area of Maverick in August and September was significantly smaller than that in June and July (Fig 4).

Fig 3

Plant height of Maverick and CUF101varieties from April to October.

The difference in plant height in Maverick is marked with *. The difference in plant height between Maverick and CUF101 is marked with #. (*, #, p < 0.05; **, ## p< 0.01).

Fig 4

Leaf area of Maverick and CUF101 varieties from April to October.

* and **indicate significant difference at p < 0.05 and p < 0.01, respectively.

mRNA abundance of four ACDEGs from April to October in the two alfalfa varieties

Among the 44 DEGs, mRNA abundance of six DEGs was significantly correlated with daylength and temperature in 2011. However, the comparison of the data between 2011 and 2016 revealed that the change in mRNA abundance of one gene in both varieties differed between the two sampling years and the change in mRNA levels of another gene followed the same trend in both varieties and thus was not significantly different between the two varieties. Therefore, these two genes were excluded from further analysis. The mRNA profiles of IAA-amino acid hydrolase ILR1-like 1(ILR-like 1), ABA receptor PYL8(PYL8), monogalactosyldiacylglycerol (MGDG) synthase-3(MGDGS-3), and Ribulose bisphosphate carboxylase/oxygenase (Rubisco) activase are discussed herein. The mRNA profiles of the four genes exhibited the same regular, specific trends from April to October in 2011 and 2016, with significant differences detected between Maverick and CUF101 (Figs 5–8). The mRNA data of other genes are not shown.

Fig 5

mRNA abundance of indoleacetic acid (IAA)-amino acid hydrolase ILR1-like 1 in Maverick and CUF101 varieties from April to October in 2011 and 2016.

* and ** indicate significant difference between August-September and June-July at p<0.05 and p<0.01, respectively.

Fig 8

mRNA abundance of Rubisco activase in Maverick and CUF101varieties from April to October in 2011 and 2016.

** indicates significant difference between September, October, July, and August at p<0.01.

The change in mRNA abundance of ILR1-like 1 and the PYL8 followed the opposite trends in the autumn. Thus, the mRNA abundance of ILR1-like 1 decreased in Maverick but increased in CUF101 from June to September, whereas the mRNA abundance of the PYL8 increased in Maverick but decreased in CUF101 in the same period. In Maverick, the mRNA abundance of ILR1-like 1 in August to October was significantly lower than that in June and July. In CUF101, the abundance of the same mRNA was significantly higher in July to September than it was in June. The mRNA abundance of the PYL8 was significantly higher in September and October in Maverick but significantly lower in CUF101 as compared to those in August(Figs 5 and 6).

Fig 6

mRNA abundance of the abscisic acid receptor PYL8 in Maverick and CUF101 varieties from April to October in 2011 and 2016.

** indicates significant difference between September-October and August at p<0.01.

The mRNA abundance of MGDGS-3 in the two varieties first increased and then decreased starting from June (in Maverick) and August (in CUF101). In Maverick, the mRNA abundance of MGDGS-3 was significantly lower from July to October than it was in June, but in CUF101, it was significantly higher from July to September compared with that in June. The mRNA abundance of MGDGS-3 from July to September was significantly lower in Maverick compared to that in CUF101 (Fig 7). Similarly, the abundance of the Rubisco activase mRNA increased initially in both varieties, which was followed by a decrease from July onwards. In both alfalfa varieties, its mRNA abundance in September and October was significantly lower than that in July and August(Fig 8).

Fig 7

mRNA abundance of monogalactosyldiacylglycerol synthase-3 in Maverick and CUF101 varieties from April to October in 2011 and 2016.

* and ** indicate significant difference between July, August, September, and June at p<0.05 and p<0.01, respectively.

mRNA profiles of four DEGs in leaves of the two alfalfa varieties under artificial growth conditions

The change in mRNA abundance of ILR1-like 1 in Maverick followed the opposite trend to that in CUF101 with increasing illumination from 8h to 16h. Thus, its mRNA abundance in Maverick gradually increased with increasing illumination from 8h to 16h and was significantly less than that in CUF101 (Fig 9A).

Fig 9

mRNA abundance of four differentially expressed genes (DEGs) in Maverick and CUF101 varieties grown under 8h to 16h daylength and artificial growth conditions.

Significant difference in the mRNA content of DEGs at 8h,12h, and 16h of illumination in each variety is marked with *. Significant difference in mRNA content between the two varieties is indicated with # (*, #, p<0.05; **, ##,p<0.01). Error bars indicate standard deviation (SD).

The mRNA abundance of the PYL8 in the two varieties gradually decreased with increasing illumination from 8h to 16h and was significantly greater than that in CUF101 (Fig 9B).

As the illumination increased from 8 h to 16 h, the mRNA abundance of MGDGS-3 increased and that of the rubisco activase mRNA was not significantly altered. The levels of both genes of mRNA in Maverick were significantly lower than those in CUF101 (Fig 9C and 9D).

The changes in mRNA abundance of the ILR1-like 1 and the PYL8 in the same variety followed the same trend as temperature increased from 16°C to 32°C.Thus, their mRNA abundance gradually increased from 16°C to 32°C in Maverick, whereas in CUF101, the abundance of the two mRNAs at 24°C was greater than that at 16°C and 32°C (Fig 10A and 10B). Similarly, the abundance of the Rubisco activase mRNA in the two varieties gradually increased from 16°C to 32°C, but it reached higher levels in Maverick than in CUF101 (Fig 10C). In contrast, the levels of the MGDGS-3 mRNA gradually increased in Maverick but decreased in CUF101 with increasing temperature from 16°C to 32°C (Fig 10D).

Fig 10

mRNA abundance of four differentially expressed genes (DEGs) in Maverick and CUF101 varieties grown under artificial growth conditions at different temperatures (from 32°C to 16°C).

Significant difference in mRNA content of DEGs between 16°C, 24°C and 32°C for each variety is marked with *. Significant difference in mRNA content between the two varieties is marked with # (*, #, p<0.05; **, ##, p<0.01). Error bars indicate standard deviation (SD).

Correlation analysis of plant growth rate, mRNA abundance of genes, and daylength or temperature under natural conditions

The growth rate of Maverick was significantly correlated with daylength, but not with temperature, whereas the growth rate of CUF101 was not significantly correlated with either daylength or temperature (Table 2).

Table 2

Correlation of the growth rate with daylength and temperature in Maverick and CUF101 varieties.

	Maverick growth rate		CUF101 growth rate
	Pearson correlation	Significance	Pearson correlation	Significance
Daylength	0.858**	0.007	0.593	0.080
Temperature	0.107	0.410	0.445	0.158

**represent significant correlation at 0.01 level.

The mRNA abundance of ILR1-like 1 in Maverick was significantly positively correlated with daylength and moderately positively correlated with temperature in both varieties in the two experimental years. In contrast, the abundance of the ILR1-like 1 mRNA in CUF101 had no significant correlation with daylength and temperature (Tables 3 and 4).

Table 3

Correlation between mRNA abundance in the leaves of Maverick and CUF101 varieties and daylength in 2011 and 2016.

	Daylength
	Maverick				CUF101
	2011		2012		2011		2012
	Pearson correlation	Significance	Pearson correlation	Significance	Pearson correlation	Significance	Pearson correlation	Significance
IAA-amino acid hydrolase ILR1-like 1	0.975**	0.000	0.700*	0.040	0.232	0.308	-0.137	0.385
Abscisic acid receptor PYL8	-0.664	0.052	-0.92**	0.001	-0.184	0.346	-0.035	0.47
MGDG synthase 3	0.610	0.073	0.318	0.244	0.012	0.49	0.054	0.454
Rubisco activase	0.294	0.261	0.438	0.163	0.579	0.086	0.604	0.075

IAA: indoleacetic acid; MGDG: monogalactosyldiacylglycerol.

*and **represent significant correlation at 0.05 and 0.01 level, respectively. IAA: indoleacetic acid; MGDG: monogalactosyldiacylglycerol.

Table 4

Correlation between mRNA abundance in the leaves of Maverick and CUF101 varieties and temperature in 2011 and 2016.

	Temperature
	Maverick				CUF101
	2011		2012		2011		2012
	Pearson correlation	Significance	Pearson correlation	Significance	Pearson correlation	Significance	Pearson correlation	Significance
IAA-amino acid hydrolase ILR1-like 1	0.491	0.131	0.627	0.066	0.545	0.103	0.599	0.078
Abscisic acid receptor PYL8	-0.726*	0.032	-0.430	0.167	0.192	0.34	0.654	0.055
MGDG synthase 3	0.665	0.051	0.713*	0.036	0.506	0.123	0.714*	0.035
Rubisco activase	0.612	0.072	0.807*	0.014	0.899**	0.003	0.926**	0.0015

IAA: indoleacetic acid; MGDG: monogalactosyldiacylglycerol.

* and ** represent significant correlation at 0.05 and 0.01 level, respectively.IAA: indoleacetic acid; MGDG: monogalactosyldiacylglycerol.

The mRNA abundance of the PYL8 was significantly negatively correlated with daylength and moderately negatively correlated with temperature in Maverick, but had no significant correlation with the two parameters in CUF101 (Tables 3 and 4).

In both varieties, the abundance of the MGDGS-3 mRNA was significantly positively correlated with temperature, and that of the Rubisco activase was moderately positively correlated with temperature; neither of the two was correlated with daylength (Tables 3 and 4).

Correlation analysis of mRNA abundance of genes and daylength or temperature under artificial growth conditions

The correlation of the mRNA abundance of ILR1-like 1 and that of Rubisco activase was in the two varieties with illumination time was opposite. The mRNA abundance of ILR1-like 1 and Rubisco activase in Maverick was significantly positively and negatively correlated with illumination time, respectively. Correlation of the mRNA abundance of PYL8 and MGDGS-3 in Maverick and CUF101 with illumination time was same. The mRNA abundance of PYL8 was significantly negatively correlated with illumination time. The mRNA abundance of MGDGS-3 was significantly positively correlated with illumination time. (Table 5).

Table 5

Correlation between mRNA abundance in the leaves of Maverick and CUF101 varieties and illumination time under artificial growth conditions.

	illumination time
	Maverick		CUF101
	Pearson correlation	Significance	Pearson correlation	Significance
IAA-amino acid hydrolase ILR1-like 1	0.812**	0.008	-0.984**	0
Abscisic acid receptor PYL8	-0.987**	0	-0.981**	0
MGDG synthase 3	0.767*	0.016	0.666*	0.05
Rubisco activase	-0.756*	0.019	0.617	0.077

* and ** represent significant correlation at 0.05 and 0.01 level, respectively.

IAA: indoleacetic acid; MGDG: monogalactosyldiacylglycerol.

Correlation of the mRNA abundance of MGDGS-3 in Maverick and CUF101 with temperature was opposite, the mRNA abundance of MGDGS-3in Maverick was significantly positively correlated with temperature. Correlation of the mRNA abundance of other three genes in Maverick and CUF101 with temperature was same, their mRNA abundance was significantly positively correlated with temperature (Table 6).

Table 6

Correlation between mRNA abundance in the leaves of Maverick and CUF101 varieties and temperature under artificial growth conditions.

	Temperature
	Maverick		CUF101
	Pearson correlation	Significance	Pearson correlation	Significance
IAA-amino acid hydrolase ILR1-like 1	0.997**	0	0.771*	0.015
Abscisic acid receptor PYL8	0.886**	0.001	0.062	0.874
MGDG synthase 3	0.947**	0	-0.842**	0.004
Rubisco activase	0.971**	0	0.996**	0

* and ** represent significant correlation at 0.05 and 0.01 level, respectively.

IAA: indoleacetic acid; MGDG: monogalactosyldiacylglycerol.

Thus, in Maverick, the correlation of the mRNA abundance of the PYL8 and Rubisco activase with illumination time and temperature followed the opposite trends, whereas that of the mRNA abundance of the ILR1-like 1 and MGDGS-3with illumination time and temperature was the same (Tables 5 and 6). In CUF101, the correlation of the mRNA abundance of MGDGS-3 with illumination time and temperature was the opposite, but that of the ILR1-like 1,PYL8, and Rubisco activase was the same (Tables 5 and 6).

Discussion

At the whole-plant level, communication between various organs in a plant is involved in the coordination of growth processes at different organizational levels. Thus, growth of individual organs is regulated by long-distance communication from other organs[75]. FD is the overall growth performance of the whole plant during special phases and environment conditions. Leaf plays a key role in the growth of the whole plant through photosynthesis, respiration, transpiration, and other basic functions involved in plant growth, and changes in those functions affect the growth of the whole alfalfa plant. Thus, the leaf plays an important role in alfalfa FD. The comparison of the leaf area in different months between fall-dormant and non-dormant alfalfa indicated that this parameter may be used as an important index of differences between various alfalfa varieties in terms of FD.

In the present study, beside a few annotated DEGs similar to dormancy-regulation genes, the GO and KEGG analysis revealed that ACDEGs play roles in basic biological processes and pathways involved in the regulation of plant growth and development, such as in the response to macro environmental factors (e.g., light, photoperiod, temperature), the photosynthesis and respiration of leaves, and the regulation of phytohormones. Therefore, these genes are expected to be involved in the regulation of FD.

Temperature and especially photoperiod are important environmental factors that regulate FD[6]. In the present study, the growth of Maverick was induced by daylength, whereas FD of Maverick was accompanied by the shortening of daylength. These findings confirm that photoperiod is a key environmental factor regulating FD in alfalfa and thus corroborate previous study[6]. Under artificial growth conditions, the expressions of the ILR1-like 1 and Rubisco activase in Maverick and CUF101 were oppositely affected by illumination, and expression of MGDGS-3 in Maverick and CUF101 was oppositely affected by temperature. Comparison of the correlations between mRNA abundance of the four genes and the two environmental factors (daylength and temperature) under natural and artificial growth conditions showed that the expression of ILR1-like 1 is upregulated by daylength and temperature, whereas in CUF101, it is mainly affected by temperature. The expression of the PYL8 in Maverick is affected by daylength, and that of MGDGS-3 is affected by daylength and temperature in Maverick or by daylength in CUF101. Temperature was the main factor affecting the levels of Rubisco activase in the two cultivars. Therefore, these results suggest that the responses of Maverick and CUF101 to changes of daylength and temperature are different.

Leaf is an important organ involved in phytohormone synthesis. Phytohormones synchronize developmental processes by adjusting plant growth in response to intrinsic and environmental cues [76, 77]. Thus, IAA and ABA play key roles in plant growth and development[78], indole-3-pyruvate monooxygenase YUCCA3, methylesterase 17, and ILR1-like 1 in ACDEGs participate in IAA synthesis and activation[33-36], and PYL8 and early responsive to dehydration 15 in ACDEGs are associated with the ABA signaling pathway[42, 79]. However, the present findings suggest that downregulation of ILR1-like 1 transcription and upregulation of PYL8 transcription were most likely the only factors involved in alfalfa FD in response to shortened photoperiod and a drop in temperature.

Auxins, such as IAA, possess various functions including induction of cell) elongation and cell division, which are important for plant growth and development[80]. ILR1-like 1 releases IAA by hydrolyzing specific amino acid conjugates of the plant growth regulator IAA [34]. The GO analysis showed that ILR1-like 1 was involved in metabolic processes resulting in cell growth(GO:0008152). Decreased expression of ILR1-like 1 in Maverick observed after reduction in daylength and drop in temperature as the summer transitioned into autumn indicates that its expression was regulated by daylength and temperature. However, in CUF101, the expression of ILR1-like 1 increased with shorter daylength in the autumn. Our previous study showed that, compared with semi- and non-dormant alfalfa, the decrease in IAA content in fall-dormant alfalfa was greater and more rapid in response to daylength shortening under artificial growth conditions[14]. In addition, plant height decreased in the fall-dormant alfalfa but remained constant in non-fall-dormant alfalfa from summer to autumn. Therefore, it is speculated that the reduction of ILR1-like 1 participates in alfalfa FD; specifically, its decrease leads to an increase in amino acid conjugates of IAA and a reduction of IAA levels with daylength shortening. Biological activity of an IAA conjugate is opposite from the function of the IAA itself[81]. A previous study showed that exogenous IAA-Ala treatment of tomato inhibits IAA-induced shoot growth and root initiation[82], which explains the observed increase in amino acid conjugates of IAA during IAA-inhibited plant growth.

ABA-mediated signaling plays a critical role in the responses of plants to environmental stresses[83]. In the present study, PYL8 was involved in plant hormone signal transduction (KO 04075). The expression of the PYL8 in Maverick was significantly negatively regulated by daylength and it increased with shortening of daylength. Furthermore, our previous study showed that, compared with semi- and non-dormant alfalfa, ABA content in fall-dormant alfalfa increased more rapidly and by a greater amount in response to daylength shortening under artificial growth conditions, reaching significantly higher levels compared to those in non-dormant alfalfa when the daylength was 13h or less[14]. In addition, plant height of the fall-dormant alfalfa decreased from summer to autumn, whereas that of the non-dormant alfalfa was not reduced. Previous studies demonstrated that ABA is associated with bud dormancy[84], inhibition of seed germination, and prevention of loss of seed dormancy[85, 86]. The PYL8, which is required for ABA-mediated responses (such as stomatal closure and germination inhibition), inhibits the activity of group-A protein phosphatases type 2C when activated by ABA, thus positively regulating the ABA signaling pathway[37-40]. Therefore, it is speculated that increased levels of the PYL8 are involved in alfalfa FD by enhancing the ABA signaling pathway.

Leaf photosynthesis plays a key role in plant growth. MGDG is the most abundant integral lipid in the thylakoid membrane and the photosystem II (PSII) complex[87, 88], which maintains both the linear electron transport process and the photostability of the PSII apparatus[89]. MGDGS-3was found to be involved in the glycolipid biosynthetic process(GO:0009247). The final step of the MGDG biosynthesis is catalyzed by the MGDG synthase[90, 91]. The mRNA abundance of MGDGS-3 in Maverick decreased with shorter illumination time and temperature drop, reaching levels that were significantly lower than those in CUF101 in the autumn; the plant height of Maverick was significantly lower than that of CUF101.Similarly, MGDG-deficient transgenic tobacco plant M18 exhibits retarded growth[89]. Therefore, reduction of MGDGS-3 levels is involved in alfalfa FD. Given that decreased MGDG content in Arabidopsis thaliana and tobacco have been associated with reduction in MGDG synthase levels[89, 92], the same is expected to occur in alfalfa. In addition, reduced MGDG levels reduce thylakoid membrane and the rate of photosynthesis[92, 93]. Therefore, it is suggested that the decrease in MGDG synthase participates in alfalfa FD by reducing leaf photosynthesis in response to temperature drop.

Rubisco(ribulose-1,5-bisphosphate carboxylase/oxygenase) is a key protein in plants. The change in its expression affects the photosynthesis and plant growth by altering the availability of N [94]. Rubisco can be activated by rubisco activase[95]. Previous studies demonstrated that the decrease in the expression of Rubisco activase may lead to the reduction of the photosynthetic rate and plant growth due to reduced activity of Rubisco. In addition, moderately high temperature was found to inhibit Rubisco activase-mediated activation of Rubisco[96]. Therefore, the reduction of Rubisco activase, which trigged by the fall in temperature, significantly reduces the activation of Rubisco or the light-saturated photosynthetic rate[97-99]. Considering that enhanced thermostability of Rubisco activase in Arabidopsis has been shown to improve CO2 assimilation rates and plant growth under heat stress[100, 101] and as our results showed that the expression of Rubisco activase was positively regulated by temperature, it was expected that the change in the expression of Rubisco activase would affect alfalfa growth. However, the changes in mRNA abundance of rubisco activase followed the same trend in Maverick and CUF101 and its abundance showed no difference between the two varieties, suggesting that rubisco activase is not involved in FD of alfalfa.

Conclusion

In the present study, 44 important candidate genes likely associated with alfalfa growth and FD were identified. These genes were mainly enriched in the following categories: transduction of light and photoperiod signals and leaf-derived signals (carbohydrates and phytohormones); the process of cell cycle, division, and growth; transcription factors, ubiquitination proteins; receptor kinases; and un-annotated genes. The present work demonstrates that the reduction of ILR1-like 1 and the increase of PYL8 and MGDGS-3 promote alfalfa FD in a response to changes in photoperiod or temperature.

qRT-PCR primers of 37 differentially expressed genes.

(DOCX)

Click here for additional data file.

Complete list of differentially expressed genes identified between samples of fall-dormant alfalfa varieties collected in May and September.

(CSV)

Click here for additional data file.

Complete list of differentially expressed genes identified between samples of the fall-dormant alfalfa variety and non-fall-dormant alfalfa variety collected in September.

(CSV)

Click here for additional data file.

Differentially expressed genes (fold change ≥ 4 and adjusted p-value ≤ 0.05)identified between fall-dormant alfalfa varieties in May and September.

(XLS)

Click here for additional data file.

Differentially expressed genes (fold change ≥ 4, adjusted p-value ≤ 0.05) identified between samples of fall-dormant alfalfa variety and non-fall-dormant alfalfa collected in September.

(XLS)

Click here for additional data file.

Common differentially expressed genes (fold change ≥ 4, adjusted p-value ≤ 0.05) in samples of fall-dormant alfalfa variety and non-fall-dormant alfalfa collected in September, and differentially expressed genes in samples of fall-dormant alfalfa variety collected in May and September.

(XLSX)

Click here for additional data file.

Candidate genes regulating alfalfa fall dormancy (ACDEGs) are identified after discarding the genes for drought resistance, insect resistance, etc.

The genes showing the same trends in the expression between fall-dormant alfalfa varieties in May and September were excluded.

(XLSX)

Click here for additional data file.

Pathways of candidate genes regulating alfalfa fall dormancy.

(DOCX)

Click here for additional data file.