Study on differentially expressed genes related to defoliation traits in two alfalfa varieties based on RNA-Seq

Background Alfalfa (Medicago sativa) is a widely cultivated, essential commercial forage crop. The protein content in its leaves is the critical factor in determining the quality of alfalfa. Thus far, the understanding of the molecular mechanism of alfalfa defoliation traits remains unclear. The transcriptome database created by RNA-Seq is used to identify critical genes related to defoliation traits. Results In this study, we sequenced the transcriptomes of the Zhungeer variety (with easy leaf abscission) and WL319HQ variety (without easy leaf abscission). Among the identified 66,734 unigenes, 706 differentially expressed genes (DEGs) upregulated, and 392 unigenes downregulated in the Zhungeer vs WL319HQ leaf. KEGG pathway annotations showed that 8,414 unigenes were annotated to 87 pathways and contained 281 DEGs. Six DEGs belonging to the “Carotenoid biosynthesis”, “Plant hormone signal transduction” and “Circadian rhythm-plant” pathways involved in defoliation traits were identified and validated by RT-qPCR analyses. Conclusions This study used RNA-Seq to discover genes associated with defoliation traits between two alfalfa varieties. Our transcriptome data dramatically enriches alfalfa functional genomic studies. In addition, these data provide theoretical guidance for field production practice and genetic breeding, as well as references for future study of defoliation traits in alfalfa. Electronic supplementary material The online version of this article (10.1186/s12864-018-5180-1) contains supplementary material, which is available to authorized users.


Background
Alfalfa (Medicago sativa L.) is a member of the plant family Fabaceae. Alfalfa is widely distributed across the world and cultivated as a commercial forage crop due to its high protein, vitamin and mineral content [1,2]. Mature alfalfa leaves are the primary crude protein storage organs (approximately 260 g kg − 1 dry basis), and protein content is a crucial indicator of the quality of alfalfa hay [3]. However, the leaves fall off very easily due to external pressure during harvest [4], and the resulting low-quality hay cannot satisfy livestock nutritional requirements.
The main factors in leaf senescence are endogenous abscisic acid (ABA) and ethylene (ETH) [5,6]. ABA is a natural plant hormone identified in the 1960s, with roles in promoting plant organ shedding and seed germination [7,8]. ABA is present in higher plant organs and tissues and especially abundant in mature and aging tissues or dormant organs. When plants are under biotic or abiotic stress, such as pathogen infection, drought or salinity, leaf senescence is induced by increasing ABA content [9][10][11][12]. ETH is another early-discovered plant hormone with function in the regulation of plant growth, such as the promotion of fruit ripening, acceleration of organ aging and shedding [13]. ETH is synthesized by conversion from methionine with enough oxygen, which occurs in various plant tissues and organs [14]. ETH synthesis is regulated by a variety of growth signals with tissue specificity at different development stages, and the accumulation of ETH can also stimulate further ETH production [15]. The content of ETH is significantly higher in the processes of fruit ripening, seed germination, and leaf and flower aging off [16,17]. In addition, the synthesis of ETH is affected by a variety of external factors, such as mechanical damage, pathogen infection, low temperature, drought and salinity [18].
Previous studies showed that both ABA and ETH have the function of accelerating plant organ shedding. However, the physiological mechanisms are still unclear and remain controversial [19]. A study by Addicott [20] showed that endogenous ABA had positive effects on plant organ shedding, whereas field experiments failed, possibly because indoleacetic acid (IAA), gibberellin (GA) and cytokinin (CTK) in leaves have an adverse effect on ABA. Furthermore, the effect of introducing exogenous ABA on plant organ abscission was weaker than the effect of exogenous ETH. Osborne [21] concluded that ABA caused the premature senescence of organ cells, but that this may be due to increasing ETH content in plant organs that leads to organ shedding.
Transcriptome refers to a collection of all the transcription products under a specific physiological condition, including messenger RNA, ribosomal RNA, transport RNA and noncoding RNA [22]. Transcriptome sequencing research is essential for novel gene discovery, functional gene annotation, gene differential expression and the development of molecular markers [23][24][25]. Next Generation Sequencing (NGS) technology enables complete and rapid access to almost all transcript sequence information of a particular tissue or organ of a species and has been widely used in modal and nonmodel plant research [26][27][28]. To our knowledge, there have been no previous reports of research on DEGs related to defoliation traits in alfalfa varieties. This study aimed to construct a transcriptome database of two alfalfa leaf varieties in the early blooming stage. Based on DEGs, we aimed to identify the critical genes related to defoliation traits in the Zhungeer variety (characterized by easy leaf abscission) and the WL319HQ variety (without easy leaf abscission). The results of this study provide a theoretical basis for future alfalfa molecular breeding and provide a reference for subsequent study of the transcription of alfalfa leaves.

De novo assembly transcriptome
Six libraries of the total RNA extracted from the leaves of Zhungeer and WL319HQ varieties at the early flowering stage were constructed for high-throughput sequencing. A total of 261,775,004 clean reads with a total of 39,202,954,927 nucleotides (nt) were obtained from the six sequencing libraries. Trinity method assembled a total of 66,734 unigenes with an N50 of 1496 nt. The maximum, minimum and average lengths were 143,789 nt, 201 nt and 869 nt, respectively (Fig. 1).

Identification and analysis of DEGs
Using the Nr annotation results, Blast2GO software was used to analyze GO functional annotations of unigenes, and WEGO software was used to perform the functional classification of unigenes. A total of 1,098 significant DEGs were assigned to one or more ontologies by the standard of |log2FC| > 1 and FDR < 0.05 (Additional file 1: Table S1; Additional file 2: Table S2; Fig. 3), 706 unigenes were upregulated, and the other 392 unigenes were downregulated (Zhungeer vs WL319HQ). The significant DEGs were used for the subsequent analysis.
In the 87 KEGG pathways, the DEGs related to direct or indirect effects on the content of ABA and ETH were predicted. In the "Plant hormone signal transduction" pathway (ko04075), 294 unigenes with 6 DEGs were annotated. In the "Circadian rhythm-plant" pathway (ko04712), 61 unigenes with 4 DEGs were annotated. In addition, the "Carotenoid biosynthesis" pathway (ko00906) contained 44 unigenes with 1 DEG crucial for ABA synthesis, which was annotated.

Real-time quantitative PCR analysis of DEGs related to defoliation
To test the reliability of transcriptome sequencing data in the three pathways, Unigene0044746 (Auxin response factor, ARF), Unigene0002039 (Phytochrome interacting factor 3, PIF3), Unigene0027311 (Ethylene receptor protein, ETR), Unigene0053251 (Phytochrome B, PHYB), Unigene0053032 (Cryptochrome, CRY) and Unigene0014585 (9-cis-epoxycarotenoid dioxygenase 3, NCED3) (Additional file 4: Table   S4) were selected for One-Step RT-qPCR analysis. In addition to DEGs, Unigene0030464 (Tubulin alpha, TUBA) and Unigene0032887 (indole-3-acetic acid-amido synthetase, GH) were selected in the RT-qPCR analysis as unaffected genes. The alfalfa GAPDH gene was used as the endogenous reference. The significant difference of RT-qPCR data between two varieties was analyzed by t-test, and the expression patterns revealed by RT-qPCR analysis were similar to those revealed by transcriptome sequencing for the same genes (Fig. 4).

Discussion
ABA, ETH, jasmonic acid, salicylic acid, CTK and auxin are the primary phytohormones that regulate plant senescence [29]. When plants are under environmental stress, such as drought, salinity or low temperature, the leaves will accelerate senescence and abscission due to the elevation of ABA concentration [11,30]. In addition, ETH stimulates the synthesis of cellulase and controls the release of cellulase from protoplasts into the cell wall, which promotes cell wall degradation and leaf abscission [5]. Many studies reported that ETH levels would rise if plant tissues were damaged by mechanical force [4,31,32]. Therefore, the abscission of alfalfa leaves may be due to mechanical damage to plant tissues during harvesting, resulting in elevation of ETH concentration in leaf tissues. The mechanism by which phytohormone interactions promote the abscission of leaves is not well understood and is difficult to study [33].
In higher plants, xanthoxin is not only the synthesis precursor of ABA but also a metabolic intermediate of carotenoid catabolism [34]. The RNA-Seq data revealed that the alfalfa 9-cis-epoxycarotenoid dioxygenase 3 (NCED3) gene differentially expressed between WL319HQ and Zhungeer varieties. The RT-qPCR analysis showed that the transcription of the Zhungeer variety was approximately 2.3 times higher than the WL319HQ variety. NCED is the enzyme that catalyzes the cleavage reaction of 9-cis-epoxycarotenoids to produce xanthoxin and is the rate-limiting enzyme of ABA biosynthesis in higher plants [35][36][37]. NCED has been identified in maize, kidney bean, tomato and Arabidopsis thaliana [4,35,36,38,39]. A. thaliana NCED and its role in ABA biosynthesis are the most well-studied [37]. Overexpression of A. thaliana NCED3 resulted in a significantly higher level of endogenous ABA in the plant [40], and similar results were obtained in transgenic tobacco to overexpress the tomato NCED1 gene [38].
The "Plant hormone signal transduction" (ko04075) pathway plays a vital role in the metabolism of phytohormones in higher plants. In the RNA-Seq data, the expression level of the ethylene receptor protein (ETR) gene in Zhungeer variety was approximately 1.4 times higher than the WL319HQ variety. ETR is a protein that initiates the ETH signal transduction pathway by binding the ETH and relaying the ETH signal to constitutive triple response 1 (CTR1) Raf-like protein kinase [41]. The high-affinity binding activity has been confirmed in A. thaliana, as well as the tomato [42][43][44]. In the present study, the ETR gene was downregulated in WL319HQ compared to the Zhungeer variety, indicating that Zhungeer is more likely to accumulate ETH than the WL319HQ variety.
Auxin response factors (ARF) are transcription factors that bind specifically to the DNA sequence 5'-TGTCTC-3′ found in auxin-responsive promoter elements (AuxREs). Although ARF genes are not involved in leaf senescence directly, they act as transcriptional activators or repressors and therefore promote flowering, stamen development, floral organ abscission and fruit dehiscence [45]. Study of the A. thaliana ARF1 and ARF2 mutant has shown that ARF2 may affect plant tissue shedding [46]. Nagpal et al. [47] found that the A. thaliana ARF6 gene promotes flower organ maturation and senescence. In addition, NPH4/ARF7 and ARF19 not only enhance the ARF2 phenotype but also induce leaf senescence [46]. The ARF gene was significantly downregulated in alfalfa leaves in the WL319HQ variety.
Plant circadian rhythm is a hereditary physiological regulatory mechanism that synchronizes with circadian alternation in the plant. The response of plants to environmental factors mainly depends on the hormone regulatory network, including CTK, IAA, ABA and GA [48]. Circadian rhythm plays an essential role in the plant hormonal regulation pathway [49,50]. At the same time, the plant hormone signaling network also acts on circadian rhythm to transmit different metabolic signals and environmental signals to the endogenous circadian rhythm system, thus forming a complex regulatory network [51,52]. In this study, three DEGs involved in the plant circadian rhythm pathway were identified that may affect the regulation of ABA and ETH. Compared with Zhungeer, the PHYB (− 1.775-fold), PIF3 (− 3.715-fold) and CRY (− 1.952) gene of WL319HQ were downregulated. PIF3 plays a critical role not only in light and temperature-mediated environmental signaling but also in the signaling of ABA and ETH. Light activates the phytochrome pathway to induce the degradation of PIF3 and promotes the accumulation of carotenoids, which may lead to an increase in ABA content.
In summary, it is likely that the Zhungeer alfalfa variety is more susceptible to abscission than the WL319HQ variety due to higher levels of plant hormone accumulation. Higher expression of NCED3 and PIF3 genes caused increasing ABA content, which may be a reason why the Zhungeer variety is more likely to defoliate. In addition, upregulation of ETR and ARF expression affects ETH; auxin accumulation is another indirect evidence as well.

Conclusion
This study is the first to use RNA-Seq to identify the defoliation traits between alfalfa varieties. It is essential to understand the physiological effects of various hormones Fig. 4 The relative expression level of DEGs involved in "Plant hormone signal transduction", "Circadian rhythm-plant" and "Carotenoid biosynthesis" pathways revealed by RT-qPCR analysis. TUBA: Tubulin alpha, no significant difference (P = 0.165); GH: indole-3-acetic acid-amido synthetase, no significant difference (P = 0.246); ARF: Auxin response factor, significant difference (P = 0.033); PIF3: Phytochrome interacting factor 3, significant difference (P = 0.043); ETR: Ethylene receptor protein, significant difference (P = 0.014); PHYB: Phytochrome B, significant difference (P = 0.047); CRY: Cryptochrome, significant difference (P = 0.049); NCED3: 9-cis-epoxycarotenoid dioxygenase 3, significant difference (P = 0.04) on plants, the interaction between hormones and their relationship with the agricultural production environment. Our transcriptome data dramatically enriches alfalfa genome research. This research will provide theoretical guidance for field production practice and genetic breeding and provides a reference for the future study of defoliation traits in alfalfa.

Sample collection
We selected two varieties of the tetraploid M. sativa: the Zhungeer variety that has easy leaf abscission and the WL319HQ variety that is not characterized by easy leaf abscission. Both varieties were grown in experimental fields located in Baotou City, Inner Mongolia, China (110°3 7′E, 40°5'N). Leaves of the Zhungeer and WL319HQ varieties were collected from two-year-old plants at the early flowering stage in June. All samples were immediately frozen in liquid nitrogen and stored at − 80°C for future RNA extraction.

RNA extraction and RNA-Seq library construction
Total RNA was extracted from three biological replicates of the leaves of each M. sativa variety using the RNeasy Plant Mini Kit (Qiagen, Germany), following the manufacturer's protocol. The integrity of RNA samples was measured by Agilent 2100 Bioanalyzer (Agilent Technologies, USA). The Poly (A) mRNA was enriched by NEBNext Oligo(dT) 25 beads (NEB, USA) from 50 μl total RNA. Then, the enriched mRNA was constructed to a cDNA library by NEBNext Ultra RNA Library Prep Kit for Illumina (NEB, USA), following the manufacturer's protocol.

Sequencing and raw reads filtering
The transcriptome sequencing was performed using an Illumina HiSeq 4000 sequencing platform with pairedend 150 bp (PE150) sequencing strategy at Gene Denovo Biotechnology Co. (Guangzhou, China). Raw reads obtained from the sequencing that may affect the following assembly and analysis were removed, including reads with adapters, reads containing more than 10% unknown nucleotides and low-quality reads containing more than 50% low-quality bases (Q-values ≤10).

De novo transcriptome assembly and unigenes detection
Since M. sativa genome information was not previously available, full-length transcriptome de novo assembly was carried out using the reference genome independent Trinity method [53]. Trinity is a software package consisting of three modules: Inchworm, Chrysalis and Butterfly. First, Inchworm assembled clean reads to a collection of linear contigs by the greedy k-mer method. Second, the related minimal overlapping contigs created by Inchworm were clustered into connected components by Chrysalis, and a de Bruijn graph was constructed for each cluster of related contigs. Finally, Butterfly analyzed the paths from the de Bruijn graphs that were taken by reads and pairs of reads, and the complete transcripts and unigenes sets were obtained from the outputs [54].

Analysis of DEGs
Gene expression analysis of unigenes was performed and normalized by RPKM (reads per kb per million reads). The edgeR package (http://www.r-project.org/) was used to identify DEGs across samples or groups. Genes with a fold change ≥2 and a false discovery rate (FDR) < 0.05 were considered significant DEGs. DEGs were further subjected to GO functional enrichment analysis and KEGG pathway analysis.

Real-time quantitative PCR analysis
One-Step Real-time quantitative PCR (RT-qPCR) was carried out to confirm the RNA-Seq data from the extracted M. sativa leaves' total RNA. The experiments were performed using One-Step TB Green PrimeScript RT-PCR Kit II (Takara, Japan) on an ABI 7500 system (Applied Biosystems, USA). Primers for One-Step RT-qPCR were designed using Primer Premier 5.0 software (Premier, Canada). The alfalfa GAPDH gene was selected as the endogenous reference [55]. All the primers were synthesized by Invitrogen (Beijing, China) and are listed in Additional file 5: Table S5. A 20 μl reaction mix was set up containing 50 ng total RNA and 0.4 μM of each primer. The One-Step RT-qPCR analysis was performed on the extracted total RNA in three biological replicates of the leaves of each M. sativa variety, and each RNA sample was carried out in technical triplicate. The program was set at 42°C for 5 min then 95°C for 10 s, followed by 40 cycles of 95°C for 5 s then 60°C for 34 s. The final melting curve analysis was set at 95°C for 15 s, 60°C for 1 min and 95°C for 15 s. The relative expression changes of the endogenous reference and tested genes were analyzed by the 2 -△△ CT method [56].