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][2][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 nondormant alfalfa, and that of non-dormant alfalfa is the weakest [6][7][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]. a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 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 Simbiot1 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).

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).

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).
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).
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).

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 IAAamino 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 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).
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).

mRNA profiles of four DEGs in leaves of the two alfalfa varieties under artificial growth conditions
Maverick gradually increased with increasing illumination from 8h to 16h and was significantly less than that in CUF101 (Fig 9A).
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).
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).
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).
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).
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).
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][34][35][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][38][39][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][98][99]. Considering that enhanced thermostability of Rubisco activase in Arabidopsis has been shown to improve CO 2 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.
Supporting information S1 Table. qRT-PCR primers of 37 differentially expressed genes.