Revealing Genetic Differences in Fiber Elongation between the Offspring of Sea Island Cotton and Upland Cotton Backcross Populations Based on Transcriptome and Weighted Gene Coexpression Networks

Abstract

Fiber length is an important indicator of cotton fiber quality, and the time and rate of cotton fiber cell elongation are key factors in determining the fiber length of mature cotton. To gain insight into the differences in fiber elongation mechanisms in the offspring of backcross populations of Sea Island cotton Xinhai 16 and land cotton Line 9, we selected two groups with significant differences in fiber length (long-fiber group L and short-fiber group S) at different fiber development stages 0, 5, 10 and 15 days post-anthesis (DPA) for transcriptome comparison. A total of 171.74 Gb of clean data was obtained by RNA-seq, and eight genes were randomly selected for qPCR validation. Data analysis identified 6055 differentially expressed genes (DEGs) between two groups of fibers, L and S, in four developmental periods, and gene ontology (GO) term analysis revealed that these DEGs were associated mainly with microtubule driving, reactive oxygen species, plant cell wall biosynthesis, and glycosyl compound hydrolase activity. Kyoto encyclopedia of genes and genomes (KEGG) pathway analysis indicated that plant hormone signaling, mitogen-activated protein kinase (MAPK) signaling, and starch and sucrose metabolism pathways were associated with fiber elongation. Subsequently, a sustained upregulation expression pattern, profile 19, was identified and analyzed using short time-series expression miner (STEM). An analysis of the weighted gene coexpression network module uncovered 21 genes closely related to fiber development, mainly involved in functions such as cell wall relaxation, microtubule formation, and cytoskeletal structure of the cell wall. This study helps to enhance the understanding of the Sea Island-Upland backcross population and identifies key genes for cotton fiber development, and these findings will provide a basis for future research on the molecular mechanisms of fiber length formation in cotton populations.

Keywords: G. hirsutum × G. barbadense; RNA-seq; STEM; WGCNA; backcross population; fiber elongation.

Figure 1

Relationship between 24 fiber samples. (A) Pearson’s correlation coefficient and clustering heatmap. The abscissa and ordinate are the number of samples; the order is determined by the correlation clustering results, and the color reflects the correlation between samples. (B) Principal component analysis of identified genes. L represents an equal mix of HL-9, HL-62, and HL-194 RNA samples; S represents an equal mix of HL-34, HL-78 and HL-159 RNA samples; 0 means 0 days after flowering; values 1, 2, 3 represent different biological replicates.

Figure 2

Confirmation of RNA-seq results using qRT–PCR experiments. (A) qRT–PCR results of 8 genes. With associated significance levels (p values) added above each histogram. ns indicates no significant effect at the p < 0.05 level, and ** and *** represent significant effect at the p < 0.01 and p < 0.001 levels, respectively. Pink bars, L; blue bars, S. (B) RNA-seq heatmap of 8 genes.

Figure 3

Statistical and enrichment analysis of DEGs. (A) Multiple comparisons of L and S in fiber development. Numbers indicate the number of DEGs. Red, upregulated; green, downregulated. (B) GO-term enrichment analysis of 6055 genes. (C) KEGG enrichment analysis of 6055 genes.

Figure 4

Gene expression patterns of L and S and GO enrichment of profile 19. (A) Gene expression patterns for the four developmental stages of L and S as inferred from STEM analysis. Each square represents the expression trend, and the text indicates the number of genes and IDs contained in the profile. (B) GO enrichment of L and S in Profile 19. The p value indicates the significance of GO terms. Gray represents no enrichment.

Figure 5

Venn diagram and KEGG pathway annotation of DEGs. (A) Number of DEGs shared by L and S in profile 19 and 5 DPA comparisons. (B) Number of DEGs shared by L and S in profile 19 and 10 DPA comparisons. (C) Number of DEGs shared by L and S in profile 19 and 15 DPA comparison. (D) KEGG pathway annotations for the 167 DEGs in A. (E) KEGG pathway annotations for the 138 DEGs in B. (F) KEGG pathway annotations for the 144 DEGs in C.

Figure 6

Weighted gene correlation network analysis. (A) Hierarchical tree diagram representing the WGCNA identification module. (B) Correlation diagram between modules and samples. The horizontal axis represents samples at different fiber stages, and the vertical axis represents modules. The numbers in the squares represent the correlation coefficients and p values between modules and traits.

Figure 7

Coexpression network analysis of the (A) magenta module, (B) black module, (C) brown module, and (D) blue module. The central hub gene is represented by an orange circle.