MAPK and hormone signaling and carotenoid biosynthesis pathways are involved in regulating flower opening of pear tree by response to temperature mediated by melatonin

Flower opening marks the transition of plants from vegetative growth to reproductive growth. The flower opening of plants are affected by various environmental factors such as temperature, light intensity, and length of daylight^[1,2]. In Arabidopsis, there are six classic flowering pathways which are: the vernalization pathway, the ambient temperature pathway, the photoperiodic pathway, the gibberellin pathway, the age pathway, and the autonomous pathway^[1,2]. In the past decade, researchers have made a clearer and deeper understanding of the vernalization pathway, the photoperiodic pathway, and the gibberellin pathway^[3−5]. The molecular mechanism of ambient temperature regulating flower opening were studied sufficiently in model plants^[6].

In Arabidopsis, high temperature (27 °C) could promote flowering^[7], however, it showed a phenotype of late flowering at low temperatures (16 °C)^[6,7]. The regulation of plant flowering by environmental temperature is affected by multiple pathways, in which SVP or FLM and their homologous genes play a key role^[8,9]. SVP is a flowering inhibitory factor, which encodes a MADS-box DNA binding protein and acts as a transcription factor. It inhibits plant flowering by binding to the promoter elements of FT and SOC1 genes. The autonomous regulatory factors FCA and FVE participate in the environmental temperature regulation pathway by inhibiting SVP^[10,11]. FLM transcribes a pair of mutually exclusive splices FLM-β and FLM-δ. This splicing form is temperature-dependent. FLM-β is the main splicing type at 16 °C, while FLM-δ is more important at 27 °C. FLM-β and FLM-δ regulate plant flowering by competing with SVP to respond to environmental temperature^[12]. In addition, multiple microRNAs (miRNAs) are also involved in the regulation of plant flowering by ambient temperature. Research has shown that miRNAs are central regulators of plant development^[13]. At 16 °C, miR156 overexpression down-regulated SPL3, thereby delaying flowering and high expression of SVP and FLM-δ directly bound to the miR172 promoter, down-regulating the expression of miR172, thereby contributing to the formation of inhibitors^[14]. PIF4 is a transcription factor of the bHLH family and it is related to the thermal morphogenesis. Reports showed that regulation of plant flowering by temperature under short-day conditions required PIF4 mediation^[15]. And PIF4 and PIF5 could initiate the expression of FT at warm nights and could initiate FT-independent mechanisms to promote flowering during warm days^[15]. PWR and HDA9 have been reported to interact with each other and modulate the flowering time by repressing AGL19 expression, whereas HOS15 and HDA9, together with the photoperiodic evening complex, regulate flower opening through repression of GI transcription^[16].

Melatonin is an active small molecule of indole-like substance. It is involved in the regulation of many physiological processes such as development, flowering, nutrient absorption and fruit ripening in plants^[17]. As for flowering, Kolár et al. found that the flowering process of Chenopodium rubrum was induced by photoperiod and the subsequent development of apical meristem was affected by the concentrations of exogenous melatonin, and found that melatonin could affect some early steps of its flowering^[18]. Murch et al. found that there was a peak of melatonin accumulation in the development stage of white datura buds and its content decreased with the maturity of buds^[19]. Park et al. found that melatonin could be induced to synthesize during rice flower bud development to regulate flowering time^[20]. The results of Shi et al. showed that melatonin at least partially inhibited the flower transition of Arabidopsis thaliana through DELLA protein^[21]. During two consecutive years of monitoring, Zhang et al. found that apple tree flowering was always associated with a decrease in melatonin content and melatonin could mediate apple tree flowering by simultaneously regulating the expressions of key genes in the photoperiodic pathway, the vernalization pathway, the gibberellin pathway and the autonomous pathway^[22].

Pear (Pyrus), a perennial deciduous fruit tree of the Rosaceae family, is one of the most important cultivated fruit tree species in the world. Influenced by global warming, the temperature of spring in northern China rises and falls sharply in recent years, especially from late March to early May^[22]. The intense temperature changes are the biggest disadvantageous to agricultural production and leads to a reduction in production. In terms of pear tree and other fruit trees that flower during this period, the adverse impact will undoubtedly double because the sudden lower or higher temperatures will lower the quality of flowers, which determines the quality and the yield of products. That is to say, the frequent occurrence of lower or higher temperatures in spring is becoming a serious threat to the pear industry. To avoid the sudden temperature changes by regulating the flower opening of pear trees is one of the effective measures for reducing losses caused by the precipitous lower or higher temperature. Ambient temperature is one of the pathways regulating plant flowering. However, the molecular mechanism of pear tree flower opening regulated by ambient temperature has not been researched in earnest, and we want to know more about it. From the previous representation, we found that melatonin plays a vital role in mediating apple tree flowering in spring, thus we indicated that it may also function in mediating pear tree flowering for the two plants have a close genetic relationship. Therefore, in this study, we treated the pear branches with flower buds with exogenous melatonin of different concentrations under lower or higher temperatures to determine whether melatonin plays a role in temperature regulating pear tree flowering. Our results showed that melatonin functions in pear tree flowering under both lower and higher temperatures. Then we sampled the pear inflorescences treated with exogenous melatonin under lower and higher temperatures respectively and sent them for transcriptome sequencing. Subsequently, we analyzed the transcriptome data and screened out the capital pathways and the key genes in each pathway of ambient temperature regulating pear flower opening mediated by melatonin. Our study could provide a theoretical basis for further researches on regulation of plant flower opening and provide some new ideas for reducing losses caused by lower or higher temperature in pear industry.

'Xin li No.7' is an early maturing variety of pear with a strong sweet flavor, high yield, and strong disease resistance. The branches with flower buds of 'Xin li No.7' pear tree were clipped from the teaching and experimental base of Hebei Agricultural University. They were then randomly divided into eight groups ensuring that the number of flower buds in each group was roughly the same. Among all the groups, four of them were treated with melatonin solution of 0 (the control), 20, 200 and 1,000 μM respectively at a higher temperature (25 °C), and the other four groups were treated with 0 (the control), 20, 200 and 1,000 μM melatonin respectively at a lower temperature (12 °C) (Table 1). During the treatments, illumination conditions and humidity were the same, illumination being 10000 Lux and humidity being 70%. When the percentage of blooming flowers in one treatment under the same temperature reached 30%, the inflorescences of the four groups were sampled, frozen with liquid nitrogen and stored at −80 °C for further use.

RNA extraction, library preparation and sequencing

The total RNA of the pear inflorescence samples was extracted by the Polysaccharide polyphenols plant total RNA extraction kit (QIAGEN, Germany), and its integrity was detected precisely through the Agilent 2100 bioanalyzer. The first strand of cDNA was synthesized in M-MuLV reverse transcriptase system using fragmented mRNA as a template and random oligonucleotides as primers. Then the RNA strands were degraded by RNaseH, and the second strand of cDNA was synthesized from dNTPs in DNA polymerase I system. Subsequently, the purified double-stranded cDNA was repaired by terminal repair, added A tail, and connected to the sequencing connector. The cDNA of about 370−420 bp was screened by AMPureXP beads, the PCR amplification was carried out, and the PCR product was purified by AMPureXP beads. Finally, the library was obtained. Transcriptome sequencing was conducted through the Illumina NovaSeq 6000 (Illumina, USA), and completed by Novogene Co., Ltd (Beijing, China).

Transcriptome data and gene annotation analysis

The raw data was obtained by transforming the image data measured by the high-throughput sequencer into sequence data (reads) via the CASAVA base recognition. To ensure the quality and reliability of data analysis, the raw data needed to be filtered. It mainly included the removal of reads with connectors (adapter), the removal of reads containing N (N means that the base information cannot be determined), and the removal of low-quality reads (reads whose base number of Q phred ≤ 20 accounts for more than 50% of the total read length). At the same time, the contents of Q20, Q30, and GC in clean data were calculated. All the subsequent analyses were based on high-quality clean data. Downloaded the reference genome and the gene model annotation files directly from the genome website (www.ncbi.nlm.nih.gov), and employed the HISAT2v2.0.5 software to build an index of the reference genome and compare the paired terminal clean reads with the reference genome.

Differential expression genes and functional enrichment analysis

DESeq2 software (1.20.0) was used to analyze the differential expression between the two comparison combinations. The filter condition is |log2(Fold Change) | ≥ 1 & padj ≤ 0.05. Cluster Profiler (3.8.1) software (https://bioconductor.org/packages/clusterProfiler/) was used to analyze the GO functional enrichment and KEGG pathway enrichment of different gene sets. The threshold of significant enrichment was padj < 0.05 for the GO functional enrichment and the KEGG pathway enrichment. The interaction relationship in STRING protein interaction database (https://cn.string-db.org/) was used to analyze the differential gene protein interaction network. The figures were drawn using cystoscope software.

Identification of transcription factor analysis

Conserved Domains in NCBI were used to predict transcription factors, the transcription factor families were counted according to the prediction results, and the related transcription factors were screened. The string diagram were generated using Origin (2022 edition).

To clarify whether melatonin treatment can affect the flower opening of pear trees at lower or higher temperatures, pear branches with flower buds were sampled and treated with melatonin of different concentrations at 12 and 25 °C respectively.

The results showed that 20 μM and 200 μM MT could delay pear flower buds opening, while 1,000 μM MT could accelerate pear flower opening under higher temperature conditions, compared with the control (Fig. 1a−d). Excitingly, under lower temperature condition, compared with the control, all the exogenous MT treatments could accelerate pear flower buds opening (Fig. 1e−h). In terms of the percentage of blooming flowers of each treatment, it showed that in higher temperature (Fig. 1i), the percentage of blooming flowers of pear flower buds of 1,000 μM MT treatment (17%) was 4.25 times higher than that of the control (4%) on the 11^th day. The percentage of blooming flowers of pear flower buds of 20 μM (9%), 200 μM (9%) and 1,000 μM (29%) MT treatment groups was 0.33, 0.33 and 1.07 times higher than that of the control (27%), respectively, on the 12^th day. The percentage of blooming flowers of pear flower buds of 20 μM (31%), 200 μM (26%) and 1,000 μM (54%) MT treatment groups was 0.65, 0.55 and 1.15 times higher than that of the control (47%), respectively, on the 13^th day. And the flowering rates of 20 μM (43%), 200 μM (49%) and 1,000 μM (81%) MT treatment groups were 0.72, 0.82 and 1.35 times that of the control (60%), respectively on the 14^th day. While at lower temperature, all the flower opening of pear flower buds was delayed compared with that under higher temperature (Fig. 1j). What's more, there were 11%, 8% and 7% blooming flowers in the 20 μM, 200 μM and 1,000 μM MT treatments respectively on the 28^th day. The percentage of blooming flowers of pear flower buds of 20 μM (32%), 200 μM (15%) and 1,000 μM (37%) MT treatment groups were 6.4, 3.0, and 7.4 times that of the control (5%), respectively, on the 29^th day. The results indicate that melatonin could mediate ambient temperature to regulate flower opening of pear trees.

Figure 1. 

Effects of exogenous melatonin on the phenotypes of pear flower opening and on the percentage of blooming flowers of pear flower buds at different temperatures. (a)−(d) Phenotypes of pear flower opening at higher temperature. (e)−(h) Phenotypes of pear flowers opening at lower temperature. (i) Percentage of blooming flowers of pear flower buds treated with melatonin of different concentrations at higher temperature. (j) Percentage of blooming flowers of pear flower buds treated with melatonin of different concentrations at lower temperature.

Transcriptome sequencing data and quality assessment

To reveal the molecular mechanism of melatonin mediating ambient temperature to regulate flower opening of pear trees, transcriptome analyses were carried out, and the sequencing data were obtained (Table 2). The results indicated that the number of raw reads obtained ranged from 40.06 million to 47.49 million, while that of the clean readings obtained by filtering and quantization ranged from 40.03 million to 44.74 million. The error rate was 0.03%, the GC content ranged from 45.73% to 47.12%, and the score of Q30 was above 90.69%. Compared with the reference genome of pear tree, there were 74.39%−76.36% of all reads totally mapped to it, and among which there were 7.25%−7.77% mapped to the multiple positions and 66.94%−68.76% mapped to the unique positions, suggesting that the transcriptome data obtained were of high quality and could be used for later analyses.

To confirm the relationships of all samples treated with melatonin under lower or higher temperature conditions, the Pearson correlation between samples and the principal component analysis (PCA) were performed respectively. As shown in Fig. 2a, under lower temperature conditions, the R² values among the three replicates of the same treatment were all larger than 0.953 and that among the four treatments were above 0.934. For PCA analysis (Fig. 2b), three replicates of the control (C1), the 20 μM MT (T1) and 200 μM MT (T2) clustered together, while that of the 1,000 μM MT (T3) were relatively dispersed, indicating that samples of the control, the 20 μM MT and 200 μM MT treatments were of high reliability and more DEGs existed among them. As for the results at higher temperature (Fig. 2c & d), it showed that the R² values among the three replicates of the same treatment were all larger than 0.949 and that among the four treatments were above 0.943 (Fig. 2c). In addition, three replicates of the control (C2), the 200 μM MT (T5) and 1,000 μM MT (T6) clustered together, while that of the 20 μM MT (T4) were relatively dispersed (Fig. 2d), indicating that samples of the control, the 200 μM MT and 1,000 μM MT treatments were of high reliability.

Figure 2. 

Overview of transcriptome samples. (a) Pearson correlation analysis heat map of samples treated at lower temperature. (b) Principal component analysis of samples treated at lower temperature. (c) Pearson correlation analysis heat map of samples treated at higher temperature. (d) Principal component analysis of samples treated at higher temperature

Identification of DEGs among different pair treatments

For the sake of analyzing the number of DEGs significantly regulated in different treatments, a volcanic map of DEGs was constructed (Fig. 3a). As shown in Fig. 3b, compared with the C1 group, there were 610 (268 up-regulated and 342 down-regulated), 587 (453 up-regulated and 134 down-regulated) and 180 (139 up-regulated and 41 down-regulated) in T1, T2, and T3 groups, respectively. Compared with the C2 group, 384 (273 up-regulated and 111 down-regulated), 218 (152 up-regulated and 66 down-regulated), 349 (257 up-regulated and 92 down-regulated) in T4, T5, and T6 groups, respectively. The proportions of DEGs in different treatments were further compared with a venn diagram (Fig. 3c & d). It was found that there was a total of 28 DEGs significantly expressed in the four treatments at lower temperature and a total of 46 DEGs significantly expressed in the four treatments at higher temperature.

Figure 3. 

Identification of the number of differentially expressed genes (DEGs) among different treatments at lower or higher temperatures. (a) The number of DEGs presented by volcanic maps of all the pair treatments. (b) The number of up-regulated and down-regulated DEGs in each pair treatment. (c) DEGs among the three pair treatments presented by a Venn diagram at lower temperature. (d) DEGs among the three pair treatments presented by a Venn diagram at higher temperature.

GO annotation analysis of DEGs among different treatments

GO (Gene Ontology) is a comprehensive database describing gene function, which can be divided into three parts: biological process (BP) and cellular composition (CC) molecular function (MF). To gain insights into the biological functions of DEGs among all the pair treatments, GO annotation analysis was performed (Fig. 4). From the results of GO enrichment analysis, 30 terms were selected from each pair treatment to draw the histogram displayed. In the T1 and C1 pair treatment, 610 differential genes were enriched to 22 GO functions, including the movement of cell or subcellular component (GO:0006928) and microtubule-based movement (GO:0007018) in BP, including copper ion binding (GO:0005507), hydrolase activity that acting on glycosyl bonds (GO:0016798), hydrolase activity that hydrolyzing O-glycosyl compounds (GO:0004553), heme binding (GO:0020037), tetrapyrrole binding (GO:0046906), oxidoreductase activity that acting on paired donors with incorporation or reduction of molecular oxygen (GO:0016705), iron ion binding (GO:0005506), etc. in MF.

Figure 4. 

GO enrichment histogram between groups. (a) Histogram of gene ontology annotation analysis between different paired groups at lower temperature. (b) Histogram of gene ontology annotation analysis between different paired groups at higher temperature.

In the T2 and C1 pair treatment, 587 differential genes were enriched to 16 GO functions, including defense response (GO:0006952) and response to stress (GO:0006950) in BP, including sequence-specific DNA binding (GO:0043565), flavin adenine dinucleotide binding (GO:0050660), oxidoreductase activity acting on paired donors with incorporation or reduction of molecular oxygen (GO:0016705), heme binding (GO:0020037), tetrapyrrole binding (GO:0046906), monooxygenase activity (GO:0004497) in MF.

In the T3 and C1 pair treatment, 171 differential genes were enriched to 5 GO functions, including sequence-specific DNA binding (GO:0043565), oxidoreductase activity acting on single donors (GO:0016703), flavin adenine dinucleotide binding (GO:0050660), monooxygenase activity (GO:0004497), metal ion transmembrane transporter activity (GO:0046873).

In the comparative combination of T4 and C2, 384 differential genes were enriched to 14 GO functions, including cell wall (GO:0005618), external encapsulating structure (GO:0030312), apoplast (GO:0048046), extracellular region (GO:0005576), and cell periphery (GO:0071944) in CC, including sequence-specific DNA binding (GO:0043565), hydrolase activity hydrolyzing O-glycosyl compounds (GO:0004553), hydrolase activity acting on glycosyl bonds (GO:0016798), copper ion binding (GO:0005507), oxidoreductase activity that acting on paired donors with incorporation or reduction of molecular oxygen (GO:0016705), iron ion binding (GO:0005506), etc. in MF.

In the comparative combination of T5 and C2, the different genes were mainly enriched in heme binding (GO:0020037) tetrapyrrole binding (GO:0046906) sequence-specific DNA binding (GO:0043565) ATP-dependent helicase activity (GO:0008026) purine NTP-dependent helicase activity (GO:0070035).

In the comparative combination of T6 and C2, the different genes were mainly enriched in hydrolase activity, acting on glycosyl bonds (GO:0016798) chitinase activity (GO:0004568) hydrolase activity, hydrolyzing O-glycosyl compounds (GO:0004553) ubiquitin-like protein transferase activity (GO:0019787) phosphatase activity (GO:0016791) phosphoric ester hydrolase activity (GO:0042578) O-methyltransferase activity (GO:0008171). And including aminoglycan metabolic process (GO:0006022) aminoglycan catabolic process (GO:0006026) chitin metabolic process (GO:0006030) chitin catabolic process (GO:0006032) amino sugar metabolic process (GO:0006040) cell wall macromolecule catabolic process (GO:0016998) cell wall macromolecule metabolic process (GO:0044036) amino sugar catabolic process (GO:0046348) drug catabolic process (GO:0042737) cell wall organization or biogenesis (GO:0071554) carbohydrate derivative catabolic process (GO:1901136) organonitrogen compound catabolic process (GO:1901565) in BP.

KEGG enrichment analysis of DEGs among different treatments

From the enrichment results of KEGG, the 20 most significant KEGG pathways under each temperature condition were selected to draw the scatter plot presented (Fig. 5). KEGG enrichment analysis fully revealed the biological pathways of DEGs among different treatments. Among all the pathways of DEGs enriched in T2 and C1 pair treatment, MAPK signaling pathway-plant (pxb04016), plant-pathogen interaction (pxb04626), plant hormone signal transduction (pxb04075), sesquiterpenoid and triterpenoid biosynthesis (pxb00909), glutathione metabolism (pxb00480), carotenoid biosynthesis (pxb00906), steroid biosynthesis (pxb00100), cysteine and methionine metabolism (pxb00270), alpha-linolenic acid metabolism (pxb00592), and fructose and mannose metabolism (pxb00051) pathways were significantly different (Fig. 5b). Among all the pathways of DEGs enriched in T3 and C1 pair treatment, phenylpropanoid biosynthesis (pxb00940) pathway was significantly different (Fig. 5c). Among all the pathways of DEGs enrichment in T5 and C2 pair treatment, carotenoid biosynthesis (pxb00906) pathway was significantly different (Fig. 5e). While the T1 and C1 pair treatment, T4 and C2 pair treatment, and the T6 and C2 pair treatment had no significantly enriched KEGG pathways.

Figure 5. 

The top 20 KEGG pathways selected in each pair treatment under different temperature conditions. The vertical axis represents the pathway name, the horizontal axis represents the rich factor, the size of dots represents the number of DEGs, and the padj is reflected by the color of the dots. (a) KEGG enrichment results of the T1 and C1 pair treatment. (b) KEGG enrichment results of the T2 and C1 pair treatment. (c) KEGG enrichment results of the T3 and C1 pair treatment. (d) KEGG enrichment results of the T4 and C2 pair treatment. (e) KEGG enrichment results of the T5 and C2 pair treatment. (f) KEGG enrichment results of the T6 and C2 pair treatment.

Key DEGs involved in MAPK signaling pathway

MAPK signaling pathways are involved in plant structure, growth and development, biological and abiotic stress responses^[23]. The MPK6 promoter of Arabidopsis thaliana has been found to strongly express YFP (yellow fluorescent protein) in flowers, highlighting the integral role of MPK6 in flower growth and development^[24]. In Arabidopsis thaliana, pollen tube growth depends on MPK3 and MAPK6 signal transductions^[25]. As shown in the (Fig. 6a), 18 DEGs were enriched on the MAPK pathway. As shown in the (Fig. 6a & Supplemental Fig. S1), compared with the control, the effects of three treatments on gene expression were generally consistent. MT treatment up-regulated the gene encoding MAPK3 (103948905) and the gene encoding MAPK18 (103943998). The expression of genes encoding type 2C serine/threonine protein phosphatases (PP2C) (103951896, 103943902, 103967224, and 103943417) were up-regulated. The transcription factor WRKY family genes 103962153, 103945640, 103958932, the gene 103936532 encoding end chitinase, the genes 103967071, 103934339, and 103934341 encoding the pathogenesis-related protein were up-regulated the gene 103964989 and 103928899 encoding calcium-binding protein and the gene 103931812 encoding ethylene transcription response factor were also up-regulated. Gene 103967068 encoding disease-related protein was down-regulated. These results indicate that genes of the transcription factor WRKY family and genes encoding MAPKs and PP2Cs may play a role in pear flower opening mediated by melatonin in lower temperature. In order to explore the interaction relationships of these genes, corresponding protein interaction prediction analysis was conducted (Fig. 6b). It showed that genes 103952175 and 103948905 may also be key nodes in regulating the genes responding to MT treatments.

Figure 6. 

Analysis of key DEGs in the MAPK signaling pathway and key DEGs in the plant hormone signal transduction pathway at lower temperature. (a) The expression of key DEGs in MAPK pathway in each MT treatment at lower temperature. Red color represents genes that were up-regulated, blue color represents genes that are down-regulated. (b) Analysis of protein interaction network in the MAPK signaling pathway. Orange color indicates the key node proteins. (c) The expression of key DEGs in plant hormone signal transduction pathway in each MT treatment at lower temperature. Red color represents genes that were up-regulated, blue color represents genes that are down-regulated. (d) Distribution map of expressed transcription factor families.

DEGs involved in plant hormone signal transduction pathway

Plant hormone signal transduction pathways also play important roles in plant growth and development^[26]. In plant hormone signaling pathway, a total of 15 DEGs in the T2 and C1 pair treatment were enriched in the six hormones (auxin, cytokinin, gibberellin, abscisic acid, ethylene) signaling pathways (Fig. 6c & Supplemental Fig. S2). Under lower temperature in 200 μM MT treatment, two auxin-responsive genes were identified, which were up-regulated gene 103961617 of GH3 gene family and up-regulated gene 103964554 encoding SAUR family protein (SAUR). In the case of cytokinin (CTK) signaling, an up-regulated gene 103952145 encoding the two-component response regulator A-ARR family (A-ARR) was identified. In the ABA signaling pathway, the expression of ABA response element binding factor (ABF) genes 103951896, 103943902, 103967224 and 103943417 were up-regulated. In the SA signaling pathway, the expressions of SA response element binding factor TGACG motif-binding protein family (TGA) genes 103947866 and 103957003 were up-regulated. In the ethylene signaling pathway, the expression of the ethylene response factor (ERF) gene 103931812 was also up-regulated. The genes 103967071, 103934339, 103934341 encoding pathogenesis-related protein (except 103967068) and the gene 103927787 encoding cyclin were also up-regulated. After 20 μM and 1,000 μM MT treatment, the expression trend of most genes was consistent with that after 200 μM treatment except for gene 103947866 and gene 103952145.

Expressed transcription factors analysis among these treatments

Numerous transcription factors (TFs) have been shown to play an important role in the stress response of plants. MYB transcription factors have been reported to play an important role in response to biological and abiotic stresses. WRKY transcription factor is involved in seed dormancy, plant growth and development, and aging process^[27]. As a group of plant-specific transcription factors, NAC transcription factors play an important role in plant growth and development, secondary metabolite synthesis and stress response^[28]. bZIP transcription factor plays an important role in plant growth and development, stress response, and regulation of biosynthesis of secondary metabolites^[29]. AP2/ERF is a plant-specific transcription factor that plays a key regulatory role in various stress response processes^[30]. In the T2 and C1 pair treatment, 9, 3, 1, 11, 12 and 18 DEGs of AP2/ERF, bZIP, bHLH, MYB, NAC, and WRKY family were found, respectively (Fig. 6d).

DEGs involved in carotenoid biosynthesis pathway

Carotenoids are a class of secondary metabolites widely distributed in nature and play an important role in the growth and development of plants^[31]. The diversity of color is caused by the differences in the types and contents of carotenoids in plant petals, fruits, and leaves^[32]. In this research, four genes were found enriched in the carotenoid synthesis pathway in pear flower buds under high temperature. Compared with the control, the expressions of three genes encoding 9-cis-epoxycarotenoid dioxygenase (NCED1) (103945979), 9-cis-epoxycarotenoid dioxygenase (NCED2) (103957334), abscisic acid 8'-hydroxylase (CYP707A2) (103955532) were up-regulated in nearly all the melatonin treatments (Fig. 7). Compared with the control, the expression of gene β-carotene hydroxylase (BCH) (103966173) was down-regulated under 20 μM and 200 μM MT treatment, while up-regulated under 1,000 μM MT treatment.

Figure 7. 

Analysis of key DEGs in the carotenoid biosynthesis pathway at higher temperature.* indicates significant difference of melatonin treatments and the control at the level of p < 0.05 .

Light and temperature are important environmental factors that regulate plant flower opening^[33]. In recent years, through the study of environmental temperature pathways, researchers have found that high temperature can promote flower opening^[34], while low temperature will inhibit flower opening. Temperature is a strong developmental signal. In the range of about 12−27 °C, there are significant differences in the growth rate and developmental response of plants, but there is no significant induction of stress response, which will not cause plants to produce stress response^[35]. Melatonin is an important small molecule, which plays an important role in plant growth and development, flowering, nutrient absorption, and fruit ripening^[36]. In addition, exogenous melatonin can also help plants overcome photoperiodic obstacles, thus changing the flower opening of plants^[18]. On apple trees, for example, melatonin applied externally before flower bud development delayed the flower opening of the apple tree in a dose-dependent manner^[22]. In this study, after the exogenous application of melatonin, the flower opening of pear was delayed by 20 μM and 200 μM melatonin at higher temperature compared with the control and the 1,000 μM melatonin treatment, while the flower opening of pear was all accelerated by the melatonin treatments at lower temperature (Fig. 1). The results indicated that melatonin could alleviate the effects of higher temperature and lower temperature on the flower opening of pear tree. However, the biochemistry and molecular mechanisms involved in this process remain to be further elucidated.

To reveal the molecular mechanism of melatonin treatment regulating flower opening of pear flower buds under different temperatures, we conducted transcriptomic analysis and found that DEGs in the T2 and C1 pair treatment in lower temperature, were mainly involved in the MAPK signaling pathway and the plant hormone signal transduction pathway. In horticultural plants, some studies have shown that MAPK-class connectivity plays an important role in coping with temperature stress^[37]. In tomato, SlMPK1 is a negative regulator of heat tolerance of tomato plants. Silencing SlMPK1 in transgenic tomatoes enhanced high-temperature tolerance, while transgenic tomatoes overexpressing SlMPK1 had lower high temperature tolerance, lower antioxidant enzyme activity, and higher H₂O₂ and MDA content. SlMPK3 is a low-temperature stress response gene. Overexpression of SlMPK3 increases the activity of antioxidant enzymes, enhances the levels of intracellular proline and soluble sugar, and enhances the resistance of plants under cold stress^[38]. Some studies have shown that MAPK genes in mulberry are also involved in the response to extreme temperatures^[39]. After cold treatment, some of the MAPK cascade genes in chrysanthemum plants were induced to express^[40]. Plant hormones play a very important role in plant flowering^[41]. The gibberellin pathway is one of the six classic pathways affecting flowering^[42]. Plant hormones such as gibberellin, jasmonic acid, abscisic acid, and auxin, have important effects on chromatin compaction mediated by DNA methylation and histone posttranslational modification, suggesting that epigenetic regulation may play a role in flowering through hormonal action^[41]. Melatonin is an important regulator of plant hormone-related genes expression. In Arabidopsis thaliana, melatonin inhibits DELLA protein and thus affects floral transformation^[21]. In the present study, 18 and 15 DEGs were found in the MAPK and the hormone signaling pathways, respectively, including up-regulated expression of WRKY, PP2C, bZIP, ERF, etc. (Fig. 6). At lower temperature, most of the genes in the MAPK and plant hormone transduction pathways of the melatonin treatment were up-regulated. Reports showed that PP2C is the core component of the ABA signaling pathway and is involved in the response to abiotic stress^[43], ERF participates in the defense response and alleviates the effects of stress^[44], and CaM4 can reduce the harm of reactive oxygen species and maintain the stability of the intracellular environment^[45]. As a strong antioxidant, melatonin can alleviate the damage caused by low temperature^[46]. Therefore, as shown in Supplemental Fig. S1, upregulation of MAPK18 and PP2C may improve the stress adaptation of pear flower buds, upregulation of ERF1 may strengthen the defense response of pear flower buds, and upregulation of CaM4 may be in involved in the maintenance of the homeostasis of reactive oxygen species of pear flower buds. Genes of GH3 and SAUR are involved in plant development. Studies have shown that they may be related to the development of floral organs^[47,48]. So as shown in Supplemental Fig. S2, upregulation of GH3 and SAUR may enhance cell enlargement and flower bud development, and upregulation of A-ARR may strengthen the cell division of pear flower buds. In a word, changes in these genes caused by melatonin treatments may alleviate the damage to pear flower buds caused by lower temperature and improve the development of pear flower buds, resulting in earlier flower opening after melatonin treatment. However, the underlying mechanism needs to be further studied.

At higher temperature, we found that DEGs in the T5 and C2 pair treatment were mainly involved in the carotenoid synthesis pathway. The carotenoids, as one of the natural pigments in plants^[49], play an important role in the formation of flower color. Studies have found that environmental conditions have a regulatory effect on the synthesis of carotenoids^[50,51]. In higher temperature, the expression levels of genes NCED1, NCED2, CYP707A2 were nearly all up-regulated in the three melatonin treatments, while the expression of gene BCH was down-regulated in the treatments of 20 μM and 200 μM MT and up-regulated in the treatment of 1,000 μM MT. β-carotene becomes zeaxanthin under the catalysis of BCH^[52,53]. The precursor of ABA de novo synthesis is zeaxanthin. Zeaxanthin is catalyzed by zeaxanthin epoxidase (ZEP) to produce antheraxanthin and violaxanthin. Violaxanthin and neoxanthin were catalyzed by NCED to produce, and xanthin finally becomes ABA after several steps of catalysis^[54]. The degradation of ABA enhances the formation of saffron acid (PA) in plants under catalysis of CYP707A^[55]. Studies have shown that carotenoid synthesis is closely related to the accumulation of ABA^[56]. Therefore, at higher temperature, the expression changes of NCED1, NCED2, CYP707A2 and BCH caused by melatonin treatment may be involved in the degradation of β-carotene and the synthesis and degradation of ABA, which changes the development of floral organs, finally resulting in the postponement of pear flower opening of the 20 μM and 200 μM treatments. But the underlying mechanism also needs to be further studied.

Combining all the data, we propose a probable model in which melatonin regulates the flower opening of pear tree through the MAPK and the hormone signaling pathways at lower temperature, and a rough model in which melatonin changes the flower opening of pear tree through the carotenoid synthesis pathway at higher temperature (Fig. 8). Our study can provide a theoretical basis for understanding and further study the molecular mechanism of exogenous melatonin regulating the flower opening of pear tree in lower or higher temperatures.

Figure 8. 

Possible models of temperature regulating flower opening time of pear trees mediated by melatonin. (a) Possible model of lower temperature regulating flower opening of pear trees mediated by melatonin. (b) Possible model of higher temperature regulating the flower opening of pear trees mediated by melatonin. The red and blue arrows in the figure represent up-regulated and down-regulated genes, respectively.

The authors confirm contribution to the paper as follows: study conception and design: Zhang H; data collection: Xu P, Wang X; analysis and interpretation of results: Xu P, Luo S, Zhang H, Xu J, Ma H, Zhang Y; draft manuscript preparation: Xu P, Cheng A, Zhang H. All authors reviewed the results and approved the final version of the manuscript.

The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.