Transcriptome sequencing and assembly
The original data obtained by sequencing with Illumina Hiseq 4000 were transformed into raw reads by base calling. The total base number obtained from nine libraries was 64.52 G and the total number of nucleotides was between 5.23 G and 7.59 G in each sample. The Q30 ratio of each sample was greater than 94% and GC content was relatively consistent, around 47% (Additional file 1). A total of 124,106 transcripts and 55,084 non-redundant unigenes were obtained from nine databases, respectively. The average length of unigenes was 661 bp and length of N50 was 1470 bp. Consequently, the sequencing data quality was high and could meet the requirements of subsequent analysis. There were 29,582 unigenes with the length of 200–500bp accounting for 53.70%, 10,103 unigenes of length 500–1000 bp accounting for 18.34%, 10,112 unigenes ranged from 1000 to 2000 bp accounting for 18.35%, and 5,287 unigenes of length more than 2000bp accounting for 9.6%, respectively (Fig. 1).
Functional annotation of unigenes
Based on Nr, Swiss-prot (a manually annotated and reviewed protein sequence database), Pfam (protein family), KOG (Clusters of Orthologous Groups of proteins), KEGG and GO databases, which were used to annotate all unigenes with comprehensive gene function information. In the present study, a total of 37,303 unigenes were successfully annotated in R. xanthina f. spontanea database, which were shown all of 67.72% (55,084). Furthermore, Fragaria_vesca presented the highest frequency in the annotation results with a total of 19,964 comments, accounting for 53.52% of all the sequences, followed by Nelumbo nucifera (9.0%), Prunus persica (3.27%), Phaseolus vulgaris (3.13%), Vitis vinifera (2.78%) and Prunus mume (2.36%)(Fig. 2 and Additional file 2).
A total of 31,258 (56.75%) unigenes were assigned to GO terms in the cellular component, molecular function and biological process categories that were further classified into 50 GO terms (Fig. 3 and Additional file 3).
Within the cellular component, a total of 22,417 DEGs were assigned under 4℃and − 20℃ stress (23℃ as the control), which indicated the union of all the DEGs were mainly related to the nucleus, cytoplasm, integral component of membrance and chloroplast. With regard to the molecular function, most of the DEGs were enriched for molecular function, protein binding and ATP binding. In the biological process, biological process, regulation of transcription and DNA–template were mainly involved.
To identify the metabolic pathways during cold stress of R. xanthina f. spontanea leaves, 21,992 DEGs were mapped to the KEGG database, and the top 19 KEGG pathways significantly (p-value < 0.05) enriched were shown in Fig. 4 (Additional file 4). Among these pathways, carbohydrate metabolism (2596, 11.80%), translation (2325, 10.57%) and folding, sorting and degradation (2036, 9.26%) were the most extensively over-presented pathways. In addition, the metabolic pathways were related to environmental adaptation, transport and catabolism. Overview, amino acid metabolism, translation and signal transduction were involved in almost every aspect of plant life. Furthermore, it was demonstrated that high-informative and wide coverage were obtained from transcriptome sequencing of leaves in R. xanthina f. spontanea, and can be used to analyze the gene products of metabolic pathways and information processing pathways at the molecular level.
DEGs in response to low-temperature stress
In our studies, the amount of unigenes expression among three libraries under low-temperature (23℃, 4℃ and − 20℃)was used to carry out comparative analysis of differences, which was shown significantly difference (Padj < 0.05, Fig. 5). 358 DEGs (3031 up-regulated and 3891 down-regulated); 42 DEGs (867 up-regulated and 1763 down-regulated); 279 DEGs (2869 up-regulated and 2872 down-regulated) responded to low-temperature stress within the 23℃, 4℃and − 20℃ were detected, respectively. And there were 1536 common genes between T1 vs T2 and T1 vs T3; 3977 common genes appeared in T2 vs T1 and T2 vs T3, and 1068 common genes existed in T3 vs T1and T3 vs T2 (Fig. 5; 23℃, 4℃and − 20℃ was labeled as T1, T2 and T3, respectively). Otherwise, 468 DEGs had been present in above three groups, which can be associated with cold-tolerance in single petal R. xanthina f. spontanea. Therefore, it was essential to utilize GO functional enrichment analysis and KEGG pathway enrichment analysis of the DEGs.
GO enrichment analysis of the common DEGs
In the study, to reveal which biological functions are significantly related to the common DEGs we obtained, a GO functional-enrichment analysis was carried out via the agriGO website with a p score cut-off of 0.05. The results indicated that the DEGs were involved in biological process, cellular component and molecular function occupied 64.84%, 9.38%, 25.78%, respectively. Consequently, most DEGs were significant correlated with some biological functions. And we found the DEGs were classified into 83 biological process, mainly focused on transcriptional regulation, DNA-template; transcriptome, response to chitin, response to abscisic acid, response to cold, response to water deprivation, ethylene-activated signaling pathway and response to wounding. As to cellular component, the DEGs were involved in chloroplast, chloroplast thylakoid membrane, integral component of plasma membrane, etc. With respect to molecular function, the DEGs mainly focused on transcription factor activity, sequence-specific DNA binding (Fig. 6 and Additional file 5).
KEGG pathway enrichment analysis of the common DEGs
In order to more precisely investigate variation of metabolic pathways in leaves during low-temperature stress, statistical pathway enrichment analysis for the DEGs was carried out based on KEGG database and 293 DEGs under low-temperature stress were assigned to 85 different KEGG pathways (p < 0.05). It consisted of four significant enriched pathways, which was plant-pathogen interaction, starch and sucrose metabolism, plant circadian rhythm and photosynthesis-antenna proteins, respectively (Table 1 and Additional file 6). When the temperature reached 4℃ and − 20℃, most DEGs in plant-pathogen interaction and starch and sucrose metabolism pathways were down-regulated. In addition, the genes in plant circadian rhythm pathway were approximately the same with down-regulated and up-regulated. However, seven DEGs in photosynthesis-antenna proteins pathway were shown up-regulated.
Table 1
Summary of KEGG pathway functional annotations for the common DEGs
Pathway
|
ID code
|
The number of DEGs in metabolic pathway
|
Number of the common DEGs
|
T2-T1
|
T3-T1
|
Up
|
down
|
Up
|
Down
|
Plant-pathogen interaction
|
ko04626
|
32
|
6
|
26
|
4
|
28
|
Starch and sucrose metabolism
|
ko00500
|
21
|
4
|
17
|
3
|
18
|
Circadian rhythm - plant
|
ko04712
|
8
|
3
|
5
|
4
|
4
|
Photosynthesis - antenna proteins
|
ko00196
|
7
|
7
|
0
|
7
|
0
|
qRT-PCR validation of the DEGs
Based on the results of transcriptome annotation under low-temperature stress, eight DEGs were screened relating to cold-resistant (Fig. 7 and Additional file 7). The expression fold change was analyzed by the use of qRT-PCR before and after low-temperature stress. The results indicated that the expression of these PGR5、CHLH、BBX24、STN7、EXPA8、LRR-RLK and CIPK12 unigenes were up-regulated, one of bZIP60 was down-regulated, which the same trend was observed in analysis of high-throughput sequencing. Consequently, these genes mentioned above may be a direct correlation with cold-resistant in R. xanthina f. spontanea, and the results further confirmed reliability of our transcriptome data.