A reference-grade genome assembly for Gossypium bickii and insights into its genome evolution and formation of pigment glands and gossypol

Abstract

The pigment gland is a morphological characteristic of Gossypium and its related genera. Gossypium bickii (G₁) is characterized by delayed pigment gland morphogenesis in the cotyledons. In this study, a reference-grade genome of G₁ was generated, and comparative genomics analysis showed that G₁ was closest to Gossypium australe (G₂), followed by A- and D-genome species. Two large fragment translocations in chromosomes 5 and 13 were detected between the G genome and other Gossypium genomes and were unique to the G₁ and G₂ genomes. Compared with the G₂ genome, two large fragment inversions in chromosomes 12 and 13 were detected in G₁. According to the phylogeny, divergence time, and similarity analysis of nuclear and chloroplast genomes, G₁ was formed by hybridization between Gossypium sturtianum (C₁) and a common ancestor of G₂ and Gossypium nelsonii (G₃). The coordinated expression patterns of pigment gland formation (GoPGF) and gossypol biosynthesis genes in G₁ were verified to be consistent with its phenotype, and nine genes that were related to the process of pigment gland formation were identified. A novel gene, GbiCYP76B6, regulated by GoPGF, was found to affect gossypol biosynthesis. These findings offer insights into the origin and evolution of G₁ and its mechanism of pigment gland formation and gossypol biosynthesis.

Keywords: Gossypium bickii; evolution; genome; gossypol; pigment glands.

Figure 1

Morphology and chromosomal and chloroplast features of G. bickii. (A) The character of pigment glands in G₁. (i) The dormant seed, (ii) the partially enlarged dormant seed, (iii) the seed after germination for 36 h, (iv) the partially enlarged cotyledon after germination for 36 h, (v) the true leaf of the seedling, (vi) the partially enlarged true leaf of the seedling, (vii) the hypocotyls of the seedling, and (viii) the partially enlarged hypocotyls of the seedling. Scale bars, 0.5 mm. (B) Chromosomal features of G₁. (i) 13 pseudochromosomes of the G₁ genome, (ii) gene density, (iii) non-coding RNA (ncRNA) density, (iv) transposable element (TE) density, (v) Gypsy density, (vi) Copia density, and (vii) GC content. The density was calculated in 500-kb intervals across the chromosomes. (C) Chloroplast genome features of G₁. The inner circle line indicates the range of the inverted repeats (IRa and IRb), which separate the genome into small single-copy and large single-copy regions. The inside of the inner circle represents the GC content. The outer circle line represents the gene location information for each class, classified by different colors. The components on the outside of the outer circle are transcribed in the clockwise direction, and on the inside, the components are transcribed in the counterclockwise direction.

Figure 2

Genome evolutionary and structural variation of Gossypium. (A) Phylogenetic tree and divergence times among Gossypium and other species in Malvaceae. (B) The divergence between G₁ and other species by K_s value for orthologous genes between G₁ and other cotton species. (C) Characterization of genomic structural variation among the cotton species. Genomic synteny blocks are connected by gray lines. The large translocations between G groups (G₁ and G₂) and other cotton species are highlighted by red links, the large inversion in chromosome 12 between G₁ and G₂ is highlighted by blue links, and the inversions in chromosome 13 between G₁ and other Gossypium species are highlighted by cyan links. (D) Hi-C heatmap verification of chromosome inversion and translocation fragments. The Hi-C-assisted assembly data for the G₁ genome were used to scan chromosomes 5 and 13 of G₁ and A₁ and chromosomes 12 and 13 of G₁ and G₂, respectively, to verify inversion and translocation fragments. (E) The top 10 enriched KEGG pathways of the translocation and inversion genes in the G₁ genome. All translocation genes in chromosomes 5 and 13 and inversion genes in chromosomes 12 and 13 were selected for KEGG pathway enrichment analysis.

Figure 3

Evolutionary analysis of nuclear and chloroplast genomes of C₁, G₁, G₂, and G₃. (A) Fiber, leaf, sepal, bract, and petal phenotypes of four species in Gossypium. (B) Nuclear and chloroplast genome phylogenetic trees and divergence time analysis of four species in Gossypium. The phylogenetic tree for the nuclear genome was constructed based on SNPs, and the chloroplast genome was based on 72 single-copy genes from chloroplast genomes. (C) Identity scores of nuclear and chloroplast genomes for four species in Gossypium. The calculation of nuclear genome similarity was based on SNPs, and the chloroplast genome was based on whole-genome pairwise comparison. (D) Model for the evolution of four species in Gossypium. Green and yellow represent nuclear and chloroplast genomes, respectively. The blue numbers represent the divergence times.

Figure 4

Morphology, gossypol contents, and differential expression analysis of A₁, A₂, and G₁. (A) Morphology of ovules at 10, 20, and 30 dpa; dormant mature seed (0 h); and seeds at 12, 24, 36, and 48 h after germination. Scale bars, 0.5 mm. (B) Glanded and glandless cotyledons in A₁, A₂, and G₁. The glandless cotyledons at 10 dpa and glanded cotyledons at 20 dpa in A₁ and A₂ and the glandless cotyledon at 24 h after germination and glanded cotyledon at 36 h after germination in G₁ are shown. (C) Gossypol contents in ovules at 10, 20, and 30 dpa and in seeds at 0, 12, 24, 36, and 48 h after germination in A₁, A₂, and G₁ (mean ± SD, n = 3). (D) Expression of differentially expressed genes (DEGs) in the gossypol biosynthesis pathway, represented by a heatmap and estimated by computing fragments per kilobase of transcript per million mapped fragments (FPKM) values. (E) Common and specific genes in five DEG lists: four upregulated groups (A₁_30 dpa versus A₁_10 dpa, A₂_30 dpa versus A₂_10 dpa, G₁_12 h versus G₁_0 h, and G₁_24 h versus G₁_0 h) and one not differentially expressed group (G₁_30 dpa versus G₁_10 dpa). (F) Common and specific genes in the five common DEGs from (D) and DEGs of CCRI12 versus CCRI12gl. (G) Correlations between GoPGF and 18 other DEGs. These genes were selected by computing the Pearson correlation coefficient (r ≥ 0.5, P ≤ 0.05) from 29 common genes in (E). (H) Expression of nine functional genes selected in (G) represented by a heatmap and estimated by computing FPKM values.

Figure 5

Function and mechanism of GbiCYP76B6 in gossypol biosynthesis. (A) Relative GbiCYP76B6 expression levels and gossypol contents in leaves and stems of seedlings after the silencing of GbiCYP76B6 in G₁. TRV:00 was used as a negative control. Independent experiments with mean ± SD, n = 3; ∗P < 0.05. (B) Relative expression levels of GbiCYP76B6 in GoPGF-silenced plants and FPKM expression values of GhCYP76B6 in CCRI12 (glanded) and CCRI12gl (glandless). ∗P < 0.05; ∗∗P < 0.01. (C) Relationship between CYP76B6 and gossypol biosynthesis in the terpene biosynthesis pathway. (D)cis-acting elements in the 2000-bp promoter upstream of the GbiCYP76B6 coding region. (E) Mechanistic model of GoPGF and CYP76B6 in gossypol biosynthesis.