Genetic basis of local adaptation in the cold-tolerant mangrove Kandelia obovata

Abstract

Understanding the genetic basis of local adaption is crucial in the context of global climate change. Mangroves, as salt-tolerant trees and shrubs in the intertidal zone of tropical and subtropical coastlines, are particularly vulnerable to climate change. Kandelia obovata, the most cold-tolerant mangrove species, has undergone ecological speciation from its cold-intolerant counterpart, Kandelia candel, with geographic separation by the South China Sea. In this study, we conducted whole-genome re-sequencing of K. obovata populations along China’s southeast coast, to elucidate the genetic basis responsible for mangrove local adaptation to climate. Our analysis revealed a strong population structure among the three K. obovata populations, with complex demographic histories involving population expansion, bottleneck, and gene flow. Genome-wide scans unveiled pronounced patterns of selective sweeps in highly differentiated regions among pairwise populations, with stronger signatures observed in the northern populations compared to the southern population. Additionally, significant genotype-environment associations for temperature-related variables were identified, while no associations were detected for precipitation. A set of 39 high-confidence candidate genes underlying local adaptation of K. obovata were identified, which are distinct from genes under selection detected by comparison between K. obovata and its cold-intolerant relative K. candel. These results significantly contribute to our understanding of the genetic underpinnings of local adaptation in K. obovata and provide valuable insights into the evolutionary processes shaping the genetic diversity of mangrove populations in response to climate change.

1 Introduction

Understanding the evolutionary processes and genetic basis of local adaptation has been a longstanding interest in evolutionary biology (Hoban et al., 2016; López-Goldar & Agrawal, 2021). Abiotic and biotic effects such as low and high temperature, drought, flooding, herbivore and pathogen stresses impose different selective pressures across habitats (Vanwallendael et al., 2019). The interactions among selection, gene flow and genetic drift shape genetic variation within and between populations, leading to evolutionary differentiation at different spatial scales (Orsini et al., 2013; von Takach, Penton, et al., 2021). Local adaptation can occur when genetic differentiation allows a single population to become better adapted in a particular set of environmental conditions within its range (Angert et al., 2020). The genomic architecture of natural populations thus provides an opportunity to reveal and link molecular mechanisms underlying species’ phenotypic diversity with diverse environments in which the species lives. This knowledge is important for informing management decisions in the context of rapid contemporary environmental changes (Isabel et al., 2020).

Climate adaptation is one of the most prevalent forms of local adaptation that has significant impact on the distribution of a species. Classical experimental approaches explore to which extent morphological, physiological and transcriptional responses are associated with changes in environmental factors, typically through transplants and common gardens (Carley et al., 2021; Oomen & Hutchings, 2022). Nevertheless, as most traits are polygenic, genes involved in ecological, physiological, or transcriptional changes may only be weakly related to fitness. With recent advance of sequencing technology, population genetic or genomic analyses are typically utilized to identify selective signatures linked to local adaptation, with allele frequencies expected to demonstrate ecological differentiation or correlation with specific environmental factors (von Takach, et al., 2021; Hu et al., 2022; Wang J. et al., 2022). As adaptation to climate change progresses, genetically distinct populations undergoing local adaptation may eventually become separate species if they no longer interbreed or exchange genes. Therefore, studying the genetic basis of climate adaptation can also provide insights into the evolutionary processes leading to the fixation of genetic differentiation and subsequent speciation (Levin, 2019; Dool et al., 2022).

Mangroves are a group of phylogenetically diverse woody plants that inhabit tropical and subtropical intertidal zones with great ecological and economical significance (Tomlinson, 1986; Sribianti et al., 2021). At the dynamic interface between sea and land, the structure and biodiversity of mangrove communities are sensitive to climate change (Jennerjahn et al., 2017). Mangrove forests have faced and survived several catastrophic climate change events since their origination during the Late Cretaceous-Early Tertiary period (Singh et al., 2022). Sea-level rise significantly influences the historical distribution of mangrove forests, leading to speciation in five common mangrove species via multiple cycles of mixing, isolation, and mixing (Lovelock et al., 2015; He et al., 2022). Meanwhile, low temperature stress governs the latitudinal range limits for contemporary mangrove flora (Lin, 1999) and plays a crucial role in mangrove afforestation and restoration (Ellison, 2000; Quisthoudt et al., 2012; Su et al., 2021; Chen et al., 2023). Given the ecological importance of mangroves and the potential impact of climate change on their habitats, further research into their adaptation to low temperatures is vital. The complex demographic histories of mangroves involving bottleneck, genetic drift, population expansion, and gene flow also provide a challenge for studying local adaptation.

Kandelia obovata, regarded as the most cold-tolerant mangrove species, is able to survive chilling temperatures as low as 4.2°C (Maxwell, 1995; Sheue et al., 2003; Zhou et al., 2007). This species can be found from the Gulf of Tonkin in the northeast to southern Japan, separated from its cold-intolerant relative, K. candel, by the South China Sea (Sheue et al., 2003), representing a good case in point for ecological speciation. Besides the clear differences in cold tolerance between Kandelia species, common garden studies have identified a greater cold tolerance in the northern populations of K. obovata compared to the southern populations (Sheue et al., 2003; Zhao et al., 2021), indicating the occurrence of local adaptation within K. obovata. Whole-genome bisulfite sequencing and RNA-seq of the K. obovata transplants suggest that modifications of DNA methylation in MADS-box genes may contribute to the adaptation to new environments, whereas the suppressed expression of lignin biosynthesis genes appears to play a role in maladaptation (Zhao et al., 2021). Additionally, the physiological and expressional analyses have highlighted several key genes and pathways that are potentially involved in cold tolerance in K. obovata, such as genes involved in calcium signaling, cell wall modification, and post-translational modifications of ubiquitination pathways (He et al., 2023). However, how genetic differentiation is maintained and whether the same genes underlying local adaptation to cold stress within species can also be responsible for between-species differences in cold tolerance remain largely unknown.

In this study, we investigated the genomic architecture of K. obovata populations along a latitudinal gradient in China using whole-genome resequencing, based on the chromosome-anchored genome assembly of K. obovata (Hu et al., 2020). Our analysis aimed to infer the demographic history and detect selective signatures within K. obovata populations, with the underlying hypothesis that local adaptation has occurred despite gene flow, driven by strong selection pressures related to temperature. We expected the signatures of local adaptation to be more pronounced in the northern population compared to the southern population. Moreover, we identified candidate genes under selection by comparing polymorphism within K. obovata with divergence between K. obovata and K. candel. We expected that the genes involved in local adaptation within K. obovata would differ from those under selection in the lineage of K. obovata since its divergence from K. candel.

2 Materials and Methods

2.1 Plant materials and population resequencing

We sampled leaves from a total of 46 K. obovata individuals, collecting plants at least 20 meters apart within each of three natural populations along the southeast coast of China (Supplementary Table S1). Specifically, we collected 14 individuals from Shacheng Harbor (Fuding, Fujian Province), 17 from Yanzao Village (Shenzhen, Guangdong Province), and 15 from Bamen Bay (Wenchang, Hainan Province). The fresh leaves were air-dried with silica gel before DNA extraction. Genomic DNA was then extracted using a modified CTAB protocol (Yang et al., 2008). Subsequently, DNA libraries were obtained using Nextera Mate Pair Sample Preparation Kit (Illumina USA), followed by whole-genome re-sequencing on an Illumina HiSeq 2500 platform (Illumina, San Diego, CA, USA) with paired-end reads of 150 bp (PE150). The average sequencing depth for each individual was approximately 40× (Supplementary Table S1).

2.2 Reads mapping, variant callings and SNP filtering

A high-quality chromosome-scaled assembly of the K. obovata genome, with a total length of 177 Mb, was published using PacBio, Illumina and Hi-C sequencing (Hu et al., 2020) and served as the reference genome for this study. Raw Illumina reads were initially processed by trimming and filtering using Trimmomatic v0.39 (Bolger et al., 2014), followed by mapping onto the K. obovata reference genome using Burrows-Wheeler-Alignment (BWA) v0.7.12-r1039 (Li & Durbin, 2009). Subsequently, sorted bam files were generated from sam files using SAMtools v1.6 (Li et al., 2009), and PCR duplicates were removed using MarkDuplicates in the Picard toolkit. Variants were called using HaplotypeCaller and genotyped using GenotypeGVCFs in Genome Analysis Tool Kit (GATK) (McKenna et al., 2010). The analysis solely focused on single nucleotide polymorphism sties (SNPs) and employed specific filtering criteria to reduce false positives: (1) SNPs with a read number less than three for each individual were removed; (2) SNPs with a minor allele frequency (MAF) less than 0.05 were discarded; (3) SNPs were further filtered using VariantFiltration with the following parameters: quality by depth (QD) < 2.0, Fisher strand (FS) > 60.0, mapping quality (MQ) < 40.0, mapping quality tank sum test (MQRankSum) > -12.5, and read pos rank sum test (ReadPosRankSum) < -8.0. Finally, SNP annotation was conducted based on the K. obovata genome using snpEff 4.3t (Cingolani et al., 2012). Following sequence alignment, removal of PCR duplicates, and SNP filtering, a set of high-quality SNPs sites were retained for the subsequent analyses.

2.3 Genetic diversity and population differentiation

The entire K. obovata genome was divided into non-overlapping 20-kb bins. Nucleotide diversity (θ_π; Tajima, 1989) for each K. obovata population and fixation index (F_ST; Weir & Cockerham, 1984) between pairwise populations were calculated within each bin using VCFtools v0.1.15 (Danecek et al., 2011). Pairwise genetic distance (D_XY; Nei & Miller, 1990) between populations were calculated within each bin using pixy v1.2.7 (Korunes & Samuk, 2021).

2.4 Population structure analysis

A principal component analysis (PCA) was performed to visualize inter-individual genetic relationships using PLINK v1.90 (Purcell et al., 2007). A phylogenetic tree based on the high-quality SNPs was constructed using IQtree (Nguyen et al., 2015) with the parameter: -m GTR+F+G4+ASC, for which the input PHYLIP file was converted using vcf2phylip.py (https://doi.org/10.5281/zenodo.1257058). The number of genetic clusters (K) was identified using ADMIXTURE v1.3.0 with default parameters (Alexander et al., 2009). Various values of K ranging from 2 to 5 were tested, and the best K was selected based on the minimum error rate of K value.

2.5 Demographic modelling

The demographic history of K. obovata was inferred based on its observed population structure through the construction of a two-dimensional joint unfolded site frequency spectrum (2D-SFS) using the easySFS tools (https://github.com/isaacovercast/easySFS#easysfs), with all SNPs and the projection number equal to the individual number (i.e. 14, 17 and 15) for each of the three populations, respectively. We considered various scenarios of divergence, bottleneck, expansion, and/or migration, which were represented by different models: (i) Models 1-4: Divergence, including one-step isolation in model 1 and three two-step isolations of the K. obovata populations in models 2-4; (ii) Model 5-7: Bottleneck, with one population experiencing the reduction of population size and then recovery; (iii) Model 8-10: Expansion, with one population experiencing exponential population changes; (iv) Model 11-14: Migration, with asymmetric gene flow and differences in occurrence between different populations; (v) Model 15: Composite model incorporating each of the estimated best scenario for divergence (Model 1), bottleneck (Model 6), expansion (Model 10), and migration (Model 14) (Supplementary Figure 1).

These demographic models were compared, and demographic parameters were inferred using a coalescent simulation-based method as implemented in fastsimcoal2.6 (Excoffier et al., 2021). The initial ranges for the parameter estimation were listed in Supplementary Table S2. The log-likelihood for a set of demographic parameters was estimated using 100,000 coalescent simulations, with 80 conditional maximum algorithm cycles in each run and global maximum likelihood estimates obtained from 100 independent runs. The maximum likelihood value of the 100 independent runs for each model was used to compare between models using the Akaike information criterion (AIC) and Akaike’s weight of evidence tests. The model with the highest Akaike’s weight value was considered as the optimal model. The parameter confidence intervals (CIs) for the optimal model were obtained from 100 parametric bootstrap samples, independently run 100 times in each bootstrap. When converting estimates to years, it was assumed that the mutation rate and the average generation interval time in K. obovata were 7.86 × 10⁻⁸ per site per generation and 20 years per generation (average time from seed germination to seed production), respectively (He et al., 2022).

2.6 Detection of positive selection

Sliding window analysis was employed to detect outlier genomic regions with strong population differentiation using pairwise fixation index (F_ST) and the cross-population composite likelihood ratio (XP-CLR) test (Chen et al., 2010) with a 20-kb window size and 5-kb step size. Pairwise F_ST estimates were calculated using VCFtools v0.1.15. The XP-CLR test was performed using the xpclr python module (https://github.com/hardingnj/xpclr) with default parameters. Outlier windows with both F_ST and normalized XP-CLR values at least 1.96 standard deviations (SD) above the mean (one-tailed p-value = 0.025) were identified as the highly differentiated regions (HDRs). Sliding window analysis of Fay and Wu’s H-statistic (Fay & Wu, 2000) in HDRs was conducted for each population separately, in comparison with the whole genome, using ANGSD v0.921 (Korneliussen et al., 2014) with a 20-kb window size and 5-kb step size.

To detect positively selected genes, the McDonald-Kreitman (MK) test (McDonald & Kreitman, 1991) was conducted using whole-genome resequencing data for all K. obovata individuals obtained in this study and the reference genome of K. candel (He et al., 2022) as an outgroup. BlastN (Altschul et al., 1990) was used to identify homologous genes between the two Kandelia species using an e-value threshold of < 1e-10, resulting in the identification of 16,536 one-to-one orthologous genes. The number of non-synonymous (D_n) and synonymous (D_s) substitutions were compared to the number of non-synonymous (P_n) and synonymous (P_s) polymorphisms within K. obovata for coding sequences of each Kandelia homologous gene using K. candel as an outgroup. Genes with p-value lower than 0.05 in the one side Fisher’s exact test with Benjamini-Hochberg multiple test (Benjamini & Hochberg, 1995) correction were classified as positively selected genes.

2.7 Genotype-environment associations

Redundancy analysis (RDA) was used to identify associations between SNP variations across diverse populations and environmental parameters (Forester et al., 2018). Climate data including the mean annual temperature, mean annual minimum temperature, and mean annual precipitation for the three sampling locations were obtained from the Central Meteorological Observatory (http://www.nmc.cn/). The RDA was conducted using the RDA function from the vegan package (Oksanen et al., 2007) with the input file of genotype matrix comprising all SNPs as transformed by PLINK v1.90. Genotypes at each SNP site was encoded as follows: 0 for homozygotes identical to the reference, 1 for heterozygotes, and 2 for homozygotes differing from the reference. The proportion of variance explained by the environmental variables was evaluated using the RsquareAdj function. The significance of the linear relationship between each constrained axis and the environmental variables was then assessed using the anova.cca function. In this study, two constrained axes (RDA1 and RDA2) were found to be significant. Outlier SNPs for each significant constrained axis were identified using a cutoff of 1.96 SD greater or less than mean (two-tailed p-value = 0.05) and were annotated with known genes based on the K. obovata genome.

2.8 K. obovata gene functional annotation and GO enrichment analysis

To annotate the K. obovata genes, BlastP (Altschul et al., 1990) was used to identify homologous genes between the K. obovata and Arabidopsis thaliana using an e-value threshold of < 1e-10. Gene Ontology (GO) functions of each K. obovata gene were annotated using eggNOG-mapper (Cantalapiedra et al., 2021). To evaluate potential over-representation of functional gene classes, we conducted GO enrichment analysis using the R package ClusterProfiler (Yu et al., 2012) by applying the annotation mentioned above. The statistical significance of over-represented GO terms within the input gene sets was assessed through Fisher’s exact tests, with a significance threshold set at False Discovery Rate (FDR) < 0.2.

3 Results

3.1 Genomic data and genetic diversity

A total of 0.82 Tb of data with an average sequencing depth of approximately 40× per individual were obtained for the 46 K. obovata individuals from the three populations, representing the natural distribution of K. obovata from north to south in China (Table 1; Figure 1A). After quality control and filtering, the clean reads were mapped to the reference genome of K. obovata, achieving an average mapping rate of 93.33%. The average depth of uniquely mapped reads per site was 40.23 (Supplementary Table S1). In total, we obtained 189,909 high-confidence single nucleotide polymorphism (SNP) sites, corresponding to an average density of 1.07 SNPs per kilobase in the K. obovata genome. Out of these SNPs, 25,988 (13.7%) were in protein-coding regions, including exons and introns, 27,551 (14.5%) were found in putative regulatory regions, including promoters (2kb upstream of TSS) and UTRs, and the majority (56.7%) were intergenic (Supplementary Figure 2).

Figure 1

When comparing between populations, the number of high-confidence SNP sites was highest for the Shenzhen population (160,726), followed by the Wenchang (140,250) and Fuding (115,239) population. The number of private or population-specific SNPs were 8,736, 32,042, and 11,683 for the Fuding, Shenzhen and Wenchang populations, respectively (Figure 1B). The site frequency spectrum of individual K. obovata populations revealed that the Fuding populations had the highest proportion (16.7%) of fixed SNPs followed by the Shenzhen population (5.5%), while only 75 SNPs (0.6%) were fixed in the Wenchang population (Figure 1C). However, the Wenchang population had a higher percentage (9.4%) of rare SNPs (frequency < 0.05), than the Fuding (5.8%) and Shenzhen population (4.4%, Figure 1C). Meanwhile, all three populations possessed considerable medium-frequency SNPs, indicating a high level of heterozygosity within individuals of K. obovata (Figure 1C).

We divided the K. obovata genome into 20-kb non-overlapping bins and used Tajima’s π (θ_π) as a measurement of nucleotide diversity at the individual SNP level. The average nucleotide diversity per site was calculated for each bin and compared between populations of K. obovata. Nucleotide diversity was unevenly distributed across the K. obovata genome, with the most pronounced diversity observed in chromosomes 11, 13, 14 and 15 (Figure 1D). Consistent with the observed number of high-quality SNPs, the Shenzhen population exhibited the highest nucleotide diversity averaged across the genome (mean ± SD, θ_π = 3.27 × 10^-4 ± 9.63 × 10^-4), followed by the Wenchang population (θ_π = 2.87 × 10^-4 ± 8.64 × 10^-4), while the Fuding population showed the lowest nucleotide diversity (θ_π = 2.32 × 10^-4 ± 8.87 × 10^-4) (Figure 1E). The differences in nucleotide diversity between pairwise populations were all significant (Mann-Whitney U test, all p-value < 0.05; Figure 1E).

3.2 Genetic differentiation and population structure

To assessed population differentiation, the fixation index F_ST (Weir and Cockerham, 1984) and the average nucleotide diversity D_XY (Nei & Miller, 1990) were calculated for each SNP and averaged within each bin for pairwise populations of K. obovata (Figure 2). F_ST values revealed substantial population differentiation, particularly evident when comparing the Fuding population with the others: Fuding vs. Shenzhen (F_ST = 0.30 ± 0.19) and Fuding vs. Wenchang (F_ST = 0.30 ± 0.16), in contrast to Shenzhen vs. Wenchang (F_ST = 0.24 ± 0.13) (Table 1 and Figure 2A). D_XY values were similar among pairs of populations (D_XY = 0.30 ± 0.15 to 0.33 ± 0.11), all supporting a high level of population differentiation in K. obovata (Figure 2B). The principal component analysis (PCA) revealed a clear clustering pattern, indicating that individuals within a population tend to group together and are distinctly separated from individuals in other populations (Figure 2C). Similarly, the structure analysis indicated that the pattern of ancestry was best represented by K = 3, which is supported by the lowest minimum K-value error rate (cross-validation error = 0.41) (Figure 2D). The three clusters that resulted from K = 3 clearly separated individuals according to their population origination. When K = 2, individuals from Fuding population formed a cohesive cluster, while individuals from Shenzhen and Wenchang clustered together. Subsequent analyses at K values of 4 and 5 revealed highly variable levels of admixture within the Shenzhen and Wenchang populations. However, there is no significant admixture within the genetic cluster in Fuding population, indicating a relatively more isolated genetic profile for this population.

Figure 2

3.3 Inference of demographic history

The optimal demographic model (model 15; log-likelihood = -1089415.63, AIC = 2178865.27, ΔAIC = 0), as shown in Figure 3 and Supplementary Table 3, supports a composite history inclusive of divergence, bottleneck, expansion, and migration of K. obovata. With a generation time of 20 years and mutation rate of 7.86 × 10⁻⁸ per site per generation (He et al., 2022), the optimal demographic model suggests that the three K. obovata populations diverged from their ancestral population with an effective population size (N_e) of 57,833 (95% CI: [33,569, 120,531]) approximately 93,080 years ago (95% CI: [43,800, 24,175,600]). This coincides the late onset of the last glacial period, marked by highly unstable sea levels. Subsequently, the Fuding population underwent a gradual population expansion with a growth rate of 1.05e-4 (95% CI: [1.58e-7, 4.41e-4]), increasing from an N_e of 4,319 (95% CI: [265, 19,075]) to 6,142 (95% CI: [874, 40,402]). In contrast, the ancestral Shenzhen population, with an N_e of 17,072 (95% CI: [265, 19,075]), experienced a ten-thousand-year bottleneck around 45,000 years ago (95% CI: [15,320, 269,360]), coinciding with the start of the last glacial maximum. The Wenchang population maintained a consistently small value of N_e of 778 (95% CI: [564, 5,884]) over the last hundred thousand years with minimal fluctuation. Migration rates per generation between populations varied substantially, with gene flow mainly occurring from Fuding to Wenchang (m_FW) at a rate of 6.35e-4 (95% CI: [3.88e-5, 1.03e-3] and from Shenzhen to Wenchang (m_SW = 2.62e-4, 95% CI: [3.41e-5, 7.70e-4]) or Fuding (m_SF = 1.48e-4, 95% CI: [1.63e-5, 2.47e-4]), while the others were relatively small (Figure 3).

Figure 3

Notably, both the Shenzhen and Fuding populations have experienced changes in population size according to the best model. Despite high genetic diversity, a model that assumed a bottleneck occurred in the Shenzhen population (model 6; log-likelihood = -1091480.43, AIC = 2182978.86, ΔAIC = 4113.594) was more likely than one assuming a bottleneck in the Fuding population (model 7; log-likelihood = -1091661.914, AIC = 2183341.828, ΔAIC = 4476.562) (Supplementary Figure 1 and Supplementary Table 3). In contrast, although the Fuding population is marginal, a model that assumed expansion in the Fuding population (model 10; log-likelihood = -1091358.265, AIC = 2182730.53, ΔAIC = 3865.26) was more likely than in the other two populations (Supplementary Figure 1 and Supplementary Table 3).

3.4 Detection of positive selection

Sliding widow analysis of population differentiation used a cutoff of at least 1.96 SD above the mean (p < 0.025) to detect highly differentiated outliers for both F_ST and XP-CLR values. For the F_ST values, 1452, 1457, and 1494 outliers were detected for three pairwise comparisons: Fuding and Wenchang (F-W), Fuding and Shenzhen (F-S), and Shenzhen and Wenchang (S-W), respectively (Figures 4A–C). For XP-CLR values, 125, 264, and 430 outliers were detected for the same pairwise comparisons (Figures 4A–C). The regions where both methods detected outliers were considered as the highly differentiated regions (HDRs) or the candidate regions for selection. Sliding window analysis of Fay and Wu’s H-statistic (Fay & Wu, 2000), designed to detect high-frequency hitchhiking alleles associated with selective sweeps, revealed significantly lower H values in HDRs than in the whole genome for each of the three populations (Mann-Witney U test, all p-value < 0.001; Supplementary Figure 3), suggesting selective sweeps occurred in HDRs of each population. Negative H values were observed in HDRs of the Fuding (mean ± SD, H-statisic= -2.43 ± 3.22) and Shenzhen populations (H-statistic = -2.43 ± 2.12), but not in the Wenchang population (H-statistic= 0.13 ± 1.03; Supplementary Figure 3), indicating that selective signatures of HDRs are stronger in the two northern populations compared to the Wenchang population.

Figure 4

A total of 4, 27, and 95 HDRs were identified for the F-W, F-S, and S-W comparisons, respectively, comprising 4, 8, and 67 genes and a total of 2,868 SNPs (Figure 4D; Supplementary Table 4). Among them, 77 genes were identified in a single population pair, while one gene, Plant U-box 45 (PUB45; Maker00002595), was shared by the F-W and S-W comparisons. Notably, no genes were shared by all three comparisons, resulting in a total of 78 genes in HDRs (Figure 4D). Furthermore, Gene Ontology (GO) analysis revealed that the identified genes are enriched in functional categories of supramolecular fiber organization (GO: 0097435) and those related to response to stimuli, including cellular response to extracellular stimulus (GO: 0031668), cellular response to external stimulus (GO: 0071496), response to extracellular stimulus (GO: 0009991), and cellular response to stress (GO: 0033554). Additionally, other enriched functional categories are involved in cell cycles and the regulation of developmental process (Figure 4E).

Using K. candel as an outgroup, we also detected genes under positive selection by the McDonald-Kreitman (MK) test (McDonald & Kreitman, 1991). Applying the G test of independence described in the original MK test, we identified 27 genes with a positive selection signature at a cutoff of an adjusted p-value < 0.05 (Supplementary Table S5). Only two genes, AP2-like ethylene-responsive transcription factor 1 (AIL1; Maker00004062) participating in the ethylene-activated signaling pathway (Kim et al., 2006) and callose synthase 5 (CALS5; Maker00007136) involved in the regulation of pollen tube growth (Dong et al., 2005) were also found in the list of HDRs-related selected genes (Figure 4F). The 27 candidate genes were enriched in housekeeping functional categories including macromolecule modification (GO: 0043412), regulation of cellular process (GO: 0050794), protein modification process (GO: 0036211), and regulation of nitrogen compound metabolic process (GO: 0034641) with a significance threshold of p-value < 0.05 (Figure 4G). However, none of these overrepresentations retained significant after FDR correction, likely due to the limited size of gene list.

3.5 Genotype-environment associations

The three populations, located at different latitudes, differ in their local environments. Redundancy analysis (RDA) detected significant associations between SNPs and two environmental variables (F-statistics, both p-value < 0.001), the mean annual temperature (MAT) and the mean annual minimum temperature (MAMT), but not between SNPs and the mean annual precipitation (MAP, Supplementary Figure 4A). This suggests that temperature is the primary environmental variable shaping population differentiation. Two constrained axes (RDA1 and RDA2) explained about 29% of the total variation. Using a cutoff of 1.96 SD greater or less than mean (two-tailed p-value = 0.05), a total of 16,256 environment-associated SNPs were identified as candidate SNPs involved in local adaptation (Figure 5; Supplementary Figure 4B). The temperature-related SNPs exhibited distinct patterns of allele frequency spectrum compared to genome-wide SNPs in all three populations (Kolmogorov-Smirnov test, p-value < 2.2 × 10^-16 for all comparisons), characterized by an elevated proportion of fixed SNPs in the Fuding and Shenzhen populations (Figure 5C). Among the 16,256 temperature-associated SNPs, 1,118 were located in the HDRs, indicating a significant enrichment of temperature-associated SNPs in HDRs compared to the whole genome (χ² test, p-value < 2.2 × 10^-16; Figure 5D). These temperature-related SNPs are within 1,074 genes which were enriched in a range of functional categories involved in sexual reproduction, DNA repair, cell death, cell cycle and immune response (Supplementary Table S6).

Figure 5

A total of 39 genes containing temperature-associated SNPs within their exons and located within HDRs were finally identified (Figure 5E). These genes, consisting of 2, 6, and 32 for the F-W, F-S, and S-W comparisons, respectively, were considered as high-confidence candidates underlying local adaptation to temperature in K. obovata (Table 2). In the F-W comparison, two genes were identified, namely Plant U-box 45 (PUB45; Maker00002595) and a gene (Maker00006413) containing eight temperature-associated SNPs but no homolog in Arabidopsis (Table 2). Plant U-box 45, known for its function in protein ubiquitination and cold response (He et al., 2023), was also identified as a high-confidence candidate in the S-W comparison. Notably, the high-confidence candidates in the F-S comparison with known functions included CRY2 interacting splicing factor 1 (CIS1; Maker00002948) involved in the regulation of flowering (Zhao et al., 2022), tRNA pseudouridine synthase A 5 (TRUA5, Maker00014866) involved in RNA modification, Pol-like 5 (PLL5; Maker00017856) participating in leaf development (Song & Clark, 2005), a 2-oxoglutarate (2OG) and Fe(II)-dependent oxygenase superfamily protein (Maker00003467), and a UDP-N-acetylglucosamine (UAA) transporter family gene (Maker00018193). Meanwhile, for the S-W comparison, candidate genes were found to be responsive to cold and other stresses, such as growth-regulating factor 5 (GRF5; Maker00004032; Lantzouni et al., 2020), ATP-binding cassette G40 (ABCG40; Maker00008456; Baron et al., 2012) and Integrin-linked kinase 1 (ILK1; Maker00007966; Brauer et al., 2016), or participate in various development processes, such as cryptochrome 1 (CRY1; Maker00010750) involved in flowering (Gao et al., 2023), cotton Golgi-related 2 (CGR2; Maker00005232) involved in leaf morphogenesis (Kim et al., 2015) and lonesome highway (LHW; Maker00018105) involved in root development (Ohashi-Ito & Bergmann, 2007) (Table 2).

4 Discussion

K. obovata offers valuable opportunities for investigating the genetic basis of local adaptation in species that have undergone ecological speciation with complex demographic histories. Here, by examining how genome-wide patterns of population structure relate to genotype-environment associations, in conjunction with inferring demographic history and detecting selection signatures, we empirically contribute to the broader understanding of how genetic variation leads to climate adaptation. Despite complex demographic history, our results demonstrate that natural selection drives local adaptation in K. obovata due to differential selective pressures related to temperature.

4.1 Complex demography and restricted gene flow

The complex demographic history of the K. obovata populations has a significant impact on their genetic diversity and population structure. The Fuding population, occurring in the northmost in China, exhibits the lowest genetic diversity and is more isolated from the other two populations. This is consistent with the expectation for marginal populations, which often show lower genetic diversity and increased genetic differentiation from more central populations (Eckert et al., 2008). However, the inference of the demographic model revealed that the Fuding population has experienced a mild population expansion, while the central population in Shenzhen underwent a population bottleneck around the start of the last glacial maximum (Figure 3). These divergent patterns of population size change may align with the cold-tolerant nature of K. obovata, with the northern population being better adapted to cold weather than the southern population (Lu et al., 2022). Additionally, these two populations exhibit distinct patterns in the site frequency spectrum compared to the Wenchang population, with the former two showing a high number of fixed SNPs while the Wenchang population has almost none (Figure 1C). These results suggest that genetic drift associated with frequent population size change in the Fuding and Shenzhen populations may have caused rapid fixation of alleles in these regions, while the southern population in Wenchang may represent the ancestor state of this species, considering that K. obovata is diverged from the cold intolerant Kandelia candel (Sheue et al., 2003). Nevertheless, it is unclear how the Wenchang population could maintain a constant population size, especially given that global cooling during the last glacial maximum is known to be a major cause of extinction for tropical woody species such as mangroves (Song et al., 2021).

It is widely acknowledged that mangroves exhibit long-distance gene flow due to the dispersal of floating propagules (Van der Stocken et al., 2022). Previous studies on K. candel have indicated substantial regional gene flow over considerable distances (Chiang et al., 2001; Geng et al., 2008). In this study, the estimated gene flow from Fuding to Wenchang (Nm_FW = 3.90) is much greater than that from Wenchang to Fuding (Nm_WF = 0.05), despite the broad confidence intervals (Figure 3). This finding is consistent with the direction of current flow in the South China Sea during early winter (Shaw & Chao, 1994), coinciding with the maturation and shedding of the hypocotyl in K. obovata in the northern regions (Khoon et al., 2004). Thus, the asymmetric migration between the Fuding and Wenchang populations provides additional evidence for the influence of ocean currents on gene flow in mangroves. However, we observed strong population structure among the three K. obovata populations in China (Figure 2D), indicating local gene flow may be insufficient to counteract genetic differentiation caused by either genetic drift or natural selection. Notably, the best-fitting demographic model (Figure 3) reveals substantial gene flow occurring only from the Fuding and Shenzhen populations to the Wenchang population (Nm_FW = 3.90 and Nm_SW = 1.13), while all other Nmestimates, ranging from 0.02 to 0.64, are relatively small. Contrastingly, a recent study on grey mangrove (Avicennia marina) populations across the Red Sea, the Arabian Sea and the Persian/Arabian Gulf revealed a moderate population genetic structure correlating with geographic distance, which supports clades both among and within the seas surrounding the Arabian Peninsula (Friis et al., 2024).

What might have caused the limited gene flow among geographically proximate K. obovata populations? One possible explanation is the impact of human activities, which have led to increased fragmentation. Human activities such as land development, urbanization, and habitat destruction can create barriers to gene flow, resulting in population isolation and reduced genetic exchange between populations (Guo et al., 2016). Moreover, the small effective population size for all populations in this study may have also played a role in limiting gene flow. With a small effective population size, genetic drift has a more pronounced effect, leading to increased genetic differentiation between populations (Song et al., 2006). This effect may be particularly significant for the Shenzhen and Fuding populations, which have experienced population bottleneck and expansion, respectively.

4.2 Local adaptation and genotype-environment associations

Local adaptation involves significant changes in allele frequency. Therefore, conducting a genome scan to identify outlier values of the allelic differentiation is the preferred method for detecting loci associated with local adaptation. It has been recognized that high neutral differentiation among populations can make it more difficult to detect high outlier loci (Pérez-Figueroa et al., 2010). This is evident in the case for K. obovata, where the average F_ST values are approximately 0.3 across the genome (Figure 2B). Utilizing the 1.96 SD cutoff of the empirical distribution, we discovered that the number of F_ST outliers was largely comparable for all three pairwise population comparisons, whereas the number of XP-CLR outliers varied substantially among comparisons (Figure 5A-C). Given that F_ST is the quotient of two variances, the large expected variability in F_ST among neutral loci, influenced by complex demography, population structure, and migration, may impact the power to discern high outlier F_ST values, potentially resulting in a high rate of false positive or false negative (Le Corre & Kremer, 2012). In contrast, the XP-CLR method, which is independent of window size and robust to uncertainty regarding demography history, identifies regions in the genome where the change in allele frequency at the locus occurred too quickly to be due to random drift (Chen et al., 2010). The different results between the two methods lies in the fact that the detection of F_ST outliers depends on neutral loci, while the XP-CLR method is sensitive to recent selective sweeps. By combining both methods, we observed an enrichment of environment-associated SNPs in highly differentiated regions (HDRs) that showed both outlier F_ST and outlier XP-CLR values, indicating the effectiveness of our analyses in identifying loci associated with local adaptation in the presence of high neutral differentiation among populations.

The three K. obovata populations, located at varying latitudes, are subject to differing intensities of selective pressure. Genotype-environment associations in K. obovata were detected exclusively for temperature-related variables, indicating that temperature is the primary factor driving differential selective pressures leading to local adaption. This supports the notion that winter temperature shapes mangrove distributions and assemblage composition in China (Wu et al., 2018). Local temperature has been consistently identified as significant drivers of varying selective pressures in many organisms, such as grey mangrove (Friis et al., 2024), lichen-forming fungi (Valim et al., 2023), and mountain pine (Méndez-Cea et al., 2023). In this study, the Fuding and Shenzhen populations showed a higher proportion of fixation in temperature-associated SNPs (Figure 5C) compared with genome-wide SNPs (Figure 1C), implying that natural selection has facilitated the rapid evolution of genes related to local adaptation. However, the evidence for local adaptation in the Wenchang population is weak, as reflected by the allele frequency spectra (Figure 5C) and results of Fay and Wu’s H (Supplementary Figure 3), likely due to the minimal selective pressure for temperature adaptation in this tropical region.

Gene flow may have also impacted local adaptation in the Wenchang population, which receives migrates from the Fuding and Shenzhen populations (Nm_FW = 3.90 and Nm_SW = 1.13; Figure 3). Contrary to expectations based on local adaptation gradients, many more HDRs were found between the Shenzhen and Wenchang populations than between the Fuding and Wenchang populations (Figure 4C). This inconsistency could stem from more substantial gene flow from Fuding to Wenchang compared to from Shenzhen to Wenchang (Figure 3). Both theoretical and empirical studies have shown that gene flow can either promote or counteract local adaptation depending on the extent of standing variation and the strength of natural selection (Tigano & Friesen, 2016). The scant number of selected candidate genes identified between the Wenchang and Fuding populations suggests that the weak selection in the Wenchang population is insufficient to overcome the homogenizing effect of gene flow. Nevertheless, we still identified a candidate gene under selection, PUB45, in both the comparison between the Fuding and Wenchang populations and the comparison between the Shenzhen and Wenchang populations (Figure 4D). This gene, which belongs to ubiquitin ligases enzymes, was reported to play a role in the cold acclimation of K. obovata seedlings in a previous study using transcriptome analysis (He et al., 2023). Collectively, our findings highlight the complexities involved in identifying loci responsible for local adaptation, especially in the presence of complex demographic histories, population structure, and gene flow.

4.3 Similarity and discrepancy between genes involved in local adaptation and those responsible for inter-species divergence

The question of whether the same genes responsible for local adaptation within a species also contribute to differences between species is a fundamental topic in evolutionary biology. K. obovata and its cold-intolerant relative K. candel are separated by the South China Sea, exhibiting genetic discontinuity and differential adaptation (Sheue et al., 2003). Among the 78 candidate genes within HDRs in K. obovata, only two genes, AIL1 and CALS5 overlap with those identified being under selection by the McDonald-Kreitman (MK) test when comparing K. obovata to K. candel (Figure 4F). This overlap may be due to both intra- and interspecific variations being influenced by similar selective pressures in cold environments. AIL1 has been identified as an ethylene response factor and a transcription factor responsive to cold in grape (Ma et al., 2022). CALS5 plays a crucial role in preserving the proper formation of callose walls during pollen development and responding to biotic stress (Dong et al., 2005). Evidence supporting the contribution of the same genetic variation to adaptive traits within and between species has been found in various organisms, such as antifreeze proteins in fish, indicating that the genetic mechanisms involved may be complex due to the polygenic nature of adaptive traits and the influence of gene-environment interactions (Berthelot et al., 2019).

The lack of overlap for the majority of candidate genes under selection, as identified by between-population comparisons with F_ST and XP-CLR and by inter-species comparison with MK test, can be explained by several factors. Firstly, intra-specific selection within K. obovata may target different genetic variations than inter-specific selection between K. obovata and K. candel. Although temperature adaptation is considered the primary force driving ecological speciation in Kandelia, the two species also differ in shoot, leaf, floral, fruit and hypocotyl characters, in addition to physiological differentiation (Sheue et al., 2003). The 27 selected genes detected by MK test are enriched in functional categories related to several protein modifications (Figure 4G). This is consistent with the notion that inter-specific selection often reflects more fundamental differences in niche occupation and overall lifestyle that have accrued over longer periods of evolutionary time (Tarjuelo et al., 2017). Second, the time scale of selection differs within and between species. Genes that were important in the initial divergence of the cold-tolerant K. obovata and cold-intolerant K. candle might not be the same genes that are currently under selection within K. obovata, which could be adapting to more recent or localized environmental changes. In line with this idea, four out of the top five enriched functional categories of the HDR candidate genes are associated with response to stimulus (Figure 4E). The discrepancies in genes identified as under selection between different methodologies may be a common issue (Balick et al., 2022; Wang X. et al., 2022; Benjelloun et al., 2023).

One limitation of this study is that only three K. obovata populations in China were surveyed. Considering the strong population structure among the surveyed populations (Figure 2D), we might capture both adaptive and neutral changes when identifying candidates under selection as high differential outliers between populations. Nevertheless, the enrichment of temperature-associated SNPs within HDRs indicates that most candidates are likely adaptive. On the other hand, the small number of populations surveyed may lead to an underestimation of genes responsible for local adaptation, as different environmental pressures might act on different gene interaction networks via gene-environment interactions (Tiffin & Ross-Ibarra, 2014). However, as temperature is the major factor shaping the evolutionary trajectories of Kandelia species, our results should capture the majority genetic variation critical for climate adaptation within and between species. A wider geographic sampling across the complete range of K. obovata and K. candel will be necessary to conduct further research on candidate functional loci associated with their ecological divergence.

5 Conclusion

Our investigation into the local adaptation of Kandelia obovata populations across China provides critical insights into the genetic dynamics and adaptive evolution of mangroves under global climate change. Through whole-genome resequencing analysis, we uncovered a strong population structure with complex demographic events such as expansion, bottlenecks, and gene flow, highlighting the intricate historical context within which local adaptation has occurred. Notably, genetic differentiation is high among the geographically close K. obovata populations, likely due to limited gene flow as a result of human activities. Genome-wide scans of population differentiation pinpointed regions under selective sweeps, with more intense signals in northern populations. Our findings emphasize the importance of temperature in driving genetic adaptations, as opposed to precipitation, which showed no discernible genotype-environment associations. The southmost Wenchang population exhibited minimal selective sweep signatures, reflecting the low selective pressure in this tropic region, potentially confounded by gene flow from other populations. We identified 39 candidate genes with high confidence responsible for local adaptation, enriched in stimulus response functions and largely different from those genes involved in the speciation of K. obovata from K. candel, which are associated with basic cellular functions. These findings set the stage for further research to explore the molecular basis of local adaptation and resilience of mangroves to environmental stress. Such knowledge is vital for conservation strategies and predictive modeling of species responses in an era of rapid climate change.

References

• 1

  AlexanderD. H.NovembreJ.LangeK. (2009). Fast model-based estimation of ancestry in unrelated individuals. Genome Res.19, 1655–1664. doi: 10.1101/gr.094052.109

  • CrossRef
  • Google Scholar
  • View reference in article

• 2

  AltschulS. F.GishW.MillerW.MyersE. W.LipmanD. J. (1990). Basic local alignment search tool. J. Mol. Biol.215, 403–410. doi: 10.1016/S0022-2836(05)80360-2

  • CrossRef
  • Google Scholar
  • View reference in article

• 3

  AngertA. L.BontragerM. G.AringgrenJ. (2020). What Do We Really Know about Adaptation at Range Edges? Annu. Rev. Ecol. Evol. Syst.51, 341–361. doi: 10.1146/annurev-ecolsys-012120-091002

  • CrossRef
  • Google Scholar
  • View reference in article

• 4

  BalickD. J.JordanD. M.SunyaevS.DoR. (2022). Overcoming constraints on the detection of recessive selection in human genes from population frequency data. Am. J. Hum. Genet.109, 33–49. doi: 10.1016/j.ajhg.2021.12.001

  • CrossRef
  • Google Scholar
  • View reference in article

• 5

  BaronK. N.SchroederD. F.StasollaC. (2012). Transcriptional response of abscisic acid (ABA) metabolism and transport to cold and heat stress applied at the reproductive stage of development in Arabidopsis thaliana. Plant Sci.188, 48–59. doi: 10.1016/j.plantsci.2012.03.001

  • CrossRef
  • Google Scholar
  • View reference in article

• 6

  BenjaminiY.HochbergY. (1995). Controlling the false discovery rate: A practical and powerful approach to multiple testing. J. R. Stat. Soc B.57, 289–300. doi: 10.1111/j.2517-6161.1995.tb02031.x

  • CrossRef
  • Google Scholar

• 7

  BenjellounB.LeempoelK.BoyerF.StuckiS.StreeterI.Orozco-terWengelP.et al. (2023). Multiple genomic solutions for local adaptation in two closely related species (sheep and goats) facing the same climatic constraints. Mol. Ecol.00, e17257. doi: 10.1111/mec.17257

  • CrossRef
  • Google Scholar
  • View reference in article

• 8

  BerthelotC.ClarkeJ.DesvignesT.William DetrichH.FlicekP.PeckL. S.et al. (2019). Adaptation of proteins to the cold in Antarctic fish: A role for methionine? Genome Bio. Evol.11, 220–231. doi: 10.1093/gbe/evy262

  • CrossRef
  • Google Scholar
  • View reference in article

• 9

  BolgerA. M.LohseM.UsadelB. (2014). Trimmomatic: A flexible trimmer for Illumina sequence data. Bioinformatics30, 2114–2120. doi: 10.1093/bioinformatics/btu170

  • CrossRef
  • Google Scholar
  • View reference in article

• 10

  BrauerE. K.AhsanN.DaleR.KatoN.ColuccioA. E.PiñerosM. A.et al. (2016). The Raf-like kinase ILK1 and the high affinity K+ transporter HAK5 are required for innate immunity and abiotic stress response. Plant Physiol.171, 1470–1484. doi: 10.1104/pp.16.00035

  • CrossRef
  • Google Scholar
  • View reference in article

• 11

  CantalapiedraC. P.Herņandez-PlazaA.LetunicI.BorkP.Huerta-CepasJ. (2021). eggNOG-mapper v2: Functional Annotation, Orthology Assignments, and Domain Prediction at the Metagenomic Scale. Mol. Biol. Evol.38, 5825–5829. doi: 10.1093/molbev/msab293

  • CrossRef
  • Google Scholar
  • View reference in article

• 12

  CarleyL. N.MojicaJ. P.WangB.ChenC. Y.LinY. P.PrasadK. V. S. K.et al. (2021). Ecological factors influence balancing selection on leaf chemical profiles of a wildflower. Nat. Ecol. Evol.5, 1135–1144. doi: 10.1038/s41559-021-01486-0

  • CrossRef
  • Google Scholar
  • View reference in article

• 13

  ChenG.GuX.CapinhaC.LeeS. Y.CuiB.YangF.et al. (2023). Large-scale changes in macrobenthic biodiversity driven by mangrove afforestation. J. Appl. Ecol.60, 2066–2078. doi: 10.1111/1365-2664.14462

  • CrossRef
  • Google Scholar
  • View reference in article

• 14

  ChenH.PattersonN.ReichD. (2010). Population differentiation as a test for selective sweeps. Genome Res.20, 393–402. doi: 10.1101/gr.100545.109

  • CrossRef
  • Google Scholar

• 15

  ChiangT. Y.ChiangY. C.ChenY. J.ChouC. H.HavanondS.HongT. N.et al. (2001). Phylogeography of Kandelia candel in East Asiatic mangroves based on nucleotide variation of chloroplast and mitochondrial DNAs. Mol. Ecol.10, 2697–2710. doi: 10.1046/j.0962-1083.2001.01399.x

  • CrossRef
  • Google Scholar
  • View reference in article

• 16

  CingolaniP.PlattsA.WangL. L.CoonM.NguyenT.WangL.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. Fly6, 80–92. doi: 10.4161/fly.19695

  • CrossRef
  • Google Scholar
  • View reference in article

• 17

  DanecekP.AutonA.AbecasisG.AlbersC. A.BanksE.DePristoM. A.et al. (2011). The variant call format and VCFtools. Bioinformatics27, 2156–2158. doi: 10.1093/bioinformatics/btr330

  • CrossRef
  • Google Scholar
  • View reference in article

• 18

  DongX.HongZ.SivaramakrishnanM.MahfouzM.VermaD. P. S. (2005). Callose synthase (CalS5) is required for exine formation during microgametogenesis and for pollen viability in Arabidopsis. Plant J.42, 315–328. doi: 10.1111/j.1365-313X.2005.02379.x

  • CrossRef
  • Google Scholar
  • View reference in article

• 19

  DoolS. E.PickerM. D.EberhardM. J. B. (2022). Limited dispersal and local adaptation promote allopatric speciation in a biodiversity hotspot. Mol. Ecol.31, 279–295. doi: 10.1111/mec.16219

  • CrossRef
  • Google Scholar

• 20

  EckertC. G.SamisK. E.LougheedS. C. (2008). Genetic variation across species’ geographical ranges: The central-marginal hypothesis and beyond. Mol. Ecol.17, 1170–1188. doi: 10.1111/j.1365-294X.2007.03659.x

  • CrossRef
  • Google Scholar
  • View reference in article

• 21

  EllisonA. M. (2000). Mangrove restoration: Do we know enough? Restor. Ecol.8, 219–229. doi: 10.1046/j.1526-100x.2000.80033.x

  • CrossRef
  • Google Scholar
  • View reference in article

• 22

  ExcoffierL.MarchiN.MarquesD. A.Matthey-DoretR.GouyA.SousaV. C. (2021). Fastsimcoal2: Demographic inference under complex evolutionary scenarios. Bioinformatics37, 4882–4885. doi: 10.1093/bioinformatics/btab468

  • CrossRef
  • Google Scholar
  • View reference in article

• 23

  FayJ. C.WuC. I. (2000). Hitchhiking under positive Darwinian selection. Genetics155, 1405–1413. doi: 10.1093/genetics/155.3.1405

  • CrossRef
  • Google Scholar
  • View reference in article

• 24

  ForesterB. R.LaskyJ. R.WagnerH. H.UrbanD. L. (2018). Comparing methods for detecting multilocus adaptation with multivariate genotype–environment associations. Mol. Ecol.27, 2215–2233. doi: 10.1111/mec.14584

  • CrossRef
  • Google Scholar
  • View reference in article

• 25

  FriisG.SmithE. G.LovelockC. E.OrtegaA.MarshellA.DuarteC. M.et al. (2024). Rapid diversification of grey mangroves (Avicennia marina) driven by geographic isolation and extreme environmental conditions in the Arabian Peninsula. Mol. Ecol.33 (4), e17260. doi: 10.1111/mec.17260

  • CrossRef
  • Google Scholar
  • View reference in article

• 26

  GaoJ.ZhangR.ZhengL.SongL.JiM.LiS.et al. (2023). Blue light receptor CRY1 regulates HSFA1d nuclear localization to promote plant thermotolerance. Cell Rep.42 (9), 113117. doi: 10.1016/j.celrep.2023.113117

  • CrossRef
  • Google Scholar
  • View reference in article

• 27

  GengQ.LianC.GotoS.TaoJ.KimuraM.IslamM. S.et al. (2008). Mating system, pollen and propagule dispersal, and spatial genetic structure in a high-density population of the mangrove tree Kandelia candel. Mol. Ecol.17, 4724–4739. doi: 10.1111/j.1365-294X.2008.03948.x

  • CrossRef
  • Google Scholar
  • View reference in article

• 28

  GuoZ.HuangY.ChenY.DukeN. C.ZhongC.ShiS. (2016). Genetic discontinuities in a dominant mangrove Rhizophora apiculata (Rhizophoraceae) in the Indo-Malesian region. J. Biogeogr.43, 1856–1868. doi: 10.1111/jbi.12770

  • CrossRef
  • Google Scholar
  • View reference in article

• 29

  HeS.WangX.DuZ.LiangP.ZhongY.WangL.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

  • CrossRef
  • Google Scholar
  • View reference in article

• 30

  HeZ.FengX.ChenQ.LiL.LiS.HanK.et al. (2022). Evolution of coastal forests based on a full set of mangrove genomes. Nat. Ecol. Evol.6, 738–749. doi: 10.1038/s41559-022-01744-9

  • CrossRef
  • Google Scholar
  • View reference in article

• 31

  HobanS.KelleyJ. L.LotterhosK. E.AntolinM. F.BradburdG.LowryD. B.et al. (2016). Finding the genomic basis of local adaptation: Pitfalls, practical solutions, and future directions. Am. Nat.188, 379–397. doi: 10.1086/688018

  • CrossRef
  • Google Scholar
  • View reference in article

• 32

  HuY.FengC.YangL.EdgerP. P.KangM. (2022). Genomic population structure and local adaptation of the wild strawberry Fragaria nilgerrensis. Hortic. Res.9, uhab059. doi: 10.1093/hr/uhab059

  • CrossRef
  • Google Scholar
  • View reference in article

• 33

  HuM. J.SunW. H.TsaiW. C.XiangS.LaiX. K.ChenD. Q.et al. (2020). Chromosome-scale assembly of the Kandelia obovata genome. Hortic. Res.7, 75. doi: 10.1038/s41438-020-0300-x

  • CrossRef
  • Google Scholar
  • View reference in article

• 34

  IsabelN.HollidayJ. A.AitkenS. N. (2020). Forest genomics: Advancing climate adaptation, forest health, productivity, and conservation. Evol. Appl.13, 3–10. doi: 10.1111/eva.12902

  • CrossRef
  • Google Scholar
  • View reference in article

• 35

  JennerjahnT. C.GilmanE.KraussK. W.LacerdaL. D.NordhausI.WolanskiE. (2017). “Mangrove ecosystems under climate change,” in Mangrove Ecosystems: A Global Biogeographic Perspective: Structure, Function, and Services, 211–244. Eds. Rivera-MonroyV. H.LeeS. Y.KristensenE.TwilleyR. R. (Cham: Springer International Publishing). doi: 10.1007/978-3-319-62206-4_7

  • CrossRef
  • Google Scholar
  • View reference in article

• 36

  KhoonG. W.AksornkoaeS.TriN. H.VongwattanaK.FanH.SantosoN.et al. (2004). Mangroves in the South China Sea (UNEP/GEF/SCS Technical Publication No. 1).

  • Google Scholar
  • View reference in article

• 37

  KimS. J.HeldM. A.ZemelisS.WilkersonC.BrandizziF. (2015). CGR2 and CGR3 have critical overlapping roles in pectin methylesterification and plant growth in Arabidopsis thaliana. Plant J.82, 208–220. doi: 10.1111/tpj.12802

  • CrossRef
  • Google Scholar
  • View reference in article

• 38

  KimS.SoltisP. S.WallK.SoltisD. E. (2006). Phylogeny and domain evolution in the APETALA2-like gene family. Mol. Biol. Evol.23, 107–120. doi: 10.1093/molbev/msj014

  • CrossRef
  • Google Scholar
  • View reference in article

• 39

  KorneliussenT. S.AlbrechtsenA.NielsenR. (2014). ANGSD: analysis of next generation sequencing data. BMC Bioinf.15, 1–13. doi: 10.1186/s12859-014-0356-4

  • CrossRef
  • Google Scholar
  • View reference in article

• 40

  KorunesK. L.SamukK. (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

  • CrossRef
  • Google Scholar
  • View reference in article

• 41

  LantzouniO.AlkoferA.Falter-BraunP.SchwechheimerC. (2020). Growth-regulating factors interact with dellas and regulate growth in cold stress. Plant Cell32, 1018–1034. doi: 10.1105/tpc.19.00784

  • CrossRef
  • Google Scholar
  • View reference in article

• 42

  Le CorreV.KremerA. (2012). The genetic differentiation at quantitative trait loci under local adaptation. Mol. Ecol.21, 1548–1566. doi: 10.1111/j.1365-294X.2012.05479.x

  • CrossRef
  • Google Scholar
  • View reference in article

• 43

  LevinD. A. (2019). Plant speciation in the age of climate change. Ann. Bot.124, 769–775. doi: 10.1093/aob/mcz108

  • CrossRef
  • Google Scholar
  • View reference in article

• 44

  LiH.DurbinR. (2009). Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics25, 1754–1760. doi: 10.1093/bioinformatics/btp324

  • CrossRef
  • Google Scholar
  • View reference in article

• 45

  LiH.HandsakerB.WysokerA.FennellT.RuanJ.HomerN.et al. (2009). The sequence alignment/map format and SAMtools. Bioinformatics25, 2078–2079. doi: 10.1093/bioinformatics/btp352

  • CrossRef
  • Google Scholar
  • View reference in article

• 46

  LinP. (1999). Mangrove ecosystem in China (Beijing, China: Science Press).

  • Google Scholar
  • View reference in article

• 47

  López-GoldarX.AgrawalA. A. (2021). Ecological interactions, environmental gradients, and gene flow in local adaptation. Trends Plant Sci.26, 796–809. doi: 10.1016/j.tplants.2021.03.006

  • CrossRef
  • Google Scholar
  • View reference in article

• 48

  LovelockC. E.CahoonD. R.FriessD. A.GuntenspergenG. R.KraussK. W.ReefR.et al. (2015). The vulnerability of Indo-Pacific mangrove forests to sea-level rise. Nature526, 559–563. doi: 10.1038/nature15538

  • CrossRef
  • Google Scholar
  • View reference in article

• 49

  LuW. X.ZhangB. H.ZhangY. Y.YangS. C. (2022). 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

  • CrossRef
  • Google Scholar
  • View reference in article

• 50

  MaX.ZhaoF.SuK.LinH.GuoY. (2022). Discovery of cold-resistance genes in Vitis amurensis using bud-based quantitative trait locus mapping and RNA-seq. BMC Genomics23, 551. doi: 10.1186/s12864-022-08788-y

  • CrossRef
  • Google Scholar
  • View reference in article

• 51

  MaxwellG. S. (1995). Ecogeographic variation in Kandelia candel from Brunei, Hong Kong and Thailand. Hydrobiologia295, 59–65. doi: 10.1007/BF00029111

  • CrossRef
  • Google Scholar
  • View reference in article

• 52

  McDonaldJ. H.KreitmanM. (1991). Adaptive protein evolution at the Adh locus in Drosophila. Nature351, 652–654. doi: 10.1038/351652a0

  • CrossRef
  • Google Scholar
  • View reference in article

• 53

  McKennaA.HannaM.BanksE.SivachenkoA.CibulskisK.KernytskyA.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

  • CrossRef
  • Google Scholar
  • View reference in article

• 54

  Méndez-CeaB.García-GarcíaI.GazolA.CamareroJ. J.de AndrésE. G.ColangeloM.et al. (2023). Weak genetic differentiation but strong climate-induced selective pressure toward the rear edge of mountain pine in north-eastern Spain. Sci. Total Environ.858, 159778. doi: 10.1016/j.scitotenv.2022.159778

  • CrossRef
  • Google Scholar
  • View reference in article

• 55

  NeiM.MillerJ. C. (1990). A simple method for estimating average number of nucleotide substitutions within and between populations from restriction data. Genetics125, 873–879. doi: 10.1093/genetics/125.4.873

  • CrossRef
  • Google Scholar
  • View reference in article

• 56

  NguyenL. T.SchmidtH. A.Von HaeselerA.MinhB. Q. (2015). IQ-TREE: A fast and effective stochastic algorithm for estimating maximum-likelihood phylogenies. Mol. Biol. Evol.32, 268–274. doi: 10.1093/molbev/msu300

  • CrossRef
  • Google Scholar
  • View reference in article

• 57

  Ohashi-ItoK.BergmannD. C. (2007). Regulation of the Arabidopsis root vascular initial population by Lonesome highway. Development134 (16), 2959–2968. doi: 10.1242/dev.006296

  • CrossRef
  • Google Scholar

• 58

  OksanenJ.LegendreP.O’HaraB.StevensM. H. H.OksanenM. J.SuggestsM. (2017). Community Ecology Package. Available at: https://CRAN.R-project.org/package=vegan.

  • Google Scholar

• 59

  OomenR. A.HutchingsJ. A. (2022). Genomic reaction norms inform predictions of plastic and adaptive responses to climate change. J. Anim. Ecol.91, 1073–1087. doi: 10.1111/1365-2656.13707

  • CrossRef
  • Google Scholar

• 60

  OrsiniL.VanoverbekeJ.SwillenI.MergeayJ.De MeesterL. (2013). Drivers of population genetic differentiation in the wild: Isolation by dispersal limitation, isolation by adaptation and isolation by colonization. Mol. Ecol.22, 5983–5999. doi: 10.1111/mec.12561

  • CrossRef
  • Google Scholar

• 61

  Pérez-FigueroaA.García-PereiraM. J.SauraM.Rolán-AlvarezE.CaballeroA. (2010). Comparing three different methods to detect selective loci using dominant markers. J. Evol. Biol.23, 2267–2276. doi: 10.1111/j.1420-9101.2010.02093.x

  • CrossRef
  • Google Scholar

• 62

  PurcellS.NealeB.Todd-BrownK.ThomasL.FerreiraM. A. R.BenderD.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

  • CrossRef
  • Google Scholar

• 63

  QuisthoudtK.SchmitzN.RandinC. F.Dahdouh-GuebasF.RobertE. M. R.KoedamN. (2012). Temperature variation among mangrove latitudinal range limits worldwide. Trees-Struct. Funct.26, 1919–1931. doi: 10.1007/s00468-012-0760-1

  • CrossRef
  • Google Scholar

• 64

  ShawP. T.ChaoS. Y. (1994). Surface circulation in the South China Sea. Deep-Sea Res. Pt. I: Oceanographic Research Papers41 (11–12), 1663–1683. doi: 10.1016/0967-0637(94)90067-1

  • CrossRef
  • Google Scholar

• 65

  SheueC. R.LiuH. Y.YongJ. W. H. (2003). Kandelia obovata (Rhizophoraceae), a new mangrove species from Eastern Asia. Taxon52, 287–294. doi: 10.2307/3647398

  • CrossRef
  • Google Scholar

• 66

  SinghS.RangarajanS.ThattaiD. (2022). Mangroves: Sea level rise and elevation - A review. AIP Conference Proceedings. 2515 (1), 020015. doi: 10.1063/5.0102952

  • CrossRef
  • Google Scholar

• 67

  SongS. K.ClarkS. E. (2005). POL and related phosphatases are dosage-sensitive regulators of meristem and organ development in Arabidopsis. Dev. Biol.285, 272–284. doi: 10.1016/j.ydbio.2005.06.020

  • CrossRef
  • Google Scholar

• 68

  SongS.DeyD. K.HolsingerK. E. (2006). Differentiation among populations with migration, mutation, and drift: Implications for genetic inference. Evolution60, 1–12. doi: 10.1554/05-315.1

  • CrossRef
  • Google Scholar

• 69

  SongH.KempD. B.TianL.ChuD.SongH.DaiX. (2021). Thresholds of temperature change for mass extinctions. Nat. Commun.12, 4694. doi: 10.1038/s41467-021-25019-2

  • CrossRef
  • Google Scholar

• 70

  SribiantiI.MuthmainnahHikmahKiswandi. (2021). Economic valuation of mangrove ecosystem environmental services based on green economy. IOP Conf. Series: Earth Environ. Sci.886 (1), 012116. doi: 10.1088/1755-1315/886/1/012116

  • CrossRef
  • Google Scholar

• 71

  SuJ.FriessD. A.GasparatosA. (2021). A meta-analysis of the ecological and economic outcomes of mangrove restoration. Nat. Commun.12, 5050. doi: 10.1038/s41467-021-25349-1

  • CrossRef
  • Google Scholar

• 72

  TajimaF. (1989). Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. Genetics123, 585–595. doi: 10.1093/genetics/123.3.585

  • CrossRef
  • Google Scholar

• 73

  TarjueloR.MoralesM. B.ArroyoB.MañosaS.BotaG.CasasF.et al. (2017). Intraspecific and interspecific competition induces density-dependent habitat niche shifts in an endangered steppe bird. Ecol. Evol.7, 9720–9730. doi: 10.1002/ece3.3444

  • CrossRef
  • Google Scholar

• 74

  TiffinP.Ross-IbarraJ. (2014). Advances and limits of using population genetics to understand local adaptation. Trends Ecol. Evol.29, 673–680. doi: 10.1016/j.tree.2014.10.004

  • CrossRef
  • Google Scholar

• 75

  TiganoA.FriesenV. L. (2016). Genomics of local adaptation with gene flow. Mol. Ecol.25, 2144–2164. doi: 10.1111/mec.13606

  • CrossRef
  • Google Scholar

• 76

  TomlinsonP. B. (1986). The botany of mangroves. (Cambridge University Press). doi: 10.2307/4109956

  • CrossRef
  • Google Scholar

• 77

  ValimH. F.GrandeF. D.WongE. L. Y.SchmittI. (2023). Circadian clock- and temperature-associated genes contribute to overall genomic differentiation along elevation in lichenized fungi. Mol. Ecol.33 (4), e17252. doi: 10.1111/mec.17252

  • CrossRef
  • Google Scholar

• 78

  Van der StockenT.VanschoenwinkelB.CarrollD.CavanaughK. C.KoedamN. (2022). Mangrove dispersal disrupted by projected changes in global seawater density. Nat. Clim. Change12, 685–691. doi: 10.1038/s41558-022-01391-9

  • CrossRef
  • Google Scholar

• 79

  VanwallendaelA.SoltaniA.EmeryN. C.PeixotoM. M.OlsenJ.LowryD. B. (2019). A molecular view of plant local adaptation: incorporating stress-response networks. Annu. Rev. Plant Biol.70, 559–583. doi: 10.1146/annurev-arplant-050718-100114

  • CrossRef
  • Google Scholar

• 80

  von TakachB.AhrensC. W.LindenmayerD. B.BanksS. C. (2021a). Scale-dependent signatures of local adaptation in a foundation tree species. Mol. Ecol.30, 2248–2261. doi: 10.1111/mec.15894

  • CrossRef
  • Google Scholar

• 81

  von TakachB.PentonC. E.MurphyB. P.RadfordI. J.DaviesH. F.HillB. M.et al. (2021b). Population genomics and conservation management of a declining tropical rodent. Heredity126, 763–775. doi: 10.1038/s41437-021-00418-9

  • CrossRef
  • Google Scholar

• 82

  WangJ.HuZ.LiaoX.WangZ.LiW.ZhangP.et al. (2022). Whole-genome resequencing reveals signature of local adaptation and divergence in wild soybean. Evol. Appl.15, 1820–1833. doi: 10.1111/eva.13480

  • CrossRef
  • Google Scholar

• 83

  WangX.HuY.HeW.YuK.ZhangC.LiY.et al. (2022). Whole-genome resequencing of the wheat A subgenome progenitor Triticum urartu provides insights into its demographic history and geographic adaptation. Plant Commun.3, 100345. doi: 10.1016/j.xplc.2022.100345

  • CrossRef
  • Google Scholar

• 84

  WeirB. S.CockerhamC. C. (1984). Estimating F-statistics for the analysis of population structure. Evolution38 (6), 1358–1370. doi: 10.1111/j.1558-5646.1984.tb05657.x

  • CrossRef
  • Google Scholar

• 85

  WuY.RicklefsR. E.HuangZ.ZanQ.YuS. (2018). Winter temperature structures mangrove species distributions and assemblage composition in China. Global Ecol. Biogeogr.27, 1492–1506. doi: 10.1111/geb.12826

  • CrossRef
  • Google Scholar

• 86

  YangG.ZhouR.TangT.ShiS. (2008). Simple and efficient isolation of high-quality total RNA from Hibiscus tiliaceus, a mangrove associate and its relatives. Prep. Biochem. Biotech.38, 257–264. doi: 10.1080/10826060802164991

  • CrossRef
  • Google Scholar

• 87

  YuG.WangL. G.HanY.HeQ. Y. (2012). ClusterProfiler: An R package for comparing biological themes among gene clusters. OMICS: A J. Integr. Biol.16, 284–287. doi: 10.1089/omi.2011.0118

  • CrossRef
  • Google Scholar

• 88

  ZhaoY.ZhongY.YeC.LiangP.PanX.ZhangY. Y.et al. (2021). Multi-omics analyses on Kandelia obovata reveal its response to transplanting and genetic differentiation among populations. BMC Plant Biol.21, 1–14. doi: 10.1186/s12870-021-03123-1

  • CrossRef
  • Google Scholar

• 89

  ZhaoZ.DentC.LiangH.LvJ.ShangG.LiuY.et al. (2022). CRY2 interacts with CIS1 to regulate thermosensory flowering via FLM alternative splicing. Nat. Commun.13, 7045. doi: 10.1038/s41467-022-34886-2

  • CrossRef
  • Google Scholar

• 90

  ZhouR.ZengK.WuW.ChenX.YangZ.ShiS.et al. (2007). Population genetics of speciation in nonmodel organisms: I. Ancestral polymorphism in mangroves. Mol. Biol. Evol.24, 2746–2754. doi: 10.1093/molbev/msm209

  • CrossRef
  • Google Scholar