Integrated Transcriptomic and Metabolomics Analyses Reveal Molecular Responses to Cold Stress in Coconut (Cocos nucifera L.) Seedlings

Coconut is an important tropical and subtropical fruit and oil crop severely affected by cold temperature, limiting its distribution and application. Thus, studying its low-temperature reaction mechanism is required to expand its cultivation range. We used growth morphology and physiological analyses to characterize the response of coconuts to 10, 20, and 30 d of low temperatures, combined with transcriptome and metabolome analysis. Low-temperature treatment significantly reduced the plant height and dry weight of coconut seedlings. The contents of soil and plant analyzer development (SPAD), soluble sugar (SS), soluble protein (SP), proline (Pro), and malondialdehyde (MDA) in leaves were significantly increased, along with the activities of superoxide dismutase (SOD), peroxidase (POD), and catalase (CAT), and the endogenous hormones abscisic acid (ABA), auxin (IAA), zeatin (ZR), and gibberellin (GA) contents. A large number of differentially expressed genes (DEGs) (9968) were detected under low-temperature conditions. Most DEGs were involved in mitogen-activated protein kinase (MAPK) signaling pathway-plant, plant hormone signal transduction, plant–pathogen interaction, biosynthesis of amino acids, amino sugar and nucleotide sugar metabolism, carbon metabolism, starch and sucrose metabolism, purine metabolism, and phenylpropanoid biosynthesis pathways. Transcription factors (TFs), including WRKY, AP2/ERF, HSF, bZIP, MYB, and bHLH families, were induced to significantly differentially express under cold stress. In addition, most genes associated with major cold-tolerance pathways, such as the ICE-CBF-COR, MAPK signaling, and endogenous hormones and their signaling pathways, were significantly up-regulated. Under low temperatures, a total of 205 differentially accumulated metabolites (DAMs) were enriched; 206 DAMs were in positive-ion mode and 97 in negative-ion mode, mainly including phenylpropanoids and polyketides, lipids and lipid-like molecules, benzenoids, organoheterocyclic compounds, organic oxygen compounds, organic acids and derivatives, nucleosides, nucleotides, and analogues. Comprehensive metabolome and transcriptome analysis revealed that the related genes and metabolites were mainly enriched in amino acid, flavonoid, carbohydrate, lipid, and nucleotide metabolism pathways under cold stress. Together, the results of this study provide important insights into the response of coconuts to cold stress, which will reveal the underlying molecular mechanisms and help in coconut screening and breeding.


Introduction
Coconut (Cocos nucifera L.) is a tropical fruit and oil cash crop that is essential on a global scale and is rich in nutrients, including fats, proteins, sugars, fatty acids, and amino acids [1,2]. Coconuts taste sweet. The flesh has the effects of tonifying deficiency, strengthening, clearing summer heat, and quenching thirst [3]. Coconut, in the palm family, is mainly planted in tropical and subtropical coastal areas [4,5]. The most suitable temperatures for growth are between 27 and 32 • C, with an average annual temperature of 29 • C. Even if the average annual temperature for one month is 18 • C, the yield decreases significantly. Average temperatures lower than 15 • C cause the fall of flowers and fruits, cracking of fruits, and yellowing of leaves [6,7]). Thus, low temperatures substantially affect the growth and yield of coconut. The demand for coconut has been increasing, prompting research into the mechanism of its cold resistance and the cultivation of a new cold-resistant variety [8].
Low temperature is a major environmental stress factor that substantially affects plant growth, development, quality, and geographical distribution [9][10][11] and can change morphological, biochemical, and physiological characteristics [10]. Cold stress impedes enzyme activity, photosynthesis, and other biochemical processes, leading to the accumulation of reactive oxygen species (ROS), which leads to oxidative damage and membrane instability [12]. Oxidative stress induced by cold stress is the product of the accumulation of ROS. ROS are produced for signaling purposes and are anoxic byproducts of aerobic metabolism. For specific responses to gene expression, ROS trigger signal transduction pathways, such as mitogen-activated protein kinase (MAPK) signal transduction pathways. Overly toxic ROS are cleared by enzymes, such as superoxide dismutase (SOD), catalase (CAT), and peroxidase (POD) [13,14]. During cold stress, the leakage of the plasma membrane destroys this enzyme system, leading to metabolic disorders. This further leads to a decrease in plant energy supply, inhibits photosynthesis and biomass production, and causes plant death [12,15,16].
To cope with cold stress, plants adopt various physiological, biochemical, and molecular regulations, including changes in membrane lipid composition, induction of enzyme activity, changes in gene expression, hormone control, signal transduction, and osmotic regulation mechanisms to maintain plant cell homeostasis [17]. Plant adaptation to cold stress varies by species [16]. Moreover, plants in tropical and subtropical regions are generally sensitive to cold stress [18,19].
Plant responses to cold stress involve complex regulatory networks that bring about a wide range of physiological and metabolic changes, such as the activation of SOD and the accumulation of soluble sugars and low-molecular-weight substances [20]. This further aids in ROS removal and relieves osmotic pressure to maintain cell homeostasis [21]. Under cold stress, a series of molecular changes occur, along with changes in metabolic processes [22]. When plants are subjected to low-temperature stress, several cold response genes participate in signal transduction, and regulatory pathways are induced, resulting in the accumulation of specific osmoregulatory substances [9]. Under low-temperature stress, a significant increase in the accumulation of protective substances, such as soluble sugar and protein, has been observed in many plants [23,24]. Studies have shown that soluble sugar levels in different plants increase significantly under low-temperature stress, such as in citrus (Citrus reticulata) [25], red spruce [23], and Pinus halepensis [26]. The concentration of proline is a critical factor in cold tolerance. As an osmoregulatory compound, proline affects wheat [27] and influences ROS homeostasis, plasma membrane integrity, and cold resistance of the crops. In addition, as a global signal transduction regulator, the plant hormone abscisic acid (ABA) stabilizes membrane structure, regulates stomatal movement, and controls osmotic stress tolerance in plants through transcriptional regulation of downstream stress-related genes [28]. The application of exogenous ABA enhanced the cold tolerance of Magnolia liliiflora [29] and tomatoes [30]. There is also increasing evidence that jasmonic acid (JA) is involved in the regulation of plant cold tolerance [31]. These reports show the complexity of plant responses to low temperatures. However, metabolic changes in response to cold stress vary by plant species [32].
Altered gene expression in response to low temperatures is a commonly used strategy for managing plant chilling injuries [33]. Many components involved in cold response signaling pathways have been isolated and characterized, including messenger molecules (e.g., Ca 2+ ), Ca 2+ -associated protein kinases, and several key transcription factors (TFs) [34]. One of the best-studied examples is the plant ICE-C-repeat/DREB binding factors (CBF)-COR signaling module [35]. C-cyclic peptide binding factor/dehydration reaction element binding protein (CBFs/DREB), a member of the AP2/ethylene reaction factor (AP2/ERF) family, plays a central role in cold acclimatization [36]. Transcription levels of CBF are substantially up-regulated by inducible factor CBF expression protein (ICE), a MYC-type alkaline helix-loop-helix family TF; CBF then activates downstream COR gene expression by binding to cis-elements in the promoter [37]. In addition, many cold-resistant proteins and protective substances (e.g., soluble sugar and proline) are synthesized in plant cells, which regulate osmotic potential and maintain membrane integrity [38,39].
Transcriptome analysis is a highly effective method to determine cold response genes in many crops, such as Brassica napus [40] and Brassica juncea [41]. With the rapid development of high-throughput sequencing and mass spectrometry technologies, advances in transcriptomics and metabolomics have enabled global molecular and physicochemical analysis of the regulatory mechanisms of plant cold tolerance [22,42]. Transcriptome changes using RNA sequencing (NA-SEQ) have revealed plant responses to cold stress in Arabidopsis [43], winter rapeseed [20], rapeseed [40], peanut [44], and Argyranthemum frutescens [45]; WRKY, NAC, MYB, AP2/ERF, and bZIP have been reported as the most abundant TF families of low-temperature stress in many species [40,46]. In addition, various metabolites, including ABA, carbohydrates, amino acids, and polyamines, have been used to investigate cold-stress responses in several species [44,[47][48][49]. Transcriptomic and metabolomic analysis and mining of these big data facilitate the exploration of the mechanism of plant cold resistance [9,50]. Comparative transcriptomic analysis of two peanut cultivars (NH5 and FH18) with contrasting cold resistance identified several cold response TFs, including bHLH, MYB, and NAC [51]. In addition, lipidomics analysis of these two varieties indicated that membrane lipid and fatty acid metabolism contributed indirectly to the cold resistance of peanut [52].
A comprehensive analysis of transcriptomics and metabolomics is essential for understanding the complex regulatory networks involved in plant cold-stress responses. Under cold stress, the combined transcriptomics and metabolomics of Xanthoceras sorbifolia showed enhanced metabolism of amino acids and sugars [32]. Low-temperature stress induced significant changes in the transcriptome and metabolome. Key pathways related to ABA/JA signaling and proline biosynthesis play an important role in regulating cold resistance in wheat [9]. Most differentially expressed genes (DEGs) and differentially accumulated metabolites (DAMs) are mainly enriched in different carbohydrates and amino acid metabolisms. Among them, starch, sucrose, and phenylalanine metabolism were significantly enriched, which played a crucial role in the adaptation of Brassica napus to cold stress [53]. Most DEGs are involved in amino acid biosynthesis, plant hormone signaling, and MAPK signaling pathways. In addition, metabolomics analysis showed that the contents of free polyamines (PAs), plant hormones, and osmotic solution mainly included increased putrescine, spermidine, spermidine, ABA, JA, raffinose, and proline when pepper was responding to cold stress. Importantly, regulation of the ICE-CBF-COR pathway through Ca 2+ , MAPK, and ROS signaling plays a crucial role in regulating the response of pepper to cold stress [54]. Moreover, comprehensive analysis of the transcriptome and metabolome showed that most of the genes and metabolites involved in carbohydrate metabolism, the TCA cycle, and flavonoid biosynthesis were up-regulated in cold-tolerant pepper cultivars under cold stress [55]. This comprehensive approach has also been widely used to study the cold-stress response of peanut [44], pumpkin [56], and tobacco [10,57] crops.
Despite the substantial progress in understanding cold-stress responses in food crops and ornamentals, knowledge of this response in palm woody oil crops is limited. In plant biology, it is crucial to have a comprehensive understanding of the cold-stress responses of various plant species. Studies have examined cold-stress responses of the woody oil crop of the palm family, oil palm [58][59][60], and the cold resistance of coconut [61,62]. However, the molecular mechanisms underlying the signaling pathways and related gene networks require further research owing to the complexity of cold resistance traits in crops.
In this study, we aimed to explore the transcriptomic and metabolomic responses of coconut seedlings under cold stress and identify the relationship between gene expression and corresponding metabolite accumulation and depletion. Moreover, we aimed to evaluate the main growth and physiological indexes of coconut seedlings, explore the transcriptome and metabolome changes in coconut under low-temperature stress, identify the signaling pathways and regulatory networks related to cold tolerance, and discuss the response mechanism of the coconut under cold stress. Identifying the underlying molecular mechanisms of the coconut hypothermia response will facilitate the development of future coconut hardy varieties.

Physiological Differences in Response to Cold Stress
The effects of cold stress on the growth of coconut seedlings were evaluated using pot experiments, in which seedlings were exposed to cold stress (5 • C) (low-temperature treatment [LT]) and room temperature (25 • C) (the control [CK]). The LT directly affected the growth of coconut seedlings ( Figure 1a) and significantly reduced their plant height and dry weight (p < 0.05): after 30 d of the LT, the plant height and dry weight of coconut seedlings were 31.08% and 49.07% less, respectively, than those of the CK (Table S1). In addition, under cold stress, the soil and plant analyzer development (SPAD) value was significantly lower (55.14%) than that of the CK (p < 0.01) (Figure 1b). SS and SP contents increased significantly under cold stress, peaked on day 20, and gradually decreased from day 20 to 30, generally increasing by 41.33% and 41.38%, respectively (p < 0.05) (Figure 1c,d). Thus, cold stress seriously affected the growth, photosynthetic index, and quality of coconut seedlings and had significant effects on the physiological characteristics of the coconut seedlings ( Figure 2). Under cold stress, malondialdehyde (MDA) and PRO contents increased significantly; they showed trends similar to those of prolonged cold stress, with an increase of 34.64% and 31.42%, respectively (p < 0.05) (Figure 2a,b). The activities of CAT, POD, and SOD in the leaves of coconut seedlings under cold stress were significantly higher (p < 0.05), 25.35%, 48.73%, and 33.08%, respectively, than those of CK, and under the prolonged cold stress, they displayed similar trends. However, the increase in SOD activity was lower than that in POD and CAT (Figure 2c-e). In addition, cold stress had significant effects on endogenous hormone levels in the leaves of coconut seedlings. Under processes of cold stress, ABA content increased slowly within 10 d of the LT, but rapidly from day 20 to 30. During the LT, ABA content increased significantly by 51.95% (p < 0.01) (Figure 2f), and the contents of IAA, ZR, and GA increased within 20 d of the LT, peaked on 20 d, and gradually decreased from days 20 to 30. The contents of IAA, ZR, and GA after LT were thus significantly higher, by 27.43%, 24.57%, and 22.02%, respectively (p < 0.05) (Figure 2g-i), than those of CK. Thus, cold stress seriously affected the physiological function indexes and hormone levels of coconut seedlings.

Evaluation of Transcriptome Sequencing Data
After sequencing quality control, 38.42 GB of clean data were obtained from six samples (i.e., three biological duplicate samples from two experimental treatments), with the percentage of Q30 bases in each sample no less than 93.48%. According to the statistical results, the efficiency of comparison between the reads of each sample and the reference genome ranged between 94.04% and 94.87% (Table S2). Fragments per kilobase of transcript per million fragments mapped (FPKM) > 1 was used as the threshold to determine gene expression, and the FPKM value in the CK 30 samples was higher than that in the LT 30 samples ( Figure S1a). The correlation analysis of FPKM values between the LT 30 and CK 30 samples is shown in Figure S1b. Principal component analysis (PCA) showed that the LT 30 and CK 30 samples were clustered separately, indicating significant differences in gene expression between sample groups. The three bioreplicated samples under the LT 30 and CK 30 conditions were strictly clustered together, indicating the high biorepeatability of the samples treated in each group ( Figure S1c). Figure S1d shows the volcano maps with significantly up-regulated and down-regulated differences between the two groups.
Int. J. Mol. Sci. 2023, 24, x FOR PEER REVIEW 5 of Figure 1. Different physiological responses of coconut seedlings to cold stress. (a) Experimen photos of coconut seedlings. (b) Changes in SPAD of coconut at 0, 10, 20, and 30 days in LT and C (c,d) SP and SS contents of coconut at 0, 10, 20, and 30 days in LT and CK. The error bar represe SD of the mean of the three biological replicates. * asterisks represent significant differences betwe the two treatments determined by Student's t-test (* p < 0.05; ** p < 0.01). (c,d) SP and SS contents of coconut at 0, 10, 20, and 30 days in LT and CK. The error bar represents SD of the mean of the three biological replicates. * asterisks represent significant differences between the two treatments determined by Student's t-test (* p < 0.05; ** p < 0.01). (c-e) The SOD, POD, and CAT activities of coconut seedlings at 0, 10, 20, and 30 days in LT and CK. (f-i) The ABA, IAA, ZR, and GA contents of coconut seedlings at 0, 10, 20, and 30 days in LT and CK. * asterisk indicates a significant difference between the two treatments determined by Student's t-test (* p < 0.05; ** p < 0.01).

Evaluation of Transcriptome Sequencing Data
After sequencing quality control, 38.42 GB of clean data were obtained from six samples (i.e., three biological duplicate samples from two experimental treatments), with the percentage of Q30 bases in each sample no less than 93.48%. According to the statistical results, the efficiency of comparison between the reads of each sample and the reference genome ranged between 94.04% and 94.87% (Table S2). Fragments per kilobase of transcript per million fragments mapped (FPKM) > 1 was used as the threshold to determine gene expression, and the FPKM value in the CK30 samples was higher than that in the LT30 samples ( Figure S1a). The correlation analysis of FPKM values between the LT30 and CK30 samples is shown in Figure S1b. Principal component analysis (PCA) showed that the LT30 and CK30 samples were clustered separately, indicating significant differences in gene expression between sample groups. The three bioreplicated samples under the LT30 and CK30 conditions were strictly clustered together, indicating the high biorepeatability of the samples treated in each group ( Figure S1c). Figure S1d shows the volcano maps with significantly up-regulated and down-regulated differences between the two groups.  (red legend) and CK (blue legend). (c-e) The SOD, POD, and CAT activities of coconut seedlings at 0, 10, 20, and 30 days in LT and CK. (f-i) The ABA, IAA, ZR, and GA contents of coconut seedlings at 0, 10, 20, and 30 days in LT and CK. * asterisk indicates a significant difference between the two treatments determined by Student's t-test (* p < 0.05; ** p < 0.01).

The Core Genes and Transcription Factors in Response to Cold Stress in Coconut
ICE-CBF-COR is the most important reaction pathway for plants at low temperature. The genes involved in cold−regulated protein, cold-responsive protein kinase, and temperature-induced lipocalin−1 were up-regulated under cold stress ( Figure 6). One TF ICE1(SCRM) gene encoding COCN_GLEAN_10023355 was up-regulated; however, another TF ICE1(SCRM) gene encoding COCN_GLEAN_10008268 was down-regulated. Six DREB/CBF transcription factors in the field of dehydration-responsive elementbinding protein (DREB1E (COCN_GLEAN_10009450), expressions of COCN_GLEAN_10008816 and COCN_GLEAN_10001620), DREB1C, DREB2C, and

qRT-PCR Validation
Nineteen DEGs involved in the MAPK signaling pathway-plant, amino sugar and nucleotide sugar metabolism, plant hormone signal transduction, carbon metabolism, starch and sucrose metabolism, and biosynthesis of amino acids were selected to verify the accuracy of the RNA-seq data using qRT-PCR. Under low temperatures, the expressions of the ETR2, MPK5, and WRKY24 genes associated with the MAPK signaling pathway-plant were significantly up-regulated, and BHLH25 was significantly down-regulated. The GAE1 and UXS2 genes related to amino sugar and nucleotide sugar metabolism were significantly up-regulated, and BHLH35 was significantly down-regulated. The XTH22, DPBF3, and SCL9 genes related to plant hormone signal transduction were significantly up-regulated, and IAA10 was significantly down-regulated. Expressions of the CYSC and SHM7 genes related to carbon metabolism were significantly up-regulated. The GLU3, GLC1, and SUS4 genes related to starch and sucrose metabolism were significantly up-regulated. The expressions of the FBA5 and CYSC genes related to the biosynthesis of amino acids were significantly up-regulated, and those of HAPS-1 were significantly down-regulated. These results suggest that a good correlation with FPKM values was obtained from the RNA-seq data (R 2 = 0.987), which confirms the reliability of gene expression data measured by RNA-seq ( Figure S6).

Metabolite Profiles of Coconut in Response to Cold Stress
Metabolomics analysis was performed to characterize metabolites of coconut in response to cold stress. A total of 423 metabolites in positive-ion mode and 225 metabolites in negative-ion mode were identified in the comparison of cold treatment and control samples (Table S7). DAMs were defined by the criteria of |log 2 FC| ≥ 1, p < 0.05, and VIP > 1. In the positive-ion model, a total of 206 DAMs were identified, while in the negative-ion mode, a total of 97 DAMs. Heat maps of DAMs between the two groups also revealed responses to cold stress at the metabolite level and successfully divided the test samples into two broad categories (Figure 7a,e). According to principal component analysis (PCA), PC1, PC2, and PC3 accounted for 64.69%, 13.32%, and 10.64% of metabolites in the positive-ion mode, respectively (Figure 7b), and for 64.40%, 15.66%, and 9.16% of metabolites in the negative-ion mode, respectively (Figure 7f), indicating a significant separation between the LT 30 and CK 30 groups. In the OPLS-DA analysis, the separation of the two groups was considerably good ( Figure S7). This indicates good repeatability between the three experimental (LT 30 ) or three control (CK 30 ) samples and a significant degree of separation between LT and CK samples. In the positive-ion mode, a total of 95 (up-regulated)/111 (down-regulated) DAMs ( Figure S7, Table S8a) were identified, including phenylpropanoids and polyketides (16%); lipids and lipidlike molecules (19%); benzenoids (9%); organoheterocyclic compounds (14%); organic oxygen compounds (8%); organic acids and derivatives (12%); nucleosides, nucleotides, and analogues (4%); organic nitrogen compounds (2%); alkaloids and derivatives (1%); and others (15%) (Figure 7d, Table S8a). A total of 38 (up-regulated)/59 (down-regulated) DAMs were identified in the negative-ion mode ( Figure S7, Table S8b), including phenylpropanoids and polyketides (25%); lipids and lipid-like molecules (8%); benzenoids (9%); organoheterocyclic compounds (12%); organic oxygen compounds (11%); organic acids and derivatives (3%); nucleosides, nucleotides, and analogues (6%); and others (14%) (Figure 7h, Table S8b). Each point represents a KEGG pathway, and the X-axis is Rich_factor, which represents the ratio of differential metabolites annotated to a pathway to the ratio of all metabolites annotated to that pathway. The Y-axis is the path name. The shade of the dot represents p value. The size of the circle indicates the number of differential metabolites enriched in the pathway. The differential metabolites in the pathway are up-regulated, represented by the lower triangle, and the differential metabolites in the pathway are down-regulated, represented by the upper triangle. There are both up-regulated and down-regulated differential metabolites in the pathway, which are represented by circles. (d,h) HMDB database classification summary of metabolites in the positiveand negative-ion modes. Items in the same box in the figure represent HMDB level classification information, corresponding to the super class and class information of the HMDB database. Column length represents the number of metabolites annotated by the classification.

Amino Acid Metabolism Considerably Contributes to Cold Stress in Coconut
In plants, some amino acids, such as proline and arginine, as osmotic substances, are critical regulatory elements in stress responses. Interestingly, KEGG pathways involved in amino acid metabolism were found to be abundant in the expression levels of both metabolites and genes in LT30 versus CK30 under cold stress. Based on the correlation analysis between the two groups (p < 0.05) and the enrichment KEGG pathway of genes and metabolites, we constructed the correlation network of DEGs and DAMs that may be involved in cold stress (Figures 9 and S10a, Table S11). In the lysine biosynthesis pathway, the genes involved in amino acid metabolism encoding AGD2, LOC105058024, COCN_GLEAN_10015782, DHDPS2, AK1, DAPF, and DHDPS2 were down-regulated under cold stress, and were significantly negatively correlated with the up-regulated metabolites, such as L-Lysine (+) and L-Saccharopine (+), but were significantly positively correlated with the down-regulated L-2-Aminoadipate adenylate (+). However, the upregulated expression of the LYSA gene was significantly positively correlated with L-Lysine (+) and L-Saccharopine (+) and significantly negatively correlated with L-2-Aminoadipate adenylate (+). In lysine degradation, 29 DEGs were detected that were involved in amino acid metabolism, and 15 DEGs were up-regulated, also showing significant positive correlation with L-Lysine (+) and L-Saccharopine (+). Among them, the genes encoding FRK1, ASHR3, LOC103697661, CCACVL1_23055, PAP4, ASHR3, and LKR/SDH were significantly up-regulated, while 14 DEGs were down-regulated and were significantly negatively correlated with L-Lysine (+) and L-Saccharopine (+), in which the genes encoding TCTP and ALDH3F1 were significantly down-regulated. In tyrosine

Discussion
Plants resist abiotic stress through various physiological and metabolic regulations [10,32,57]. In this study, we performed physiological, transcriptomic, and metabolomic analyses to uncover the regulatory network of coconut's response to cold stress. The physiological changes in coconut primarily occurred at lower temperatures; significant changes were observed in metabolite accumulation and gene expression using metabolomics and RNA-Seq analysis.

Physiological Response of Coconut under Cold Stress
Intracellular ROS was believed to play a decisive role in regulating signal transduction events. Abiotic stress led to the accumulation of ROS [63]. Excessive ROS enhances membrane lipid peroxidation and destroys membrane fluidity, resulting in electrolyte leakage (EL) [64], thus destabilizing the plasma membrane of plant cells [65]. Malondialdehyde (MDA) is a lipid peroxide derivative, production of which can aggravate membrane damage [66]. The increase in MDA level indicates the degree of damage to the cell membranes. For example, low-temperature treatment increased the MDA concentration and EL value of yellow horn (Xanthoceras sorbifolia), indicating that cold stress may lead to plasma membrane destruction [32]. However, under abiotic stress conditions, ROS homeostasis largely depends on the ROS clearance system [67]. The activity of several antioxidant enzymes (including SOD, CAT, and POD) has been found as a trigger to reduce ROS damage to plant cells, with SOD being an important ROS scavenging enzyme [68]. In this study, we observed an increase in MDA content (Figure 2a), indicating that low temperatures have a considerable effect on coconut seedlings. The activity of reactive oxygen scavenging enzymes (including SOD, CAT, and POD) was reportedly elevated (Figure 2c-e), indicating that they remove ROS produced by cold stress. SOD is a scavenger of peroxide anions, which can disproportionate the anions into H 2 O 2 and O 2 [63]. At the same time, POD and CAT catalyze the conversion of H 2 O 2 to oxygen and water [69]. The increase in SOD activity was found to be lower than that of POD and CAT, which may be due to the higher decomposition ability of POD and CAT to H 2 O 2 produced by SOD. In conclusion, the enhancement of cold-induced enzyme activity expression may improve the detoxification ability of reactive oxygen species to oxidative stress in coconut. Soluble sugars and proteins are considered to be key osmoregulatory substances in plants, and their accumulation in cytosol can prevent protoplast dehydration and improve the plant's cold resistance [24,25]. In this study, the contents of soluble sugar and protein reportedly increased significantly under cold stress, reaching a peak on day 20, and decreased gradually from days 20 to 30 (Figure 1c,d), indicating that soluble sugar and protein positively regulate the damage to the coconut cell membrane at low temperatures and reduce it with prolonged cold stress. This decreased content may be attributed to the inhibition of photosynthesis in the late stage of cold stress [16], which substantiates the significant decline in SPAD (Figure 1b). Proline (PRO) is an important amino acid in osmoregulation, which also helps to stabilize subcellular structures, scour free radicals, and buffer cell reductant-oxidant potential under stress conditions [70]. Rapid decompression of PRO after decompression can provide numerous reducing agents, which helps to support mitochondrial oxidative phosphorylation and ATP production. Recovery from stress and repair of stress-induced damage are also important indicators in the study of plant cold resistance [71]. Similar to previous studies, PRO levels found in cold-stressed coconuts were higher than the control (Figure 2b). This suggests that the high level of PRO in coconut can provide a large amount of PRO reducing agents, thus contributing to energy production, which is beneficial for plants to recover from adversity and repair damage caused by stress.

Transcription Factor in Response to Cold Stress in Coconut
Transcription factors play an important role in plant growth and development, and also participate in the plant transcriptional network in response to abiotic stresses [72]. The TF families related to plant stress resistance mainly include AP2/EREBP, MYB, WRKY, and bHLH [73]. Consistent with findings from other plants, our transcriptomic data showed that most differentially expressed TFs belonged to the WRKY, AP2/ERF, HSF, bZIP, MYB, and bHLH families. Previous studies have shown the diversity and complexity of regulatory pathways of TFs in response to cold stress [74]. The AP2/ERF family has been divided into the AP2/ERF, RAV, and DREB subfamilies, which play an important role in multiple stress responses [75]. For example, overexpression of the rice OsDRE1F gene can improve the tolerance of Arabidopsis and rice to low temperature [76], while the birch BpERF13 gene can improve the tolerance of Arabidopsis and rice to low temperature by up-regulating the CBF gene and reducing the number of ROS, thus demonstrating the positive regulation of cold stress [77]. In this study, AP2/ERF transcription factors were up-regulated in 74% of AP2/ERF TFs and down-regulated in 26% of AP2/ERF TFs under cold stress ( Figure 6, Table S6). These results are consistent with previous reports of cold-stress regulation caused by the presence of the AP2/ERF family [78,79]. WRKY-TFs represented a valuable family for resistance to abiotic stresses, such as cold, heat, and salt. For example, a total of 17 members of the family were induced by refrigeration in japonica rice [78]. In addition, a group of WRKY genes were also found to be cold-regulated in indica rice, and 18 WRKY genes were fully up-regulated in cold-tolerant varieties genotypes [80]. In addition, for the VbWRKY32 in the seed stage of Verbena bonariensis, the transcription level of cold response genes were up-regulated, thereby increasing the activity of antioxidant enzymes and the content of osmoregulatory subunits, thus improving the survival ability under cold stress [81]. In this study, 67% of WRKY TFs were up-regulated and 33% (14 DEGs) were down-regulated after low-temperature treatment ( Figure 6, Table S6). The expressions of 43 WRKY DEGs under cold stress may be helpful to further study the mechanism of their cold resistance. In addition, MYB family members have been shown to be key factors in regulating responses to abiotic stress. Some MYB TFs are involved in the regulation of cold stress. For example, LcMYB4 expression in Leymus chinensis was rapidly induced by cold treatment and actively regulated the cold tolerance of Arabidopsis [82]. Similar to previous studies, MYB genes in coconut were induced by the down-regulated and up-regulated mode under cold stress. For example, 38% of MYB TFs were up-regulated and 62% (eight DEGs) were down-regulated ( Figure 6, Table S6). In this study, the transcription of heat shock transcription factors (HSFs) in coconut was induced by cold stress, and 89% of HSF TFs were up-regulated ( Figure 6, Table S6). Some reports have shown that the abundance of HSF TFs (HSFA4A, HSFA6B, HSFA8, and HSFC1) in Arabidopsis Thaliana was enhanced by cold stress, while induction was diminished in ice1 mutants, suggesting that HSF was involved in the cold acclimation pathway [83]. In addition, overexpression of HSF TFs can stimulate the synthesis of protective metabolites, such as galactosol-affinity sugars, to improve plant tolerance to abiotic stresses [84]. Therefore, it was proposed that the relatively high transcription level of HSF TFs in coconut under cold stress may be more conducive to its response to cold stress. In addition, bZIP and bHLH TFs have also been reported to be involved in stress response [40,85]. In this study, 80% of bZIP TFs were upregulated. In the bHLH family, two DEGs encoding BHLH41 (COCN_GLEAN_10013682, COCN_GLEAN_10011130) were significantly up-regulated. However, one DEG encoding BHLH41 (COCN_GLEAN_10013779) was significantly down-regulated ( Figure 6, Table S6). In our results, major TF families, such as AP2/ERF, WRKY, MYB, HSF, bZIP, and bHLH, identified a large number of DEGs that induced expressions under cold stress. Therefore, these TFs have a significant effect on the cold resistance of coconut and may be involved in complex mechanisms related to low-temperature regulation.

The Core Regulator Genes in Response to Cold Stress in Coconut
In many plants, the ICE-CBF-COR signaling pathway is the most important signaling pathway in response to cold stress. This pathway is regulated by CBF/DREB transcription factors to induce low-temperature tolerance [86]. Numerous genetic and molecular analyses have identified C-repeat/DREB binding factors (CBFs) as key transcription factors that play a role in cold domestication and are essential in higher plants [87]. Overexpression of CBF could induce COR expression and improves frost resistance [88]. In this study, the expression of cold-regulated protein, cold-responsive protein kinase, and temperatureinduced lipocalin-1 genes were up-regulated under low-temperature stress. A TF ICE1 (SCRM) encoding COCN_GLEAN_10023355 was up-regulated; however, another TF ICE1 (SCRM) encoding COCN_GLEAN_10008268 was down-regulated. Six DREB/CBF transcription factors (DREB1E (COCN_GLEAN_10009450, COCN_GLEAN_10008816, and COCN_GLEAN_10001620), DREB1C, DREB2C, and DREB1G were significantly up-regulated, while two DREB/CBF transcription factors (DREB3 and DREB1F) were significantly downregulated ( Figure 6, Table S6); in addition, some genes of aquaporins (TIP1-1, PIP2-4, SIP2-1, PIP2-2, NIP1-1, PIP1-2) that were bona fide participants in the cold response at the molecular levels [89,90] were significantly down-regulated, and the gene of aquaporin (SIP1-2) was significantly up-regulated (Table S6), which not only confirmed the universality of the ICE-CBF-COR signaling pathway in response to low temperatures, but also proved that it is helpful to further study the molecular mechanism of this pathway mediating the cold tolerance of coconut.
Analysis of transcriptome data and KEGG of DEGs showed that many genes are involved in signal perception (MAPK signaling pathway-plant), transduction, and regulation under low-temperature stress. After the plant received the cold signal, the membrane permeability of the Ca 2+ receptor was enhanced, allowing more Ca 2+ to enter the plasma, resulting in rapid accumulation of Ca 2+ in the cytosol. Thus, Ca 2+ is an important second messenger that plays a key role in response to cold stress [91]. For example, in mesophyll cells of Arabidopsis, cold stress can lead to an immediate increase in cytoplasmic free Ca 2+ concentration and activation of Ca 2+ permeable channels [92]. Ca 2+ sensors sense changes in intracellular Ca 2+ levels through phosphorylation and then transduce signals to turn on downstream signaling processes for cold-specific gene expression, which helps plants adapt to cold stress [93]. Ca 2+ acts as a secondary messenger in response to cold stress and is recognized by Ca-binding proteins. The main intracellular Ca 2+ sensors in plants are calmodulin (or calmodulin-binding protein)/calcium-binding protein (or calmodulin-like protein)/calmodulin-binding protein (CAM/CML/CBP), calciumdependent protein kinases (CPK), and calcineurin B-like proteins/CBL-interacting protein kinases (CBL/CIPK) [9,94]. During this study, in the calcium-binding proteins (CMLs), 85% of DEGs were up-regulated. In calmodulin-related proteins (CAM/CML/CBP), 75% of DEGs were up-regulated. However, in calcium-dependent protein kinases (CPKs), 70% of DEGs were down-regulated. Moreover, in calcineurin B-like protein (CBL), two DEGs encoding CBL9 (CUFF17.251.1) and CBL7 were up-regulated, but three encoding CBL1 (CUFF34.779.2, CUFF5.369.1) and CBL9 (COCN_GLEAN_10019851) were downregulated. In CBL-interacting protein kinase (CIPK), 75% of DEGs were up-regulated ( Figure 6, Table S6). Similarly, some studies also suggested that the primary genes involved in intracellular Ca 2+ sensors are involved in Ca 2+ signaling and promote an increase in intracellular Ca 2+ concentration to activate various transcriptional cascades and potential candidate genes that activate downstream COR gene expression [95].
Endogenous hormones and their signaling pathways play a key role in regulating plant defense mechanisms against various biological stresses [100,101]. ABA accumulation is a key event in abiotic stress response [102]. ABA had been reported to accumulate in many species, including wheat, pepper, and bananas, all of which exhibit increased expression of ABA biosynthesis pathway genes under low-temperature stress [9,103,104]. Abiotic stress induces ABA accumulation and combines with PYP/PYL receptors to inhibit PP2C activity, thereby activating protein kinases (SnRK2s) and causing them to phosphorylate downstream transcription factors [105,106]. In this study, ABA content was significantly increased under cold stress ( Figure 2f). Most ABA signaling pathway genes encoding PYL3, PYL8, PP2C, RK2, and DPBF3 were up-regulated; however, a few ABA signaling pathway genes encoding PYL10, ABI5, DPBF3, and PYL1 were significantly down-regulated (Tables S3 and S6), suggesting that cold-treatment-induced endogenous ABA accumulation directly triggers the expression of some important ABA signaling pathway genes. This further promoted the transcription of the COR gene in response to coconut cold resistance.
Auxin is an endogenous small molecule that works synergistically with other hormonal pathways and has a significant effect on plant growth and stress response [107,108]. Auxin response factor and auxin-responsive protein are key components of the auxin signaling pathway and are regulators of plant growth and stress response [109]. The auxin signal of strawberry was severely blocked under low temperature, indicating that it plays a very important role in cold stress [110]. In this study, under cold stress, IAA content increased significantly, reaching a peak at day 20, slowly decreasing on day 30, and then stabilizing (Figure 2g). In auxin-responsive proteins and the factors, 24 DEGs including the auxin-responsive proteins (SAUR32, ARGOS, SAUR71, Os09g0247700, IAA10, SAUR36) and the auxin response factors (ARF1, ARF8, ARF12, ARF6, ARF9, ARF18) were significantly up-regulated (p < 0.05). However, 21 DEGs were down-regulated, in which the auxin-responsive proteins (IAA10, IAA16, IAA6, SAUR71, PIN1C, PIN3A, AUX28, LAX2, SAUR50, SAUR36) and the auxin response factors (ARF19, ARF17, ARF2) were significantly down-regulated (Tables S3 and S6). These genes, which were related to auxin-responsive proteins and factors, and the auxin-responsive proteins were induced under cold stress. Therefore, the auxin signaling pathway and its related genes played a very important role in response to cold stress in coconut.

The Major Significant Enriched KEGG Pathways in Transcription and Metabolism in Response to Cold Stress in Coconut
Low temperatures resulted in significant differences in gene and metabolite expression patterns in coconut seedlings. Large amounts of DEGs were detected from coconut seedlings exposed to low temperatures (Table S3). According to KEGG enrichment analysis, most DEGs-enriched KEGG pathways were "plant hormone signal transduction", "MAPK signaling pathway-plant", "plant-pathogen interaction", "biosynthesis of amino acids", "amino sugar and nucleotide sugar metabolism", "glycerophospholipid metabolism", "circadian rhythm-plant", "carbon metabolism", "starch and sucrose metabolism", "glycolysis/Gluconeogenesis", "purine metabolism", "phenylpropanoid biosynthesis", "steroid biosynthesis", and "lipoic acid metabolism" ( Figure 5, Table S5). Similar results were found in the cold-stress response of Xanthoceras sorbifolia and Brassica napus [32,53].
Many DAMs, including phenylpropanoids and polyketides, lipids and lipid-like molecules, benzenoids, organoheterocyclic compounds, organic oxygen compounds, organic acids and derivatives, nucleosides, nucleotides, and analogues, were identified from coconut seedlings exposed to low temperatures (Figure 7d,h, Table S8a). Similar results were found in cold-stress responses of pepper and peanut [44,55]. In positiveand negative-ion mode, most DAM co-enriched KEGG pathways were "amino acid metabolism", "biosynthesis of other secondary metabolites", "membrane transport", and "metabolism of cofactors and vitamins". In addition, the positive-ion DAMs significantly enriched KEGG pathway included "metabolism of other amino acids", "lipid metabolism", and "translation", while the negative-ion DAMs significantly enriched KEGG pathways contained "Carbohydrate metabolism" and "Nucleotide metabolism" (Figures 7c,g and S8, Table S9). The KEGG enrichment pathways of these DAMs have been confirmed in other plants in response to cold stress [9,111].

Combined Analysis of Transcription and Metabolism in Response to Cold Stress in Coconut
According to the combination analysis of transcription and metabolism in response to cold stress in coconut, KEGG enrichment analysis showed that DEGs and DAMs in positive-and negative-ions modes were significantly enriched in the KEGG pathways, including "amino acid metabolism", "metabolism of other amino acids", "biosynthesis of other secondary metabolites", "metabolism of cofactors and vitamins", "carbohydrate metabolism", "lipid metabolism", "nucleotide metabolism", "translation", and "membrane transport" (Figure 8, Table S10).

Amino Acid Metabolism in Response to Cold Stress in Coconut
Previous reports have shown that various stress response metabolites in plants were synthesized by the amino acid metabolism pathway, suggesting that amino acid metabolism played an important role in plant response to stress [112]. In this study, many genes and metabolites were induced from the "amino acid metabolism" pathway, including "lysine biosynthesis", "tryptophan metabolism", "cysteine and methionine metabolism", "lysine degradation", "histidine metabolism", "phenylalanine, tyrosine and tryptophan biosynthesis", "tyrosine metabolism", "2-Oxocarboxylic acid metabolism", "biosynthesis of amino acids", "cyanoamino acid metabolism", and "glutathione metabolism". Moreover, the DEGs and DAMs of these pathways were significantly correlated (p < 0.05) (Figure 9, Tables S10 and S11). For example, in the lysine biosynthesis pathway, the seven genes encoding AGD2, LOC105058024, COCN_GLEAN_10015782, DHDPS2, AK1, DAPF, and DHDPS2 were down-regulated and were significantly positively correlated with the downregulated L-2-Aminoadipate adenylate (+), while being significantly negatively correlated with the up-regulated metabolites, such as L-Lysine (+) and L-Saccharopine (+); however, on the contrary, it was in terms of the up-regulated LYSA gene correlation with metabolites. In lysine degradation (ko00310), 15 DEGs were up-regulated and significantly positively correlated with L-Lysine (+) and L-Saccharopine (+), while 14 DEGs were down-regulated. In the biosynthesis of amino acids (ko01230) pathway, a total of 169 DEGs were identified to be involved in amino acid metabolism; 46% of DEGs were up-regulated and were significantly positively correlated with shikimate (+, −), L-Isoleucine (+), L-Lysine (+), L-Saccharopine (+), phenylpyruvate (−), S-Adenosylhomocysteine (−), and citric acid (−), but significantly negatively correlated with S-Adenosylhomocysteine (+), L-Glutamic acid (+), and L-Histidine (+); however, 54% of DEGs were down-regulated (Figures 9 and S10a, Table S11). The study showed that the genes involved in the amino acid metabolic pathway were up-regulated or down-regulated in response to cold response, which was closely related to or consistent with the metabolites of this pathway. In short, the changes observed at both the transcriptome and metabolic levels strongly suggest that accumulation of these cold response metabolites via synthetic or metabolic pathways may contribute to cold-stress defense, thereby minimizing cold damage.

Carbohydrate Metabolism in Response to Cold Stress in Coconut
Carbohydrates are the main products of photosynthesis. They are the nutrients and important regulators required for plant growth, energy metabolism, and stress response [21]. Soluble sugars are known to function not only as osmoprotectants [113,114], but also to provide cold resistance in plants as a ROS scavenger [115]. In this study, carbohydrate metabolism mainly included starch and sucrose, carbon, and galactose metabolisms. Metabolome analysis showed that D-FRUCTOSE (−), citric acid (−), 2-Dehydro-3-deoxy-D-gluconate (−), and raffinose were up-regulated; however, sucrose (−) and L-Glutamic acid (+) were down-regulated. In the starch and sucrose metabolism pathway, 54 DEGs were up-regulated and significantly positively correlated with D-FRUCTOSE (−) (p < 0.05), but were significantly negatively correlated with sucrose (−); however, 66 DEGs were down-regulated (Figures 10 and S10b , Table S11). Similarly, previous studies have also confirmed that most DEGs and DAMs are mainly enriched in different pathways involved in carbohydrate metabolism, among which starch and sucrose metabolism were significantly enriched, which played a crucial role in the adaptation of Brassica napus to cold stress [53]. In carbon metabolism, 105 DEGs were up-regulated and significantly positively correlated with citric acid (−) and 2-Dehydro-3-deoxy-D-gluconate (−), but significantly negatively correlated with L-Glutamic acid (+); however, 60 DEGs were down-regulated. In the galactose metabolism pathway, 31 DEGs were up-regulated and significantly positively correlated with D-FRUCTOSE (−) and raffinose (−), but significantly negatively correlated with sucrose (−), whereas 17 DEGs were down-regulated (Figures 10 and S10b, Table S11). Several studies have highlighted the role of these genes in cold stress [113,116]. For example, overexpression of galactinol synthase in cold-induced wheat increases the levels of galactinol and raffinose and confers a higher tolerance to low-temperature stress [117]. Similarly, results obtained also suggest that genes and metabolites involved in carbohydrate metabolism contribute significantly to cold resistance in coconuts.

Flavonoid Metabolism in Response to Cold Stress in Coconut
A subset of plant-specific accumulative metabolites, and the genes involved in secondary metabolism, including flavonoid biosynthesis, flavone and flavonol biosynthesis, and phenylpropanoid biosynthesis, may play a key role in plant response to cold stress [117]. In this study, metabolome detection and analysis identified many flavonoid and secondary metabolites (Table S8). In the flavone and flavonol biosynthesis pathway, three DEGs encoding 3MAT, UGT73E1, and CYP75B137 were up-regulated and were significantly positively correlated with astragalin (+), apiin (+, −), and apigenin 7-O-neohesperidoside (−), but significantly negatively correlated with apigenin 7-O-neohesperidoside (+), astragalin (−), and quercitrin (−); however, DEGs encoding 3MAT was down-regulated. In the flavonoid biosynthesis pathway, 17 DEGs were up-regulated and were significantly positively correlated with (−)-Epicatechin (+), daidzin (+), chlorogenic acid (+), neohesperidin (+), and naringin (−), but significantly negatively correlated with chlorogenate (+), whereas 19 DEGs were down-regulated (Figures S10 and S11, Table S11). It can be seen that cold stress in coconut seedlings has significant effects on flavone and flavonol biosynthesis and flavonoid biosynthesis genes, and related metabolites. Previous studies have also shown that the gene expression of flavonoid biosynthesis is closely related to enhancement of low-temperature tolerance [118]. Flavonoids also played an important role in alleviating cold-stress responses by eliminating ROS produced by cold stress [119]. Flavonoids rich accumulation can reduce cold-stress damage [55,111]. Therefore, it can be inferred that flavone and flavonol biosynthesis and flavonoid biosynthesis genes and related metabolites are very important for the cold resistance of coconut seedlings. Phenylpropanoid biosynthesis is one of the most important metabolites in plants, producing a large number of secondary metabolites, such as lignin and flavonoids [120]. Moreover, in the phenylpropanoid biosynthesis pathway, 66 DEGs were detected, involved in the biosynthesis of other secondary metabolites; 21 were up-regulated and were significantly positively correlated with chlorogenic acid (+) and 3,5-Dimethoxy-4-hydroxycinnamic acid (−), but significantly negatively correlated with 4-Hydroxycinnamic acid (+), coniferol (+), 4-Hydroxy-3-methoxycinnamaldehyde (−), and chlorogenate (−), while 45 DEGs were down-regulated (Figures S10 and S11, Table S11). Recent studies have also shown that cold stress induced the expression of structural genes in the phenylpropanoid pathway, including chalcone synthetase (CHS) and 4-coumarin-CoA ligase (4CL). Flavonoids and lignin in loquat (Eriobotrya japonica) [121] and apple (M. domestica) [122] have promoted the adaptation to the low-temperature environment. In this study, it was also found that cold stress induced the expressions of these genes and associated metabolites in the phenylpropanoid biosynthesis pathway of coconut seedlings ( Figure S10, Table S11). Therefore, the development of gene and metabolite differences in phenylpropanoid biosynthesis was a positive response to cold stress in coconut seedlings. These findings suggest that advances in genetic and metabolic pathways are highly effective in reducing cold-stress damage. At the same time, the aggregation of several genes involved in the metabolism and synthesis of amino acids, sugars, and flavonoids and specific cold-reaction metabolites revealed the cold-stress response and helped to elucidate the molecular mechanism of cold-stress tolerance.

Lipid Metabolism in Response to Cold Stress in Coconut
Lipid metabolism is a biomarker of lipid damages in the plant's response to cold stress [15]. For the glycerophospholipid metabolism pathway in this study, 50% of DEGs were up-regulated and significantly negatively correlated with choline (+) and phosphocholine (+); however, 50% of DEGs were down-regulated. In the arachidonic acid metabolism pathway, six DEGs was up-regulated and significantly positively correlated with 15-keto-Prostaglandin E2 (+) and delta-12-PGJ2 (+), whereas four DEGs were downregulated ( Figure 10, Table S11). This study demonstrated that the expression of genes and metabolites in glycerophospholipid metabolism and arachidonic acid metabolism pathways involved in the lipid metabolism of coconut were induced by cold stress. Moreover, significant correlation has been identified between their DEGs and DAMs, which suggested that lipid metabolism was a very important metabolic pathway in the active response of coconut cold stress.

Nucleotide Metabolism in Response to Cold Stress in Coconut
Nucleotide metabolism can be attributed to cold-induced damage. The study found that nucleotide metabolites, such as cyclic AMP (adenosine monophosphate), adenosine, 2-hydroxyadenosine, cytidine, adenine, β-nicotinamide mononucleotide, guanosine monophosphate, guanosine, and 2-hydroxy-6-aminopurine, increased when A. arguta was exposed to cold, as cold stress can lead to nucleotide damage and stimulate the increase of nucleotide metabolism [15]. For the pyrimidine metabolism pathway in this study, related to nucleotide metabolism, 9 were up-regulated and significantly negatively correlated with 5-Methylcytosine (+), deoxycytidine (+, −), pseudouridine (−), and uridine (−); however, 32 DEGs were down-regulated ( Figure 10, Table S11). The results showed that cold-stressed coconut seedlings induced the expression of genes and associated metabolites in the pyrimidine metabolism pathway involved in nucleotide metabolism, suggesting that nucleotide metabolism was very important for the positive response to cold stress in coconut.

Plant Materials and Cold Treatments
The coconut fruits (Yellow dwarf coconut) (Wenye No. 2) were grown at the Institute of Coconut Research, Chinese Academy of Tropical Agricultural Sciences (Haikou, China). They were planted in plastic culture bags (25 cm in diameter and 35 cm in depth) filled with a mixture of soil, organic fertilizer, and coconut bran (10:2:1). The plants were grown in a culture room for 2 months, with a relative humidity of 80-90%, a photoperiod of 16 h/8 h (light/dark), and a temperature cycle of 25 • C/25 • C (day/night). Next, the highly consistent coconut seedlings (60 cm) were selected and divided into two groups for the treatment experiment. One group (20 plants) was stored in the greenhouse for 30 d. The coconut leaves were harvested at 0 (CK 0 ), 10 (CK 10 ), 20 (CK 20 ), and 30 d (CK 30 ). The other group of coconut seedlings (15 plants) were transferred to the controlled greenhouse for cold treatment at 5 • C, with the same photoperiod. The coconut leaves were harvested at 0 (LT 0 , represented by CK 0 samples; thus, no collection was required), 10 (LT 10 ), 20 (LT 20 ), and 30 d (LT 30 ). Leaves treated with CK 0 , CK 10 , CK 20 , CK 30 , LT 10 , LT 20 , and LT 30 were selected for growth index and physiological index analysis; in addition, leaves treated with CK 30 and LT 30 were subjected to transcriptome sequencing and metabolite analysis, with each treatment being conducted in three biological replicates. At the end of each sample collection, the leaves were immediately frozen in liquid nitrogen and stored at −80 • C until further analysis.

Measurements of Plant Height, Dry Weight, and Soil and SPAD Values
Plant height (accuracy: 1 mm) was measured using a tape measure. The SPAD value of coconut leaves was measured using the SPAD chlorophyll analyzer (SPAD-502 Plus, Konica Minolta, Tokyo, Japan). To assess the dry weight of coconut seedlings, fresh stems, leaves, and roots of coconut seedlings were dried at 105 • C for 15 min and 70 • C for 72 h, and then weighed using an electronic balance (Labpro, Shanghai, China). The dry weight of the whole plant (stem + leaf + root) was calculated.

Physiological Index Measurements
Coconut seedlings were treated at 25 • C and 5 • C for 0 (CK), 10 (CK 10 , LT 10 ), 20 (CK 20 , LT 20 ), and 30 days (CK 30 , LT 30 ), and leaves of each treatment were collected for physiological analysis, all of which were performed with three biological replicates. ABA, IAA, ZR, GA, SS, SP, Pro, and MDA contents and SOD, POD, and CAT activities were determined. First, 0.1000 g coconut leaf sample was accurately weighed and was mixed with pre-cooled PBS in a ratio of 1:10 by weight (g) to volume (mL). These mixed samples were ground at high speed and centrifuged at 2500 rpm for 10 min. Then, 50 µL of supernatant was measured for determination. The IAA, GA, ABA, ZR, MDA, SP, Pro, SOD (A001-3-2), CAT(A007-1-1), and POD (A084-3-1) kits and standards were obtained from the Nanjing Institute of Biological Engineering, and the measurements were carried out in strict accordance with the manufacturer's instructions and according to the methods reported by Li (2000) [123]. We used a 1 cm light path colorimetric tube and a blank colorimetric tube to set the baseline. Absorbance was measured with an enzyme labeling meter (DG5033A, Nanjing Huadong Electronics Group Medical Equipment, Nanjing, China). Wavelengths were set to 595 nm (SP), 620 nm (SS), 532 nm (MDA), 520 nm (Pro), 550 nm (SOD), 405 nm (CAT), 420 nm (POD), and 450 nm (IAA, ABA, GA, ZR). All measurements were taken within 10 min of adding the termination solution. The concentration/activity was calculated from the absorbance value combined with the formula for each physiological indicator.

RNA Extraction and RNA-Sequencing (RNA-Seq)
Total RNA was extracted from frozen samples using the improved cetyltrimethyl ammonium bromide (CTAB) method. The purity and integrity of RNA were evaluated visually by agarose gel electrophoresis. RNA concentrations were measured using a Nan-oDrop 2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). The Agilent 2100 Bioanalyzer system (Agilent Technologies, Palo Alto, CA, USA) was used to quantify the integrity of RNA. Library assembly and RNA-seq analysis were performed at Beijing Biomarker Biotechnology Company (Beijing, China) and Beijing Biomarker Cloud Technology Company (Beijing, China) [124]. Sequencing was performed by using the NEBNext ® Ultra™ II RNA Library Prep Kit (New England Biolabs, Ipswich, MA, USA) and adding an index code to each sample. These libraries were sequenced on the Illumina ® HiSeq2500 platform (Illumina, San Diego, CA, USA). Each sample was sequenced thrice. Raw reads were filtered by removing low-quality reads and adapters. Using the transcript splicing Alignment Hierarchical Index (HISAT 2) program [125], clean sequences were located to the reference coconut genome (http://creativecommons.org/licenses/by/4.0/, accessed on 1 August 2022) [126]. Gene functions were annotated using the following databases: NCBI Non-redundant Protein Sequence (Nr), Homologous Protein Cluster (COG/KOG), Swiss-PROT protein sequence Database, Kyoto Encyclopedia of Genes and Genomes (KEGG), Homologous Protein Families (Pfam), and Gene Ontology (GO) [127,128]. To verify transcriptional expression levels for all samples of each transcription region, RESM software (3.8.6) was used to calculate the fragments per kilobase of transcript per million fragments mapped (FPKM) [129]. DESeq software (1.6.3) was used to analyze the differential gene expression among the samples, and the Benjamini-Hochberg method was used to determine its significance. The definition of DEGs is based on |fold change (FC)| ≥ 2 and FDR < 0.01(log 2 |fold change (FC)| ≥ 1). The False Discovery Rate (FDR) was obtained by correcting the difference significance p-value [130]. The GOseq R software package (2.18.0) was used for GO enrichment analysis of DEGs [131][132][133]. KEGG pathway enrichment analysis of DEGs was performed using KEGG Orthology Based Annotation System (KOBAS) software (3.0) [134].

Metabolite Analysis
Sample preparation, metabolomics, and data analysis was conducted by Beijing Biomarkers Biological Technology Co., Ltd. (Beijing, China) (http://www.biomarker. com.cn/, accessed on 15 August 2022). The frozen coconut leaves were ground into powder in liquid nitrogen, and 100 mg of the powder was added to a 1.5 mL Eppendorf tube. Then, 1.0 mL 70% methanol solution was extracted at 4 • C for 24 h and centrifuged at 10,000× g and 4 • C for 10 min. The extract was filtered through a 0.22 µm nylon membrane and analyzed by liquid chromatography-mass spectrometry (LC-MS). Extract from each treatment and mixture of three duplicate samples was used to prepare quality control samples. During the analysis, each quality control sample was measured together with the corresponding three experimental samples to check the stability of the analysis conditions. The metabolites of leaf extract (10 µL) were analyzed by ultra-performance LC-electrospray ionization mass spectrometry (UPLC-ESI-MS/MS). Aqueous chromatographic separation was performed at 40 • C using a UPLC HSS T3 C18 column (2.1 mm × 100 mm, i.e., 1.8 µm) (Waters, Milford, MA, USA). The mobile phase consisted of water containing 0.04% acetic acid (mobile phase A) and acetonitrile containing 0.04% acetic acid (mobile phase B). The linear gradient procedure for elution was set to 0-11.0 min from 5% to 95% B, 11.0-12.0 min from 95% to 5%, and 12.0-15.0 min from 5% to 5%. The flow rate of the mobile phase was 0.40 mL/min. The API 4500 QTRAP LC-MS/MS system (AB SCIEX, Framingham, MA, USA) was used for LC-MS/MS analysis. ESI source parameters were as follows: turbo spray, ion source, source temperature 550 • C, ion spray voltage 5.5 kV. The curtain gas pressure was 25 pounds per square inch (psi), the ion source gas I pressure was 55 psi, and the gas II pressure was 60 psi. A multi-reaction monitoring experiment was carried out with 5 psi nitrogen colliding gas, and the results of quadrupole scanning were obtained.
The identification of metabolites was based on the public database and cloud technology database of Beijing Biomarker Biotechnology Co., Ltd. (Beijing, China) [124]. The raw data collected using MassLynx V4.2 was processed using Progenesis QI 3.0 software for peak extraction, peak alignment, and other data processing operations; the data were identified based on Progenesis QI software, the online METLIN (2019) database, and the self-built database of Baimaike (http://www.biomarker.com.cn/, accessed on 1 September 2022); and theoretical fragment identification was performed. The quality deviation was within 100 ppm. The structures of the metabolites were analyzed using standard metabolic procedures. The metabolites were quantitatively analyzed by the multiple reaction monitoring method. All the identified metabolites were analyzed by partial least squares discriminant analysis (PLS-DA). Principal component analysis (PCA) and orthogonal PLS-DA (OPLS-DA) were used to identify potential biomarkers. For the choice of biomarkers, with projected variable importance (VIP) ≥ 1, folding change (FC) ≥ 2 (up-regulated), or ≤0.5 (down-regulated), p < 0.05 was used as the criterion for screening significant differentially accumulated metabolites (DAMs).

Integrated Metabolome and Transcriptome Analyses
Pearson correlation tests were performed to determine the correlation between DEGs and DAMs. The Pearson correlation coefficient (PCCs) between DEGs and DAMs was calculated using the "corrplot" in the R package. Genes and metabolites with PCC > |0.8| and p < 0.05 were used to establish correlation heat maps and networks. In the correlation analysis, we compared the KEGG pathway where DEGs and DAMs were enriched (p < 0.05). DEGs and DAMs with high correlation in KEGG pathways (p < 0.05) were selected for analysis and network mapping [135][136][137].

Quantitative Real-Time PCR (qRT-PCR) Analysis
The DEGs of coconut seedlings identified by RNA-seq were verified by qRT-PCR. A gene-specific RT-qPCR primer (Table S12) was designed. The qPCR was run on the Photocycling ® 480II real-time system (Roche, Carlsbad, CA, USA) on a 96-well plate, using HiefeqPCR SYBR green main mixture (Not-Rox) (Yeasen Biotech, Shanghai, China) according to manufacturer's instructions. The thermal cycle steps consisted of denaturation at 95 • C for 5 min, followed by 40 cycles at 95 • C for 10 s and 60 • C for 30 s. All qRT-PCR analyses were performed using three bio-replicates, each of which contained three technical replicates. The internal reference gene (β-actin) was used for normalization. The 2 −∆∆CT method was used to calculate the expression levels of differential target genes in the reference control group [138].

Statistical Analysis
Data were analyzed by Excel (2010) (Microsoft Corporation, Redmond, WA, USA) and expressed as mean ± standard deviation (SD). For statistical analysis, SPSS software was used (Version 20.0; SPSS, Chicago, IL, USA) to conduct one-way analysis of variance (ANOVA) and Student's t-test, and the results were used to determine the difference and statistical significance at p < 0.05.

Conclusions
We investigated the response of coconuts to cold stress at the physiological, transcriptional, and metabolomic levels. Our results suggest that coconuts reduce cold damage by enhancing CAT, POD, and SOD activity and increasing the accumulation of soluble sugars and proteins in cells. The MAPK signaling pathway-plant, plant hormone signal transduction, plant-pathogen interaction, biosynthesis of amino acids, amino sugar and nucleotide sugar metabolism, carbon metabolism, and starch and sucrose metabolism pathways played an important role in the response of coconut to cold stress. TFs, including WRKY, AP2/ERF, HSF, bZIP, MYB, and bHLH, potentially played a critical role in regulating cold-stress-related genes to improve cold resistance in coconut. These results provide a comprehensive understanding of the mechanism of the low-temperature stress response in coconut and guidance for the development of cold-tolerant coconut varieties by identifying candidate genes and incorporating them into molecular breeding programs. This study also deepens the understanding of the complex regulatory mechanisms in plants subjected to cold stress. However, further research should explore the functions of differential genes and metabolites involved in these key regulatory pathways to further understand the mechanisms of the coconut's low-temperature response.