Transcriptome profiling of Malus sieversii under freezing stress after being cold-acclimated

Background Freezing temperatures are an abiotic stress that has a serious impact on plant growth and development in temperate regions and even threatens plant survival. The wild apple tree (Malus sieversii) needs to undergo a cold acclimation process to enhance its freezing tolerance in winter. Changes that occur at the molecular level in response to low temperatures are poorly understood in wild apple trees. Results Phytohormone and physiology profiles and transcriptome analysis were used to elaborate on the dynamic response mechanism. We determined that JA, IAA, and ABA accumulated in the cold acclimation stage and decreased during freezing stress in response to freezing stress. To elucidate the molecular mechanisms of freezing stress after cold acclimation, we employed single molecular real-time (SMRT) and RNA-seq technologies to study genome-wide expression profiles in wild apple. Using the PacBio and Illumina platform, we obtained 20.79G subreads. These reads were assembled into 61,908 transcripts, and 24,716 differentially expressed transcripts were obtained. Among them, 4410 transcripts were differentially expressed during the whole process of freezing stress, and these were examined for enrichment via GO and KEGG analyses. Pathway analysis indicated that “plant hormone signal transduction”, “starch and sucrose metabolism”, “peroxisome” and “photosynthesis” might play a vital role in wild apple responses to freezing stress. Furthermore, the transcription factors DREB1/CBF, MYC2, WRKY70, WRKY71, MYB4 and MYB88 were strongly induced during the whole stress period. Conclusions Our study presents a global survey of the transcriptome profiles of wild apple trees in dynamic response to freezing stress after two days cold acclimation and provides insights into the molecular mechanisms of freezing adaptation of wild apple plants for the first time. The study also provides valuable information for further research on the antifreezing reaction mechanism and genetic improvement of M. sieversii after cold acclimation. Supplementary Information The online version contains supplementary material available at 10.1186/s12864-021-07998-0.


Background
Central Asian mountains are particularly important and have been listed as a global biodiversity hotspot. These hotspots contain ancestors of some domestic fruit, including apples. Malus sieversii (also called "Xinjiang wild apple"), the ancestor of domesticated apple [1], is mainly distributed in the Tianshan Mountains of China, Kazakhstan, Kyrgyzstan, and Uzbekistan. M. sieversii is a germplasm resource with excellent cold resistance traits and is usually used as a popular rootstock for breeding cold-resistant domesticated apple in China [2,3]. Therefore, it is of great significance to understand the cold resistance mechanism of wild apple for its further protection and utilization.
As an abiotic stress, low temperature not only affects the geographical distribution of crops, but also has a negative impact on crop yield and quality. Cold stresses are divided into chilling (0-15°C) and freezing (< 0°C) damage, can lead to tissue damage and even plant death, and often bring large economic losses to agricultural production [4]. When plants are exposed to low temperatures, the cell membrane is the main site of cold-induced damage, which triggers a series of physiological and biochemical reactions [5]. After exposure to low temperature for a period of time, many plants can improve frost resistance, a process called cold acclimation (CA) [6]. During cold acclimation, plants are able to sense environmental temperature changes, and then trigger a series of protective mechanisms, including changing the composition, structure and function of membrane; synthesis of cryoprotectant molecules, such as soluble sugar, protein and proline; increase the activity of protective enzymes [7]. CA can also induce the expression of some cold-induced genes during the stress period, and these genes likely play a role in stabilizing the cell membrane against freezing injury [8]. For example, CBF/DREB1 protein from Arabidopsis transcription factor family has been identified to control the expression of a regulon of cold-induced genes that increase plant freezing tolerance [9]. In apple (Malus domestica), homologue MdCIbHLH1 of Arabidopsis ICE1 contributes to improved cold tolerance by directly activating CBF genes [10]. Furthermore, it is known that the MYB-bHLH-WDR complex, especially the bHLH and MYB components, is involved in modulating anthocyanin accumulation at low temperatures [11]. Additionally, the accumulation of sucrose and other monosaccharides during cold acclimation also contributes to membrane stability because these molecules can protect the membrane from freezing damage in vitro [12]. The entire CA process in nature typically requires much time and is easily affected by climate fluctuations. However, with global warming, the winters and early springs in temperate regions are becoming progressively milder, and temperature patterns are becoming increasingly irregular. This increases the frequency of temperature changes that may cause premature subzero temperatures, thereby subsequent increasing the risk of freezing injury [13]. Additionally, changing phenological patterns, such as an early start of growing season and early flowering time [14], consistent with climate warming, it may increase the risk of tissue damage caused by frost. The likelihood of such scenarios is typically high during early spring. Therefore, it is very important to understand the mechanism of freezing stress after cold acclimation to protect and utilize plants. However, little is known about the molecular mechanism of freezing stress after a short period of cold acclimation in M. sieversii.
According to previous research, such as on Hordeum vulgare [15], Arabidopsis thaliana [16], Glycine max [17], and Camellia sinensis [18], the molecular mechanism of plants responding to freezing stress has been explained. Transcriptome sequencing was used to obtain differentially expressed genes (DEGs) and further explore the role of related signaling pathways and candidate genes to explain the changes in the plant response to freezing stress. For instance, strawberry (Fragaria×ananassa) responds to freezing stress through several pathways: flavonoid biosynthesis, plant hormone signal transduction, MAPK (mitogen-activated protein kinase) signaling, starch and sucrose metabolism, circadian rhythm [19]. A previous transcriptome study identified 4173 and 7734 DEGs in two apple cultivars, respectively, during the chilling and freezing treatments [20]. Oil rapeseed (Brassica napus) responds to low-temperature stress through two conserved (primary metabolism and plant hormone signal transduction pathway) and two novel (plant-pathogen interaction and circadian rhythm pathway) signaling pathways [21]. Transcriptomic analysis of Magnolia wufengensis under cold stress showed that the response mechanism was related to photosynthesis, and plant hormone signal transduction, primary and secondary metabolism pathways [22].
RNA-seq based on the Illumina platform is a powerful method for investigating gene expression profiles and revealing signal transduction pathways in a wide range of biological systems [23]. In recent years, RNA-seq technology has been widely used to clarify the response of various plants to low temperature stress, including C. sinensis [24], and A. thaliana [25]. In apple (Malus domestica), comparative transcriptome analysis revealed 377 commonly upregulated and 211 commonly downregulated DEGs involved in the crosstalk between drought, cold and high-salinity stress [26]. Single molecule real-time (SMRT) sequencing based on PacBio (Pacific Biosciences) platform can capture the full length of transcripts, which provides a simpler and more accurate method for gene annotation and isoform identification. The third generation sequencing technology represented by PacBio has the advantage of long read length [27]. PacBio and RNA-seq are highly complementary to each other. To obtain an overall view of the molecular regulation during freezing stress after the CA process, we combined PacBio sequencing and RNA-seq to investigate the transcriptome of M. sieversii seedlings at low temperatures. This study will provide a theoretical basis for freezing stress resistance and lay a foundation for understanding the cold response molecular mechanism in M. sieversii, which provides a rich genetic resource for further research investigating the freezing response in apple and other woody plants.

Results
Physiological and plant hormone changes of M. sieversii under low temperature stress In this study, 8-week-old apple plantlets were used as experimental materials to study the transcriptomic and physiological changes of M. sieversii under freezing stress (Fig. 1A). The well-developed wild apple seedlings were cold acclimated at 4°C in a chamber for two days. Subsequently, the temperature was gradually lowered from 4°C to − 4°C (2°C/h) for 4 h, and samples were collected at 0, 1, 6, and 10 h intervals at − 4°C (Fig. 1B). To evaluate the physiological damage caused by low temperature, we measured relative electrolyte leakage (REC) values, maximum photosynthetic efficiency (Fv/ Fm) values, malondialdehyde (MDA) content, and sucrose content. The plants treated at − 4°C showed a relatively higher sucrose content, which was approximately 200% higher during the 1 h interval (Fig. 1C). Furthermore, under low-temperature treatment, the Fv/ Fm values were significantly reduced compared with those of the control (Fig. 1E). In addition, the REC values and MDA production in the treated group were relatively higher than those in the control (Fig. 1D, F), indicating that plasma membranes may be damaged Bars represent the mean ± SD (n = 3). Accumulation levels were analyzed using a T-test. Within each figure, asterisks above the bars indicate statistical significance (*p < 0.05; **p < 0.01) under freezing stress. Low-temperature treatment significantly affected the endogenous hormone content of Xinjiang wild apple leaves. Compared with the control, the contents of endogenous jasmonic acid (JA) and abscisic acid (ABA) after cold acclimation (0 h) increased significantly. Compared with cold acclimation (0 h), the contents of endogenous JA and ABA in the treated leaves showed a downward trend (0-10 h), but these values were still higher than those of the control (Fig.  1G,H), except for the ABA content at 10 h. Moreover, we also determined the content of IAA (Fig. 1I). Surprisingly, the content of IAA in the cold acclimation stage (0 h) did not decrease, but under frozen stress (0 h-10 h), the content decreased significantly. In general, cold and freezing environments produce different hormonal responses.

Analysis of sample Iso-Seq results
To identify and characterize the transcript dynamics of M. sieversii under freezing treatment after cold acclimation, we employed PacBio SMRT and RNA-seq technologies for transcriptome sequencing. In this study, the sample materials use for transcriptomic analysis and the treatment methods were consistent with the physiological measurements in Fig. 1. Subsequently, the cDNA libraries constructed from treated and CK samples were sequenced. After correcting the data, a total of 11,781,904 (20.79G) subreads were obtained (Table 1). In addition, a total of 611,407 circular consensus sequence (CCS) reads were obtained after filtering. Finally, 485,006 full-length non-chimeric reads (FLNC) were obtained and poly-A tails in which an average FLNC read length was 1765 bp. In total, 432,508,752 nucleotides were obtained, and 191,961 reads were generated. The minimum length was 193 bp, and the maximum length was 14,168 bp. The N50 value was 2613 bp, and the N90 value was 1344 bp. The length of each subread sequence was counted after quality control of the raw data, and the length distribution map was generated (Additional file 1).
GMAP software was used to match the consensus sequences to the reference genome sequences, and a total of 176,234 (91.81%) reads were mapped to the Malus domestica genome sequences (GDDH13 v1.0); according to the results, the sequences can be divided into four types: unmapped, 15,727 (8.19%); multiple mapped, 10,438 (5.44%); reads mapped to the forward strand (+), 95,421 (49.71%); and reads mapped to the antisense strand (−), 70,375 (36.66%) ( Fig. 2A). To obtain a comprehensive genetic annotation, the functions of novel genes were annotated based on seven databases. Of the 3917 novel genes, 3808 were significantly matched in at least one of the databases. Among them, 3359, 2242, 3216, 1625, 2038, 3708, and 2038 novel genes were found in the NR (NCBI nonredundant protein sequence), SwissProt, KEGG (Kyoto Encyclopedia of Genes and Genomes), KOG/ COG (Clusters of Orthologous Groups of proteins), GO (Gene Ontology), NT (NCBI nucleotide sequences by BLAST, E-value 1e-10), and Pfam (Protein family) databases, respectively (Additional file 1). The histogram demonstrated that 1097 genes were simultaneously annotated across all seven databases. We compared the 61,908 transcripts with the M. domestica genome gene set, and they were classified into three groups, with following results are shown in Fig.  2B: isoforms of known genes (13,069, 21.11%), novel isoforms of known genes (43,794,70.74%), and isoforms of novel genes (5045, 8.15%). Alternative splicing (AS) events were analyzed with SUPPA software. We detected 22,316 transcripts with alternative splicing (AS) events, including skipped exons (SE), mutually exclusive exons (MX), alternative 5′splice site (A5), alternative 3′splice site (A3), retained intron (RI), alternative first exon (AF) and alternative last exon (AL) (Fig. 2C, Additional file 6). The transcripts were analyzed by principal component analysis (PCA) (Additional file 2). The first two axes of the PCA represent the differentiation between the biological replicates of the different treatments. In addition, 1710 fusion transcripts were identified (Additional files 3, 7).

Identification and analysis of DETs
In total, we obtained 24,716 DETs in different comparisons (Additional file 8). To better understand the gene expression changes between the CK and cold-treated apple leaf samples in response to freezing after 2 d of  3B, Additional files 9, 10). We conducted a GO enrichment and KEGG pathway analysis on all upregulated and downregulated DETs to explore their biological functions (Additional files 11,12,13). In addition, KEGG pathway analyses were also conducted to identify several important pathways ( Fig. 3C-E, Additional file 4,). In the pairwise, DETs were significantly enriched in the "biosynthesis of secondary metabolites", "photosynthesis", "peroxisome", "starch and sucrose metabolites" and "plant hormone signal transduction" categories. All the pathways worked together to form a complex and efficient signaling network that reduced or even eliminated cold stress damage. The pathway enrichment results indicated that the major pathways included plant hormone signal transduction, starch and sucrose metabolism, peroxisomes, transcription factors, and photosynthesis.

Responses of the plant hormone signaling network to cold stress
In the ABA signal transduction pathway (Fig. 5A), 2 PYL genes (MD08G1043500, MD12G1178800, MD12G1178800_novel01) were upregulated under freezing stress conditions, whereas another PYL gene (MD05G1300200), a protein phosphatase 2C

Transcription factors (TF)
TFs are known to play essential roles in plant responses to abiotic stress. The iTAK and HMMER analysis of apple transcripts allowed the identification of nonredundant TFs. In our work, we analyzed the expression patterns of 1801 transcripts of 676 genes in the top 15 transcription factor families (Fig. 6, Additional file 14). Heatmaps of abundance-altered TFs were plotted to depict the transcript abundance profiles of each TF family (Fig. 6B, C). The MYB superfamily was the largest, including 154 MYB and 178 MYB-related members. The second-largest family was the bHLH family, with 161 members, followed by the bZIP (143)

Validation of DETs by qRT-PCR analysis
Sixteen DETs related to cold adaptation in different signaling pathways and TFs were randomly selected to validate the gene expression levels using qRT-PCR. They

Discussion
Apple is one of the most important fruit crops worldwide. M. sieversii, a progenitor of domesticated apple, represents a reservoir of genetic diversity and is widely used as a rootstock in China [28]. The selection of rootstocks is usually based on stress tolerance [29]. Cold stress is one of the most severe abiotic stresses in apple plants and is one of the many limiting factors for global apple production. Through cold acclimation, plants can alleviate the harm of low temperatures and enhance their tolerance to cold stress [30]. In tea plants, differential changes in transcription and metabolism levels occurred in response to CA and cold stress conditions [31]. Since the cold-stress response network is very complex and contains an array of biochemical and physiological modifications, many defense mechanisms are activated during cold acclimation, thus the transcriptome changes will be less intense even when the plant is exposed to lower temperatures. To discover the key functional and regulatory genes of cold resistance in M. sieversii, we examined seedlings under cold stress after 2 d of cold acclimation and performed fullsequence transcriptomic analysis. Full-sequence transcriptomic data have great potential for discovering new or previously unrecognized isoforms, meanwhile, different expression patterns appearing between isoforms of the same gene. For example, in present work, isoforms MD05G1126500_novel01, MD08G1162500_novel01, MD1 2G1040700_novel07, MD12G1126100_novel02, MD05G 1126500_novel01, MD06G1228000 are up-regulated, another isoforms of the corresponding gene MD05G1126500_ novel02, MD08G1162500, MD12G1040700_novel06, MD12G1126100_novel08, MD05G1126500_novel02, MD06G1228800_novel01 are down-regulated. In addition, the transcriptome results revealed that isoforms involved in transcription regulation, signal transduction of hormones, transcription factors, and oxidative stress were differentially expressed under cold stress, we also compared the activities of MDA and Fv/Fm and the contents of ABA and JA in seedlings of M. sieversii. To further explore the response mechanism of apple under freezing stress, we identified a series of genes from different pathways. This study provides a novel understanding of the roles of CA in freezing stress in M. sieversii.

Sugar and starch metabolism
The metabolism of sugar and starch plays a vital role in the cold resistance of plants. Sugar can be used as a signaling molecule in hormone synthesis, plant growth and development, and various antistress reactions, as a nutrient to supply energy and as a protective agent to avoid cold damage [32,33]. SUS is a key enzyme involved in sucrose metabolism and is part of a family of genes encoding different sucrose isozymes, including SUS1 and SUS3 (MD11G1307000, MD02G1100500) were upregulated. Here, the induction of two genes encoding SUS may point towards the degradation of sucrose to UDP-glucose and fructose. UGDH (MD16G1062900, MD13G1064400) was induced. Under low-temperature stress, the task of increasing glycolysis can be completed [34]. Starch hydrolysis was associated with increased maltose accumulation and increased expression of BAM, which encodes the beta-amylase necessary for starch mobilization and raises the levels of soluble sugars capable of acting as osmolytes or antioxidants [35][36][37]. Under cold stress conditions, the activity of BAM1 was largely unaffected, but BAM3 expression increased [38]. In this study, two related SS genes (MD09G1052000 and MD00G1130900) were repressed, and another SS1 (MD01G1218000) was induced. Notably, the gene expression of BAM family members, such as MD16G1231500, MD06G1044200, MD06G1043900, and MD06G1044100, was significantly upregulated under freezing stress (Fig. 4A). This phenomenon may indicate that Xinjiang wild apples respond to freezing stress by accelerating the decomposition of starch to provide energy and achieve heat preservation.

Antioxidant defense system
In the antioxidant defense system, various nonenzymatic and enzymatic ROS detoxification mechanisms have developed in plant cells. Major enzymatic systems include superoxide dismutase (SOD), ascorbate peroxidase (APX), peroxidase (POD), and catalase (CAT) [39].
MDA levels and electrolyte permeability are commonly used as indicators of membrane structural damage and plant resilience [40,41]. Under abiotic stress, antioxidant enzymes can eliminate harmful substances to the greatest extent and effectively reduce membrane damage. Plants with suppressed APX production induce SOD and CAT to compensate for the loss of APX, whereas plants with suppressed CAT production induce APX and GPX. In addition, POD and CAT are hydrogen peroxide scavengers and catalyze the conversion of hydrogen peroxide (H 2 O 2 ) into water. In this study, genes encoding CAT (MD06G1008700), GST (MD10G 1172300), GPX6 (MD06G1081300), APXT (MD04G11 04900), SPS1 (MD10G1258500) and RBOHD (MD1 3G1134500, MD13G1134500_novel02, MD13G1134500_ novel10) were upregulated under freezing stress conditions.
However, the genes encoding POD (MD03G1013200, MD10G1206100, MD04G1171700, MD03G1059200, and MD01G1162100) were significantly downregulated (Fig. 4B). Overall, the antioxidant defense system-related genes also rapidly responded to freezing stress in Xinjiang wild apples, enabling plants to cope with freezing stress. Antioxidant defense systems involve many kinds of antioxidant genes, which have unique functions under stress conditions, forming a complex antioxidant defense system in vivo.

Responses of the plant hormone signaling network to cold stress
We measured ABA accumulation at − 4°C for different times after cold acclimation. When cold acclimation ended, a large amount of ABA accumulated. ABA levels were increased in plants that were stressed, such as those under drought, cold, or high temperature stresses, and ABA played a crucial role in plant stress resistance [42][43][44]. During cold acclimation, ABA accumulation can induce the expression of some cold-resistant genes to improve freezing resistance [45]. Numerous studies found a direct link between the accumulation of ABA in woody plants and the increase in cold tolerance [46]. Important components of ABA include PP2C and PYL, which is closely related to the level of plant stress resistance [47,48]. PYL (MD05G1300200) and PP2C (MD01G1139200, MD01G1139200_novel04) were induced in response to cold stress. ABF2 regulates various aspects of the ABA response by controlling the expression of a subset of ABAresponsive genes [49,50]. In this study, ABF (MD08G1099600, MD08G1099600_novel01, MD15G1 081800, MD15G1081800_novel02, MD15G1081800_ novel03) was upregulated under freezing stress conditions. These results indicated that ABA stimulated this signaling pathway and that DEGs were involved in the ABA signaling response to freezing stress.
We also measured the content of JA, which show a consistent trend as that for ABA content. After cold acclimation (0 h), the content of JA increased significantly, while the content of JA decreased during the freezing treatment (0 h-10 h). However, the JA content was still higher than that of the control. Similar to ABA, cold stress induces a number of JA biosynthetic and signaling genes that ultimately trigger the accumulation of JA and increase cold tolerance by regulating the ICE-CBF/ DREB1 cascade [51]. Consistent with the positive role of JA in freezing stress, the production of endogenous JA was triggered by cold treatment. Exogenous application of MeJA significantly improved plant freezing tolerance with or without cold acclimation [52]. In JA signal transduction pathway, JAZ genes (MD02G1096100, MD06G1228900, MD09G1178600, MD10G1260700, MD15G1116100, MD15G1116100_novel01, MD16G112740, and MD17G1164400) was upregulated under freezing stress in this study. The bHLH transcription factor MYC2 (MD01G1086900, MD06G1119900, and MD11G1254100) acted as a direct target of the JAZ protein to prevent JA signaling. In this study, JAZ protein was induced under freezing stress, whereas the transcription factor MYC2 was repressed. Taken together, these results indicated that the JA signaling pathway regulated Xinjiang wild apple freezing tolerance after cold acclimation.
BR and ET also play major roles in plant defense responses to cold stress [53]. It is worth noting that in BR signal transduction pathways, BAK1 (MD04G1054400), a transmembrane receptor kinase; BZR1 (MD15G1031900, MD08G1035400), a core member of the BR signaling pathway; and CYCD3 (MD15G1288100) were repressed under cold stress (Fig. 5). ET levels increase in response to cold stress and can induce the accumulation of antifreeze proteins [54]. In this study, genes encoding CTR1 (MD09G1236700, MD10G1168800), EIN3 (MD02G1266200), and ERF1 (MD04G1228800, MD04G1228800_novel01) were upregulated. In the ET signaling pathway, CTR1 is a positive regulator in the ethylene signaling pathway. EIN3 is controlled by a ubiquitin proteasome system, and ERF1 is an immediate target for EIN3 [55]. The results indicated that these genes in the ethylene signaling transduction pathway responded to freezing stress. These results indicate that both the BR and ET signaling pathways may function in M. sieversii cold response processes after acclimation.
In particular, the IAA content did not decrease significantly during cold acclimation but decreased significantly during freezing stress. We demonstrated that the endogenous hormone IAA may have a special role in cold acclimation. Moreover, the decreasing trend of the IAA content under freezing stress is similar to that of ABA and JA, which may involve coordination with ABA and JA. Although the role of auxin after cold stress is not clear, genes related to the IAA signal transduction pathway were differentially expressed in this study. IAA (MD10G1176400, MD13 G1204700) and SAUR (MD12G1027600, MD07G1297400, MD07G1297400_novel01) were induced in response to low temperature (Fig. 5E), whereas ARF (MD01G1083400, MD14G1131900) and SAUR (MD10G1059800, MD10G 1061300) were repressed. Analysis of genes related to the IAA signal transduction pathway showed that the TIR1 (Novelgene0397) transcript was highly expressed. TIR1 can cause AUX/IAA ubiquitination, relieving the inhibition of ARF activity by AUX/IAA and hence initiating the transcription of IAA-induced genes. Recent studies have shown that plant changes induced by cold stress are closely related to the auxin response [56][57][58]. These genes were differentially expressed in response to low temperature, indicating that the auxin signaling pathway may participate in the M. sieversii cold response. In conclusion, the M. sieversii freezing response after CA involves a complicated regulatory network containing multiple types of plant hormone signaling pathways.

Photosynthesis
The direct impact of abiotic stress on the activity of photosynthesis was the disruption of all photosynthetic components, such as photosystems I and II, electron transport, carbon fixation, ATP generation, and stomatal conductance [59], which resulted in serious damage to the photosynthetic machinery of plants [60]. In this study, DEGs associated with photosynthesis showed that most of these genes were downregulated, such as PsaN (MD01G1004300), PsaL (MD04G1244400, Novel-gene0781), PsaO (Novelgene3166), PsbP (MD17G1 132900, MD09G1146300), PetC (MD09G1098300, MD09G1098300_novel02), and ATPc (MD05G1347300). These results suggest that photosynthesis may be inhibited by freezing conditions. It is worth noting that two genes were upregulated, PSII core complex gene psbY (MD06G1160500) and PSII reaction center gene PSB28 (MD10G1176000). The evidence of genes that encode proteins linked to photosynthesis in the previous study as psbY is indicating that the slowdown in the degradation of chlorophylls by 1-MCP would also affect other photosystem II-related proteins [61]. PsbY gene upregulation expression may play a critical role in photosynthetic systems under freezing stress. Low temperatures also cause a steady-state light imbalance, which reduces the quantum efficiency of PSII. This implies that freezing stress could induce a series of changes in genes related to photosynthetic pathways, thus decreasing Fv/ Fm values.

Transcription factor dynamics during freezing responses
In this study, 1801 transcripts of 676 TFs from the top 15 TF families responding to freezing stress were examined (Fig. 6A), and the expression patterns of the top 15 TF families were analyzed (Fig. 6B,C). Many of these families of TFs were reported to play central roles in plant responses to abiotic stress, including cold stress [62][63][64][65][66]. A total of 49 AP2/ERF members were detected under freezing stress. DREB1/CBF (MD04G1067800, MD04G1067800_novel01, MD04G1067800_novel02, MD04G1165400, MD04G1165400_novel01), a subfamily member of AP2/ERF, was identified to be induced at least tenfold and twofold under freezing stress, respectively. Similarly, in previous studies, the expression of the CBF gene was upregulated in apple bark at different freezing temperatures [20]. Similarly, MYB113 (MD09G126510), MYB1 (MD09G1278600, MD09G1278600_novel02), MYB62 (MD13G1039900), MYB108 (MD13G1148200, MD13G1148200_novel01) and MYB4 (MD16G1218900) strongly responded to freezing stress after cold acclimation, and this response continued throughout the whole freezing stress period. In previous studies, three MYB transcription factors (MdMYB12, MdMYB22, and MdMYB114), which had several CBF/DREB response elements in their promoters, were significantly induced by low-temperature exposure and their expression also correlated highly with anthocyanin accumulation [11]. In apple, MYB transcription factors, such as MYB1 and MYB113, MYB114 are involved in the regulation of anthocyanin synthesis pathways [67,68]. During freezing stress, low temperature induces the expression of MYB transcription factors and enhances anthocyanin synthesis to enhance freezing resistance. In addition, recent investigations in apple show that MdMYB4-transgenic apple calli exhibited enhanced tolerance to cold stress, as indicated by reduced electrolyte leakage value, and MDA content than the control calli [69]. The upregulation of MYB4 may increase cold resistance in wild apples. It also appears likely that WRKY TFs play important roles in M. sieversii cold regulation. Under freezing stress, 118 WRKY members were differentially expressed, among which WRKY71 (MD10G1266400, MD10G1266400_novel01), WRKY6 (MD10G1324500), and WRKY70 (MD12G1189900) were induced at least 5fold. In Arabidopsis, WRKY70 in promoting BRmediated gene expression to regulated plant growth and drought response, WRKY71 activity hastens flowering, thereby providing a means for the plant to complete its life cycle in the presence of salt stress [70,71]. In wild apple, the high expression level of WRKY70 and WRKY71 may affect plant growth to cope with low temperature environment. Our transcriptome data also showed that other transcription factors, such as bHLH

Conclusion
In this study, we first present the transcriptome of M. sieversii using PacBio single-molecule long-read sequencing to analysis of the responses to freezing temperature after cold acclimation. The expression patterns of 4410 DETs were determined, and identify several important pathways genes such as genes encoding superoxide dismutase, sucrose synthase, plant hormone-signaling transduction kinase, and peroxidase. In addition, a total of 1801 transcripts were confirmed as TFs. Some different members of the AP2/ERF, WRKY, MYB, and bHLH gene families were significantly up/downregulated. Furthermore, the transcription factors DREB1/CBF, MYC2, WRKY71, WRKY6, WRKY70, MYB4 and MYB88, were strongly induced during the whole stress period. In short, cold acclimation can induce the continuous expression of some antifreeze genes during freezing stress. These findings provide a characterization of gene transcription and facilitate the understanding of the mechanisms of tolerance to low temperature after cold acclimation in M. sieversii.

Plant materials
Seeds of M. sieversii were collected in September 2016 from one area (43°23′2.20″N; 83°35′43.48″E) in a natural wild reserve forest in Yili, Xinjiang. With the permission of the Xinyuan County Forestry Bureau, seeds were obtained from the field and identified [72]. The field investigation and sample collection experiments were performed under institutional guidelines in accordance with local legislation. The germplasm of M. sieversii was identified by Ph.D. Wenjun Li, who worked at the Xinjiang Institute of Ecology and Geography, Chinese Academy of Sciences. A voucher specimen was deposited in the Natural Museum, Xinjiang Institute of Ecology and Geography, Chinese Academy of Sciences. (MS-20160919-015). M. sieversii seedlings were planted in 1/ 2 Murashige and Skoog (MS) medium with 30 g·L − 1 sucrose and 6 g·L − 1 agar in a greenhouse with 32 μmol·m − 2 s − 1 light intensity, 40-50% relative humidity, and a 16 h light/8 h dark photoperiod at 24°C. Eightweek-old well-developed wild apples grew in medium at 4°C under a 16 h photoperiod for two days (cold acclimation). Then, the chamber was cooled gradually (2°C/ h) from 4°C to − 4°C for 4 h (Boxun, LRX-580C-LED). Remove upper and lower three leaves and the middle leaves were separately harvested at the time points of CK, 0, 1, 6 and 10 h and each sample contained three biological replicates. Leaf samples were frozen in liquid nitrogen and stored at − 80°C before physiological change analysis, hormone concentration determination, and RNA extraction. The RNA-seq was conducted using fifteen samples (CK, 0, 1, 6, 10 h) and the PacBio sequencing was implemented using the mixture of the samples.

Measurement of physiological changes and endogenous hormone contents
Before sequencing the transcriptome, we determined the physiological indicators. MDA was measured according to a previously published method [73]. Chlorophyll fluorescence was recorded with a PAM2500 fluorometer (Walz, Germany), and the maximal quantum efficiency of PSII was calculated as Fv/Fm (Fv = Fm-F0). REL was conducted using the procedure described by this method [74], and it was calculated based on the following equation: The determination of endogenous ABA and JA content was performed using HPLC (Comin Suzhou, China), with three replicates for each analysis.

RNA quantification and transcriptome sequencing
Total RNA was isolated from the leaves of freeze-treated and untreated M. sieversii samples using the E.Z.N.A. Plant RNA Kit (Omega Bio-Tek, Norcross, GA, USA). Equal quantities of RNA from three biological replications at each stage were pooled to construct a complementary cDNA library. RNA concentration and quality were measured using a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, Waltham, Massachusetts, USA). cDNA library construction was performed by high-throughput sequencing using the PacBio and Illumina platforms at Novogene (Beijing, China). Additional nucleotide errors in consensus reads were corrected using the RNA-seq data with the software LoRDEC [75]. Index of the reference genome sequences was built using HISAT2 (v2.1.0) and paired-end clean reads were aligned to the reference genome sequences using HISA T2. Consensus reads were aligned to the reference genomic sequence (GDDH13, Version 1.0) using GMAP with parameters-no-chimeras-cross-species-expand-off sets1-B5-K50000-f samse-n1 against the reference genome sequences [76]. Gene structure analysis was performed using the TAPIS pipeline [77]. The GMAP output bam format file and gff/gtf format genome sequences annotation file were used for gene and transcript determination. Alternative splicing events and alternative polyadenylation events were then analyzed.
Fusion transcripts were determined as transcripts mapped to two or more long-distance genes and were validated by at least two Illumina reads.

Iso-Seq data analysis
Raw reads were processed into error-corrected reads of interest (ROIs) according to the conditions full passes ≥0 and sequence accuracy> 0.75. By detecting whether each ROI sequence contains a 5′ primer, 3′ primer, and poly-A tail, the sequences can be divided into full-length sequences and non-full-length sequences. Based on the poly-A tail signal and the 5′ and 3′ cDNA primers in the ROIs, FLNC transcripts were identified [78]. The ICE (Isoform-clustering) algorithm in SMRT Analysis (v2.3.0) software was used to obtain a consistent sequence isoform, and the full-length consensus sequences were polished using Quiver. Full-length transcripts with a post correction accuracy above 99% were generated for further analysis. The redundant sequences were removed from the high-quality and corrected low-quality transcript sequences of each sample using CD-HIT (v4.6.8) software [79].

Differential expression analyses
Analysis of transcript expression levels used a fragments per kilobase of transcript per million mapped reads (FPKM) method. Differential expression analysis was performed using the DESeq R package (Version 1.18.0). The P values were adjusted using the Benjamini & Hochberg method. A corrected P-value of 0.05 and log2 (fold change) of 1 were set as the thresholds for significant differential expression.
GO enrichment analysis of DETs or lncRNA target genes was performed by the GO seq R package [80], in which gene length bias was corrected. GO terms with corrected P-values less than 0.05 were considered significantly enriched for differentially expressed genes. KEGG is a database resource for understanding the high-level functions and utilities of biological systems, such as cells, organisms, and ecosystems [81], from molecular level information, especially large-scale molecular datasets generated by genome sequencing and other highthroughput experimental technologies (http://www. genome.jp/kegg/). KOBAS (v.3.0) software was used to test the statistical enrichment of DETs in KEGG pathways [82]. Plant transcription factors were predicted using iTAK (v1.7A) software [83].

Quantitative real-time PCR analysis of DETs
The qRT-PCR assays were performed to validate the consistency of RNA-Seq analysis. QRT-PCR analysis was performed using a Light Cycler 480 machine with SYBR Green I Master Mix (TaKaRa, China). Three biological replicates were set for each treatment. The relative expression of each transcript was calculated after normalization to the reference gene. The relative quantification from three biological replications was normalized to the reference gene and was calculated by the 2 -ΔΔCt method. All primer sequences are shown in Additional file 15. Three technical replicates were performed based on the 2 -ΔΔCt algorithm using EF1α as a quantitative control.

Statistical analysis
Statistical analysis was performed using ANOVA, and then a new multirange test was performed using Duncan's range test in SPSS version 16.0. A significance level of P < 0.05 was applied to indicate significant values. Three independent biological replicates were used for RNA-seq and physiological and hormonal assays. Graphs were created using Origin 8.0 (Version 9.1; Origin Lab, Northampton, MA, USA).