Propagules of K. obovata were picked from Zhangjiangkou Mangrove National Nature Reserve in March 2017 (Fujian, China; 23°53'-23°56'N, 117°24'-117°30'E), the natural habitats of K. obovata (Sheue et al., 2003; Su et al., 2019). Propagules were potted in the sand, and the base of propagules was immersed with Hogland’s solution and 10‰ artificial seawater. Propagules grew into seedlings in a plant culture room, with 27/22°C (day/night) and 16/8 h (day/night) photoperiod, the photosynthetic photon flux density (PPFD) was 130μmol m^-2 s^-1. All plants used in our experiments were cultivated under the same conditions. The K. obovata seedlings were cultivated from March 2017 to February 2018 and were then used in the experiment.

The exact cold acclimation temperatures are species-dependent (Kazemi-Shahandashti and Maali-Amiri, 2018). An appropriate cold acclimation temperature should ensure the survival of plants during the experiment exposure. To find the appropriate temperature for cold treatment, we performed a test experiment, where three groups of seedlings were put into 0 ± 1°C, 5 ± 1°C and 10 ± 1°C low-temperature incubators for 24 hours. In this test experiment, we found that seedings died within a week under 0 ± 1°C, but survived under 5 ± 1°C and 10 ± 1°C. We therefore adopted 5°C in the cold acclimation experiment.

The cold wave was manipulated as a four-day cold treatment with 5 °C followed by a 15-day recovery treatment at 27°C ( Figure 1 ). The experiment included four cold treatment (COLD-1/2/3/4) and three recovery treatments (REC-1/2/3).

In the experiment, six propagules of K. obovata were planted in three pots (individuals 1 and 4 in pot I, individuals 2 and 5 in pot II, and individuals 3 and 6 in pot III) and grew into seedlings in a warm plant culture room (27°C). Besides samples for each treatment (COLD-1/2/3/4, REC-1/2/3), we collected samples from individuals 1-6 before any treatment was given as controls (CTR). At the end of each treatment, leaf samples for physiological analysis were collected from individuals 1-3 and leaf samples for transcriptomic analysis were collected from individuals 4-6. Samples for both analyses were collected from the same pots to reduce the variation in the physiological and transcriptomic analyses. The experiment was conducted in two growth rooms. Each time during cold treatment (COLD-1/2/3/4), the plants were moved to the cold culture room (5°C) for 4 days. During recovery phases (REC-1/2/3) and CTR, the seedlings grew in the warm culture room (27°C). In addition to the difference in temperature between our cold culture room and warm culture room, other culture conditions were identical. The photoperiod in both cold and warm treatment periods was changed to 14/10 h (day/night) to simulate the shortened duration of illumination in autumn and winter.

Plants exposure to cold stress will impair the cellular metabolism and lead to accumulated reactive oxygen species (ROS), including superoxide radical anions (O^2-) and hydrogen peroxide (H₂O₂) (Novaković et al., 2018). Malondialdehyde (MDA) is the final product of lipid peroxidation, which can be used as a measurement of oxidative damage under cold stress, research in Magnolia denudata found that the MDA concentration increased with the enhancement of the cold acclimation (Wu et al., 2021). Antioxidant enzymes, including superoxide dismutase (SOD) and peroxidase (POD), can protect proteins, metabolites, and cell structures against oxidation and peroxidation (Hincha and Zuther, 2020). SODs catalyze oxidation and reduction of superoxide anions and generate H₂O₂, a stable type of ROS, and PODs are efficient in reducing H₂O₂ to H₂O. In this experiment, content of MDA, activities of SOD and POD, and content of chlorophyll of leaves were measured to evaluate the physiological states of K. obovata seedlings.

Leaves of K. obovata seedling individuals 1, 2, and 3 were stored in the icebox and tested immediately by using the reagent kits according to the manufacturer’s instructions (Jiancheng Co., Ltd., Nanjing, A003-1, A001-1, A084-3, A147-1).

The difference of physiological indicators between phases was compared using one-way ANOVA with p values of < 0.05 considered to be significant.

Leaves of K. obovata seedling individuals 4, 5, and 6 were rapidly frozen in liquid nitrogen and then moved to the -80°C refrigerator and then sent to the Experimental Department of Novogene Co. (Beijing, China) for total RNA extraction, library construction, and RNA-sequencing. Initial quality control was done using FastQC (V 0.11.9) followed by trimming with Trimmomatic (V 0.38) (Bolger et al., 2014) to remove adapters and low-quality base pairs. The selection of trimming steps and their associated parameters are supplied on the command line (ILLUMINACLIP:/Trimmomatic-0.38/adapters/TruSeq3-PE-2. fa: 2: 30: 10: 1: true HEADCROP: 12 LEADING: 3 TRAILING: 3 SLIDINGWINDOW: 4: 20 MINLEN: 50). These clean reads were subjected to downstream analyses.

Paired-end clean reads were aligned to the reference K. obovata genome (Qiao et al., 2020) using HISAT2 (V 2.1.0) (Kim et al., 2019). FeatureCounts (Liao et al., 2013) was used to count the reads numbers mapped to each gene. Before differentially expressed gene analysis, genes with a read count smaller than ten in all samples combined were removed from further analyses. The R package “DESeq 2” (Anders and Huber, 2010) was used to analyze the differentially expressed genes (DEGs). Genes with at least |log2FoldChange| > 2 and padj < 0.01 were designated as differentially expressed. All genes with expression at any phase in the experiment were used for Principal Component Analysis (PCA). PCA was performed using the packages DESeq2 and ggplot2 after variance stabilizing transformation of the count data by the function of DESeq2.

To know the cold acclimation mechanism of K. obovata seedlings, we defined cold acclimation-related genes (CARGs) based on the analysis results of DEGs. We set CTR, REC-1, and REC-2 as ‘Warm’ status, genes differentially expressed (|log2FoldChange| > 2 and padj < 0.01) between COLD-1/Warm versus COLD-2/Warm or COLD-4/Warm were named CARGs.

Gene ontology enrichment analysis (Ashburner et al., 2000) was performed using the OmicShare tools, a free online platform for data analysis ² , FDR was calculated using the hypergeometric test. All CARGs were submitted to do the GO enrichment analysis, with FDR <0.05 as the threshold, GO terms satisfying this condition were defined as GO terms that were significantly enriched.

To investigate the possible koCBFs target genes (Shi et al., 2018), a Perl script (Supplementary Information 1) was used to extract the upstream 2000 bp putative gene promoter regions before the start codon of all K. obovata genes to find the existence of the CRT/DRE motif.

WGCNA R package (Version 1.69) (Zhao et al., 2010; Liang et al., 2020) was adopted for the co-expression network analysis. All 2,681 DEGs were used to construct a signed co-expression network. The soft-thresholding power value was chosen as empiric value as 18 dues to no suitable estimated soft power value.

The biological meanings of the eigengenes can be further analyzed by calculating the Pearson correlation of physiological indices data with the eigengene of each module. We used the absolute value of the Pearson’s correlation of Module Membership (MM) and Gene Significance (GS) to select hub genes (Zhao et al., 2010; Liang et al., 2020). Considering the potential cold acclimation-related biological functions of genes in Modules turquoise, blue and yellow, we identified hub genes within-module by combining GS of POD and SOD and MM of Modules turquoise, blue and yellow. Hub gene sets in different modules were identified with the criteria: | MM | > 0.8 and | GS | > 0.3. Hub-transcription factors were manually screened based on gene annotations. In addition, it has been reported that there are a series of genes in Arabidopsis that are co-regulated with ATCBFs and co-regulate the cold acclimation process, they are called “First Wave” genes (Park et al., 2018). Based on this report, we also found transcription factors that were expressed in parallel with koCBF4 in our data.

We found that cold waves significantly changed seedlings’ malondialdehyde (MDA) content (F = 3.17, p = 0.029), superoxide dismutase (SOD) enzyme activity (F = 4.768, p = 0.005), and peroxidase (POD) enzyme activity (F = 3.241, p = 0.033) ( Figures 2A-C , see Supplementary Figure 1 for the scatterplots showing the specifics of each sample). The SOD enzyme activity decreased in the first two cold waves but increased slightly in the latter waves. The POD enzyme activity showed a minor increase in the first cold wave but a greater increase in the later cold waves. Compared to the control, the MDA content only decreased in the second cold wave and appeared to recover in subsequent cold waves. Chlorophyll contents also tended to differ between the first two and the subsequent cold waves ( Supplementary Figure 2 ).

The transcriptome sequencing generated an average of 26 million paired-end reads from each of the 18 samples ( Supplementary Table 1 ). These samples were clustered by the principal component analysis (PCA) into two groups ( Figure 2D ). One group includes samples from the first cold treatment (COLD-1) and warm conditions (CTR and recovery phases), and the other included samples from the later cold waves (COLD-2 and COLD-4). Consistent with the PCA results, the k-means analysis also showed that COLD-1 can be separated from COLD-2 and COLD-4 (k=2, Supplementary Figure 3 ). Also, the differentially expressed genes (DEGs) between COLD-1 and CTR were less than the other two comparisons (COLD-2 vs. REC-1 and COLD-4 vs. REC-2, Supplementary Figure 4 ). These results suggested that seedling responses differed between the first and later cold waves.

The first cold experience initiated the cold acclimation process, and the later cold waves manipulated the scenario to show the cold acclimation effects. Both PCA and k-means ( Figure 2D ; Supplementary Figure 3 ) results showed that CTR, REC-1 and REC-2 samples were clustered with their replicates overlapping. Thus we merged CTR, REC-1 and REC-2 as ‘Warm’, and defined genes whose expression was not affected by the first cold (not DEGs in COLD-1 vs. Warm) but significantly affected by the later cold (DEGs in both COLD-2 vs. Warm and COLD-4 vs. Warm) as cold acclimation-related genes (CARGs) ( Figure 2E , see details in Method ). This defined analysis revealed 1,135 CARGs ( Supplementary Table 2 ), and 3/4 CARGs were upregulated in COLD-2 and COLD-4 treatments, referred to as upregulated CARGs ( Supplementary Figure 5 ). GO enrichment analyses showed that CARGs were significantly enriched (FDR < 0.05, Supplementary Figure 6 ) in pathways related to transcription factor activity, oxidation-reduction process, calcium ion binding, protein ubiquitination, and cell wall.

In the calcium ion binding pathway, we found that two Vacuolar cation/proton exchanger 2 (CAX2) genes were upregulated CARGs ( Figure 3 ), which are Ca²⁺/H⁺ exchangers (Kazemi-Shahandashti and Maali-Amiri, 2018). Genes encoding calcium sensors, including CPKs, CIPKs, and CMLs, were also upregulated CARGs ( Figure 3 ), but CML18 was downregulated CARGs. In the cell wall pathway, we found that previously reported cell wall integrity sensors are CARGs in our list, including the arabinogalactan-proteins (AGPs), Fasciclin-like AGPs (FLAs), one Rapid alkalinization factor (RALFL33), and two LRR receptor-like serine/threonine-protein kinases (LRR-RLKs), MIK2, and THE1 ( Figure 3 ) (Johnson et al., 2003; Cheung and Wu, 2011; Novaković et al., 2018). In particular, RALFL33 belonged to both in calcium-dependent signaling network and cell wall integrity, suggesting cold acclimation of K. obovata involves signaling between calcium and cell wall integrity sensors.

Consistent with the physiological responses of peroxidase ( Figure 2B ), we found that five class III peroxidase superfamily protein genes (PERs) were upregulated CARGs ( Figure 3 ). CARGs also included four cellulose synthases like (CSLs) genes, five xyloglucan endotransglucosylase/hydrolase (XTHs) genes, four pectin methylesterase (PMEs) genes, four EXORDIUM (EXOs) genes, and two Expansin-like genes, EXLA1 and EXPA8. The products encoded by these genes (Novaković et al., 2018) are regulators in cell wall modification affecting the hardness or softness of cell wall.

Most CARGs enriched in the protein ubiquitination process encoded E3 ligases ( Figure 3 ), including 11 Arabidopsis thaliana Toxicos en Levadura (ATLs), seven F-box protein genes, and 12 Plant U-box type E3 ubiquitin ligases (PUBs). Among these E3 ligases, ATL54 was previously reported to function in secondary cell wall formation (Noda et al., 2013). These results thus suggest a connection between ubiquitination regulation and cell wall modification underlying cold acclimation mechanisms of K. obovata.

Co-expression network analysis with all 2,681 DEGs ( Figure 2E ) and the physiological indices data revealed six modules ( Figures 4A ; Supplementary Figure 7 ), including Module turquoise (1,504 genes), blue (563 genes), brown (306 genes), yellow (101 genes), green (77 genes), and grey (130 genes). The eigengene of a module represents a consensus expression pattern of genes in this module (Zhao et al., 2010; Liang et al., 2020). The Eigengenes of Modules turquoise, blue and yellow are similar to the expression pattern of CARGs ( Figure 4B ; Supplementary Figure 8 ). Eigengenes in Module turquoise and yellow were positively correlated with POD, eigengene in Module blue was positively correlated with SOD ( Figure 4C ; Supplementary Figure 9 ), which means that the gene expression patterns and physiological parameters covaried for the three modules.

From these three modules (turquoise, blue and yellow), we identified 1,387 hub genes ( Figure 4 ; Supplementary Figure 10 ). 777 (56%) hub genes were also CARGs ( Supplementary Figure 10B , Supplementary Table 2 ), and 75 were hub-transcriptional factors ( Supplementary Table 2 ). These hub-transcriptional factors included C-repeat/DREB binding factors (CBFs) and were classified into seven families: bHLHs (4), DREBs (3), NACs (5), ERFs (16), MYBs (10), WRKYs (16), and ZATs (8).

CBFs were previously demonstrated to be crucial transcription factors that enhanced the freezing tolerance of plants (Thomashow, 1998). We identified four copies of CBFs in the K. obovata genome, among which koCBF4 was a hub-transcription factor, koCBF3 was CARGs and the other two koCBFs were also upregulated in the repeated cold treatments ( Supplementary Figure 11 , Supplementary Table 2 ). To identify the genes potentially regulated by koCBFs, we screened the 1,135 CARGs and found 491 (43%) genes have at least one CBF binding motif ( Supplementary Table 3 , Supplementary Figure 12 ). Thus, koCBFs may play a critical role in regulating the cold acclimation co-expression network.

We also found 14 hub-transcription factors expressed in parallel with koCBF4 ( Figure 4D , see details in Method), including MYB73, WRKY22, ZAT12, two CZF1s, two WRKY40s, three WRKY33s, and four ZAT10s. CZF1, ZAT10, and ZAT12, they were previously recognized in Arabidopsis thaliana as the CBF-independent transcription factor in cold responses (Park et al., 2018). The parallel expression of these hub-transcription factors and koCBF4 suggests that the transcriptional regulation underlying K. obovata cold acclimation involved both the CBF-dependent and CBF-independent pathways.

We found that both the physiological and transcriptional responses of K. obovata seedlings differed between the first and subsequent cold waves. The MDA content dropped in the second cold wave but increased to a higher level in the third and fourth cold treatments. Another article exploring the mechanism of short-term cold acclimation of K. obovata found that, when seedlings were subjected to cold stimulus for the second time, the MDA content of seedlings have undergone a single cold treatment (no recovery) is significantly lower than that of seedlings without cold memory (Liu et al., 2018). Previous knowledge (Novaković et al., 2018) shows that the influx of Ca²⁺ can produce Reactive Oxygen Species (ROS), and the accumulation of superoxide can oxidize polyunsaturated fatty acids in the cell plasma membranes to produce MDA. Thus, MDA content is an important indicator of membrane lipid peroxidation, and SOD and POD are important enzymes in detoxifying ROS. In this study, the changing of MDA content and SOD enzyme activity showed different expression trends in the first half and the second half of our long-term cold waves experiment, the levels of MDA and SOD tend to increase after an initial decrease. It was worth mentioning that the changes in the expression of PERs (class III peroxidase superfamily protein genes) were very consistent with the changes in POD enzyme activity that we detected. PERs were significantly up-regulated ( Figure 3 ), and the activity of the POD enzyme was also increased ( Figure 2D ). This consistency showed that the POD enzyme does play a role in the cold acclimation of K. obovata. However, despite the elevated tendency of SOD and POD enzyme activities, a high level of MDA indicated that the ROS effect had not been completely removed. MDA content presents a stable level in COLD-3 and COLD-4 might be explained by the retention of ROS, which can be used as a necessary signal transduction mechanism to promote cold acclimation ( Figure 5 ).

Temperature fluctuations usually affect the photosynthetic efficiency of plants (Yuan et al., 2021). In our experiment, although the chlorophyll contents did not show strong changes ( Supplementary Figure 2 ), we found that several genes related to photosystem were significantly downregulated, especially in COLD-4 ( Supplementary Figure 13A ). The proteins encoded by these genes are important components of each subunit in the photosynthetic system, among them, PSAs and PSBs are important subunits of photosystem II complex (PSII) and photosystem I complex (PSI), respectivily. Through the schematic diagram of plant photosystem ( Supplementary Figure 13C ), we found that almost one third of subunits of PSI and PSII were affected ( Supplementary Figure 13B ). This result showed that even though the temperature of each cold treatment was the same, the response of K. obovata seedlings to the regular cold waves was different in the early and late stages. The differential expression of genes related to the photosystem was also found in the study of Brassica campestris’s response to temperature fluctuations (Yuan et al., 2021). A previous study (Liu et al., 2018) also found that cold stress can disturb leaf photosynthesis of K. obovata seedlings. Reducing the photosynthesis rate may help K. obovata seedlings survive better, reflecting the necessity of cold acclimation for mangrove plants to survive the winter. However, this conclusion warrant further investigation.

In the gene expression part, the expressional profile did not show significant changes under the first cold wave; but about ten times more genes were differentially expressed in the later cold stress. Different responses between first and latter cold waves suggested that the first cold experience may alter the gene expression in later cold stress. This phenomenon is consistent with the concept of cold acclimation that the first cold experience can reprogram the response to subsequent cold events (Preston and Sandve, 2013), which is considered an adaption to cold climate.

We found that cold acclimation genes (CARGs) enriched in the calcium signaling pathway included two cation/proton exchangers (CAX2) and 21 Calcium sensors (i.e. CPKs, CIPKs, and CMLs). This observation is consistent with previous studies showing that cytosolic Ca²⁺ concentration was elevated under emergent cold stress and captured by calcium sensors (Yuenyong et al., 2018), which were upstream to CBFs in cold response signaling (Hiraki et al., 2018). As CARGs represented genes whose expression remained unchanged in the first cold wave but changed significantly during later waves, these results indicated that the Ca²⁺-mediated signaling pathway was not only involved in the first cold exposure but also in generating strong signals in later cold waves ( Figure 5 ).

Besides the calcium signaling pathway, we found CARGs also enriched in the cell wall modification pathways. These genes include nine cell wall integration (CWI) sensors (e.g. AGPs, FLAs, THE1, RALFL33, and MIK2) and genes regulating the cell wall components (CSLs, XTHs, PMEs, EXOs, Expansins, and PERs). The primary cell walls of plants are composed of three main components: cellulose, hemicellulose, and pectin. Cellulose microfibrils are embedded in a matrix of hemicelluloses and pectin, and changing the relative composition and connection between them gives the cell wall flexibility (Novaković et al., 2018). Cell wall status is monitored by CWI sensors, which are located in the plasma membrane. In this study, the enrichment of CWI sensors is consistent with the previous studies that cold treatment destroyed the mechanical properties of the cell wall, which was captured by CWI sensors (Johnson et al., 2003; Seifert and Roberts, 2007). THE1 and RALFL33 are known to detect the fragments of the cell wall, MIK2 senses the cellulose integrity, and AGPs and FLAs involve in stress-related wall remodeling (Novaković et al., 2018). In addition, a previous study suggested that ROS could cause cell wall damage (Novaković et al., 2018) and trigger the cell wall remodeling process initiated by CWI sensors, but details of this pathway remained to be explored.

After CWI sensors perceived the cell wall damage, the cell wall modification process may be initiated (Novaković et al., 2018). Our findings were consistent with previous studies that show alterations in these enzymes’ activities enhanced cold tolerance by affecting wall properties (Novaković et al., 2018). All these genes (CSLs, XTHs, PMEs, EXOs, Expansins, and PERs) were demonstrated to be involved in cell wall flexibility (stiffening/loosening). CSLs and XTHs are found to influence the cell wall ingredient of hemicellulose, PMEs can regulate pectin stiffness/viscosity, EXOs and Expansins can control cell expansion, and PERs can affect the chemical bonds that connect components of the cell wall (Novaković et al., 2018; Cao et al., 2019; Takahashi et al., 2020). However, as these genes are involved in both cell wall stiffening and loosening processes (Novaković et al., 2018), whether the cell wall gets stiffer or looser in cold acclimation remains unknown ( Figure 5 ).

Transcriptional binding sites analysis suggested that ten cell wall stiffening/loosening regulating genes (e.g. PERs and XTHs) had ‘CCGAC’ binding motif and maybe under the transcriptional regulation of koCBFs ( Supplementary Table 2 ). This finding is consistent with previous studies on Arabidopsis thaliana and tea plant, Camellia sinensis, that cell wall modification genes were targets of CBFs during cold acclimation (Zhao et al., 2016; Wang et al., 2019). The results indicated that the cold acclimation of K. obovata seedlings may involve CBF-related signal transmission from calcium sensors to cell walls ( Figure 5 ).

We found that CARGs were also enriched in pathways relating to Ubiquitination ( Supplementary Figure 6 ). The genes enriched in these pathways included three kinds of E3 ligases, ATLs, F-boxes, and U-boxes. Previous knowledge shows that these E3 ligases are important in substrate identifying and binding, in mediating cell division, nutrition balance, and plant growth (Song et al., 2017; Miricescu et al., 2018). These studies have emphasized the key role of ubiquitination in regulating plant development and acclimating to environmental changes, through mediating the proteolysis of key regulators in phytohormone pathways. However, there is no report on the role of E3 ubiquitin ligase in cold acclimation. The enrichment of several E3 ligases may indicate that the ubiquitination process was beneficial to the cold acclimation of K. obovata seedlings, but the substrates of these E3 ligases were not clear and required further investigation.

We found that four koCBFs were upregulated in response to cold waves ( Supplementary Figure 11 ), with koCBF3 and koCBF4 being CARGs ( Supplementary Table 2 ). CBFs are transcription factors that mediate plant cold tolerance by binding to CRT/DRE cis-elements containing the conserved ‘CCGAC’ sequence (Shi et al., 2018). About 43% of CARGs had at least one ‘CCGAC’ conserved binding motif of CBFs in their promoter region ( Supplementary Figure 12 ). We found a previously validated CBF-regulating gene, COR47 (Gilmour et al., 1992), containing four ‘CCGAC’ motifs ( Supplementary Table 2 ), suggesting that our motif prediction approach yielded robust results. koCBF4 was a hub transcription factor in the Module turquoise. These results indicated that the classic CBF-dependent pathway may function in the K. obovata seedlings’ cold acclimation process ( Figure 5 ).

Many CBF-dependent genes known to play important roles in cold acclimation did not show significant changes in this study e.g. ICE1, MPK4, MPK6, or OST1 (Shi et al., 2018). This is possible because this study adopted long-term cold waves, whereas the above genes reflect the immediate response to cold stress. These results suggest that cold acclimation and transient cold response might be underlined by different molecular mechanisms.

Besides CBFs, we found several other hub-transcription factors, such as CZF1s, ZATs in regulating the co-expression during cold acclimation ( Figure 4D , 5 ). The findings of these transcription factors are consistent with previous studies that cold acclimation can be achieved by CBF-independent pathways (Park et al., 2018). The upregulated expression of CZF1s and ZATs under cold acclimation is consistent with reports in A. thaliana and other woody plants (Park et al., 2018; Liang et al., 2020). The study of A. thaliana showed that the CBF-independent pathway may co-regulate the downstream cold response genes (Park et al., 2018). As the CBF motif is only observed on about half of CARGs ( Figure S6 , Supplementary Table 3 ), CBF-independent transcription factors are anticipated to play an equally important role in regulating the cold acclimation of K. obovata seedlings ( Figure 5 ).