Comparative Transcriptomic Analysis of Gossypium hirsutum Fiber Development in Mutant Materials (xin w 139) Provides New Insights into Cotton Fiber Development

Abstract

Cotton is the most widely planted fiber crop in the world, and improving cotton fiber quality has long been a research hotspot. The development of cotton fibers is a complex process that includes four consecutive and overlap** stages, and although many studies on cotton fiber development have been reported, most of the studies have been based on cultivars that are promoted in production or based on lines that are used in breeding. Here, we report a phenotypic evaluation of Gossypium hirsutum based on immature fiber mutant (xin w 139) and wild-type (**n W 139) lines and a comparative transcriptomic study at seven time points during fiber development. The results of the two-year study showed that the fiber length, fiber strength, single-boll weight and lint percentage of xin w 139 were significantly lower than those of **n W 139, and there were no significant differences in the other traits. Principal component analysis (PCA) and cluster analysis of the RNA-sequencing (RNA-seq) data revealed that these seven time points could be clearly divided into three different groups corresponding to the initiation, elongation and secondary cell wall (SCW) synthesis stages of fiber development, and the differences in fiber development between the two lines were mainly due to developmental differences after twenty days post anthesis (DPA). Differential expression analysis revealed a total of 5131 unique differentially expressed genes (DEGs), including 290 transcription factors (TFs), between the 2 lines. These DEGs were divided into five clusters. Each cluster functional category was annotated based on the KEGG database, and different clusters could describe different stages of fiber development. In addition, we constructed a gene regulatory network by weighted correlation network analysis (WGCNA) and identified 15 key genes that determined the differences in fiber development between the 2 lines. We also screened seven candidate genes related to cotton fiber development through comparative sequence analysis and qRT–PCR; these genes included three TFs (GH_A08G1821 (bHLH), GH_D05G3074 (Dof), and GH_D13G0161 (C3H)). These results provide a theoretical basis for obtaining an in-depth understanding of the molecular mechanism of cotton fiber development and provide new genetic resources for cotton fiber research.

1. Introduction

Cotton is a highly significant fiber crop that is extensively cultivated worldwide. Over the years, there has been a strong focus on enhancing the quality of cotton fibers, which has emerged as a prominent area of research. Cotton breeders have been actively engaged in develo** and cultivating germplasm resources that not only yield high quantities but also possess exceptional fiber quality [1,2]. Cotton fiber development is a complex process that occurs in four consecutive and overlap** stages. These stages include fiber initiation, elongation, thickening of the secondary cell wall, and maturation. Fiber initiation occurs from three days before flowering to three days post anthesis (DPA). Elongation occurs from 3 to 16 DPA. During this period, the cotton fibers undergo significant elongation as they grow in length. The secondary cell wall (SCW) stage typically spans from 16 to 40 DPA. During this phase, the cotton fibers undergo SCW deposition, resulting in increased fiber diameter and strength. The thickening of the SCW is an important process in enhancing the overall quality and strength of cotton fibers. The final stage of cotton fiber development is maturation, which occurs from 40 to 50 DPA [3]. The number of fibroblast-differentiated protrusions affects the number of mature fibers to a certain extent, and the progression of rapid elongation through primary wall synthesis determines the length of the mature fibers. Cellulose is deposited in the SCW, and the cell wall is thickened during the thickening period of the SCW stage [4,5,6]. Scientists are working to elucidate the key regulatory mechanisms involved in fiber development, including fibroblast differentiation, cell wall synthesis, and cellulose biosynthesis [7,8,9]. Cotton fibers develop from ovule epidermal cells, and improving cotton fiber quality has long been a research hotspot. The creation and breeding of germplasm resources with high yields and excellent fiber quality has also become a long-term research direction and goal for cotton breeders [1,2]. Studying the regulatory mechanism of cotton fiber cell elongation and secondary wall development is highly important for improving the cotton yield and fiber quality.

In recent years, there has been significant progress in sequencing technologies, leading to their widespread use and continuous optimization. This has resulted in a growing body of evidence highlighting the crucial role of RNA transcriptional regulation in plant growth and development [10]. Transcriptome sequencing (RNA-seq), one of the most commonly used second-generation high-throughput sequencing methods, has been widely used in the study of cotton fibers [11,12]. Through transcriptomics, many key regulatory pathways and gene expression patterns associated with fiber development have been revealed [13,14,15]. These findings have helped identify the key genes that control fiber development and have provided important insights for improving cotton fibers. Through RNA-seq of 0–35 DPA fibers of PimaS-7 and 5917, 4 candidate genes related to fiber strength were identified [16]. RNA-seq analysis at different stages of fiber development (7, 14, and 26 DPA) in the Coker 312 cotton variety allowed the identification of transcription factors and functional genes associated with this process. These genes encode proteins involved in various functional and metabolic pathways, including those involved in catalytic activity, carbohydrate metabolism, cell membrane and organelle functions, and signal transduction [17]. Analysis of the RNA-seq data at 10 and 20 DPA from fibers from 4 wild cotton species and 5 domesticated cotton species showed that wild cotton plants allocate more resources to the stress response pathway and that acclimation may lead to a reprogramming of resource allocation in the direction of increasing fiber growth via regulation of the stress response network [18]. Through RNA-seq analysis of the offspring of the ** stages. Most related studies have been based on cultivars that are promoted in production or based on lines that are used in breeding; especially, there are few reports on natural mutants [25,26,27]. Here, we performed phenotypic evaluation and comparative transcriptomic studies at seven time points during cotton fiber development in immature fiber mutant (xin w 139) and wild-type (** of the im mutant revealed that the gene encoding the PPR protein (Gh_A03G0489) is related to the immature fiber phenotype of the im mutant [37]. Although the expression levels of Gh_A03G0489 were significantly different between TM-1 and IM during the SCW period, the expression levels were relatively low [37]. The PPR gene family is the largest gene family identified to date [38]. They mainly act on mitochondria or chloroplasts to regulate organelle genes at the post-transcriptional level, thereby affecting plant growth and development [38]. We identified 523 genes encoding PPR proteins among the DEGs, and the expression levels in both lines were generally low (FPKM < 2). Notably, we found that two genes encoding PPR proteins (GH_D11G0868 and GH_A05G0653) exhibited a greater than four-fold change in expression during the SCW period (Table S3). Mei et al. reported that GH_A05G0653 was expressed mainly in the A subgenome of 4 G. hirsutum varieties at 20 DPA but showed no difference among the 3 G. barbadense varieties [39]. We also found that the expression of GH_A05G0653 in xin w 139 at 20 DPA was greater than that in ** signaling pathway, the pentose phosphate pathway. The pentose phosphate pathway produces NADPH, and NADPH is involved in fatty acid synthesis [47]. The expression of the ethylene synthesis gene GhACO1 was upregulated after the addition of very long-chain fatty acids (VLCFAs) in vitro, which increased the ethylene content in cotton fibers, thereby promoting the elongation of cotton fibers [48]. A recent study showed that BRs regulate the synthesis of GhKCS-mediated VLCFAs through GhBES1, promoting fiber elongation [49]. We found that the DEGs in xin w 139 were annotated to the biosynthesis of the unsaturated fatty acids pathway, suggesting that the difference in fiber length between the two lines may be regulated by the expression of genes involved in the metabolic pathways of unsaturated fatty acids. VLCFAs, via enhanced UDP-l-rhamnose and UDP-d-galacturonic acid biosynthesis, can also increase the elongation of cotton fibers, which may also regulate fiber elongation through the sphingolipid biosynthesis pathway [1]. However, the contribution of VLCFAs to fiber strength has not been reported. Due to the significant difference in fiber length and fiber strength between the two lines, we suspect that the pentose phosphate pathway contributes to both the lightness and length of the fibers; this hypothesis still needs to be verified, but our results can provide important information for subsequent studies.

4. Materials and Methods

4.1. Plant Materials

In this study, a fiber development mutant was found for the line **n W 139, which was selected for breeding at the Institute of Economic Crops of **njiang Academy of Agricultural Sciences in 2018 and named xin w 139 (see Figure 1a,b for the phenotype and fiber analysis of **n W 139 and xin w 139 after flocculation in 2022). The agronomic traits (plant height and number of fruiting branches), yield traits (number of bolls per plant, weight per boll and coating) and fiber quality traits (fiber length and fiber strength) of the **n W 139 and xin w 139 lines were determined in Manas County, Changji city, **njiang, in 2022 and in Toutai township, Wusu city, **njiang, in 2023. For the field experiment, 1 row of each cultivar was planted; each plot was 5 m long, the row spacing was 0.35 m, the plant spacing was 0.10 m, and the management measures were the same as those used in the local conventional field. The cotton bolls on the day of flowering were labeled 0 DPA and they were sampled at 11 a.m. at 0, 5, 10, 15, 20, 25 and 30 DPA. During the sampling, the cotton husk was quickly peeled off within 1 min of peaching, and the fibers were removed with tweezers and immediately placed in liquid nitrogen for cryopreservation (6 replicates of each sample, 3 for RNA-seq and 3 for qRT–PCR).

4.2. RNA Extraction, cDNA Library Preparation, and Sequencing

RNA extraction was performed using the TRIzol method, and the quality of the extracted RNA was assessed via 1% agarose gel electrophoresis [50]. The extracted total RNA was subsequently transported to the Maiwei Metabolism Company (Wuhan, China) on dry ice for sequencing, after which the extracted RNA was fragmented using a PCR plate with a magnetic plate holder. Reverse transcription of the fragmented mRNA to cDNA was performed using Superscript II and random primers (Invitrogen, Carlsbad, CA, USA). The RNA-seq library preparations were sequenced on an Illumina (San Diego, CA, USA) HiSeq 2500/X platform, and 150 bp paired-end reads were generated. The library fragments were purified with an AMPure XP system (Beckman Coulter, Beverly, MA, USA). Then, 3 µL of USER Enzyme (NEB, Ipswich, MA, USA) was incubated with size-selected, adaptor-ligated cDNA at 37 °C for 15 min, followed by 5 min at 95 °C, before PCR. Then, PCR was performed with Phusion High-Fidelity DNA polymerase, universal PCR primers and Index (X) Primer. Finally, the PCR products were purified (AMPure XP system), and the library quality was assessed on an Agilent Bioanalyzer 2100 system. Fastp software (version 0.23.4) was used to remove the adapter sequences and filter out low-mass reads and reads with more than 5% poly-N sequences to obtain clean reads that could be used for the subsequent analysis [51]. The genome (https://www.cottongen.org/species/Gossypium_hirsutum/ZJU-AD1_v2.1, accessed on 3 January 2024) of upland cotton TM-1 was used as a reference, HISAT2 was used for the read alignment, and String Tie was used to quantify the reads in the alignment [52,53].

4.3. Analysis of DEGs

The fragments per kilobase of exon per million fragments mapped (FPKM) is a measure that quantifies gene expression levels. The number of reads per million that aligned to exonic regions was calculated and normalized by the length of the exonic regions and the total number of mapped reads. The gene expression levels in this study were determined using the FPKM method. The fold change in gene expression was calculated using EdgeR 4.0 software based on the number of clean reads obtained from the gene alignment [54]. An FDR ≦ 0.01 and a |log2-fold change| ≧ 1 were used as the standards for screening differentially expressed genes (DEGs) [55]. KEGG is a database resource for understanding high-level functions and utilities of biological systems, such as cells, organisms and ecosystems, from molecular-level information, especially large-scale molecular datasets generated by genome sequencing and other high-throughput experimental technologies (http://www.genome.jp/kegg/, accessed on 5 January 2024). We used KOBAS software 3.0.3 to test the statistical enrichment of differentially expressed genes in the KEGG pathways. The DEG sequences were submitted to PlantTFDB (http://planttfdb.cbi.pku.edu.cn/, accessed on 5 January 2024) for TF prediction.

4.4. Construction of Co-Expression Networks

The expression profile of the DEGs was determined by dynamic branching cleavage using the R language WGCNA package; the weighting coefficient was close to 0.8 β, the correlation coefficient requirement was met, and β = 9 was selected as the weighting coefficient in this study [56]. The automatic network builder Blockwise Modules was used to construct the network to obtain gene co-expression modules. The number of genes contained in each module was unequal, and the modules with a similarity of 0.75 were combined with minModuleSize = 30 and Merge Cut Height = 0.25 as the standards. The correlation coefficient between the characteristic vector module eigengene (ME) of the module and the different time points of fiber development of the two lines was calculated. Visualization of the co-expression networks was performed using Cytoscape (version 3.10.1) software [57].

4.5. SNP/Indel Analysis

Based on the HISAT2 alignment of the reads of each sample to the reference genome sequence, the SNPs/indels were identified using GATK (version 3.2-2) software [58]. The GATK identification criteria were as follows: (1) no more than 3 consecutive single-base mismatches were present in the 35 bp range and (2) the SNP quality values that were standardized by a sequence depth greater than 2.0. SNP/Indel loci were annotated using SnpEff software (version 3.6) to annotate regions of the genome (gene upstream, downstream, exon, or intron regions) based on the location of the variant locus on the reference genome as well as the gene location information on the reference genome [59].

4.6. qRT–PCR

The total RNA was extracted using the RNAprep Pure Polysaccharide Polyphenol Plant Total RNA Isolation Kit from Tiangen (Bei**g, China). The concentration of each RNA sample was measured using a NanoDrop 2000 spectrophotometer from Thermo Fisher Scientific (Waltham, MA, USA). cDNA was synthesized through reverse transcription of the RNA using the M-MLV RTase cDNA Synthesis Kit from TaKaRa (Kyoto, Japan). qRT–PCR analysis was conducted on a Bio-Rad CFX96 real-time PCR system from Mannheim Roche Diagnostics GmbH (Mannheim, Germany), with iTaq Universal SYBR Green Supermix from Takara Bio, Inc., in a reaction volume of 20 μL. For the qRT–PCR analysis, the following reaction program was used: an initial pre-denaturation step at 95 °C for 30 s, followed by 40 cycles of denaturation at 95 °C for 5 s, annealing at 60 °C for 5 s, and extension at 72 °C for 34 s. The relative quantitative analysis of the qRT–PCR results was performed using the 2^−ΔΔCt method [60]. The Ct value of the internal reference gene was subtracted from the Ct value of the target gene to obtain the ΔCt. The mean ΔCt value of the control group (0 DPA) was subtracted from each ΔCt of the treatment group to obtain the ΔΔCt. The final expression level was calculated via the 2^−ΔΔCt method. The internal reference gene used was GhUBQ7, and each experiment was conducted with three biological replicates. All the primers utilized in this study can be found in Table S1.

5. Conclusions

In conclusion, we used RNA-seq data from seven time points for natural mutant and wild-type lines to provide a reliable dataset for studying cotton fiber development. We not only identified 20 DPA as the key period for FS and length development but also identified several important regulatory pathways involved in fiber development by identifying the DEGs and TFs between lines. Cluster analysis was performed to divide the DEGs into clusters that could describe different stages of fiber development. In addition, seven candidate genes related to cotton fiber development, including three TFs, were screened by WGCNA, sequence comparative analysis and qRT–PCR. However, the exact role of these genes and markers in upland cotton fiber development has yet to be determined. Our results provide a theoretical basis for obtaining an in-depth understanding of the molecular mechanism of cotton fiber development and provide new genetic resources for cotton fiber research.