Population genomics and comparative transcriptomics analysis provide insights into evolution and cold adaptation of mangrove Kandelia obovata

Figure 1. Genomic features of K. obovata and K. candel. (A) Lengths of 18 chromosomes. The left and right part of the circle indicates the chromosomes of K. obovata and K. candel, respectively. (B) Repeat sequence density with window size of 300 kb. (C) GC content density with window size of 300 kb. (D) Gene number density with window size of 300 kb. (E) Inter-genome collinear blocks between K. obovata and K. candel.

Figure 2. Genome comparison between K. obovata and K. candel. (A) Syntenic analysis of K. obovata (KO) and K. candel (KC) genomes. (B) Histogram of duplication depths of syntenic regions between K. obovata and K. candel. (C) Structural variations (SV) distribution including syntenic, inversion, translocation and duplication between K. obovata and K. candel genomes. K. obovata and K. candel genomes are the reference genome and query genome, respectively. (D) Ks distribution between homologous gene pairs in the syntenic blocks of KC versus KC, KO versus KC, and KO versus KO, respectively. (E) GO enrichment analysis of genes involved in structural variation and genes from syntenic gene pairs with Ks < 1.

Figure 3. Gene families in K. obovata and K. candel. (A) Venn diagram of specific and common gene families among 9 species (Arabidopsis thaliana, Avicennia marina, Kandelia candel, Kandelia obovata, Oryza sativa, Populus trichocarpa, Ricinus communis, Sonneratia alba, and Sonneratia caseolaris). (B) The upper panel indicates the global average temperature (GAT) curve in the past 200 million years (Scotese et al., 2021). The lower panel represents a phylogenetic tree of gene families including the number of expanded and contracted gene families among 9 species. WGD represents whole genome duplication. (C) GO enrichment analysis of specific, expanded, and contracted gene families of K. obovata (KO) and K. candel (KC), respectively. (D) The number of TFs from expanded and contracted gene families of KO and KC. (E) Heatmap shows the number of differentially expressed TFs in C1-C8 clusters. The number of DEGs, clustering genes and total TFs are shown inside the parenthesis.

Figure 4. Population structure and genetic diversity analysis. (A) PCA analysis of 53 individuals of K. obovata (KO) and K. candel (KC) populations. (B) A phylogenetic tree of 53 resequencing samples based on SNPs. (C) Rate of linkage disequilibrium (LD) decay for all KO populations. (D) Population structure based on SNPs (K = 2 to 5) using admixture analysis. The colors in each column indicate the proportion of resequencing samples for each population. (E) Distribution of Fst values between KO and KC from 53 resequencing samples.

Figure 5. Genomic variants related to bioclimatic variables. (A) Heatmap shows Pearson’s correlation of 19 environmental variables (bio1-19). (B) Scree plot shows the percentage of explained variances by the top 10 principal components of 19 environmental variables. (C) The scatter plot shows the correlation between latitude and the first principal component (PC1) based on the Pearson’s correlation coefficients (PCC). (D) The number of SNPs associated with latitude, temperature (bio4 and bio7), and precipitation (bio14, bio17, and bio19) by genotype-environment association analysis. (E) The histogram shows the numbers of candidate genes regulated by environment-associated SNPs. (F) Venn diagram shows the numbers of overlapped genes in TFs, temperature-associated genes, precipitation-associated genes, and latitude-associated genes. (G) GO and KEGG enrichment analysis of TFs associated with temperature, precipitation, and latitude, respectively.

Figure 6. Identification of candidate genes for cold adaptation in K. obovata. (A) Venn diagram shows the numbers of overlapped temperature-associated genes identified in this study and previous study (Zou et al., 2024). (B) The numbers of temperature-associated TFs detected by different strategies. (C) KEGG enrichment analysis of 108 temperature-associated TFs. The length of the box indicates -log10 P-value. (D) KEGG enrichment analysis of DEGs identified during K. obovata cold acclimation. (E) Heatmap shows the gene expression pattern of 18 circadian rhythm-related TFs during K. obovata cold acclimation. (F) Construction of regulatory network mediated by 18 circadian rhythm-related TFs. (G) GO and KEGG enrichment analysis of target genes regulated by 18 circadian rhythm-related TFs. The red, green, and blue colors represent KEGG pathway, biological process and molecular function, respectively. (H) Volcano plot shows differentially abundant metabolites in K. obovata under chilling stress, including chilling-responsive flavonoids and phenylpropanoids metabolites.

Figure S1. Quality assessments of K. obovata and K. candel assemblies and gene functional annotation. (A) BUSCO evaluation of K. obovata and K. candel assemblies with genome, transcriptome, and protein modes. (B) Evaluation of K. obovata and K. candel genome quality based on LAI scores. Gold: LAI ≥ 20; reference: 10 ≤ LAI < 20; draft: 0 < LAI < 10. (C) Transposable elements (TE) content of K. obovata and K. candel genomes including five types (DNA, LINEs, LTR, SINEs and unclassified type). (D) Functional annotation of predicted protein-coding genes in K. obovata and K. candel based on COG, eggNOG, GO, KEGG, KOG, NR, Pfam, Swiss-Prot, and TrEMBL databases.

Figure S2. Statistics of genomic structural variations. (A) Distribution of genomic structural variations including insertion, deletion, repeat expansion, repeat contraction, tandem expansion, and tandem contraction on 18 chromosomes between K. obovata and K. candel genomes, with K. obovata serving as the reference genome. (B) The number of sequence length of structural variations between K. obovata and K. candel genomes. The Y axis indicates log10 (count + 1).

Figure S3. Identification and KEGG enrichment analysis of cold-responsive genes in K. obovata. (A) The number of upregulated and downregulated DEGs identified during K. obovata cold acclimation. (B) Venn diagram of DEGs and clustering genes. (C) Selection of clustering number based on the minimum centroid distance. (D) The pathway enrichment analysis of clustering genes from C1, C2, C4, C5, C6, and C8 clusters.

Figure S4. Population structure and enrichment analysis. (A) PCA analysis of K. obovata (KO) including Yunxiao KO population. (B) Rate of linkage disequilibrium (LD) decay for all KO and Yunxiao KO populations. (C) KEGG enrichment analysis of 116 TFs in the selective sweep regions. (D) GO enrichment analysis of 116 TFs in the selective sweep regions. The blue and yellow colors indicate biological process and molecular function, respectively.

Figure S5. The correlation between latitude and 19 bioclimatic variables based on the Pearson’s correlation coefficients (PCC). A P-value < 0.05 was considered significant correlation.

Figure S6. Manhattan plot of genomic loci associated with latitude, bio4, bio7, bio14, bio17, and bio19 along the 18 K. obovata chromosomes, respectively. The Y axis indicates -log10 P-value.

Table S1. The detailed information of genus Kandelia samples collection sites.

Table S2. Comparison of the new and old versions of K. candel genome.

Table S3. Lists of collinear gene pairs in WGD event.

Table S4. Orthologous gene families among 9 plant species.

Table S5. A summary of public sequences data used in this study.

Table S6. Lists of differentially expressed genes during cold acclimation and de-acclimation.

Table S7. Description of global climate data from WorldClim database including 19 bioclimatic variables.

Table S8. Identification of candidate temperature-associated transcription factors using different strategies.

Table S9. The functional annotation of 18 circadian rhythm-related TFs in K. obovata based on UniProt database.

Table S10. The functional enrichment analysis of target genes regulated by 18 circadian rhythm-related TFs.

Table S11. Statistics of flavonoid and phenylpropanoid biosynthetic pathway encoding genes.