De novo transcriptome assembly
The total RNA extracted from the root crown of Medicago sativa during the wintering stage were constructed six libraries for high-throughput sequencing. In order to ensure the quality of information analysis, the original data was filtered. Among the raw reads, low quality reads, those containing adapter sequences, or low quality bases were discarded, and 459,865,666 clean reads with a total of 68,515,729,607 nucleotides (nt) were obtained from the six sequencing libraries. A total of 114,567 unigenes with an N50 of 1092 nt were assembled. The maximum length, minimum length and average length were 15,687 nt, 201 nt and 740 nt, respectively (Fig. 1).
Function annotation
A total of 57,712 unigenes (50.37% of total 114,563 unigenes) were annotated against Nr, SwissProt, Pfam, COG, GO and KEGG databases using BLASTX (E-value < 1 ´ 10-5). Among them, 56,044 (48.92%), 34,764 (30.34%), 30,449 (26.58%), 6,670(5.82%), 41,457(36.19%), and 21,483(18.75%) unigenes were annotated to the Nr, SwissProt, Pfam, COG, GO and KEGG databases, respectively (Fig. 2a). According to the Nr database, 10.59% (5,933) of unigenes showed homology (1 ´ E-10 < E-value ≤ 1 ´ E-5), 14.37% (8,051) of unigenes showed strong homology (1 ´ E-20 < E-value ≤ 1E-10), and the remaining of 62.68% (35,133) unigenes had very strong homology (E-value ≤ 1E-30) (Fig. 2b). For the species distribution of the top BLASTX hits, 78.28% (43,813) unigenes matched to the homologous sequences of Medicago truncatula, while 3,482, 2380, 1,140, 548 and 388 unigenes matched to the homologous sequences of Trifolium pretense, Trifolium subterraneum, Cicer arietinum, Cajanus cajan and Glycine max, respectively (Fig. 2c).
Identification and analysis of differentially expressed genes (DEGs)
Using the Nr annotation results, BLAST2GO 2.5.0 software was used to perform GO functional annotations of unigenes. A total of 4,442 significant DEGs were assigned to one or more ontologies by the standard of |log2 fold change| >1 and P-value < 0.05 (Additional file 1: Table S1; Fig. 3), there were 2,299 unigenes up-regulated, and the other 2,143 unigenes were down-regulated (Longmu 806 vs Sardi). Those DEGs were used for the next analysis.
In order to better understand the function of genes differentially expressed during the wintering period of alfalfa, GO functions were used to classify the function of DEGs (Additional file 2: Table S2). A total of 4,442 unigenes were summarized in three main functional categories “biological processes”, “cellular component”, and “molecular function” (Fig. 4). In the biological processes category, “cellular process” (GO:0009987), “metabolic process” (GO:0008152), “single-organism process” (GO:0044699) and “biological regulation” (GO:0065007) were the most frequent terms and contained 728, 702, 470 and 260 unigenes, respectively. In the molecular function category, genes were focused on subcategory including catalytic activity (GO:0003824) and binding (GO:0005488).
KEGG pathway enrichment analysis
KEGG analysis of differentially expressed genes using KOBAS (v2.1), a total of 1,159 unigenes were assigned to 116 KEGG pathways (Additional file 3: Table S3). The ten top KEGG pathways with the highest representation of the DEGs were: “Starch and sucrose metabolism” (ko00500), “Protein processing in endoplasmic reticulum” (ko04141), “Plant-pathogen interaction” (ko04626), “Phenylpropanoid biosynthesis” (ko00940), “Plant hormone signal transduction” (ko04075), “Ribosome” (ko03010), “MAPK signaling pathway-plant” (ko04016), “Endocytosis” (ko04144), “Glycolysis/Gluconeogenesis” (ko00010), “Spliceosome” (ko03040).
In the “Plant hormone signal transduction” pathway (ko04075), DEGs associated with ABA and ethylene biosynthetic pathways (including SnRK2 and EIN3) were up-regulated, respectively. In the “peroxidase” pathway (ko04146), DEGs associated with superoxide dismutase X1 (SOD1) was up-regulated. These annotations provided valuable resources for studying the specific functions and pathways of the alfalfa gene.
Stress response of transcription factors (TFs) to low temperature stress
Transcription factors (TFs) play important roles in the response to abiotic (low temperature, salt, and drought) stresses and directly control the expression of specific sets of stress-responsive genes [43, 44]. A total of 34 TFs families containing 1,364 unigenes were identified (Fig. 5). According to our data, these TFs families responding to low temperature stress included MYB, B3, AP2/ERF, C2C2, WRKY, NAC, FAR1, bHLH, LBD, C3H, bZIP, GRAS, C2H2, MADS, LOB, HSF, SBP, TCP, GRF, ZF-HD family (Additional file 4:Table S4). Many of these TFs families have been reported to play an important role in plant responses to abiotic stresses including low temperature stress [45], and have been utilized to improve plant abiotic stress tolerance by gene transfer technology [46].
Among these TFs, the top six TFs families with the highest representation were MYB (186), B3 (145), AP2-ERF (106), C2C2 (92), WRKY (88), and NAC (82). The MYB family members detected in our data 86 were up-regulated by low temperature stress. The AP2/ERF family was very important TFs family that regulates gene expression to cold and freezing stress in a variety of plants [47]. The AP2/ERF family members detected in our data 43 were up-regulated by low temperature stress. Moreover, most of the unigenes expression was up-regulated in the WRKY family under freezing stress, and 45 unigenes in the B3 TFs family were up-regulated under low temperature stress. In the C2C2, bHLH, bZIP, LOB, SRS, E2F/DP, GeBP and BBR-BPC TFs families, the number of TFs with up-regulated expression was less than the number of TFs with down-regulated expression. In addition to the above TFs, five unigenes of the EIL family were also identified as up-regulated under low temperature stress. Furthermore, the number of TFs with up-regulated and down-regulated expression was the same in the AP2/ERF, NAC, NF-Y, Nin-like, and NF-X1 families.
Validation of DEGs data by Real-Time Quantitative PCR (RT-qPCR) analysis
To verify the reliability of transcriptome sequencing data, differences in expression of ten genes were detected using RT-qPCR, and these genes were well characterized by the NCBI Nr database (Fig. 6a). These genes included ethylene (EIN3), ABA (SnRK2), auxin (ARF), and jasmonic acid (JAZ) response regulatory genes in the plant hormone signaling pathway, and MYB, B3, AP2/ERF, WRKY TFs family, and detection of reactive oxygen species (ROS) transcript superoxide dismutase X1(SOD1) associated with low temperature stress response (Additional file 5: Table S5). The expression of transcription factor genes of Longmu 806 verity was up-regulated, and the expression of seven transcription factor genes was significantly higher than that of sardi verity under low temperature stress. By analyzing the expression profiles of selected genes during wintering, correlation analysis was performed with RT-qPCR analysis results and transcriptome sequencing results (Fig. 6b). Moreover, RT-qPCR showed a high correlation (R2=0.8016, P<0.05) of fold change between RNA sequencing analysis and RT-qPCR. The expression patterns in PCR assay were generally in agreement with the results of the RNA-seq.