Phenotypic adaptation and genomic variation of Kandelia obovata associated with its northern introduction along southeastern coast of China

Introduction: Mangroves play a crucial role within coastal wetland ecosystems, with Kandelia obovata frequently utilized for introduction studies and cultivation research. Investigating the rapid adaptability of K. obovata across diverse environmental conditions offers valuable insights into how mangroves can effectively acclimate to global climate fluctuations.

Methods: In this study, following a common gardenexperiment, we investigated variations in morphological traits between twodistinct populations of K. obovata, Quanzhou (QZ) and Wenzhou (WZ),originating from the same introduction site Zhangzhou (ZZ). Then we performed the whole-genome resequencing on multiple populations along the southern coast of China to assess genetic divergence and diversity patterns in response to environmental factors.

Results: Our findings have uncovered divergent growth-defense trade-off mechanisms employed by these two populations when exposed to varying minimal temperatures in the coldest month within their respective habitats. Moreover, our observations have revealed discernible genetic divergence during the process of environmental acclimatization. Subsequent whole-genome re-sequencing have unveiled a significant decrease in genetic diversity within the northernmost population, suggesting that temperature plays a primary role in shaping genetic variability within the K. obovata species.

Discussion: These findings present new evidence for the rapid adaptation of K. obovata and contributes to our understanding of environmental adaptation characteristics during its introduction to northern regions, which holds significant implications for the conservation and sustainable development of mangroves.

1 Introduction

Given the escalating severity of global climate change, understanding the mechanisms by which organisms adapt to environmental fluctuations has emerged as a pivotal area of scientific investigation in evolutionary biology (Bjorkman et al., 2017). The adaptability of plants stands as a critical determinant influencing their survival capacity within an ever-changing climatic milieu (Foyer and Kranner, 2023). Exploring the adaptive responses exhibited by plants towards diverse environmental alterations can provide valuable insights into their potential coping strategies amidst global climate change. Mangroves represent arboreal ecosystems thriving in tropical and subtropical intertidal zones, characterized by remarkable biodiversity and recognized as one of the foremost carbon-sequestering habitats worldwide. They play indispensable roles in storm resistance, coastal protection, and preservation of coastal ecological equilibrium (Li et al., 2024). Temperature serves as the primary environmental factor delimiting the geographical distribution of mangroves, making them among the most vulnerable groups affected by global climate change (Chen et al., 2017). Consequently, mangrove plants offer ideal experimental subjects for investigating rapid plant responses to global climate change. However, due to intensified human activities and global climate change, mangrove forests are experiencing considerable decline worldwide. Consequently, artificial plantation of mangroves has emerged as an essential strategy for ecosystem restoration (Ellison et al., 2020). There is abundant evidence suggesting that introduced mangrove forests in new environments often undergo adaptive evolution to cope with selection pressures specific to these novel habitats (Miryeganeh, 2022).

Plants adopt phenotypic plasticity and genetic differentiation as two distinct strategies to adapt to their environment (Boquete et al., 2021; Bastias et al., 2024). During the course of evolution, plants must balance these two approaches. Initially, a plant species may utilize plasticity to mitigate natural selection pressure when colonizing a new environment. However, with time and expansion of its geographic range, plant populations may undergo genetic differentiation. In this process, plants can develop various ecotypic traits to adapt to diverse climatic and geographical conditions (Baythavong and Stanton, 2010). Studies have demonstrated that successful invasion by Brachypodium sylvatium is driven by genetic variation rather than phenotypic plasticity induced by the environment (Marchini et al., 2019). A common garden experiment (CGE) can eliminate the influence of environmental factors from provenances on results and assess whether plant populations from different locations have genetic differentiation (Zhao et al., 2021).

Research has revealed that K. obovata, a member of the Rhizophoraceae family, is a widely distributed mangrove species in China renowned for its ability to thrive in high latitudes and withstand cold temperatures (Su et al., 2019). Previous study has indicated the presence of significant genetic differentiation and variation within K. obovata populations along the southeastern coast of China (Zhao et al., 2024). These populations exhibit substantial genetic diversity and may serve as key centers of diversity for Asian mangroves. The genetic diversity of mangroves has been extensively investigated using molecular marker techniques such as Sequence-Related Amplified Polymorphism (SRAP) and Simple Sequence Repeats (SSR) in previous studies (Lu et al., 2021). However, with the rapid advancement of sequencing technology, a chromosome-level reference genome for K. obovata was published in 2020 (Hu et al., 2020). By employing resequencing and bioinformatics methods, studying the genomic level differences and adaptive evolution between populations enables us to comprehend the genetic mechanisms underlying organisms’ adaptation to selection and identify crucial candidate genes. This provides a theoretical foundation for the rational utilization of germplasm resources (Lu et al., 2023).

As part of the exploration work on mangrove protection and introduction in China, K. obovata was introduced from the Jiulong River Estuary Mangrove Provincial Nature Reserve in Zhangzhou city (ZZ), Fujian Province (117°92’E, 24°46’N) to the Luoyang River Mangrove Nature Reserve in Quanzhou Bay, Quanzhou city (QZ), Fujian Province (118°59’E, 24°59’N) in 2003. Subsequently, it was further introduced to Longgang Aojiang Estuary Mangrove, Wenzhou city (WZ), Zhejiang Province (120°96’E, 28°12’N) in 2005. After approximately two decades since its original habitat to the introduction sites, significant phenotypic variations have been observed in these two populations of K. obovata. This study aims to investigate whether these differences are associated with genetic background through CGE. Additionally, using whole-genome re-sequencing (WGRS) technology, this study analyzes and compares the levels of genetic diversity between populations of K. obovata from its original habitat and introduction sites to unravel their phylogenetic relationships and genetic structure. These findings will provide genetic data support for identifying and utilizing germplasm resources of K. obovata. Furthermore, this study explores candidate genes under strong selective pressure across different geographical environments and predicts cis-regulatory elements within their promoter regions to elucidate the impact of environmental adaptation on population differentiation of K. obovata. This provides a theoretical basis for understanding rapid ecological adaptation mechanisms of K. obovata towards environmental changes and strengthens ecological research on global climate response and adaptability of mangroves. It is also significant for mangrove conservation, new variety breeding, and serves as a reference for studying mechanisms by which plants rapidly adapt to environmental changes.

2 Materials and methods

2.1 Common garden experiment

In May 2019, we acquired approximately 20,000 hypocotyls of K. obovata from two distinct locations: the Luoyang River Mangrove Nature Reserve in Quanzhou Bay, Quanzhou City (QZ), Fujian Province, and the Longgang Aojiang Estuary Mangrove in Wenzhou City (WZ), Zhejiang Province. These individuals were then transplanted into plastic pots at the Mangrove Base of Nanhui Dongtan, Shanghai (SH) (121°97’E,30°90’N) for cultivation under CGE conditions. The objective was to assess the overwintering survival rate of K. obovata seedlings under natural conditions (Supplementary Figure 1). In May 2021, we additionally acquired mature hypocotyls of K. obovata from both QZ and WZ. We conducted a random sampling of 30 healthy hypocotyls from each population and recorded their hypocotyl length (HL), hypocotyl weight (HW), and hypocotyl diameter (HD). Subsequently, we carefully selected 1,000 individuals with comparable weights, lengths, and sizes from each population. These individuals were also transplanted into plastic pots in SH for CGE cultivation, with the purpose of subsequent phenotypic observations. The Nanhui Dongtan, located in the Yangtze estuary, is the largest and most extensive coastal wetland in the region. It experiences a subtropical monsoon climate characterized by moderate temperatures and high humidity levels. This area exhibits features of both monsoon and maritime climates. With an average annual temperature ranging from 15 to 16 °C, it has recorded its highest temperature at 37.3 °C and lowest temperature at -7.9 °C.

2.2 Determination of morphological traits

After 18 months of growth in CGE, a total of 30 robust and healthy plants were randomly selected from both the QZ and WZ populations for comprehensive measurement of plant morphological traits. The assessment covered the following parameters: plant height (PH), basal diameter (BD), crown width (W), leaf number (LN), and branch number (BN). Thirty intact leaves were carefully collected from the upper branches of both QZ and WZ populations, thoroughly rinsed with tap water, and then precisely weighed to determine their fresh weight (LFW). After being soaked in distilled water for 12 hours followed by thorough drying to eliminate surface moisture, the leaves were reweighed to obtain their saturated fresh weight (LSFW). Subsequently, scanned images of the collected leaves were analyzed using ImageJ software (1.53t) (Schneider et al., 2012) to calculate various parameters such as leaf length (LL), leaf width (LW), leaf perimeter (LP), leaf area (LA), leaf shape index (LSI = LL/LW). The collected samples were placed in trays and subjected to constant drying at 60°C for 48 hours until reaching a stable weight. The dry weight (DW) was measured to compute relative water content (RWC=(LFW-DW)/(LSFW-DW)) and leaf dry matter content (LDMC=DW/LSFW). Specific leaf area (SLA) is calculated as LA/DW.

2.3 Leaf anatomical structure

We employed the Safranin-O/Fast green staining method to prepare K. obovata samples. Paraffin sections were made with the second leaf from the top down, and ImageJ software (1.53t) was utilized for quantifying the thickness of various leaf components: leaf thickness (Lt), upper cuticle layer thickness (UCu), upper epidermis thickness (UEp), upper hypodermis thickness (UHy), upper palisade tissue thickness (UPt), spongy tissue thickness (St), lower cuticle layer thickness (LCu), lower epidermis thickness (LEp), lower hypodermis thickness (LHy) as well as lower palisade tissue thickness (LPt). The ratio of palisade to spongy tissues (P/S=Pt/St) was calculated along with cell tense ratio (CTR=Pt/LT) and spongy ratio (SR=St/LT).

2.4 Assessment of stomatal characteristics

30 leaf samples were selected from both the QZ and WZ populations in CGE. Temporary slides of the lower epidermis stomata of K. obovata leaves were prepared using the nail polish imprinting method (Pathoumthong et al., 2023), observed, and photographed under a 10x optical microscope. ImageJ software (1.53t) was utilized to quantify the stomata number (SN), stomatal area (SA), and stomatal density (SD) per view field. 10-15 view fields were examined in each sample, and the results were averaged.

2.5 Determination of cold resistance

We measured the daily minimum temperature changes at Nanhui Dongtan Mangrove Base in Shanghai during the coldest months of 2020 and 2021 as shown in Supplementary Figure 2. Despite a warm winter in Shanghai in 2020, an extended period of extreme low temperatures (-7°C) occurred from December 29, 2020 to January 2, 2021. The overwinter survival rates of mangrove plants in their natural environment at the base were separately recorded on January 21, 2020 (without extreme low temperatures) and January 5, 2021 (during the cold wave with extreme low temperatures). The data was analyzed using Microsoft Excel (2021) software. Non-paired Student’s t-test was employed to compare significant differences among different populations in various indicator values. Additionally, a Chi-square test with Yates’ continuity correction was used to assess the differences in overwinter survival rates between the WZ and QZ populations under different temperature conditions. A significance level of p>0.05 indicated non-significance; a range of 0.01<p<0.05 denoted significance; while p<0.01 represented high significance levels. R (4.2.2), Rstudio (2023.03.0 + 386) ggplot2 package and GraphPad Prism (8) software were utilized for graphical representation purposes.

2.6 Whole-genome re-sequencing

In mid to late July 2023, mature K. obovata leaves samples were collected from four locations: Dongzhaigang Mangrove Nature Reserve (DZG) in Haikou City, Hainan Province (110°58’E,19°95’N); Jiulong River Estuary Mangrove Provincial Nature Reserve in Zhangzhou City (ZZ), Fujian Province (117°92’E,24°46 N); Luoyang River Mangrove Nature Reserve in Quanzhou Bay, Quanzhou City (QZ), Fujian Province (118°59’E,24°91’N); and Longgang Aojiang Estuary Mangrove in Wenzhou City (WZ), Zhejiang Province (120°96’E,28°12’N) (Supplementary Figure 1). A total of 40 samples were collected, comprising 10 mature K. obovata leaves from each location, with a minimum distance of over 50 meters maintained between the collection points.

The CTAB method was employed to extract DNA from the leaf samples. Only high-quality DNA samples (OD₂₆₀/₂₈₀ = 1.8~2.0, OD₂₆₀/₂₃₀ = 2.0) were utilized for constructing the sequencing library. A total of 0.5 μg of DNA per sample served as input material for the preparation of the DNA library. The Truseq Nano DNA HT Sample Prep Kit (Illumina USA) was used to generate the sequencing library in accordance with the manufacturer’s recommendations, and index codes were assigned to each sample. In brief, genomic DNA samples were sonicated to obtain fragments with a size of 350 bp, followed by end-polishing, A-tailing, and ligation with full-length adapters suitable for Illumina sequencing technology; this was succeeded by additional PCR amplification steps. After purification of PCR products using the AMPure XP system, libraries underwent size distribution analysis via Agilent 2100 Bioanalyzer and quantification through real-time PCR (3nM). Finally, paired-end DNA-seq sequencing libraries were sequenced on an Illumina NovaSeq system at Shanghai Majorbio Bio-pharm Technology Co., Ltd.

2.7 Variant discovery

The raw reads of low quality (mean phred score < 20), which included reads containing adapter contamination and unrecognizable nucleotides (N base > 10), were trimmed or discarded using Fastp software (0.23.2) (Chen et al., 2018). After trimming, the reads were mapped to the reference genome (https://bigd.big.ac.cn/gwh/Assembly/990/show) using BWAMEME software (1.0.5) (Jung and Han, 2022) with default mapping parameters. The alignment bam files were sorted by SAMtools (1.15.1) (Li et al., 2009) and PCR duplicates were marked using MarkDuplicated as part of the modified GATK Best Practice pipeline (4.3.0.0) (Mckenna et al., 2010). Base quality recalibration was performed, followed by germline variant calling for Single Nucleotide Polymorphisms (SNPs) across all samples using Haplotyper and Gvcftyper programs in Sentieon genomics tools (202112.07) (Freed et al., 2017). Variants were filtered according to standard hard filtering parameters based on GATK Best Practices pipeline (4.3.0.0), and annotated using SnpEff (5.1d) (Cingolani et al., 2012). Subsequently, several filtering steps were applied to reduce false positives for SNPs and genotype calling using VCFtools software (0.1.16) (Danecek et al., 2011): (if) SNPs with more than two alleles were removed, (ii) SNPs with mean depth values less than 4 across all samples were removed, (iii) SNPs with minor allele frequency < 0.05 were removed, (iv) Only SNPs that could be genotyped in at least 70% of the samples were retained, (v) For population structure analysis, SNPs showing linkage disequilibrium patterns were pruned using Plink software (1.90b6.20) (Purcell et al., 2007).

2.8 Genetic diversity analysis

Based on filtered vs. files, we calculated observed heterozygosity (Ho), expected heterozygosity (He) and nucleotide diversity (π) using the populations module in the software Stacks (2.64) (Catchen et al., 2013). Additionally, GenAlEx (6.5) (Peakall and Smouse, 2012) was employed to determine genetic diversity parameters such as Polymorphic Information Content (PIC) and Shannon’s Information Index. These parameters were utilized to evaluate the level of genetic diversity in four populations of K. obovata studied here. A higher value for these genetic diversity parameters indicates a greater level of genetic diversity within the population.

The Genetic Differentiation Index (F_ST) between populations was calculated using the populations module in the Stacks software (2.64) (Catchen et al., 2013). F_ST values were obtained through pairwise comparisons among all populations to evaluate the extent of genetic differentiation between two populations. The level of genetic differentiation can be determined based on a range of F_ST values: 0~0.05 indicates an extremely low degree of genetic differentiation between populations; 0.05~0.15 suggests a moderate level of genetic differentiation; 0.15~0.25 signifies a considerable degree of genetic differentiation between populations; and when F_ST > 0.25, it indicates a high level of genetic differentiation among populations (Wright, 1978).

2.9 Phylogenetic analyses

The Maximum Likelihood phylogenetic tree was constructed using IQ-TREE2 software (2.1.2) (Minh et al., 2020). The ML analyses were performed on the pruned SNP sites employing IQ-TREE2 with GTR+I+G4 model and 1000 bootstraps.

The unsupervised maximum-likelihood clustering algorithm implemented in ADMIXTURE (1.3.0) (Alexander et al., 2009) was used to cluster each genome in the investigated populations of K. obovata. Initial clustering was performed for K = 1 to K = 20 ancestral clusters using default settings. The optimal value of K was determined based on the cross-validation error (CV error), and the genetic structure corresponding to this optimal value of K was outputted as the final result. To enhance the accuracy of initial clustering, pruned SNPs were utilized in the structural analysis.

To visualize the genetic relationships among samples, we conducted principal component analysis (PCA) based on pruned SNPs using Plink (1.90b6.20) (Purcell et al., 2007).

2.10 Screening for selective sweep and identification of positively selected genes

Following the sliding window strategy, we utilized PIXY software (1.2.7.beta1) (Korunes and Samuk, 2021) to perform segmentation and calculation of F_ST as well as π between QZ and WZ populations. The window size was set at 10 Mb with a step size of 10 kb. To strengthen the analysis of positive selection, we employed a permutation test to account for random differentiation that may occur between populations in the absence of selection pressure. Specifically, we performed 1,000 permutation tests by randomly shuffling the population labels and recalculating the F_ST distribution for each 10 kb window. This approach allowed us to simulate different population structures under the null hypothesis of no selection. The 95^th percentile of the permutation distribution was used as a baseline to determine whether the observed F_ST values were significantly higher than expected by chance. Finally, we employed a screening criterion where regions with F_ST values in the top 5% were identified as potential candidate regions distinguishing QZ from WZ populations. Furthermore, our selective sweep analysis combined both F_ST and π values. Specifically, regions of the genome exceeding the top 5% threshold, along with regions exhibiting extremely high π ratios (π_QZ/π_WZ), were considered as potential areas displaying strong signals of selection scanning for WZ population. Sequence Toolkit in TBtools (2.112) (Chen et al., 2020) was used to associate the selected region with candidate genes to obtain the positively selected genes of WZ population.

The TBtools software (2.112) was also employed to conduct Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis on the identified regions from the QZ and WZ populations. Significance of enrichment was determined when p-value < 0.05 for both GO and KEGG pathways. Visualization was performed using the qqman and ggplot2 packages in R (4.2.2) and Rstudio (2023.03.0 + 386) software.

2.11 Analysis of cis-element in promoters

After extracting the 2,000 bp upstream sequence of the positive selection gene in the WZ population using TBtools software (2.112), we conducted promoter cis-element analysis through the PlantCare database (Lescot et al., 2002).

2.12 Integrated analysis of transcriptome and genome data

We retrieved previously generated transcriptome data (Zhang et al., 2024), and intersected the differentially expressed genes (DEGs) under cold stress from both populations with 40 candidate positively selected genes identified from the genome based on F_ST and π values. These genes were further analyzed for transcriptome expression profiles and promoter elements. Heatmap visualization was performed using ggplot2 and pheatmap packages in R (4.2.2) and RStudio (2023.03.0 + 386). Transcriptional levels and cis-elements were visualized using TBtools Simple BioSequence Viewer (2.112) and GraphPad Prism (8).

3 Results

3.1 Morphological differences between two populations of K. obovata

In comparison to the QZ population, the WZ population’s hypocotyls exhibit relatively lower values for physiological parameters (Supplementary Table 1). There are statistically significant differences (p<0.05, non-paired Student’s t-test) in HL, HW, HD, PH, BD, LN, and BN, as well as W between these two populations. Moreover, significant differences between the two populations (p<0.05, non-paired Student’s t-test) are also observed in SLA, and DW (Figure 1). However, no notable differences are found in LW, LSI, LA, RWC or LDMC.

Figure 1

Figure 1. Differences in morphology and leaf traits between two populations of K. obovata. (A) Comparison of morphological traits between two populations of K. obovata. HL, hypocotyl length; HW, hypocotyl weight; HD, hypocotyl diameter; PH, plant height; BD, basal diameter; BN, branch number; W, crown width; LN, leaf number; SLA, specific leaf area; DW, dry weight; Lt, leaf thickness; St, spongy tissue thickness; UPt, upper palisade tissue; UHy, upper hypodermis; P/S, palisade-spongy tissue ratio; CTR, cell tense ratio; SR, spongy ratio; SA, stomatal area; SD, stomatal density. * represents p<0.05; ** represents p<0.01, *** p represents <0.001. (B) Comparison of leaf anatomy and lower epidermal stomata between two populations of K. obovata. UCu and LCu, upper and lower cuticles; UEp and LEp, upper and lower epidermis; UHy and LHy, upper and lower hypodermis; UPt and LPt, upper and lower palisade tissue; St, spongy tissue; LT, leaf thickness.

3.2 Disparities in leaf anatomical structure and stomatal characteristics between two K. obovata Populations

The cross-sectional images of K. obovata reveal a leaf structure comprising upper and lower cuticles, epidermis, hypodermis, palisade tissue, and spongy tissue. Both the upper and lower epidermis consist of a single layer of cells, while the hypodermis consists of two cell layers—the first layer being smaller without tannins and the second layer larger containing tannins. The densely arranged palisade tissue has a thicker upper layer than the lower layer. The spongy tissue, which accommodates vascular bundles, displays a looser arrangement (Figures 1A, B). Furthermore, the Lt, St, and SR are significantly greater in the QZ population compared to those in the WZ population (p<0.05, non-paired Student’s t-test). Conversely, the WZ population exhibits significantly higher thickness in UHy, UPt, P/S, and CTR compared to the QZ population (p<0.05, non-paired Student’s t-test). Stomata of K. obovata are predominantly distributed in the lower epidermis (Figure 1B). The SA is significantly smaller in QZ population than that in WZ population (p<0.05, non-paired Student’s t-test), while SD is significantly larger than that in WZ population (p<0.05, non-paired Student’s t-test) (Figure 1A).

3.3 Variation in Overwintering Survival Rate between Two K. obovata Populations under Natural Conditions

One of the primary constraints on the northward expansion of K. obovata is its susceptibility to cold temperatures. We investigated the overwinter survival rates of both populations under natural conditions following the CGE (Table 1). In January 2020 with moderate winter temperatures, the survival rate of K. obovata was notably high. Specifically, the natural overwinter survival rate of WZ population was 84.5%, surpassing that of the QZ population (67.3%). Both populations exhibited reduced survival rates under extreme low temperatures in 2021; however, even under such conditions, the WZ population demonstrated a higher winter survival rate (18.6%) compared to that of the QZ population (6.8%). Our findings indicate a significant difference in overwinter survival rates between the two populations both under extreme low temperatures and moderate winter conditions (p< 2.2e-16, Chi-square test). The WZ population exhibits superior cold tolerance and better adaptation to the northern CGE environment.

Table 1

Table 1. Natural overwinter survival rates of K. obovata in 2020 and 2021.

3.4 WGRS, SNP variation detection and annotation

In this study, we selected samples from four population: DZG as the outgroup, ZZ as the sample from original location, QZ, and WZ. 10 samples were collected from each population, resulting in a total of 40 samples. WGRS was performed on these 40 samples. Following rigorous quality control and raw data filtering, we obtained a total of 92.98 Gb of clean read data with Q30 sequences exceeding 93.36% and GC content ranging from 36.35% to 39.64% (Supplementary Table 2). Alignment of the clean data with the reference genome yielded 45.82 Gb of clean data, indicating an average alignment rate of 92.62% between the samples and the reference genome. The genome coverage rate was calculated at 90.14%, with an average sequencing depth of 12.20× and an average coverage rate (≥4×) (Supplementary Table 3). The uniform coverage across the entire genome indicated excellent sequencing randomness (Supplementary Figure 3).

After filtration, a total of 117,774 high-quality SNPs were identified in the population. Subsequent analysis revealed 69,692 transitions (ts) and 48,082 transversions (tv), resulting in a ts/tv ratio of 1.45 (Supplementary Table 4). Statistical analysis based on the annotation information of the K. obovata genome indicated that the majority of SNPs were located in intergenic regions, with only approximately 20% found within genic regions (Supplementary Table 4). Specifically, whole-genome assessment of WZ and QZ identified a total of 50,330 and 44,418 SNPs respectively, with approximately 20.68% and 19.40% located within genic regions (Supplementary Figure 4).

3.5 Analysis of genetic diversity and relationship among different populations of K. obovata

SNP analysis was employed to investigate the genetic relationships among 40 samples. The DZG population was utilized as an outgroup, and a phylogenetic tree was constructed using the maximum likelihood method. The results revealed challenges in distinguishing between the ZZ and QZ populations, indicating low genetic differentiation and close genetic relationship. However, the WZ population exhibited clear differentiation from the ZZ and QZ populations, suggesting substantial genetic divergence between WZ and ZZ/QZ populations (Figure 2A). PCA revealed that PC1 effectively distinguished the DZG outgroup from other populations, while PC2 differentiated the WZ population from ZZ and QZ populations. Only a slight difference was observed between ZZ and QZ populations on this principal component, indicating their close genetic relationship (Figure 2B). STRUCTURE analysis revealed that a K value of 3 resulted in the lowest error rate of variation coefficient (Supplementary Figure 5), leading to the division of the 40 samples into three clusters. This grouping combined QZ and ZZ into one cluster while separating DZG and WZ (Figure 2C).

Figure 2

Figure 2. Kinship analysis in different populations of K. obovata. (A) Phylogenetic tree constructed using the maximum likelihood method; (B) PCA of the four populations; (C) Genetic structure of the four populations based on the WGRS. (D) Results of genetic diversity for four populations. The numbers in brackets represent the nucleotide diversity index (π) and the numbers in the lines represent the fixation index (F_ST) between the two populations.

In addition, we computed π, Ho, He, PIC, and Shannon index for each population of K. obovata to evaluate their genetic diversity levels. The findings revealed that the DZG population exhibited the highest genetic diversity level among the four K. obovata populations, followed by QZ and ZZ; whereas the WZ population displayed the lowest genetic diversity level (Table 2). Furthermore, we assessed inter-population genetic differentiation using genetic differentiation index F_ST. The results indicated minimal genetic differentiation between the QZ and ZZ populations, while moderate genetic differentiation was observed between the WZ and QZ as well as ZZ populations (Figure 2D).

Table 2

Table 2. Genetic diversity index statistics among K. obovata populations.

3.6 Detection of selection sweep signals and analysis of cis-elements in candidate genes

Using the 1,000 permutation test, we observed that the 95^th percentile of the F_ST distribution was 0.243 in the absence of selection pressure. In contrast, the real population’s top 5% F_ST threshold was 0.534 (Supplementary Figure 6). The fact that the observed F_ST value exceeds the 95^th percentile of the permutation distribution indicates that these high F_ST values are not due to random fluctuations and are more likely to reflect genuine selection pressures. Thus, the top 5% regions of F_ST in K. obovata populations from QZ and WZ were identified based on F_ST scans, resulting in the detection of 785 windows through Manhattan plots. Following alignment with the reference genome annotation file, a total of 1,159 genes were identified (Figure 3A, Supplementary Tables 5-7). KEGG pathway enrichment analysis of these 1,159 candidate genes revealed significant enrichment in 9 pathways (Supplementary Table 8). Specifically, 15 genes showed enrichment in cysteine and methionine metabolism, while 9 genes exhibited enrichment in glutathione metabolism, along with 8 genes involved in alanine, aspartate, and glutamate metabolism; additionally, 7 genes were associated with ubiquinone and other terpenoid-quinone biosynthesis (Figure 3B).

Figure 3

Figure 3. Analysis of selective sweep signals and cis-elements of candidate genes between QZ and WZ populations. (A) Manhattan plot of F_ST distribution with chromosomes. The horizontal axis represents chromosomes, the vertical axis represents F_ST between QZ and WZ populations, and the dashed line is the threshold line, with the default value of top5%; red loci beyond the threshold line are those with significant selective sweep effects. (B) KEGG functional enrichment of genes. The horizontal axis represents rich factors, the vertical axis represents functionally enriched pathways, bubble colors indicate enrichment significance, and bubble sizes indicate the number of genes in the gene set for the candidate gene. (C) Selected regions of WZ populations based on the combination of F_ST and π ratio (QZ/WZ) screened at the top5% level. The horizontal axis represents the ratio of nucleotides, and the vertical axis represents F_ST, with WZ selected regions in red. (D) GO functional enrichment analysis under selection signaling of WZ, with the horizontal axis representing significance, the vertical axis representing functionally enriched terms, and the number in the circle denoting the number of selected genes in the gene set. (E) Prediction of cis-acting elements in the promoter region of some WZ selected genes. The cis-acting elements were classified into three categories according to their functions, with green indicating abiotic and biotic stresses, red indicating phytohormone responsive, and blue indicating plant growth and development. The numbers in the squares on the left represent the number of cis-elements per gene, and the numbers on the right represent the number of cis-elements with different functions per gene.

Furthermore, employing a method that integrates genetic differentiation coefficient (F_ST) and nucleotide diversity (π), we identified 13 regions exhibiting positive selection in WZ, falling within the top 5% of the region between F_ST and log₂(π_QZ/π_WZ) in WZ (Figure 3C). Upon aligning these selected windows with the reference genome annotation file, we determined a total of 40 candidate genes undergoing positive selection in WZ (Table 3; Supplementary Table 9). GO functional enrichment analysis revealed significant enrichment of biological processes related to stress response, anatomical structure development, sexual reproduction, reproductive development process, and post-embryonic development among the positively selected genes in WZ (Figure 3D, Supplementary Table 10).

Table 3

Table 3. Candidate genes identified by the top5% F_ST and π ratio value between WZ and QZ populations (40 genes).

We utilized the 2,000 bp upstream sequence as the promoter region for investigating 40 genes exhibiting positive selection in the WZ population. Our analysis revealed an abundance of cis-elements associated with non-biological stress responses (Supplementary Figure 7). The findings demonstrate widespread presence of stress-responsive elements (STR), low temperature-responsive elements (LTR), MYB-binding sites (MBS) induced by drought, and wound-responsive elements (WUN-motif) within the promoters of positively selected genes in the WZ population. Furthermore, we have observed a wide distribution of diverse cis-elements, including abscisic acid responsive element (ABRE), salicylic acid responsive element (as-1), ethylene responsive element (ERE), jasmonic acid responsive elements (CGTCA-motif and TGACG-motif), as well as numerous light-responsive elements such as Box4, across these promoters. (Figure 3E).

3.7 Analysis of the relationship between selected genes and cold stress response in the WZ population

In our previous study, we examined the transcriptome profile of QZ and WZ during a simulated cold wave (Supplementary Figure 8), identifying 3,810 DEGs (Zhang et al., 2024). In this study, we conducted an analysis of 1,159 genes exhibiting positive selection, as identified from the F_ST index scan (within the top 5%) of the ZZ and WZ populations based on resequencing data. Upon taking the intersection of these gene sets, we obtained 215 genes (Figures 4A, B; Supplementary Table 11), which have shown differential transcription changes in the cold wave (Figure 4C). These genes were subsequently subjected to SNP mining and functional annotation analysis. Our findings revealed a total of 592 variants, including 206 variations within 1 kb region upstream of transcription start site (TSS), 195 variations within 1 kb region downstream of transcription end site, 96 variations within introns, as well as 39 missense mutations and 18 synonymous mutations (Supplementary Table 12).

Figure 4

Figure 4. Integrated analysis of transcriptome and genome. (A) A Venn diagram illustrating the overlap of genes identified through F_ST-based genomic positive selection analysis and transcriptome DEGs under cold wave, resulting in 215 genes. (B) Visualization of SNP variation types and their chromosomal locations within the 215 genes, with each row representing a distinct chromosome and different colors denoting various variation types. (C) A heatmap displaying the transcription level changes of the 215 genes under cold wave. (D) The impact of upstream base mutations on cis-elements of K. obovata genes in two populations, indicated by red arrows and text depicting cis-elements before and after mutations. (E) The differential transcription levels of selected genes in the cold wave (* p<0.05).

The differential expression of genes can be attributed to the presence of distinct cis-acting elements within various gene promoters. Variations in upstream sequences have the potential to modify these cis-acting elements and binding sites for transcription factors, thereby influencing their binding capacity and efficiency, consequently impacting gene expression levels. Genes exhibiting variations in their upstream sequences were then screened, and promoter analysis was conducted for these varied upstream sequences. Substantial disparities in cis-regulatory elements between the WZ and QZ populations have been found. For instance, within geneMaker00002377, substitution of A for G at nucleotide position 1656 in the QZ population results in an additional stress response element (STRE) present exclusively in the WZ population. Similarly, within the geneMaker00002517 encoding drought-induced protein, there is an extra abscisic acid response element (ABRE) present solely in the WZ population compared to QZ; this element can be recognized by bZIP transcription factors leading to increased expression under cold conditions and thus enhancing cold tolerance of K. obovata (Sun et al., 2022). Furthermore, there is an elevated C allele frequency at nucleotide position 1473 within geneMaker00005748 in WZ population which encodes K. obovata homologue of E3 ligase SHOOT GRAVITROPISM9 (SGR9) (Nakamura et al., 2011), resulting in an additional AE-box response element which may enhance cold tolerance. Additionally, the substitution from G to A at nucleotide position 1333 within extracellular signal-regulated kinase (ERK) activator CEP14-encoding gene (geneMaker00015529) (Delay et al., 2013) leads to an inclusion of MYB element exclusively found in WZ that could potentially bolster plant recovery ability under low temperature stress during later stages (Figures 4D, E; Table 4). We propose that these differences may account for alterations in target gene transcription levels between the two populations and potentially contribute to changes in plant cold tolerance.

Table 4

Table 4. Effects of upstream variation on cis-elements and functions of the K. obovata..

We also identified a significant number of mutations in the downstream gene sequences of DEGs. For example, the CBF3 gene (geneMaker00008995), which constitutes the core of the cold stress pathway and plays an essential upstream initiation role in K. obovata cold resistance (Peng et al., 2020), exhibits an AA base deletion at locus 938050 on chromosome 9 in the WZ population. Furthermore, the ABI5 gene (geneMaker00005324), a negative feedback factor in abscisic acid (ABA)-signaling involved in regulating ABA signaling pathway and ROS levels (Collin et al., 2021), shows a C to A base replacement at locus 2261919 on chromosome 5 within the WZ population. The potential relationship between these variants and changes in transcription levels of corresponding genes remains uncertain yet.

4 Discussion

4.1 Variation in phenotypic and leaf functional traits of K. obovata from different populations

Under the influence of natural and artificial selection, plants will inevitably differentiate in terms of phenotype and ecological traits as they adapt to new environments, gradually forming distinct geographical sources and leading to significant differences in morphology and functional leaf traits. Morphological indices such as leaf thickness, leaf dry matter content, leaf water content, and specific leaf area can reflect the plant’s adaptive features and ability to acquire resources in different habitats. Specific leaf area is typically closely linked to the plant’s growth strategies for survival. Plants with low specific leaf area often thrive in harsher or less fertile environments, while those with high specific leaf area can effectively maintain their nutrient content (Meziane and Shipley, 1999). Furthermore, as an index of ecological adaptation, leaf anatomical structure is closely related to a plant’s cold tolerance; it is frequently used as a key metric for evaluating plant cold tolerance. Generally speaking, the higher the ratio of palisade to mesophyll tissue density within leaves (and thus tighter structural density), the greater their cold tolerance becomes. In response to adverse conditions such as high temperatures, small dense stomata are advantageous for avoiding rapid water loss due to transpiration.

QZ and WZ populations have been introduced into Fujian Quanzhou and Zhejiang Wenzhou within the last 20 years. Despite this study’s elimination of native habitat effects through a CGE, these populations have already displayed distinctive phenotypic variations and leaf traits. The research findings suggest that the QZ population at lower latitudes exhibits greater PH, BD, W, LN and BN. Moreover, there are notable disparities in leaf functional traits and anatomical structures between the two populations. The WZ population demonstrates a larger SLA, potentially attributed to its robust photosynthetic capacity facilitating rapid resource acquisition and internal nutrient maintenance. Its palisade tissue is thicker with a denser cell structure, resulting in a higher P/S that enhances cold resistance. Conversely, the QZ population features smaller SA and higher SD possibly linked to relatively elevated local temperatures and strong transpiration for favorable water balance. In contrast, the WZ population displays larger but sparser stomata which may enhance its resilience to adverse conditions by reducing transpiration rate (Figures 1A, B).

It is widely accepted that there exists a trade-off between plant growth and defense. When plants are frequently exposed to pathogens or adverse environmental conditions, they prioritize defense mechanisms, leading to a reduction in growth and reproductive capacity (Zhou et al., 2022). We posit that the seedlings of WZ adopted a distinct energy allocation strategy compared to those of QZ in order to more effectively acclimate to the high-latitude winter environment. Additionally, the functional characteristics of the leaves underwent changes aimed at enhancing their adaptability to the environment, with various leaf indices demonstrating a synergistic and trade-off relationship.

Consequently, we postulate that the WZ population of K. obovata inhabiting relatively cold environments may exhibit adaptive responses by modifying its growth-defense trade-off strategy. This could involve allocating more resources to defense mechanisms, slowing down the growth rate, increasing leaf tissue density and structural compactness, reducing SLA, and regulating stomatal characteristics as an adaptive response to the challenging habitat. These adaptive adjustments are likely to enhance the plant’s resistance to low temperatures and improve its environmental adaptability and resource utilization. Conversely, the QZ population in milder habitats may allocate greater energy towards growth and reproduction. In conclusion, this study unveils significant phenotypic disparities between WZ and QZ, suggesting the adaptive capacity of K. obovata to environmental fluctuations. The heterogeneity of the environment over the approximately 20-year adaptation period may underlie the observed differentiation in phenotypic traits, leaf functionality, and physiological characteristics between QZ and WZ populations.

4.2 Genetic diversity and kinship analysis of different populations in the same introduction

The research conducted by Yang et al. demonstrates a substantial variation in both hypocotyl and seedling growth traits of K. obovata across diverse collection sites, showing a strong correlation with geographical and climatic factors (Yang et al., 2020). Prolonged acclimatization to its indigenous habitat can result in genetic divergence concerning plant morphology and leaf functional traits, thereby ensuring stable inheritance of these variations across successive generations. Our findings demonstrate that the mangrove species K. obovata, which was introduced from ZZ to WZ and QZ within a brief two-decade period, rapidly acclimated to the local environment through the development of specific traits. Significant phenotypic differences were already evident in the hypocotyls (Figure 1; Supplementary Figure 1).

The CGE further confirms that K. obovata seedlings also demonstrate significant variability in phenotype, leaf functional traits, and cold tolerance between the QZ and WZ populations (Figure 1, Table 1). The results of PCA and phylogenetic tree analysis based on WGRS indicate that ZZ and QZ are intertwined and clustered into the same branch, while WZ forms a separate cluster (Figures 2A, B). Additionally, there exists moderate genetic differentiation between the WZ and QZ populations, with experimental findings consistent with the geographic distribution of K. obovata. These results suggest that the observed phenotypic differences between the two K. obovata populations are primarily attributed to genetic differentiation rather than phenotypic plasticity.

This study utilized WGRS technology to assess the genetic diversity levels in various K. obovata populations. The analysis involved calculating heterozygosity (He, Ho), nucleotide diversity, polymorphic information content, and Shannon index. The findings revealed a decrease in genetic diversity of K. obovata populations in China with increasing latitude. At the whole-genome level, DZG exhibited the highest genetic diversity at lower latitudes, while ZZ and QZ displayed similar levels of genetic diversity and WZ showed the lowest (Table 2). Furthermore, observed heterozygosity (Ho) for ZZ, QZ, and WZ populations exceeded expected heterozygosity (He) (Table 2), suggesting higher frequency fluctuations possibly due to founder effect and bottleneck effect during the northward expansion of K. obovata driven by human selection of breeding individuals. The initial genetic bottleneck during the establishment of the WZ population likely resulted in limited genetic variation, leading to random phenotypic changes, some of which may be attributed to the founder effect. However, as the WZ population adapted to low-temperature environment, natural selection likely accelerated the optimization of these variants, driving the positive selection of cold-tolerance-related genes. For instance, we identified significant enrichment in pathways such as “related to stress response,” and “glutathione metabolism” (Figure 3D), along with selection signals linked to ABA/JA and low-temperature response elements, which are crucial for cold stress adaptation (Figure 3E). The WZ population has developed robust cold tolerance under prolonged exposure to low temperatures, and this trait has remained genetically stable across multiple generations, as confirmed by CGE in Shanghai (Table 1; Figure 1). The WZ population consistently exhibited enhanced overwinter survival rates (Table 1) and stronger cold-resistant phenotypes, accompanied by adaptive leaf anatomical changes, including thicker palisade tissue and a higher palisade-to-spongy ratio (Figure 1). These traits showed consistent inheritance, supporting the role of adaptive evolution in shaping this phenotype. Thus, while the founder effect may have initially constrained genetic diversity in the WZ population, the enhanced cold tolerance is more likely the result of adaptive evolution driven by natural selection. In conclusion, the phenotypic changes observed in the WZ population reflect a combination of adaptive evolution and the founder effect. Further studies, such as whole-genome association studies (GWAS) or quantitative trait locus (QTL) analysis, are needed to explore other phenotypic differences beyond cold tolerance.

Transcription factors regulate the expression of a wide range of biological and abiotic stress response genes by interacting with cis-acting elements in the promoter region, enabling plants to adapt to diverse stresses. To further explore potential stress response elements in the promoters of robust WZ populations, we utilized the Plantcare software for cis-element prediction in their promoter regions. The analysis revealed a substantial presence of cis-elements associated with abiotic stress response, such as LTR and hormone response elements (Figure 3E), suggesting their significant role in plant growth, development, low-temperature stress responses, and modulation of hormone signaling molecules. Endogenous plant hormones play a pivotal role in abiotic stress (Waadt et al., 2022). This investigation identified ABRE within almost all positively selected genes. Furthermore, most genes contain jasmonic acid (JA) response elements CGTCA-motif and TGACG-motif within their promoter sequences which facilitate ABA- and JA-dependent signaling pathways through interaction with upstream transcription factors. ABA and JA, essential endogenous plant hormones, are crucial in mediating abiotic stress responses and have been shown to significantly impact on low temperature stress in plants (Huang et al., 2017; Hu et al., 2017). Studies have shown that exogenous ABA application can significantly improve cold tolerance in many plants, such as Capsicum annuum (Guo et al., 2012) and Vitis vinifera (Wang et al., 2020). Similarly, JA positively regulates the ICE-CBF pathway to enhance cold tolerance in Arabidopsis thaliana, and exogenous application of methyl jasmonate (MeJA) significantly improves cold tolerance in plants (Hu et al., 2013, Hu et al., 2017). Notably, exogenous ABA can alleviate cold stress in K. obovata by activating antioxidant enzyme activities and promoting the accumulation of osmotic regulators, thereby mitigating the negative effects of cold stress (Liu et al., 2022). Furthermore, in our previous research, we found that the enhanced cold tolerance in the WZ population is linked to JA signaling molecules and exogenous application of MeJA reduced cold-induced damage in the QZ population (Zhang et al., 2024). These findings suggest that these genes may be regulated by ABA and JA pathways thereby promoting K. obovata’s response to low temperature stress.

4.3 Low Temperature Maybe a Major Factor in the Genetic Differentiation of K. obovata

Climate exerts significant selective pressure on plant traits, with low temperature playing a crucial role in influencing the distribution of mangrove plants at high latitudes. Previous studies have shown that temperature is a major environmental factor limiting the distribution of mangroves (Cavanaugh et al., 2014). For instance, extreme freezing events are considered as a primary environmental factor determining the latitudinal distribution of K. obovata (Su et al., 2019; He et al., 2023) or other species such as Avicennia germinans (Pickens and Hester, 2011; Osland et al., 2014). A related study integrated geographical distribution data for K. obovata from both natural and introduced populations along China’s southeastern coast and combined it with climate and hydrological data to quantitatively analyze the relationship between geographical distribution patterns and key environmental factors (Zhu et al., 2021). The results indicated that the top five climate factors influencing K. obovata’s distribution were annual mean temperature, mean temperature of the coldest quarter, extreme minimum temperature, temperature seasonality, and mean temperature of the driest quarter. These findings are consistent with our observations. On the one hand, phenotypic and genomic analyses revealed temperature-associated variations, including selection signatures at cold-responsive elements and promoter polymorphisms in cold adaptation genes (Table 1; Figures 3, 4). On the other hand, based on data from the China Meteorological Science Data Sharing Service Network (http://data.cma.cn), the average temperature and lowest temperature of the coldest month in the WZ and QZ regions over the past 75 years are as follows: 8.1°C and 12.0°C for the average temperatures of the coldest month in the WZ and QZ regions respectively; -0.9°C and 5.5°C for their respective average extreme low temperatures in the coldest month (Supplementary Figure 9; Supplementary Table 13). These findings indicate that temperature likely serves as a primary driver of genetic differentiation between populations of K. obovata introduced from related source sites (WZ and QZ) due to substantial differences in latitude, average temperature, and other environmental conditions. However, it should be noted that while temperature gradients strongly correlate with latitudinal variation, other environmental parameters may synergistically shape genetic differentiation. Although WZ is located at a higher latitude and experiences colder temperatures, the observed genetic differentiation may result from a complex interplay of multiple environmental pressures, rather than temperature alone. Future studies that collect additional long-term environmental data—such as salinity, precipitation, and soil parameters—from both QZ and WZ will provide a more comprehensive analysis of the multi-layered environmental factors influencing the genetic structure and distribution of K. obovata.

Following transplantation to higher latitudes, acclimatization processes occurred within the ZZ population, leading to elimination of mangrove plants ill-suited for low temperatures while retaining those adapted to such conditions, ultimately giving rise to today’s WZ population which demonstrates enhanced adaptability to low temperatures at high latitudes (Table 1). The molecular mechanisms underlying adaptation to low temperatures by the WZ population have been partially elucidated (Su et al., 2019). In summary, both natural selection and genetic drift play crucial roles in shaping the adaptation of K. obovata to its environment. After approximately 20 years, phenotypic and nucleotide-level genetic differentiation between the WZ and ZZ populations indicates that K. obovata possesses potential for rapid adaptation, rendering it a promising candidate species for studying plant adaptability.

5 Conclusion

In conclusion, this study investigated the phenotypic and functional trait variations of WZ and QZ populations in common garden experience. The findings revealed that different populations exhibited distinct growth-defense trade-off strategies to adapt to their respective environments. Within the same introduction site, moderate genetic differentiation was observed among K. obovata populations, indicating the rapid adaptation ability of K. obovata to environmental changes. WGRS further elucidated the genetic structure of K. obovata populations within the same introduction site, with results showing that the northernmost artificially introduced population (WZ) displayed the lowest genetic diversity. Additionally, its gene promoter regions affected by strong selection contained a significant number of stress response elements related to low temperature and hormones, suggesting that temperature may be a primary driver of genetic differentiation in K. obovata. These research findings contribute to our understanding of environmental adaptation characteristics and growth/survival strategies of K. obovata in introduction areas, providing theoretical support for its cultivation as well as for mangrove forest protection and development efforts.

Data availability statement

The WGRS data presented in the study are deposited in the NCBI repository, accession number PRJNA1145277 (https://www.ncbi.nlm.nih.gov/sra/PRJNA1145277). The transcriptome data cited in this study are deposited in the NCBI repository, accession number PRJNA1093421 (https://www.ncbi.nlm.nih.gov/sra/PRJNA1093421).

Author contributions

JZ: Writing – original draft. SO: Writing – original draft. XC: Writing – original draft. SY: Writing – original draft. QC: Writing – original draft. JY: Writing – original draft. ZS: Writing – original draft. WZ: Writing – original draft. YW: Writing – original draft. YZ: Writing – original draft. PN: Writing – original draft.

Funding

The author(s) declare that financial support was received for the research and/or publication of this article. The authors declare that this study received funding from Morgan Stanley Securities (China) Co., Ltd. and Shanghai Genshi Medical Technology Development Co., Ltd. and the National Natural Science Foundation of China (32270629 and 32070201). The funders were not involved in the study design, collection, analysis, interpretation of data, the writing of this article or the decision to submit it for publication.

Acknowledgments

The present manuscript commemorates professor Yang Zhong’s enduring vision and exceptional contributions to the mangroves in China.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The author(s) declare that no Generative AI was used in the creation of this manuscript.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2025.1512620/full#supplementary-material

References

Alexander, D. H., Novembre, J., Lange, K. (2009). Fast model-based estimation of ancestry in unrelated individuals. Genome Res. 19, 1655–1664. doi: 10.1101/gr.094052.109

PubMed Abstract | Crossref Full Text | Google Scholar

Bastias, C. C., Estarague, A., Vile, D., Gaignon, E., Lee, C. R., Exposito-Alonso, M., et al. (2024). Ecological trade-offs drive phenotypic and genetic differentiation of Arabidopsis thaliana in Europe. Nat. Commun. 15, 5185. doi: 10.1038/s41467-024-49267-0

PubMed Abstract | Crossref Full Text | Google Scholar

Baythavong, B. S., Stanton, M. L. (2010). Characterizing selection on phenotypic plasticity in response to natural environmental heterogeneity. Evolution 64, 2904–2920. doi: 10.1111/j.1558-5646.2010.01057.x

PubMed Abstract | Crossref Full Text | Google Scholar

Bjorkman, A. D., Vellend, M., Frei, E. R., Henry, G. H. (2017). Climate adaptation is not enough: warming does not facilitate success of southern tundra plant populations in the high Arctic. Glob Chang Biol. 23, 1540–1551. doi: 10.1111/gcb.13417

PubMed Abstract | Crossref Full Text | Google Scholar

Boquete, M. T., Muyle, A., Alonso, C. (2021). Plant epigenetics: phenotypic and functional diversity beyond the DNA sequence. Am. J. Bot. 108, 553–558. doi: 10.1002/ajb2.1645

PubMed Abstract | Crossref Full Text | Google Scholar

Catchen, J., Hohenlohe, P. A., Bassham, S., Amores, A., Cresko, W. A. (2013). Stacks: an analysis tool set for population genomics. Mol. Ecol. 22, 3124–3140. doi: 10.1111/mec.12354

PubMed Abstract | Crossref Full Text | Google Scholar

Cavanaugh, K. C., Kellner, J. R., Forde, A. J., Gruner, D. S., Parker, J. D., Rodriguez, W., et al. (2014). Poleward expansion of mangroves is a threshold response to decreased frequency of extreme cold events. Proc. Natl. Acad. Sci. U.S.A. 2, 723–727. doi: 10.1073/pnas.1315800111

PubMed Abstract | Crossref Full Text | Google Scholar

Chen, C., Chen, H., Zhang, Y., Thomas, H. R., Frank, M. H., He, Y., et al. (2020). TBtools: an integrative toolkit developed for interactive analyses of big biological data. Mol. Plant 13, 1194–1202. doi: 10.1016/j.molp.2020.06.009

PubMed Abstract | Crossref Full Text | Google Scholar

Chen, L., Wang, W., Li, Q. Q., Zhang, Y., Yang, S., Osland, M. J., et al. (2017). Mangrove species’ responses to winter air temperature extremes in China. Ecosphere 8. doi: 10.1002/ecs2.1865

Crossref Full Text | Google Scholar

Chen, S., Zhou, Y., Chen, Y., Gu, J. (2018). fastp: an ultra-fast all-in-one FASTQ preprocessor. Bioinformatics 34, i884–i890. doi: 10.1093/bioinformatics/bty560

PubMed Abstract | Crossref Full Text | Google Scholar

Cingolani, P., Platts, A., Wang Le, L., Coon, M., Nguyen, T., Wang, L., et al. (2012). A program for annotating and predicting the effects of single nucleotide polymorphisms, SnpEff: SNPs in the genome of Drosophila melanogaster strain w1118; iso-2; iso-3. Fly (Austin) 6, 80–92. doi: 10.4161/fly.19695

PubMed Abstract | Crossref Full Text | Google Scholar

Collin, A., Daszkowska-Golec, A., Szarejko, I. (2021). Updates on the role of abscisic acid insensitive 5 (abi5) and abscisic acid-responsive element binding factors (ABFs) in ABA signaling in different developmental stages in plants. Cells 10. doi: 10.3390/cells10081996

PubMed Abstract | Crossref Full Text | Google Scholar

Danecek, P., Auton, A., Abecasis, G., Albers, C. A., Banks, E., Depristo, M. A., et al. (2011). The variant call format and VCFtools. Bioinformatics 27, 2156–2158. doi: 10.1093/bioinformatics/btr330

PubMed Abstract | Crossref Full Text | Google Scholar

Delay, C., Imin, N., Djordjevic, M. A. (2013). CEP genes regulate root and shoot development in response to environmental cues and are specific to seed plants. J. Exp. Bot. 64, 5383–5394. doi: 10.1093/jxb/ert332

PubMed Abstract | Crossref Full Text | Google Scholar

Ellison, A. M., Felson, A. J., Friess, D. A. (2020). Mangrove rehabilitation and restoration as experimental adaptive management. Front. Mar. Sci. 7. doi: 10.3389/fmars.2020.00327

Crossref Full Text | Google Scholar

Foyer, C. H., Kranner, I. (2023). Plant adaptation to climate change. Biochem. J. 480, 1865–1869. doi: 10.1042/bcj20220580

PubMed Abstract | Crossref Full Text | Google Scholar

Freed, D., Aldana, R., Weber, J., Edwards, J. (2017). The Sentieon Genomics Tools - A fast and accurate solution to variant calling from next-generation sequence data. doi: 10.1101/115717

Crossref Full Text | Google Scholar

Guo, W. L., Chen, R. G., Gong, Z. H., Yin, Y. X., Ahmed, S. S., He, Y. M. (2012). Exogenous abscisic acid increases antioxidant enzymes and related gene expression in pepper (Capsicum annuum) leaves subjected to chilling stress. Genet. Mol. Res. 11, 4063–4080. doi: 10.4238/2012.September.10.5

PubMed Abstract | Crossref Full Text | Google Scholar

He, S., Wang, X., Du, Z., Liang, P., Zhong, Y., Wang, L., et al. (2023). Physiological and transcriptomic responses to cold waves of the most cold-tolerant mangrove, Kandelia obovata. Front. Plant Sci. 14. doi: 10.3389/fpls.2023.1069055

PubMed Abstract | Crossref Full Text | Google Scholar

Hu, M. J., Sun, W. H., Tsai, W. C., Xiang, S., Lai, X. K., Chen, D. Q., et al. (2020). Chromosome-scale assembly of the Kandelia obovata genome. Hortic. Res. 7, 75. doi: 10.1038/s41438-020-0300-x

PubMed Abstract | Crossref Full Text | Google Scholar

Hu, Y., Jiang, L., Wang, F., Yu, D. (2013). Jasmonate regulates the inducer of cbf expression-C-repeat binding factor/DRE binding factor1 cascade and freezing tolerance in Arabidopsis. Plant Cell 25, 2907–2924. doi: 10.1105/tpc.113.112631

PubMed Abstract | Crossref Full Text | Google Scholar

Hu, Y., Jiang, Y., Han, X., Wang, H., Pan, J., Yu, D. (2017). Jasmonate regulates leaf senescence and tolerance to cold stress: crosstalk with other phytohormones. J. Exp. Bot. 68, 1361–1369. doi: 10.1093/jxb/erx004

PubMed Abstract | Crossref Full Text | Google Scholar

Huang, X., Shi, H., Hu, Z., Liu, A., Amombo, E., Chen, L., et al. (2017). ABA is involved in regulation of cold stress response in Bermudagrass. Front. Plant Sci. 8. doi: 10.3389/fpls.2017.01613

PubMed Abstract | Crossref Full Text | Google Scholar

Jung, Y., Han, D. (2022). BWA-MEME: BWA-MEM emulated with a machine learning approach. Bioinformatics 38, 2404–2413. doi: 10.1093/bioinformatics/btac137

PubMed Abstract | Crossref Full Text | Google Scholar

Korunes, K. L., Samuk, K. (2021). pixy: Unbiased estimation of nucleotide diversity and divergence in the presence of missing data. Mol. Ecol. Resour 21, 1359–1368. doi: 10.1111/1755-0998.13326

PubMed Abstract | Crossref Full Text | Google Scholar

Lescot, M., Déhais, P., Thijs, G., Marchal, K., Moreau, Y., Van De Peer, Y., et al. (2002). PlantCARE, a database of plant cis-acting regulatory elements and a portal to tools for in silico analysis of promoter sequences. Nucleic Acids Res. 30, 325–327. doi: 10.1093/nar/30.1.325

PubMed Abstract | Crossref Full Text | Google Scholar

Li, H., Handsaker, B., Wysoker, A., Fennell, T., Ruan, J., Homer, N., et al. (2009). The sequence alignment/map format and SAMtools. Bioinformatics 25, 2078–2079. doi: 10.1093/bioinformatics/btp352

PubMed Abstract | Crossref Full Text | Google Scholar

Li, C., Wang, F., Yang, P., Wang, F. C., Hu, Y. Z., Zhao, Y. L., et al. (2024). Mangrove wetlands distribution status identification, changing trend analyzation and carbon storage assessment of China. China Geology 7, 1–11. doi: 10.31035/cg2023049

Crossref Full Text | Google Scholar

Liu, X., Lu, X., Yang, S., Liu, Y., Wang, W., Wei, X., et al. (2022). Role of exogenous abscisic acid in freezing tolerance of mangrove Kandelia obovata under natural frost condition at near 32^°N. BMC Plant Biol. 22, 593. doi: 10.1186/s12870-022-03990-2

PubMed Abstract | Crossref Full Text | Google Scholar

Lu, W. X., Zhang, B. H., Yang, S. C. (2023). Survive the north: transplantation for conservation of mangrove forests requires consideration of influences of low temperature, mating system and their joint effects on effective size of the reforested populations. Front. Ecol. Evol. 11. doi: 10.3389/fevo.2023.1160468

Crossref Full Text | Google Scholar

Lu, W. X., Zhang, B. H., Zhang, Y. Y., Yang, S. C. (2021). Differentiation of cold tolerance in an artificial population of a mangrove species, Kandelia obovata, is associated with geographic origins. Front. Plant Sci. 12. doi: 10.3389/fpls.2021.695746

PubMed Abstract | Crossref Full Text | Google Scholar

Marchini, G. L., Maraist, C. A., Cruzan, M. B. (2019). Trait divergence, not plasticity, determines the success of a newly invasive plant. Ann. Bot. 123, 667–679. doi: 10.1093/aob/mcy200

PubMed Abstract | Crossref Full Text | Google Scholar

Mckenna, A., Hanna, M., Banks, E., Sivachenko, A., Cibulskis, K., Kernytsky, A., et al. (2010). The Genome Analysis Toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data. Genome Res. 20, 1297–1303. doi: 10.1101/gr.107524.110

PubMed Abstract | Crossref Full Text | Google Scholar

Meziane, D., Shipley, B. (1999). Interacting determinants of specific leaf area in 22 herbaceous species: effects of irradiance and nutrient availability. Plant Cell Environ. 22, 447–459. doi: 10.1046/j.1365-3040.1999.00423.x

Crossref Full Text | Google Scholar

Minh, B. Q., Schmidt, H. A., Chernomor, O., Schrempf, D., Woodhams, M. D., Von Haeseler, A., et al. (2020). IQ-TREE 2: new models and efficient methods for phylogenetic inference in the genomic era. Mol. Biol. Evol. 37, 1530–1534. doi: 10.1093/molbev/msaa015

PubMed Abstract | Crossref Full Text | Google Scholar

Miryeganeh, M. (2022). Mangrove forests: natural laboratories for studying epigenetic and climate changes. Front. Plant Sci. 13. doi: 10.3389/fpls.2022.851518

PubMed Abstract | Crossref Full Text | Google Scholar

Nakamura, M., Toyota, M., Tasaka, M., Morita, M. T. (2011). An Arabidopsis E3 ligase, SHOOT GRAVITROPISM9, modulates the interaction between statoliths and F-actin in gravity sensing. Plant Cell 23, 1830–1848. doi: 10.1105/tpc.110.079442

PubMed Abstract | Crossref Full Text | Google Scholar

Osland, M. J., Day, R. H., Larriviere, J. C., From, A. S. (2014). Aboveground allometric models for freeze-affected black mangroves (Avicennia germinans): equations for a climate sensitive mangrove-marsh ecotone. PLoS One 9, e99604. doi: 10.1371/journal.pone.0099604

PubMed Abstract | Crossref Full Text | Google Scholar

Pathoumthong, P., Zhang, Z., Roy, S. J., El Habti, A. (2023). Rapid non-destructive method to phenotype stomatal traits. Plant Methods 19, 36. doi: 10.1186/s13007-023-01016-y

PubMed Abstract | Crossref Full Text | Google Scholar

Peakall, R., Smouse, P. E. (2012). GenAlEx 6.5: genetic analysis in Excel. Population genetic software for teaching and research–an update. Bioinformatics 28, 2537–2539. doi: 10.1093/bioinformatics/bts460

PubMed Abstract | Crossref Full Text | Google Scholar

Peng, Y. L., Wang, Y. S., Fei, J., Sun, C. C. (2020). Isolation and expression analysis of two novel C-repeat binding factor (CBF) genes involved in plant growth and abiotic stress response in mangrove Kandelia obovata. Ecotoxicology 29, 718–725. doi: 10.1007/s10646-020-02219-y

PubMed Abstract | Crossref Full Text | Google Scholar

Pickens, C. N., Hester, M. W. (2011). Temperature tolerance of early life history stages of black mangrove avicennia germinans: implications for range expansion. Estuaries Coasts 34, 824–830. doi: 10.1007/s12237-010-9358-2

Crossref Full Text | Google Scholar

Purcell, S., Neale, B., Todd-Brown, K., Thomas, L., Ferreira, M. A., Bender, D., et al. (2007). PLINK: a tool set for whole-genome association and population-based linkage analyses. Am. J. Hum. Genet. 81, 559–575. doi: 10.1086/519795

PubMed Abstract | Crossref Full Text | Google Scholar

Schneider, C. A., Rasband, W. S., Eliceiri, K. W. (2012). NIH Image to ImageJ: 25 years of image analysis. Nat. Methods 9, 671–675. doi: 10.1038/nmeth.2089

PubMed Abstract | Crossref Full Text | Google Scholar

Su, W., Ye, C., Zhang, Y., Hao, S., Li, Q. Q. (2019). Identification of putative key genes for coastal environments and cold adaptation in mangrove Kandelia obovata through transcriptome analysis. Sci. Total Environ. 681, 191–201. doi: 10.1016/j.scitotenv.2019.05.127

PubMed Abstract | Crossref Full Text | Google Scholar

Sun, M. M., Liu, X., Huang, X. J., Yang, J. J., Qin, P. T., Zhou, H., et al. (2022). Genome-wide identification and expression analysis of the NAC gene family in Kandelia obovata, a typical mangrove plant. Curr. Issues Mol. Biol. 44, 5622–5637. doi: 10.3390/cimb44110381

PubMed Abstract | Crossref Full Text | Google Scholar

Waadt, R., Seller, C. A., Hsu, P. K., Takahashi, Y., Munemasa, S., Schroeder, J. I. (2022). Plant hormone regulation of abiotic stress responses. Nat. Rev. Mol. Cell Biol. 23, 680–694. doi: 10.1038/s41580-022-00479-6

PubMed Abstract | Crossref Full Text | Google Scholar

Wang, H., Blakeslee, J. J., Jones, M. L., Chapin, L. J., Dami, I. E. (2020). Exogenous abscisic acid enhances physiological, metabolic, and transcriptional cold acclimation responses in greenhouse-grown grapevines. Plant Sci. 293, 110437. doi: 10.1016/j.plantsci.2020.110437

PubMed Abstract | Crossref Full Text | Google Scholar

Wright, S. (1978). Evolution and the genetic of population, variability within and among natural populations. Chicago: Univ. Chicago Press 4, 213–220.

Google Scholar

Yang, S., Liu, X., Deng, R.-J., Chen, Q.-X., Wang, J.-W., Lu, X. (2020). Geographic variations of hypocotyl and seedling growth traits for Kandelia obovata with different provenances. Chin. J. Ecol. 39, 1769–1777. doi: 10.13292/j.1000-4890.202006.003

Crossref Full Text | Google Scholar

Zhang, J., Fan, T., Cai, X., Ouyang, S., Yang, J., Song, Z., et al. (2024). Comparative transcriptomic and metabolomic analysis of two related Kandelia obovata populations in response to cold wave. Phenomics accepted. doi: 10.1007/s43657-024-00204-7

Crossref Full Text | Google Scholar

Zhao, C. P., Jia, M. M., Zhang, R., Wang, Z. M., Mao, D. H., Zhong, C. R., et al. (2024). Distribution of mangrove species Kandelia obovata in China using time-series sentinel-2 imagery for sustainable mangrove management. J. Remote Sens. 4, 15. doi: 10.34133/remotesensing.0143

Crossref Full Text | Google Scholar

Zhao, Y., Zhong, Y., Ye, C., Liang, P., Pan, X., Zhang, Y. Y., et al. (2021). Multi-omics analyses on Kandelia obovata reveal its response to transplanting and genetic differentiation among populations. BMC Plant Biol. 21, 341. doi: 10.1186/s12870-021-03123-1

PubMed Abstract | Crossref Full Text | Google Scholar

Zhou, H., Hua, J., Zhang, J., Luo, S. (2022). Negative interactions balance growth and defense in plants confronted with herbivores or pathogens. J. Agric. Food Chem. 70, 12723–12732. doi: 10.1021/acs.jafc.2c04218

PubMed Abstract | Crossref Full Text | Google Scholar

Zhu, H., Lin, H. J., Le, Y., Li, H. P., Yue, C. L., Jiang, B. (2021). Geographical distribution pattern and environmental explanation of Kandelia obovata Sheue, H.Y. Liu & J. Yong populations along the Southeast coast of China. Plant Sci. J. 39, 476–487. doi: 10.11913/PSJ.2095-0837.2021.50476

Crossref Full Text | Google Scholar