Analysis of the expression and function of the CBL-CIPK network and MAPK cascade genes in Kandelia obovata seedlings under cold stress

Mangroves are an important component of coastal wetland ecosystems, and low temperature is the main factor that limits their extension to higher latitudes. Kandelia obovata as one of the most cold-tolerant species in mangrove ecosystems can provide basis for the northward migration of mangrove ecosystems. We took K. obovata seedlings from Zhoushan (Zhoushan, Zhejiang, China) as the research object in this study. Transcriptome sequencing based on the Illumina HiSeqTM 2500 platform was performed to compare the transcriptome changes of roots, stems, and leaves before and after freezing and to reveal the molecular mechanisms of frost resistance. A total of 1560, 370, and 416 genes were differentially expressed in the roots, stems, and leaves before and after cold snaps, respectively. Among these differentially expressed genes, 13 positive and negative regulators were attributed to the CBL-CIPK signaling network and MAPK cascade, which might be related to the frost resistance mechanism of K. obovata Transcription factors such as AP2/EREBP and bHLH were involved in regulating the synthesis pathways of ethylene, cytokinin, growth hormone, and flavonoids. Results provide new insights into the frost resistance mechanism of K. obovata seedlings.

environment (Wang et al., 2011). Mangrove ecosystems play an irreplaceable role in protecting coastal wetland ecosystems and biodiversity and maintaining coastal ecological balance. They are also important ecological security systems in the coastal region of southeastern China (Wang et al., 2011;Su et al., 2019). Among mangrove plants, Kandelia obovata has the highest natural distribution latitude in China and is one of the most coldresistant mangrove species (Wang et al., 2011). In the 1950s, this species was successfully transplanted from the northernmost part of its natural distribution in Fujian Province (27°17′ N) to further up north in Yueqing in the Zhejiang Province (28°20′ N) (Su et al., 2019;Fei et al., 2021b). However, the persistent phenomenon of climate change has led to a tendency for mangrove species to expand to higher latitudes. The introduction of mangrove K. obovata has been recorded in the Zhoushan of Zhejiang Province (29°30′N) (Osland et al., 2013;. During their northward journey, mangroves are susceptible to extreme weather conditions. In the context of global climate change, the frequency and intensity of extreme weather, particularly in the tropical-temperate transition zone, limits plant productivity and affects species distribution, ecosystem structure, and ecosystem function (Chen et al., 2017). The tendency of North Atlantic and East Asian mangroves to spread northward is significantly influenced by the frequency and severity of cold snaps (Cavanaugh et al., 2014;Cavanaugh et al., 2015;Osland et al., 2017). Cold stress is divided into chilling stress (0°C-15°C) or freezing stress (<0°C). Chilling stress may impede the growth and development of plants, and extended exposure to temperatures below 0°C can cause ice crystals to accumulate in the cytoplasm of plant cells. This phenomenon can dehydrate cells or rupture cell membranes or walls, causing irreparable harm to the plant (Shi et al., 2018). The various physiological and metabolic processes that occur in plants at low temperatures are very intricate regulatory processes (Guy, 1990). Numerous investigations into the physiological and molecular processes underlying the hypothermic response of K. obovata have been conducted and yielded insightful results (Peng et al., 2015;Liu et al., 2019;Wang et al., 2019;Fei et al., 2021a). For example, the KoWRKY40 gene in WRKY transcription factors plays key roles in responses to biotic and abiotic stresses (Liu et al., 2019;Du Z. K. et al., 2022). Among the C-repeat binding factor (CBF), KoCBF1 and KoCBF3 may have functions in plant growth and cold resistance and are involved in responses to salinity and heavy metal stress (Peng et al., 2020). Ko-Aquaporins (KoAQPs) aid in the transport of water and solutes for adaptation to coastal intertidal habitats and have complex responses to environmental factors, such as salt and temperature (Guo et al., 2022). KoOsmotin, a class of proteins with osmotic functions, is highly induced in K. obovata leaves during cold stress (Fei et al., 2021a). KoHsp70, stress-inducible heat shock protein gene, plays a role in protective response to cold stress in K. obovata . Moreover, the photosynthetic system, enzymes, and non-enzymatic antioxidants of K. obovata are extensively involved in stress response at low temperatures to increase the plant's potential for cold adaptation (Peng et al., 2015;Chen et al., 2017). Gene expression regulation tends to respond rapidly to external environmental changes by regulating protein abundance. The six levels of plant gene expression control include mRNA degradation control, protein activity control, RNA processing control, RNA transport, localization control, transcriptional control, and translational control. As transcriptional regulation is a multi-faceted, deep-level process, further research is needed in this field (Alberts, 2015). Although the gene expression profiles of many tissues have been disclosed by sequencing of critical genes for cold adaptation in adult K. obovata, information on the two important cold stress nodes of the B-like protein-interacting protein kinase (CBL-CIPK) network and mitogen-activated protein kinase (MAPK) cascade has not been completely reported (Su et al., 2019).
The CBL-CIPK network plays a key role in cold stress response, where CBL proteins interact with CBL-interacting protein kinases (CIPK) or SOS Ras/Rho Guanine Nucleotide Exchange Factor 2 (SOS2) protein kinase to form different signaling cascades. A complex CBL-CIPK network can regulate Ca 2+ transport in response to low temperatures and freezing injury to plant cells Kolukisaoglu et al., 2004;Buchwal et al., 2020). Ca 2+ affects the cold tolerance of plant cells by improving the stability of the cell wall and promoting the secretion of peroxidase Kolukisaoglu et al., 2004). The typical MAPK cascade signaling network is composed of three specific protein kinases, namely, MAPK kinase kinases (MKKKs), MAPK kinases (MKKs), and MAPKs. These protein kinases are activated by sequential phosphorylation at conserved activation sites, and the interaction between SaMKK2 and SaMAPK4/7 increases the plant's resistance to cold temperatures (Teige et al., 2004;Chen et al., 2022). Currently, the CBL-CIPK network and the MARK cascade have been isolated from many plants, and their regulatory pathways have been well studied. However, the regulated transcriptional signaling pathways of K. obovata are still poorly understood, especially during its seedling growth phase.
Due to large-scale ecological restoration operations along the Chinese coast, K. obovata has been introduced as a non-native species in Zhoushan beyond the limits of their natural range. However, the area is more susceptible to extreme cold snaps. Extreme cold events have a significant effect on the abundance and productivity of mangrove plants. However, as plants grow and mature, adults have stronger resistance and recovery from cold events than seedlings. Given that the survival and development of seedlings often determine the rate at which the population grows, additional research is required on mechanisms by which K. obovata withstand cold events during northern migration (Coldren and Proffitt, 2017). In this context, field-based experiments will be more helpful in understanding how K. obovata seedlings react to cold events or freezing damage because extreme cold snap events are relatively rare. It is difficult to replicate a natural environment similar to freezing stress on K. obovata seedlings in an artificial environment. (Osland et al., 2015). In this study, transcriptome changes were examined in the root, stem, and leaf tissues of K. obovata seedlings transplanted in the Zhoushan intertidal zone before and after a cold snap. Given the relevance of CBL-CIPK network and MAPK signaling pathways on the cold tolerance of plants, we then identified the genes involved in these pathways in K. obovata and compared in detail the gene expression of these two pathways before and after the cold snap to provide insights into the freezing resistance mechanisms of K. obovata seedlings under extremely cold events.
2 Materials and methods 2.1 Material collection K. obovata embryos were collected from Lei Zou of Guangdong in March 2020 and temporarily nurtured in an outdoor nursery in the Zhoushan area for 6 months before being moved to the intertidal zone of the Zhoushan main island in September 2020 ( Figure S1). We transplanted 300 seedlings on about 600 m 2 of tidal flat. We determined the sampling time range according to Zhoushan's local historical meteorological data. To predict the best sampling times, we followed the weather forecast to know when a cold snap would occur. Three tissues including the root, stem, and leaf were freshly sampled from different healthy K. obovata seedlings on December 29 and December 31 before and after the cold snaps ( Figure S2). The lowest temperature registered during the cold snap was -2.10°C, measured with an automatic temperature recorder (HOBO MX2201) placed on the sampling sites for the period (2020.12.16-2021.01.13). After collection, samples were washed with water, dried with paper, and stored in liquid nitrogen in an ultra-low temperature refrigerator at -80°C. At least three biological replicates were collected from the field at each sampling time, with each plant representing an independent biological replicate. In total, we analyzed 18 samples (3 biological replicates × 2 sampling times × 3 tissues).

RNA extraction, library preparation, and RNA-Seq
Extraction was performed according to the manufacturer's instructions using the RNeasy Plant RNA Extraction Kit (Qiagen, Dusseldorf, Germany) to extract RNA from each of the six unique tissues mentioned above. RNA concentration was measured using an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). NanoPhotometer ® spectrophotometer (IMPLEN, CA, USA) was used to check RNA purity.
We constructed cDNA library and conducted RNA sequences at Biomarker Biotechnology Corporation in Beijing, China. About 3 mg of RNA per sample was used for RNA preparation. A poly-T oligo-attached magnetic bead method was used to purify mRNA from total RNA. Dividing the mRNA with divalent cations was conducted at elevated temperatures using NEBNext First-Strand Synthesis Reaction Buffer (5×). The first-strand cDNA was synthesized using random hexamer primers and M-MuLV reverse transcriptase (RNase H-). DNA polymerase I and RNase H were used to synthesize cDNA on the second strand while exonucleases and polymerases were used to turn the remaining overhangs into blunt ends. To prepare samples for hybridization, the 3′ ends of the DNA fragments were adenylated and ligated to NEBNext adaptors with hairpin loop structures. In the following step, 150-200 bp-long cDNA fragments were selected from the library and purified using a Beckman Coulter Agencourt AMPure XP system (Brea, California, USA). The size-selected, adaptor-ligated cDNA was added with 3 mL of USER enzyme (New England Biolabs, Ipswich, MA, USA) at 37°C for 15 min followed by 5 min at 95°C. PCR analysis was performed using Phusion High-Fidelity DNA polymerase (Thermo Fisher, Waltham, MA, USA), universal PCR primers, and index (X) primer. A Bioanalyzer 2100 was used to evaluate the quality of the library constructed from the PCR products purified with the AMPure XP system. Clustering of index-coded samples was performed using a cBot Cluster Generation System and a TruSeq PE Cluster Kit v3-cBot-HS (Illumina, San Diego, CA, USA). The manufacturer's instructions were followed for all experimental procedures. Using the Illumina Hiseq 2000 platform, paired-end reads were generated from the library preparations.

Differential expression analysis
Raw data in fastq format were processed using in-house Perl scripts. Clean data were obtained by removing low-quality reads and those containing the adapters or poly-N sequences. We checked the quality of the unassembled read dataset by examining Q20, Q30, GC-content, and sequence duplication. All the downstream analyses were performed using high-quality clean data.
Genes from 18 transcriptomes were previously annotated in seven public databases including NR (NCBI, Non-Redundant Protein Database), Swiss-Prot, GO (Gene Ontology), Pfam database, KOG database, COG (Clusters of Orthologous Groups), and KEGG (Kyoto Encyclopedia of Genes and Genomes). Based on the annotation results of 18 transcriptomes, the unigene annotation information of the K. obovata seedling transcriptome was screened and analyzed. FPKM method was used to calculate gene expression based on the annotation results. The relative measure of transcript abundance was fragmented per kilobase of transcript per million mapped reads (FPKM) (Florea et al., 2013). Correlation and principal component (PCA) analyses were also performed for the 18 transcriptomes of genes by using the latic, cluster package.
Differentially expressed genes (DEGs) were analyzed using the DESeq 2 package to visualize the results, with FDR < 0.01 and |log2FoldChange| ≥ 2.00 as the screening criteria to obtain the set of DEGs among samples. DEGs were classified as up-or down-regulated genes based on significant positive or negative log changes in their values. The screened DEGs were subjected to GO and KEGG pathway enrichment analyses, where q-value < 0.05 indicated the significant enrichment of the DEGs.

Identification of transcription factor
The TF family of K. obovata seedlings was predicted in different tissue samples by using the Plant Transcription Factor Database (version 3.0). The DEGs were submitted to iTAK software for prediction and further classification. The dataset of the TFs of K. obovata seedlings was used as reference, and the default values of the parameters "substitution matrix" and "E-value" were used .

Identification of CBL-CIPK network and MAPK cascade genes in K. obovata seedlings
The NGDC website (https://ngdc.cncb.ac.cn) contains information on the sequences of K. obovata genes. The CBL-CIPK network and MAPK cascade genes of K. obovata can be identified by comparing the protein sequences generated by the genes of K. obovata seedlings in the database with those in the experiment. When the similarity of gene sequences was higher than 80% compared with the database, it was identified as the corresponding genes with the close homology relationship in the k. obovata.

Gene expression analysis of CBL-CIPK network and MAPK cascade in K. obovata seedlings
The expression of 13 key genes in the CBL-CIPK network and MAPK cascade in K. obovata seedling was studied under control and cold stress conditions. The fpkm values of gene expression were normalized by z-score, and the expression was plotted using the stats package for heat mapping. We performed Pearson correlation analysis to obtain correlation matrices between gene pairs. The network relationship graph between 13 genes was determined using the igraph package in R language.
3 Results 3.1 Creation of databases and functional annotation of the transcriptome of K. obovata seedlings We generated 18 cDNA libraries in 18 samples of various tissues of K. obovata seedlings at different temperatures by using RNA-seq technique. A total of 117.37 Gb of raw data were obtained to gain insights into the response network of genes in the root, stem, and leaf of K. obovata seedlings before and after cold stress. After removing low-quality sequences and adapter sequences from the raw data obtained by sequencing, the GC rates of the 18 samples were 45.58%-46.5%, and all Q30s were greater than 90%, indicating the sufficient accuracy and quality of the sequencing data for further analysis. The general sequencing statistics are shown in Table S1. On the PCA of the gene expression profiles of the 18 samples, three biological replicates of each tissue clustered together ( Figure S3), indicating the reliability of all mRNA sequence library replicates. The Pearson correlation analysis of the dataset revealed gene expression patterns and significant tissue specificity, with identical tissues clustered together ( Figure S4), consistent with the phenomenon in Figure S3. To gain insight into single genes, GO was classified using WEGO, resulting in 53 GO terms in a threelevel ontology ( Figure 1A) (Ye et al., 2018). The category of biological process was dominated by "cellular process" (GO:0009987, 6249 unigenes), "metabolic process" (GO:0008152, 5723 unigenes), "single-organism process" (GO:0044699, 4353 unigenes) dominated. In the category of cellular component, "cell" (GO:0005623, 5901 unigenes), "cellular part" (GO:0044464, 5901 unigenes), "membrane" (GO:0016020, 5088 unigenes), "organelle" (GO:0043226, 4583 unigenes), and "membrane part" (GO: 0044425, 4445 unigenes) had important roles. In the GO molecular function, a large number of single genes were annotated as "binding" (GO:0005488, 8107 unigenes) and "catalytic activity" (GO:0 003824, 6915 unigenes). The KEGG pathway annotations were used to identify biochemical pathways in the K. obovata seedling transcriptome. A total of 6,902 single genes were matched to 136 pathways. The major pathways included "plant hormone signal transduction" (524 single genes), "plant-pathogen interaction" (494 single genes), "MAPK signaling pathway-plant" (318 single genes), and "carbon metabolism" (311 single genes) pathways (Table S2).

Identification and functional enrichment analysis of DEGs in different tissues under cold stress
RNAseq was used to identify DEGs in different tissues of K. obovata seedlings before and after the cold snaps. According to the RNA-seq results, gene expression in the seedlings changed significantly before and after the cold snap, with a total of 1560, 370, and 416 DEGs in the root, stem, and leaf respectively. The Volcano map and Differential genetic Wayne plots showed the aggregation of DEGs in the roots, stems, and leaves, respectively (Figures 2 and S5). The DEGs included 346 up-regulated genes and 1214 down-regulated genes in the roots, 183 up-regulated genes and 187 down-regulated genes in the stems, and 281 up-regulated genes and 135 down-regulated genes in the leaves (Figure 2).
The DEGs of each tissue were subjected to GO and KEGG pathway enrichment analyses to determine the possible functions of the cold response genes. In GO enrichment, the leaf response annotations were mainly focused on biological processes, with "system development" and "regulation of secondary metabolic processes" as the most abundant. The most abundant GO terms in molecular functions were related to "DNA binding." The root and stem response annotations tended to be consistent, and the molecular function of "transcription factor activity" dominated. The cellular component was mainly annotated as "nuclear," and the "management of transcription factors" was highly represented in biological processes. The terms "sequence-specific DNA binding" and "response to acidic chemicals" were unique to the root (Figures 2B-D). Among the KEGG pathway enrichments, circadian-vegetative (ko04712), taurine, and hypotaurine metabolism (ko00430) were significantly enriched in the leaves. Plant-pathogen interactions (ko04626), MAPK signaling pathway (ko04016), biosynthesis of scopolamine, piperidine, and pyridine alkaloids (ko00960), beta-Alanine metabolism (ko00410), and amino and nucleotide sugar metabolism (ko00520) were significantly enriched in the roots. Starch and sucrose metabolism (ko00500), phenylpropanoid biosynthesis (ko00940), and carotenoid biosynthesis (ko00906) were enriched in the stems. Phytohormone signaling (ko04075) was enriched in all tissues ( Figures 3A-C).

Transcription factors in response to cold stress in K. obovata seedlings
Transcription factors regulate gene expression and are thus regulators of cellular processes and environmental responses. Transcription factors can control the expression of many target genes by specifically binding to their various promoters (Nakashima et al., 2009). To clarify the transcriptional regulation of cold stressresponsive genes, we identified DEGs encoding TFs under cold stress. Among transcription factors, AP2/ERF, WRKY, C2H2, and NAC were concentrated at the root. And the enrichment of these transcription factors in the root was much higher than in the stem and leaves. The major transcription factors in the stems were the same as in the roots, except for the increase in HD-ZIP transcription factors. WRKY did not appear in the stems. The transcription factors NAC, AP2/ERF, bHLH, and GARP were mainly present in the leaves ( Figures 4A-C).

FIGURE 2
Volcano plot showing the numbers of DEGs among the three K.obovata seedlings organizations. The x-axis represents the fold-change (FC) of gene expression among samples. The more different genes are at either end. The y-axis represents the significance of the shift in gene expression, which is indicated by -log10. The y is the P-value (t-test) negative log, as adjusted by the false discovery rate (FDR). Genes that were significantly differentially expressed are represented by red dots (up-regulated) and green dots (down-regulated), while insignificantly differentially expressed genes are represented by black dots. Leaf: 281 up-regulated genes and 135 down-regulated genes. Root: 346 up-regulated genes and 1214 downregulated genes. Stem: 183 up-regulated genes and 187 down-regulated genes.

Identification and functional analysis of CBL-CIPK network and MAPK cascade genes
Research has concentrated on the CBL-CIPK network and MAPK cascade of Arabidopsis thaliana and some model plants.
No study has been conducted on thorough identification and functional analysis of the expression and regulation of the CBL-CIPK network and MAPK cascade-related genes in K. obovata seedlings under cold stress. In the present work, we found that K. obovata seedlings contain 13 CBL-CIPK networks and MAPK cascade-related genes ( Figure 5A). The heat map shows that CBL4, CBL3, and CIPK2 are grouped together, but they are separated from CBF1, meanwhile, CBL2, CBL3, and MAPK6 are grouped together. The gene functions of the CBL-CIPK network and MAPK cascade are connected in the K. obovata seedlings. The co-expression regulatory network of the 13 genes in K. obovata was established using this dataset, and the mutual regulatory relationships between among the genes were visualized. The gray line represents the negative regulatory relationship between genes, the blue line represents the positive regulatory relationship between genes, and the increase in line thickness represents a stronger prediction relationship ( Figure 5B). Among the coexpression regulatory networks, genes in the CBL-CIPK network and MAPK cascade, for instance CRLK1 with MEKK1 and MEKK1 with MAPK3/6, had strong correlations. The coexpression relationship among the 13 genes could provide insights into the correlation between the CBL-CIPK network and MAPK cascade genes in K. obovata seedlings under cold stress. Significant transcription factor (TF) classification in three tissues of K. obovata seedlings, leaves (A) and roots (B) and stems (C), before and after chilling, with vertical coordinates representing the number of transcription factors (P < 0.05 and |LOG2FC| > 2).

FIGURE 3
Results of gene and genome enrichment analysis of DEGs in leaves (A) and roots (B) and stems (C) of K. obovata seedlings before and after the cold snaps and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway. The horizontal coordinates represent the degree of enrichment, the different colors represent the size of the p-value, and the size of the dots represents the number of genes.

Signal-mediated cold stress response
Under cold stress, plants can trigger the expression of genes involved in multiple signal transduction pathways and further activate downstream regulatory processes related to physiological adaptation. Cold stress generates downstream responses in plants by modifying membrane flow in plant cells (Orvar et al., 2000). Important processes, such as nutrient and ion transport, solute uptake, and osmotic pressure, are affected by the properties of cell membranes, particularly their structure and composition, proper membrane fluidity is necessary to maintain the structure and function of the membrane (Valerica Raicu, 2008). According to our studies, there were 463 DEGs related to membrane composition altered in various tissues (101, 367, and 64 in the leaves, roots, and stems, respectively), among which 69 were up-regulated in the leaves, 46 in the roots, and 27 in the stems. This finding suggests that different tissues of K. obovata seedlings respond differently to cold stress. The membrane fluidity increases in the leaves but decreases in the roots and stems. Changes in membrane fatty acid composition and content can affect membrane fluidity, thereby influencing the adaptation of K. obovata seedlings to cold stress (Zhao et al., 2012). Membrane fluidity initiates cold signaling by affecting the transport of ions and metabolites (Jeon and Kim, 2013). Ion transport mainly affects Ca 2+ ions released from apoplastic, endomembrane systems, and organelles, and CIPK acts as a Ca 2+ ions sensor in response to elevated cellular solute Ca 2+ levels (Tähtiharju et al., 1997). Under cold stress, four genes associated with Ca 2+ signaling were found in K. obovata seedlings, consistent with the expected role of Ca 2+ in early cold signaling. The results are in line with Su's findings in his study of adult hardiness (Su et al., 2019). Ca 2+ induces intracellular biochemical responses through sensing and transduction of downstream CDPKs, thereby regulating plant responses to various abiotic stress signals (Boudsocq and Sheen, 2013). CBL is a specific calcium sensor with non-classical EF-hand domains that can capture calcium signals in different plant species, including Arabidopsis . These CBLs/SCaBPs do not have enzymatic activity of their own, but CBLs are triggered by changes in Ca 2+ characteristics and interact precisely with CIPK (CBL-interacting protein kinase) or SOS2 protein kinase . Increasing lines of evidence indicate that the CBL-CIPK network plays a key role in cold stress response. In the present study, one CDPK, twoCMLs, three CBLs, and three CIPKswere found to be up-regulated, indicating the critical function of Ca 2+ -mediated signaling pathways in the response of K. obovata seedlings to cold stress.

Metabolic pathways involved in cold stress
Plants produce metabolites to adapt to a stressful, changing growing environment, which may involve interactions through signaling processes and pathways (Edreva et al., 2008). In the present study, the term "regulation of secondary metabolic processes" was a rich category in GO classification. The results of the KEGG pathway enrichment analysis showed that starch and sucrose metabolism, amino saccharide, and nucleotide sugar metabolism were enriched in different tissues. Plant responses to cold stress include sucrose accumulation due to the induction of sucrose phosphate synthase and sucrose catabolic enzymes (convertases and sucrose synthases) (Coleman et al., 2009). Moreover, nucleotide sugars are key intermediates in the starch and sucrose biosynthesis pathways (Coleman et al., 2009;Figueroa et al., 2021). At low temperatures, starch hydrolyzes into soluble sugars (sucrose), improving the osmotic pressure in cells and preventing protoplasm from dehydration and solidification. Dissolution of sucrose increases the concentration of the cell fluid, reducing the freezing point, avoiding the possibility of protein precipitation and solidification and salting, and simultaneously increasing the cell's structural stability (Sun et al.,  2015). We found the enrichment of metabolic pathways for starch and sucrose metabolism as well as amino and nucleotide sugars in the stems and roots. These metabolic pathways are significantly enriched in cold-domesticated K. obovata and play an important role in its resistance to cold stress.

Transcription factors involved in the cold stress response
To resist and adapt to cold stress, plants form complex and efficient regulatory networks, with transcriptional regulation playing a key role. Transcription factors regulate gene expression by binding to cis-acting elements in the promoter region and have an important function in the plant abiotic stress response network . Furthermore, TFs highlight the role of multiple transcriptional mechanisms in stress signaling pathways and are linked to the regulation of various abiotic stress signals and gene transcription (Shanker and Venkateswarlu, 2011). Many TF families, including AP2/EREBP, bHLH, WRKY, MYB, NAC, and MYC, coordinate the signals transduced by plants in response to different environmental stresses (Liu et al., 2014). According to our results, the largest TF families induced by cold stress are AP2-ERF and bHLH, which consist of 20 members, respectively. AP2-ERF and bHLH-TFs are involved in a variety of physiological processes, and overexpression of AP2/ERF in red pine and PtrbHLH in lemon and tobacco can enhance the cold tolerance of plants at freezing temperatures (Huang et al., 2013;Wang et al., 2020). According to Hajar's research results, the expression level of plant AP2/ERF transcription factors in roots was higher than that in stems and leaves under cold stress, which is similar to our research results (Alberts, 2015). WRKY is the largest TF family in plants, and it controls plant growth and development and plays a critical part in many stress responses. NAC is a class of transcription factors that are unique to plants.It is involved in various cold stress regulatory pathways, particularly ABA-dependent signaling pathways, which interact with the gibberellin and Auxin signaling pathways and control lateral root development (Nakashima et al., 2012;Wang et al., 2018). Since the samples we collected were seedlings, the overexpression of WRKY and NAC transcription factors in the roots might improve the plant's cold tolerance and promote root development Fei et al., 2022;Sun et al., 2022). Cold stress in plants first damages the structure of the cell membrane, thus affecting its function of the cell membrane. C2H2-ZF transcription factors can enhance the ability of plants to resist cold stress by maintaining the stability of cell membranes . At the same time, C2H2-ZF transcription factors contain ciselements related to plant hormones or abiotic stress, which are involved in tissue and organ development, especially root and flower . Therefore, C2H2-ZF transcription factors were highly enriched in roots. Other TFs, including DEGs of the GRAS, bZIP, MYB, and WHSF families, also respond to cold stress in K. obovata seedlings. Hence, these TFs may be essential regulators that drive downstream gene expression cascades.

CBL-CIPK network and MAPK cascade involved in cold stress response
The MAPK cascade is a crucial plant response to abiotic stress. The MAPK1-MKK4/5-MAPK3/6 module of Arabidopsis thaliana was the first MAPK signaling module identified in a plant. It mediates the transcriptional up-regulation of genes encoding transcription factors WRKY22 and WRKY29 in response to fungal and bacterial pathogen assault (Asai et al., 2002;Galletti et al., 2011). The MKK4/5-MAPK3/6, MKK1-MAPK3/6, and MKK6-MAPK4 modules are key players in plant response to cold stress (Wang et al., 2021). We identified 13 genes in the CBL-CIPK network and MAPK cascade signaling module of K. obovata seedlings and investigated their relevance under cold stress. When subjected to cold stress conditions, MEKK1-MKK2-MPK4 confers cold sensitivity and freezing tolerance in K. obovata seedlings (Zhao et al., 2017). Low temperatures change the fluidity of the cell membrane, which Ca 2+ channels and proteins can sense on the membrane. The calcium response protein and MAPK cascade genes are then activated, leading to downstream modulation of genes that regulate the low-temperature response (Zhao et al., 2017). MAPK contains a group of tertiary phosphorylation-dependent kinases: MEKK or MKKK, MEK or Mkk, and MPK. The MAPK cascade activates and mediates signals from the cell membrane to the nucleus, regulating physiological activity (Furuya et al., 2013). When K. obovata seedlings are exposed to cold stress conditions, ca2+/ cams regulate the activity of receptor-like kinases (CRLK1 and CRLK2) and phosphorylates MEKK1, thereby activating the MeKK1-MKK2-MPK4 cascade. Activated MEKK1-MKK2-MPK4 inhibits the protein abundance and kinase activity of Mpk3 and MPK6, actively modulating the hypothermia response, inhibiting cold-mediated activation, and enhancing the cold resistance of plants. (Furuya et al., 2013;Zhao et al., 2017).

Conclusion
Mangroves are an important component of the wetland ecosystem. Kandelia. obovata, as one of the cold-resistant mangrove species, will experience more severe extreme low temperatures as it moves north. We compared gene expression differences in three tissues of K. obovata seedlings before and after cold snaps and found 1,560, 370, and 416 DEGs in the roots, stems, and leaves, respectively. The gene expression strategies of K. obovata seedlings in response to winter cold snaps are summarized as follows: (1)In regulating transcription factors in various tissues, AP2-ERF, NAC, and C2H2 are the most involved transcription factor families in the cold stress response. In addition, bHLH, WRKY, GRAS, bZIP, MYB, and WHSF also participate in the cold stress response. (2)The fluidity of cell membranes in different seedling tissues varies in response to cold pressure, being faster in leaves and slower in roots and stems. Changes in the membrane composition affect Ca 2+ ion transport and, therefore, the CBL-CIPK network. Moreover, (3) the upregulation of secondary metabolism-related genes resulted in the accumulation of starch and sucrose in a seedling's root and stem tissues. In K. obovata seedlings, 13 cold stress-responsive CBL-CIPK network and MAPK cascade genes were identified and used to construct a paired coexpression network model. The findings provide insights into molecular regulatory mechanisms governing the cold stress response of K. obovata seedlings. This study can be used as an empirical reference for further research into the CBL-CIPK network and MAPK cascade in plants.

Data availability statement
The data presented in the study are deposited in the NCBI repository, accession number PRJNA927281. https://www.ncbi.nlm.nih.gov/ sra/PRJNA927281.

Author contributions
KT, and QL: conceptualization, investigation, validation, data analysis, methodology, and visualization, writing-original draft. XZ, PC, SX, HG, and YW: conceptualization, experiment design, reviewing and revising the manuscript. WL: collection of samples, funding acquisition. All authors contributed to the article and approved the submitted version.