Comparative Transcriptomic Analysis to Identify the Genes Related to Delayed Gland Morphogenesis in Gossypium bickii

Abstract

Cotton is one of the major industrial crops that supply natural fibers and oil for industries. This study was conducted to understand the mechanism of delayed gland morphogenesis in seeds of Gossypium bickii. In this study, we compared glandless seeds of G. bickii with glanded seeds of Gossypium arboreum. High-throughput sequencing technology was used to explore and classify the expression patterns of gland-related genes in seeds and seedlings of cotton plants. Approximately 131.33 Gigabases of raw data from 12 RNA sequencing samples with three biological replicates were generated. A total of 7196 differentially-expressed genes (DEGs) were identified in all transcriptome data. Among them, 3396 genes were found up-regulated and 3480 genes were down-regulated. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) annotations were performed to identify different functions between genes unique to glandless imbibed seeds and glanded seedlings. Co-expression network analysis revealed four modules that were identified as highly associated with the development of glandless seeds. Here the hub genes in each module were identified by weighted gene co-expression network analysis (WGCNA). In total, we have selected 13 genes involved in transcription factors, protein and MYB-related functions, that were differentially expressed in transcriptomic data and validated by quantitative reverse-transcription polymerase chain reaction (qRT-PCR). These selected genes may play an important role for delayed gland morphogenesis. Our study provides comprehensive insight into the key genes related to glandless traits of seeds and plants, and can be further exploited by functional and molecular studies.

Keywords: DEGs; Gossypium bickii; RNA-seq; WGCNA.

Figure 1

Pictorial description of delayed gland morphogenesis in Gossypium bickii and regular gland formation in Gossypium arboreum showing glands on seeds and germination stages. (a) Imbibed seed image of G. bickii showing no glands; (b) seed germination stage of G. bickii showing glands on cotyledons and hypocotyl; (c,d) image of G. arboreum imbibed seed and germination showing glands; (e,f) gland formation in cotyledon and hypocotyl of G. bickii; and (g,h) gland formation in cotyledon and hypocotyl of G. arboreum.

Figure 2

Principal component analysis (PCA) of genes identified from 12 samples with three biological replicates. Gbdd, Ga24h, Gbgl and Ga48h represent glanded seedlings of G. bickii, glanded imbibed seeds of G. arboreum, glandless imbibed seeds of G. bickii and glanded seedlings of G. arboretum, respectively.

Figure 3

Expression dynamics changes and comparative analysis of differentially-expressed genes (DEGs) between Gbgl, Ga24h, Gbdd and Ga48h following delayed gland morphogenesis at imbibed seed and seedling stages. (a) Number of DEGs showing the up-regulated and down-regulated genes; (b) Venn diagram showing common genes, all differentially-expressed genes in different stages; (c) the number of transcripts demonstrating changes in expression in Gbgl, Ga24h, Gbdd and Ga24h following pattern fold change (FC) in expression, calculated as the log2 ratio of gene expression in glandless imbibed seed with glanded imbibed seed and seedling. Gbdd, Ga24h, Gbgl and Ga48h represent G. bickii glanded seedlings, G. arboreum glanded imbibed seeds, G. bickii glandless imbibed seeds and G. arboreum glanded seedlings, respectively.

Figure 4

The expression pattern of DEGs from different biological replicates. A heat map represents the relative expression levels of genes based on fragment per kb per million of the mapped reads (FPKM) values using RNA sequencing (RNA-seq) data. Gbdd, Ga24h, Gbgl and Ga48h represent G. bickii glanded seedlings, G. arboreum glanded imbibed seeds, G. bickii glandless imbibed seeds and G. arboreum glanded seedlings, respectively.

Figure 5

Gene ontology functional classification of DEGs.

Figure 6

KEGG pathway enrichment analysis of DEGs.

Figure 7

Determination of soft-thresholding power in the gene weighted co-expression network analysis. (a) Analysis of the scale-free fit index for various soft-thresholding powers (β). (b) Analysis of the mean connectivity for various soft-thresholding powers. (c) Histogram of connectivity distribution when β = 8. (d) Checking the scale-free topology when β = 8.

Figure 8

(a) Weighted gene co-expression network analysis (WGCNA), all expressed genes in Gbgl, Gbdd, Ga24h and Ga48h RNA-seq samples, showing hierarchical dendrogram co-expression of identified modules by WGCNA; each leaf represents a single gene. (b) The heat map plot of the gene network analysis. Dark red color shows overlap of highly expressed genes which pair from respective datasets, and the lighter color represents the low overlap genes. Different modules are indicated by the diagonals.

Figure 9

Validation of 13 DEGs related to delayed gland morphogenesis identified from transcriptome analysis data by quantitative reverse-transcription polymerase chain reaction (qRT-PCR). (a) RNA-seq-based log2-fold change expression; (b) qRT-PCR-based relative expression profile.