The genome of okra (Abelmoschus esculentus) provides insights into its genome evolution and high nutrient content

Abstract

Okra (Abelmoschus esculentus) is an important vegetable crop with high nutritional value. However, the mechanism underlying its high nutrient content remains poorly understood. Here, we present a chromosome-scale genome of okra with a size of 1.19 Gb. Comparative genomics analysis revealed the phylogenetic status of A. esculentus, as well as whole-genome duplication (WGD) events that have occurred widely across the Malvaceae species. We found that okra has experienced three additional WGDs compared with the diploid cotton Gossypium raimondii, resulting in a large chromosome number (2n = 130). After three WGDs, okra has undergone extensive genomic deletions and retained substantial numbers of genes related to secondary metabolite biosynthesis and environmental adaptation, resulting in significant differences between okra and G. raimondii in the gene families related to cellulose synthesis. Combining transcriptomic and metabolomic analysis, we revealed the relationship between gene expression and metabolite content change across different okra developmental stages. Furthermore, the sinapic acid/S-lignin biosynthesis-related gene families have experienced remarkable expansion in okra, and the expression of key enzymes involved in the sinapic acid/S-lignin biosynthesis pathway vary greatly across developmental periods, which partially explains the differences in metabolite content across the different stages. Our study gains insights into the comprehensive evolutionary history of Malvaceae species and the genetic basis that underlies the nutrient content changes in okra, which will facilitate the functional study and genetic improvement of okra varieties.

Figure 1

Overview of the Abelmoschus esculentus genome assembly. Circos representation of A. esculentus genome features. The outermost layer represents the 65 pseudo-chromosomes (scale mark = 1 Mb) (I). The density of protein-coding genes (II), GC content (III), GC skew (IV), LTR/Gypsy (V), and LTR/Copia (VI) were computed using a 200 kb non-overlapping window. The innermost track (VII) shows the synteny relationship between chromosomes, with colored curves displaying the inter-subgenome syntenic regions.

Figure 2

Comparative genomics analysis between Abelmoschus esculentus and neighboring species. a Phylogenomic analysis of A. esculentus with 13 other plant species. The species divergence times (marked in green at each node) were estimated using r8s. The gain and loss of gene families was analysed using CAFE. The number of gene families experiencing expansion and contraction are enumerated below the species names in orange and blue, respectively. Different categories of orthologous genes across all species are displayed as stacked bar charts on the right. b Synteny analysis between the two subgenomes of A. esculentus and Hibiscus cannabinus. Collinear analysis reflects the correspondence of chromosomes between the A and B subgenomes and chromosomal rearrangement events. For the chromosome pairs undergoing rearrangements, different colors were assigned for the syntenic blocks. c Venn diagram showing the numbers of common and specific gene families in each species. d KEGG enrichment analysis of gene families showing significant expansion in A. esculentus after divergence from the most recent common ancestor. e Phylogenetic tree of HD-Zip I subfamily members in A. esculentus, G. raimondii, and Arabidopsis thaliana.

Figure 3

Widespread whole-genome duplication (WGD) and karyotype evolution in Malvaceae. a The histogram of inter- or intra-species synonymous substitution rates (Ks) indicates whole-genome duplication (WGD) events and species divergence time. The four dashed lines indicate the positions of Ks peaks in Malvaceae species. b Synteny depth ratio between Abelmoschus esculentus and Gossypium raimondii. c Synteny depth ratio between A. esculentus and Hibiscus cannabinus. d Synteny depth ratio between A. esculentus and Hibiscus mutabilis. e Synteny depth ratio between H. cannabinus and H. mutabilis. f Dot plot showing the syntenic relationship between A. esculentus and H. cannabinus (4:1, left panel), and H. mutabilis (4:3, right panel). g Chromosome evolution of A. esculentus and closely related species in the family Malvaceae based on the karyotypes inferred using WGDI. The inferred WGD events are indicated on the corresponding branches. The divergence times are displayed on the left side. WGD I was shared by the five Malvaceae species. WGD II occurred individually in Hibiscus syriacus. WGD III was shared by H. cannabinus, H. mutabilis, and A. esculentus. WGD IV and V represent two WGD events that were unique to A. esculentus. WGT VI was a triplication event that occurred individually in H. mutabilis. Blue ellipses indicate WGD events, and the red star represents the WGT event. The karyotype for H. syriacus is unknown (shown in grey) due to the lack of chromosome-level genome assembly.

Figure 4

Huge deletion and functional differences revealed by comparison between Abelmoschus esculentus and Gossypium raimondii.a Dot plot showing the synteny comparison between A. esculentus and G. raimondii. The colored blocks represent the most significant 8:1 syntenic relationship. b Syntenic comparison between A. esculentus and G. raimondii revealing the different syntenic pattern. The red lines represent a 1:8 syntenic relationship, the green lines represent a 1:16 syntenic relationship. c Bar plot showing the distribution of different gene ratios between A. esculentus and G. raimondii. The colored parts of the pie chart correspond to A. esculentus genes with different ratios, while the white part represents unique genes in A. esculentus. d KEGG enrichment analysis of A. esculentus genes showing different ratios (ranging from 1:2 to 1:8) against G. raimondii. e Microsynteny analysis of G. raimondii D9 and A. esculentus A27, A28, A29, A30, B27, B28, B29, and B30. Colored lines represent different A. esculentus chromosomes (top). Microsynteny analysis of G. raimondii D9 chromosome and eight A. esculentus chromosomes, and the red lines indicate the PPR-rich genomic regions (bottom). f Phylogenetic tree of the CesA and Csl subfamily genes in G. raimondii and A. esculentus. The CesA subfamily was divided into six clades, and the Csl subfamily was grouped into eight clades. g Expression pattern of CesA and Csl genes across five stages of A. esculentus. The colored rectangles on the left represent different clades in the phylogeny. The members of the CesA clades and the CslA clade are displayed in the heatmap because the members of other Csl clades had only weak or undetectable expression levels. To achieve a better heatmap effect, the FPKM value of genes with FPKM higher than 30 was uniformly set to 30 during the generation of the heatmap.

Figure 5

Insights into the biosynthesis pathway of sinapic acid and S-lignin in okra based on transcriptome and metabolome data. a Heatmap displaying the content of secondary metabolites identified in different okra samples by untargeted metabolite profiling. b Gene co-expression subnetwork module significantly related to sinapic acid content. The orange nodes represent five key enzymes participating in S-lignin biosynthesis in this subnetwork, and the green lines represent genes with significant associations with these enzyme genes. c The content of metabolites in the S-lignin biosynthesis pathway detected in different tissues. Different letters above columns represent a significant difference between the samples examined. Statistical test was carried out using one-way ANOVA (analysis of variance). d The inferred S-lignin biosynthesis pathway in Abelmoschus esculentus. The numbers above each gene symbol represent the copy number in the corresponding species. The heatmap displays the transcript levels of enzyme-encoding genes in this pathway across the five stages. e Phylogenetic tree of CCR enzyme genes in A. esculentus, Gossypium raimondii, Hibiscus syriacus, Hibiscus cannabinus, Hibiscus mutabilis, and Theobroma cacao.