Transcriptome analyses revealed chilling response genes in mango (Mangifera indica L. cv. Keitt) leaf

ABSTRACT Mangifera indica L. cv. Keitt is a cold-stress fruit plant native to China's drought river valley. Chilling stress affects productivity. Understanding the mechanisms of chilling stress is important to increasing chilling resistance in mango. Leaves of Keitti were subjected to 4 °C for 0, 3, 6, and 9 h for RNA-Seq-based transcriptome analysis, respectively. The chlorophyll content, carotenoid content, catalase, and peroxidase activities significantly increased during 9 h. The leaves responded to the stress by enhancing photosynthetic pigment content and antioxidant enzyme activity. After 3 h of chilling, 410 genes were differentially expressed. WRKY70 and PLD1 were significantly up-regulated after 9 h. Compared to 9 and 0 h, there were 1123 DEG. The DEGs are enriched in hormonal signal transduction, secondary metabolites, and the abiotic stress response. Similarly, the transcriptional factor families including NCED2, MYB73, and HLH162 up-regulated. The study will promote research on the development of chill-resistant mangoes.


Introduction
Increasing demand for fruit crops in cold regions is on the rise (Dinesh et al. 2012).The lower temperature is a key environmental factor that drastically impairs plant growth and significantly restricts productivity (Chaudhry and Sidhu 2021;Wongkaew et al. 2021).Low temperature limits the ability of chlorophyll production to function, and the plant suffers as a result to produce.There are two types of low temperature that affect plant growth and development; the chilling injury from 0°C to 20°C and the freezing injury, which falls below 0°C (Ding et al. 2019).In the environmental field, the low temperature restricts plant vegetative growth and delays floristic tissue production (Robinson et al. 2003).It causes leaf yellowing (Luo et al. 2019), leaf tip death (Martínez et al. 2019), and leaf curling (Tsai et al. 2019).In the case of reproduction, it caused the abscission of juvenile buds, flowers abortion, and plant cell death (Kiran et al. 2019;Li et al. 2017).To limit this possible radical damage, the plants develop several stressors that increased their ability to sustain their growth.The most common response mechanism is acclimation or adaptation, which mainly occurs through enhanced tolerance to environmental stress after preparation for treatment.
The tree crop mango originated in tropical and subtropical regions have lower low-temperature tolerance and thus are more susceptible to chilling stress (Giné-Bordonaba et al. 2020;Wei et al. 2017).As plants respond to stress their molecular, biochemical, and physiological pathways change (Xie et al. 2019).Transcriptional factors (TFs), induced or activated by the plant's perception of environmental signals, serve as central mediators of transcriptional reprogramming, which leads to the adaptation process (Tsai et al. 2019).Several TF families under plant hormonal signaling influence the regulation of the expression of chilling-responsive genes.WRKY genes are typically induced by salicylic acid and in turn control the expression of defense-related genes (Yu et al. 2001).In Arabidopsis thaliana and rice, WRKY TFs play important roles in transcriptional reprogramming during abiotic stresses such as chilling and freezing stress (Cardoso Filho et al. 2017;Li et al. 2013).For example, AtWRKY40 inhibited directly the expression of hormone-responsive genes that function as a negative regulator of Abscisic acid (ABA) signaling in a low-temperature stress (Chen et al. 2012).In Indica rice, an enhanced accumulation of ABA was associated with an enhanced expression of 9-cis-epoxycarotenoid dioxygenase (NCED) gene to function (Verma et al. 2019).The basic helix-loop-helix (bHLH) family also played an important role in multiple stress reactions.For instance, the expression of NtbHLH123 conferred tolerance to low temperatures in tobacco (Zhao et al. 2018).Molecular evidence indicated that the progressive expression of WRKY70 and Phospholipase D1 (PLD1) was responsible for low-temperature tolerance, where PLD1 mediated the process by initiating the lowtemperature shock response in plants (Zheng et al. 2016).In rice, osWRKY70 mediated continued low-temperature adaptation (Nordey et al. 2014;Yu et al. 2001).
Mango (Mangifera indica L. cv.Keitt) originated over 4000 years ago from tropical trees belonging to the flowering plant genus Mangifera in tropical regions sensitive to low temperatures, especially at germination, seedling stages, and before fruiting (Jagarlamudi et al. 2011).Keitt is a chilling-resistant mango cultivar that is widely grown in China's dry and hot Jingsha River valley, where the latitude reaches 26 (Dautt-Castro et al. 2015).It is important to provide a comprehensive interpretation of the transcriptomic changes in response to the chilling stress of Keitt.In this study, RNA-Seq was used to analyze the transcriptome patterns of 'Keitt' leaves under 0, 3, 6, and 9 h of 4 °C chilling stress.The analyzed results provide an overall view of the molecular regulatory network response to chilling stress in Keitt leaves.The differentially expressed genes could be candidates for further functional validation and application in breeding for versatile, low-temperature tolerant mango genotypes.

Plant material and chilling treatment
Mature functional leaves of Mangifera indica L. cv.Keitt were incubated in a 4 °C chilling chamber.Leaves were harvested at 0, 3, 6, and 9 h after 4 °C treatment and immediately immersed in liquid nitrogen before being stored at -80 °C for RNA extraction and biochemical analysis.Three replicates were used at each time point.

Measurement of biochemical parameters
On three biological experimental repeats, the peroxidase (POD) and catalase (CAT) activities, carotenoid (Caro), and chlorophyll (Chl) content of the leaves of Keitt after 0, 3, 6, and 9 h of 4°C chilling stress were quantified (Abdelaal et al. 2020).

RNA extraction, library construction, and sequencing
Total RNA was extracted using the Trizol reagent kit (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's protocol.The RNA quality was assessed using the Agilent 2100 Bioanalyzer (Agilent Technologies) and checked using RNase-free agarose gel electrophoresis.mRNA was enriched by Oligo (dT) beads (Epicentre, Madison).Then the enriched mRNA was fragmented into short fragments using fragmentation buffer and reverse transcript into cDNA with random primers.Second-strand cDNA was synthesized using DNA polymerase I, RNase H, dNTP, and buffer.Then the cDNA fragments were purified with the QiaQuick PCR extraction kit (Qiagen, Venlo), the ends repaired, an A base added, and ligated to Illumina sequencing adapters.The ligation products were size selected by agarose gel electrophoresis, PCR amplified, and sequenced using Illumina Novaseq6000 by Gene Denovo Biotechnology Co. (Guangzhou, China).The STEM tool was used for the analysis of a series of gene expression data (Ernst and Bar-Joseph 2006).

Bioinformatics analysis of RNA-Seq data
Reads were obtained from sequencing, including raw reads containing adapters or low-quality bases following assembly and analysis.Thus, to get high-quality clean reads, the reads were further filtered by fastp (version 0.18.0).The parameters read were obtained by removing reads containing adapters, reads containing more than 10% of unknown nucleotides (N), and low-quality reads containing more than 50% of low-quality (Q-value ≤ 20) bases.The short read alignment tool Bowtie2 (version 2.2.8) was used for mapping reads to the ribosome RNA (rRNA) database.The rRNA-mapped reads were then removed.The remaining clean reads were further used in assembly and gene abundance calculation.An index of the reference genome was built, and paired-end clean reads were mapped to the reference genome using HISAT 2.4 with '-RNA-Strandness RF' and other parameters set as a default.The annotated differentially expressed gene bases were obtained from the Gene Ontology (GO), Pathway Enrichment Analysis (Kyoto Encyclopedia of Genes and Genomes, KEGG), and Gene Set Enrichment Analysis (GESA) databases.The transcripts with a false discovery rate (FDR) below 0.05 and an absolute fold change of less than 2 were considered differentially expressed genes or transcripts.

Validation of RNA-Seq data by qPCR
A quantitative PCR (qPCR) was performed with six genes selected from DEGs of the significantly enriched pathways under chilling treatment.Specific primers were designed by primer 5, and alpha-tubulin was used as an internal control for analysis.The qPCR was performed with three replicates for each sample using a Bio-Rad CFX96 Real-Time instrument and the DNA Master SYBR Green kit.

Biochemical analysis of chilling exposed 'Keitt' leaf
The biochemical analysis revealed that the photosynthetic pigment content and antioxidant activities of the leaves were enhanced apparently after 3 h of the chilling stress, and they did not increase apparently from 6 h to 9 h.The increase of chlorophyll (Chl) and carotenoids (Car) content may increase the photosynthetic capacity of the leaves, resulting in the enhancement of the resistance.Compared to Chl, the content of the Car increased more (about 4 times) after 3 h of chilling treatment.It was suggested that the carotenoid responded more sensitively to chilling stress in Keitt.Carotenoids acted as a photo-protective of the photosynthesis apparatus (Solovchenko and Neverov 2017), and as a transfer agent of the light energy chlorophylls absorbed during the chilling time for photosynthesis activities (Xu and Harvey 2019).The apparent increase in the carotenoid content under the chilling time may protect the photosynthesis system and increase their resistance to chilling stress ability.
Chilling caused higher carotenoids content may result in a faster rate of carotenoid biosynthesis and increased carotenogenic gene expression (Gao et al. 2012).The enzyme activity of POD and CAT increased apparently after 3 h of chilling stress and reached its peak value at 6 h.The enzyme activity of POD responded more sensitively to chilling stress; it increased about four times after 3 h of chilling treatment.Chilling stress caused excessive tetrapyrrole intermediate accumulation and increased chlorophyll absorbance types affected the chlorophyll content in Zea mays (López-Martínez et al. 2015;Rodríguez et al. 2013).Early low-temperature stress increased peroxidase (POD) activity in tobacco seedlings (Xu et al. 2010).

RNA-Seq quality
RNA-Seq analysis of leaves after 0, 3, 6, and 9 h of 4 °C chilling treatments was carried out using Illumina-based read sequencing.The numbers of clean reads and highquality reads obtained from 12 libraries were listed in Table 1.Clean reads were identified in 99.42-99.57% of the raw data, and 89.24-90.26% of the reads were mapped to the reference genome (Table 2).In total, 72375 (K0), 73074 (K3), 73138 (K6), and 74148 (K9) genes were detected in our study with an FPKM value ≥1 (Figures 1  and 2A).
In addition, more than 66.62% of the genes were detected in all four group libraries.According to the reference genome, 26,510 known transcripts representing 76.79% and 1404 new transcripts were identified from all the samples.In the RNA-Seq analysis, it is important to note that the relative expression of a transcript is proportional to the number of cDNA fragments it originates.The transcription group data detection had a high sensitivity to chilling as a result of this.This means that their protein-encoding genes' expression level FPKM (fragments per kilobase million) has significant values across the mapped reads in order of the magnitude of the reads.
The dispersion degree of gene expression level distribution was average, and the overall gene abundance in different samples was good (Figure 2A), and more reads were produced from the sample sequenced at a larger depth.On the other hand, reducing the large data set to a small, informed, meaningful value set was also key.A violin plot was used to visualize the distribution of the data samples' probability densities.The median (white dot), interquartile range (black bar in the center of the violin), and lower/ upper adjacent values (the black lines stretched from the bar) showed higher probability density within the samples (Figure 2B).The Principal Component Analysis (PCA) was used to study the idea of a dimensionality reduction distance relationship between samples; the more similar the sample composition, the closer the distance reflected in the PCA.According to the PCA analysis, the three biological replicates of the chilling-exposed groups were clustered together, respectively (Figure 2C), showing the strength and direction of a linear relationship between the samples' variables with a correlation coefficient above 0.95, as shown in Figure 2(D).Based on the FPKM analysis of each sample, Figure 2 showed the distribution of different transcripts by expression maps.Generally, the maps were used to evaluate the difference between samples in database analysis, sequencing, and comparison coefficient.

Identification of DEGs in response to chilling stress
The assembled genes were analyzed with a False Discovery Rate (FDR) below 0.05 and absolute fold change |logFC| ≥ 2 as a cutoff considered as the differentially expressed genes/transcripts.
The distribution of the DEGs is shown in Figure 3.After 3 h of chilling stress, 253 up-regulated and 157 down-regulated genes were detected, while 94 up-regulated and 20 down-regulated genes were detected between KT3 and KT6, and 53 up-regulated and 4 down-regulated genes were detected between KT6 and KT9.After 9 h of chilling treatment, 875 and 878 genes were up-regulated, and 248 and 326 genes were down-regulated compared to 0 and 3 h, respectively.The distribution of the DEGs showed that the gene expression level responded to chilling stress quickly after 3 h of treatment, and the gene expression pattern after 9 h was quite different from that after 0 and 3 h.It was suggested that during the first 3 h of chilling treatment, some early resistance genes were differentially expressed, and after 9 h of chilling treatment, some late resistance genes were differentially expressed.

Functional enrichment analysis of the DGEs along chilling stress
The enriched pathways of DEGs in six data sets were presented in (Figure 4A).The enrichment degree of pathways was analyzed by enrichment factor.On integrating the numbers of DEGs in all the data sets, five major pathways related to the chilling acclimation mechanism were obtained.They were 'Plant-pathogen interaction' (ko04626), 'Plant hormone signal transduction' (ko04075), 'MAPK signaling pathway-plant' (ko04016), 'Phenylalanine metabolism' (ko00940), 'Carotenoid biosynthesis' (ko00906), 'Protein processing in endoplasmic reticulum' (ko04141), and 'Phosphatidylinositol signaling system' (ko04070).

Expression trend analysis of DEGs
To better understand the molecular mechanism response to chilling stress, expression trend analysis of DEGs was performed at 4 °C for 0, 3, 6, and 9 h.The gene expression trend could be divided into 20 profiles.The genes in profile 17 were early-stage chilling response genes; their expression levels increased significantly after 3 h of chilling stress and were then maintained at high levels.The genes in profile 10 were late-stage chilling response genes; their expression levels did not change significantly during the first 6 h of chilling stress and increased significantly after 9 h of chilling treatment.The genes in profile 19 were continuously up-regulated during the 9 h process, which may be important genes in the response to chilling stress (Figure 5).
Figure 6 depicts the KEGG and GO enrichment of DEGs in the 20 expression trend profiles.These continuously up-regulated genes in profile 19 were mainly enriched in plant hormone signal transduction, carotenoid biosynthesis, lipid metabolism, and then environmental adaptation.Within the DEGs of profiles, the 10 and 17 enrichment analyses showed that the majority of the genes in these trended areas responded to chemical stress, stimuli such as the GO factor and hormonal signal transduction, biosynthesis of secondary metabolism, and signaling pathways in the KEGG pathways.It was indicated that there were complex signal transduction changes in mango leaves, transmitting stress signals to enhance the biosynthesis of amino acids, lipids, and secondary metabolites to increase stress resistance.

Hormone signal transduction chilling-response genes
Plant hormones are not only involved in plant development; they also function critically in responding to stress.Several genes involved in hormone metabolism or the signal transduction process were identified as up-regulated in Keitt leaves in response to chilling stress.These genes were identified to be involved in the abscisic acid receptor, auxin response factor, and ethylene-response factor as a substantial direct role in leaves responses to 4 °C stress (see Supplementary Table 3).Invariably, these plant hormones are identified to respond to stress by regulating specific gene expression (Bari and Jones 2009).During the 9 h of chilling stress, Dehydration-responsive element-binding protein 2A (DREB2A), Auxin-induced protein 22B (AUX22B), Auxin response factor 16 (ARF16), Ethylene-responsive transcription factor 1B (ERF1B), PYL4, and PYL1 were up-regulated significantly compared to 0 h.This means that, as chilling stress period extended, the expression levels of plant hormones-related genes were increased.Auxins were an endogenous small molecule with an incredibly large impact on plant growth and response to stress in concert with other hormonal pathways (Liu et al. 2014).The Auxin Response Factors, AUX22B, ARF16 were key components of the auxin signaling pathway, and functioned as regulators response to chilling stress (Mahato et al. 2016;Wojcikowska and Gaj 2013).Auxin signal in a strawberry was severely induced after low-temperature treatment, indicating that auxin signal played important role in chilling stress (Zhang et al. 2019b).In Arabidopsis thaliana, the DNAbinding specificity of ERF/AP2 domain of DREBs was involved in low temperature-inducing gene expression (Sakuma et al. 2002).ERF1B acted as a transcriptional activator identified as a binding agent to GCC-box promoter elements in low-temperature-responsive genes regulating the expression of the gene during the chilling stress.The ERF1B integrated with ethylene signal transduction regulated ethylene/jasmonate-dependent defenses genes to function (Müller and Munné-Bosch 2015).JAs positively control chilling stress responses because it has been shown that JAs application increase the chilling tolerance of many plant species, including maize and Cucumis sativus (Chen et al. 2000;Wang et al. 2009).ABA signaling genes under normal conditions, had relatively lower expression levels, including PYL4 and PYL1 ABA receptor, after 1 h of cold stress,  expressed significantly higher (Deng et al. 2018).Previous studies showed that the biological function of OsPYLs, genes substantially improved low-temperature stress tolerance in rice (Kim et al. 2014).Moreover, the comprehensive evolution of OsPYL7, OsPYL9, and OsPYL gene functions in rice revealed significantly higher responses to low-temperature stress (Di et al. 2018).The overexpression of OsPYL3 substantially improved low-temperature stress tolerance in rice (Verma et al. 2019).

Transcriptional factors respond to chilling stress
Transcription factors (TFs) played key roles in the regulation of gene expression responses to low-temperature stress in plants (Shinozaki and Yamaguchi-Shinozaki 2000).The expression of TFs induced or activated by plant perception of environmental signals were central mediators of transcriptional reprogramming which leads to plant adaptation (Saibo et al. 2009).The TF families were differentially expressed during the 9 h chilling stress of 'Keitt' (see Supplementary S1, Table 4).Notably, large numbers of these TFs were found to be up-regulated significantly after 3, 6, and 9 h of chilling treatment.Under a low temperature of 4 °C, NCED2 was up-regulated rapidly within 3 h, which catalyzes ABA biosynthesis in the carotenoid biosynthesis pathway.In transgenic Indica rice, enhanced accumulation of ABA was associated with an enhanced expression of genes for ABA biosynthesis viz., ZEP1, NCED1, NCED2, NCED3, and NCED4 (Verma et al. 2019).NCEDs genes were involved in isoprenoid and chlorophyll metabolism, the up-regulation of NCED2 under chilling stress may increase chlorophyll accumulation as it enhanced chloroplast enlargement in Arabidopsis thaliana (Meier et al. 2011).The expression of ABA receptor genes PYL4 and PYL-1 was up-regulated continuously during the chilling treatment.PYLs were reported to function in the enhancement of the perception and transmission of the ABA signal.However, in the late stage of chilling treatment at 9 h, the expression of CYP707A2 was up-regulated, which was involved in the degradation of ABA.The expression of Ethylene Response Factors, including Ethylene-responsive transcription factor 1B (ERF1B), Endoplasmic Reticulum Nucleus Signaling 1 (ERN1), Ethylene-responsive transcription factor 61 (ERF061), Ethylene-responsive transcription factor 22 (ERF022), and Ethylene-responsive transcription factor (ERF054), was significantly up-regulated during the 9 h of chilling treatment.WRKY51 and WRKY70 were exclusively up-regulated in response to chilling treatment.WRKY70 was involved in adaptation to low-temperature stress and positively regulated the salicylic acid (SA)-mediated signal pathway but negatively regulated the jasmonic acid (JA)-mediated signal pathway in Arabidopsis thaliana (Divi et al. 2010).WRKY70 requires phytohormone signaling pathways activation via modulating the target gene to induce their expression (Mahato et al. 2016;Zhao et al. 2020).The role of VaWRKY12 from Vitis amurensis in low-temperature stress tolerance was confirmed through an ectopic overexpression in Arabidopsis thaliana and grapevine (Zhang et al. 2019a).At the seedling stage of Verbena bonariensis, VbWRKY32 was up-regulated under low-temperature stress, which antioxidant enzyme activities increased thereby improving the survival ability under chilling stress (Wang et al. 2020).In addition, members of the MYB family are pivotal factors in regulating responses to abiotic stress.The expression of MYB4 was rapidly induced by low-temperature  treatment in Leymus chinensis and positively modulates lowtemperature tolerance in Arabidopsis thaliana (Li et al. 2020).However, MYB4a of Vaccinium corymbosum transcript levels decreased gradually during the low-temperature treatment (Zhang et al. 2020).Moreover, MYB102 in tomato enhances antioxidant enzyme activity and promoted the expression of cold resistance under low-temperature stress conditions (Wang et al. 2020).The expressions of MYB genes (MYB61, MYB102, MYB73, MYB73, and PP2) in 'Keitt' were up-regulated continuously under the 9 h chilling stress.The bHLH family also plays important role in multiple stress responses.The bHLH162 can bind to the DRE/CRT regulatory element present in the promoters of target genes demonstrating a dependency of low-temperature acclimation on the photosynthetic process in pumpkin (Luo et al. 2021).The expression of NtbHLH123 confers tolerance to low temperatures in tobacco (Zhao et al. 2018).NtbHLH123 increased stress tolerance by improving the expression of several abiotic stress-responsive genes to mediate the ROS scavenging ability and other stress tolerance pathways (Zhao et al. 2018).The expression levels of PLD1 was significantly up-regulated under the 9 h of chilling stress.PLD1 gene was originally identified and was used in Arabidopsis thaliana for abiotic and biotic stress adaptation (Zhao et al. 2020).The function of AtPLD1 was to activate the expression of Hydrolyzed glycerol-phospholipids to generate phosphatidic acids to enhance the low-temperature resistance of Arabidopsis thaliana (Bezvoda 2014).In the early stage of the chilling response, PLD1 induced freezing tolerance in Arabidopsis thaliana and allows osmolyte accumulation until temperature was elevated to allow regrowth (Zhao et al. 2020).

qPCR validation of gene expression
To ensure the reliability of the RNA-Seq data, the expression patterns of six randomly selected DEGs were evaluated by qPCR (Figure 7).The expression profiles of all 6 genes showed the same trend between the qPCR and the RNA-Seq results, demonstrating that the RNA-Seq data was highly reliable for further analysis (see Supplementary 2, Table 5).MYB family transcription factor PHL7-like, fasciclin-like arabinogalactan protein 12, Peptide-N4-(N-acetyl-beta-glucosaminyl) asparagine amidase A-like, 4-coumarate-CoA ligaselike 5, probable N-acetyltransferase HLS1, and protein aspartic protease in guard cell 2.

Conclusion
Chilling stress has an important influence on plant growth and development, leading to physiological, biochemical, metabolic, and molecular changes.In the present study, we comprehensively analyzed the changes in Keitt leaves during 9 h of 4°C chilling stress.The results revealed that Keitt responded to chilling stress by increasing photosynthetic pigments (Chl and Car) and antioxidant enzyme activities at the physiological level.Under chilling stress, the phytohormones biosynthesis and plant hormone signal transduction pathways were activated immediately to regulate the amino acid metabolism, lipid metabolism, and secondary metabolism pathways to enhance resistance.Members of the ERF, WRKY, MYB, and bZIP gene families were significantly differently expressed during chilling stress, suggesting that these genes play important roles in chilling stress resistance in Keitt.These results could provide valuable information to reveal the molecular mechanism of the chilling stress resistance of Keitt.

Figure 1 .
Figure 1. Biochemical analysis of Mangifera indica L. cv.Keitt response to chilling at different times.(A), chlorophyll content, (B) carotenoid content, (C) CAT enzyme activity, (D) POD enzyme activity.The columns with the different letters mean significant difference (P < 0.05) between different times of exposure.Bar indicates standard error.

Figure 2 .
Figure 2. A representation of RNA-Seq gene expression distribution.(A) Fragments per kilobase million (FPKM) values across the orders of magnitude in the mapped reads.(B), Violin plot in the expression of the sample, (C) the principal component analysis (PCA) (D), the value of the relationship between the sample correlation coefficient.

Figure 3 .
Figure 3.The distribution of the differentially expressed genes in 'Keitt' leaves exposed to chilling stress.

Figure 4 .
Figure 4.The KEGG (A) and GO (B) enrichment of the DEGs between different chilling treatment time point.

Figure 6 .
Figure 6.KEGG and GO enrichment analysis of the DEGs in the 20 expression trend profiles.

Table 1 .
Statistics of RNA-seq reads obtained from 12 libraries.

Table 2 .
Statistics of RNA-seq reads obtained from 12 libraries compared to the reference genome.