Expression profile analysis of cotton fiber secondary cell wall thickening stage

Plant materials

Cotton cultivars, Xinhai 32 is a high-generation inbred line of Sea Island cotton and 17–24 and 62–33 are high-generation inbred lines of upland cotton. All three cultivars were bred by the Cotton Institute, Xinjiang Academy of Agricultural and Reclamation Science. Xinhai 32, 17–24, and 62–33 are referred to herein as HL, L, and S, respectively. These three cotton cultivars were planted in the field of the Xinjiang Academy of Agricultural and Reclamation Science. Ovules (seeds) were harvested at the indicated days post-anthesis (DPA). Cotton bolls of 20, 25, and 30 DPA from the three different cotton cultivars were sampled at 10:00 AM. Five cotton bolls were selected from each plant, and the fibers were isolated from the ovules by scratching the ovule with a metal strainer in liquid nitrogen. The fiber samples were then quickly ground into powder and stored in the ultra-low temperature refrigerator at −80 °C.

Fiber traits and phenotypic evaluation

Boll weight, 100 seed weight, lint index, and seed index were weighed by an analytical balance (0.0001 g, BSA224S, Sartorius, Göttingen, Germany). Fiber quality traits, including the fiber length (mm), fiber uniformity ratio (%), fiber strength (cN/tex), fiber elongation, and micronaire, were measured with an HVI 900 instrument (USTER HVISPECTRUM, SpinLab, Leipzig, Germany) at the Cotton Fiber Quality Inspection and Test Center of Ministry of Agriculture (Anyang, China) (Shang et al., 2015). Bolls were collected for seed index analysis. For this purpose, one hundred seeds from each cultivar were randomly selected and weighed as seed index (SI, g) (Liu et al., 2023; Shang et al., 2016). To measure oil content, seeds were delinted with concentrated sulphuric acid. The oil contents of the three cotton cultivars were measured at different stages of ovular development using the Soxhlet extraction method (García-Ayuso et al., 2000). To record grain weight, one hundred cotton ovules were randomly weighed.

RNA extraction, library construction, and sequencing

Total RNA from 20, 25, and 30 DPA cotton fiber from each of the three different cotton cultivars was extracted using the RNAprep pure plant kit (Tiangen, Beijing, China) according to the manufacturer’s instructions. A total of 1.0 μg purified mRNA was selected for cDNA library construction following a previous report (Chen et al., 2021). Briefly, mRNA was purified from total RNA using poly-T oligo-attached magnetic beads. The cDNA fragments of 240 bp length were selected, and the library fragments were purified using the AMPure XP system (Beckman Coulter, Beverly Hills, CA, USA). The PCR products were purified using the AMPure XP system, and the library quality was assessed on an Agilent Bioanalyzer 2100 system. After generating clusters, the libraries were sequenced using the Illumina HiSeqTM 2500 platform as 150 bp paired-end reads. Three biological replicates were performed for the nine samples.

RNA-seq reads quality control, mapping, and differentially expressed genes (DEGs) analysis

The FASTX-Toolkit (http://hannonlab.cshl.edu/fastx_toolkit/) was used to process the raw reads in fastq format according to a previous study (Chen et al., 2021). Clean data (clean reads) were obtained by removing reads containing adapters, reads containing poly-N sequences, and low-quality reads from the raw data. All subsequent analyses were performed based on clean high-quality data. Among the nine fiber samples from G. hirsutum and G. barbadense (Hu et al., 2019; Wang et al., 2019), RNA-seq data were mapped to their respective reference genomes using HISAT2 software (Kim, Langmead & Salzberg, 2015; Pertea et al., 2016). Reads with at most one mismatch were used to calculate the expression levels of genes. Gene expression values were calculated following the method of the previous study (Chen et al., 2022). Differential expression analysis of the two groups was performed using DESeq2 and presented using fragments per kilobase of transcript per million fragments mapped (FPKM) (Love, Huber & Anders, 2014). The resulting p values were adjusted using Benjamini and Hochberg’s approach for controlling the false discovery rate (FDR). Genes with an adjusted p value < 0.01 and two-fold change (up and down) were defined as differentially expressed. TBtools (Chen et al., 2020) was used to display the gene expression patterns of the FPKM values. Clean data were available from the Genome Sequence Archive in the BIG Data Center of Sciences (https://bigd.big.ac.cn/) under accession number CRA009299. The statistical power of this experimental design, calculated in RNASeqPower is 0.84.

Gene functional annotation and enrichment analyses

The functions of differentially expressed genes were annotated using the following databases: Nr (NCBI nonredundant protein sequences, ftp://ftp.ncbi.nih.gov/blast/db/), Gene Ontology (GO) (Gene Ontology, http://www.geneontology.org/), and Kyoto Encyclopedia of Genes and Genomes (KEGG) (http://www.genome.jp/kegg/). To analyze the enriched GO of the DEGs, we used the GOseq R package based on Wallenius noncentral hypergeometric distribution (Young et al., 2010). KOBAS was employed to assess the statistical enrichment of the DEGs in the KEGG pathways (Mao et al., 2005).

RNA extraction, cDNA synthesis, and RT-qPCR expression analyses

These experiments were conducted according to the methods reported previously (Cao et al., 2020a, 2022; Cui et al., 2022; Zhang et al., 2021a). In brief, total RNA from 20, 25, and 30 DPA fiber samples of the three different cotton cultivars were extracted using the RNAprep pure plant kit (Tiangen, Beijing, China). DNase I treatment was applied to the RNA samples before synthesizing cDNA using TransScript^® First-Strand cDNA Synthesis SuperMix from TransGen Biotech, China, and the resulting products were diluted fivefold before use. Specific forward and reverse gene primers (Table S1) were designed using Primer v5.0 software for real-time quantitative PCR, which was performed using SYBR-Green PCR Mastermix (TaKaRa, Shiga, Japan) on a cycler (Mastercycler RealPlex; Eppendorf Ltd., China). The G. hirsutum and G. barbadense histone-3 (GhHIS3 and GbHIS3) genes were used as internal references, and the relative amount of amplified product was calculated following the 2^−∆∆Ct method (Livak & Schmittgen, 2001).

Statistical analysis

The R package available at https://www.r-project.org/ was utilized for the analysis of variance and Student’s t-test. The Shapiro-Wilk test was utilized to test for normality, confirming that the data followed a Gaussian distribution. The least significant difference (LSD) was used to test for significance at either the 1% or 5% levels. The analysis included at least three biological replicates for each sample.

Physiological traits differences among the G. barbadense (Xinhai 32, HL) and G. hirsutum (17-24, L and 62-33, S) cotton cultivars

Physiological traits of the three cotton cultivars, all high-generation inbred lines developed by the Cotton Institute, Xinjiang Academy Agricultural and Reclamation Science, were characterized. Table 1 shows the main properties of these cultivars, which we have designated as S (cv 62–33), L (cv 17–24), and HL (Xinhai 32). In comparison to the upland cotton cultivars (L and S), the Sea Island cotton (HL) exhibited longer fiber, greater fiber strength, increased oil content, a higher seed index, greater fiber uniformity, and greater fiber elongation. In contrast, the upland cotton cultivars (L and S) had significantly higher boll weight, lint percentage, lint index, and micronaire than did the HL cultivar. Notably, these properties, although greater in the G. hirsutum cultivars, are outside of the optimal range. In particular, micronaire (one of the most important measures of cotton fiber quality) was considered grade A (micronaire range: 3.7–4.2) in HL, compared the B grade (micronaire range: 3.5–3.6 and 4.3–4.9) observed in the upland cotton cultivars, L and S. Between the two cultivars of upland cotton, fiber length, fiber strength, and oil content also varied, with cv 17–24 (L) exhibiting significantly higher values than cv 62–33 (S). These results indicate that major differences in agronomically important fiber properties exist between Sea Island cotton and upland cotton, as expected, but also between the two upland cotton cultivars, with cv 17–24 (L) exhibiting better agricultural performance.

Transcriptome data generation of 20, 25, and 30 DPA fibers of Xinhai 32 (HL), 17–24 (L), and 62–33 (S) cotton cultivars

As these three cultivars showed great variance in fiber length and strength, especially fiber strength (Table 1), we evaluated gene expression at three timepoints during SCW synthesis (i.e., 20, 25, and 30 DPA fiber) from Xinhai 32 (HL), 17–24 (L), and 62–33 (S). RNA from at least three biological replicates was pooled at each timepoint for each accession, hereafter referred to as HL20, HL25, HL30, L20, L25, L30, S20, S25, and S30, and over 108 million high-quality reads were generated using the Illumina HiSeqTM 2500 sequencing platform (Table 2). The HISAT2 software was used to align these clean reads to the reference cotton genomes (Hu et al., 2019; Wang et al., 2019) with at most one base mismatch. The ratio of mapped reads ranged from 72% in sample L30 to 77% in sample HL20 (Table 2), and the number of uniquely mapped clean reads ranged from 53% in sample L30 to 70% in sample L20. These data indicated that the RNA-seq data in this study were reliable for the subsequent analyses.

Identification of differentially expressed genes

Differential gene expression was surveyed for the nine cotton fiber samples. Eighteen comparisons among these nine samples were performed to capture expression differences between samples at an assessed timepoint and within a single sample at each timepoint (Fig. 1A). In each comparison, the number of differentially expressed genes (DEGs) varied between 761 genes in HL25-VS-L25 and 10,464 genes in HL20-VS-HL25 (Fig. 1A). It is worth noting that upregulation was more frequent than downregulation in most comparisons (Fig. 1A). Because these cultivars produce fibers with different properties, we considered the overlap in gene expression at each assessed timepoint to identify genes important for SCW biosynthesis that differed in expression among the three species/cultivars. Venn diagrams were constructed for the DEGs from each species/cultivar comparison at each surveyed stage (i.e., 20, 25, and 30 DPA stages; Fig. 1). As expected, the intraspecies cultivar comparisons (i.e., L vs S) typically exhibited fewer uniquely DEGs than the interspecies comparisons (i.e., HL vs S or L), with the exception of the 25 DPA stage of development (Fig. 1). Although at the 25 DPA timepoint fewer DEGs were identified (Fig. 1), the intraspecies (S vs L) comparison resulted in 25–40% more DEGs than either interspecific comparison (Fig. 1); the fewest DEGs at this stage were between the HL and L cultivars, perhaps indicating that underlying fiber length expression differences are largely responsible for the differences in gene expression at this stage. The results showed that 423 DEGs were shared by all species/cultivar comparisons at the 20 DPA stage (Fig. 1B), possibly indicating genes that underlie the differences among cultivars. Far fewer DEGs were observed at 25 DPA (39 DEGs), although this number is commensurate with the general reduction in DEGs among cultivars at this stage (Fig. 1C). At 30 DPA, 361 DEGs were shared among the cultivar comparisons (Fig. 1D).

GO analysis for DEGs

We considered the possible biological functions of these DEGs between Sea Island and/or the upland cottons using GO category enrichment (Fig. S1). Results from the three categories, including biological process, cellular component, and molecular function, suggest enrichment of 20, 11, and 11 functional categories in the 20, 25, and 30 DPA fiber comparisons, respectively. GO terms associated with important biological processes included metabolic, cellular, developmental, and single-organism processes, biological regulation, response to stimulus, and signaling. Cellular components, such as cell, cell part, membrane, membrane part, organelle and organelle parts were enriched. Molecular function enrichment consisted of catalytic activity, transporter activity, binding, nucleic acid-binding transcription factor activity, antioxidant activity, and receptor activity.

To compare the difference between the three cultivars, the GO enrichment of functional categories in the 25 DPA fiber was analyzed (Fig. 2). The results showed that the top three GO category enrichments in biological processes among the three groups of comparisons (HL-25-vs-L-25, HL-25-vs-S-25, L-25-vs-S-25) were cellular process, metabolic process, and localization, respectively (Figs. 2A–2C). Additionally, the cellular components, cell, cell part and membrane accounted for the largest proportion in all three groups (Figs. 2A–2C). The most enriched molecular functions all consisted of catalytic activity, and transporter activity between the comparisons of the three cultivars (Figs. 2A–2C). These results indicated that the GO enrichment analysis was not specific to the three cultivars.

Pathway enrichment of DEGs from the three cotton cultivars at different stages of fiber development (20, 25, and 30 DPA)

To investigate the biological functions of these DEGs during fiber development in the three cotton cultivars, we performed KEGG pathway enrichment analysis for the DEGs. At the 20 DPA fiber stage, the DEGs were assigned to 126 KEGG pathways according to the functional categorization. For the interspecific comparison HL-20-vs-L-20, the top 20 KEGG pathways of enriched DEGs were categorized into the following functional pathways (Fig. 3A). While the greatest number of DEGs mostly belonged to the plant hormone signal transduction pathways, there was also involvement of several other pathways whose RichFactor was closer to 50%, such as fatty acid elongation. Other pathways with notable enrichment are: i) carbohydrate metabolism: glycolysis/gluconeogenesis, pyruvate metabolism, galactose metabolism, fructose and mannose metabolism; ii) fatty acid metabolism: fatty acid elongation, and biosynthesis of unsaturated fatty acids; and iii) amino acid metabolism: valine, leucine and isoleucine degradation, lysine biosynthesis, and cysteine and methionine metabolism. For the other interspecific comparison, HL-20-vs-S-20, the top 20 KEGG pathways of enriched DEGs were involved in glycolysis/gluconeogenesis, galactose metabolism, fructose and mannose metabolism, fatty acid elongation, biosynthesis of unsaturated fatty acids, and leucine and isoleucine degradation (Fig. 3B); however, the greatest representation of genes was for RNA transport and the ribosome. In this comparison, the RichFactor varied more narrowly (~12%, vs 25% in HL20 vs L20). For the intraspecific comparison L-20-vs-S-20, the top 20 KEGG pathways of enriched DEGs mainly included the glycolysis/gluconeogenesis, galactose metabolism, fatty acid elongation, biosynthesis of unsaturated fatty acids, plant hormone signal transduction, and biosynthesis of secondary metabolites (Fig. 3C). The greatest numbers of genes were found in the metabolic pathways and biosynthesis of secondary metabolites. While the RichFactor for these was low, it is worth noting that the RichFactor was generally low for all pathways, possibly indicating the recruitment of few genes/pathways from these somewhat broad categories. These results suggested that the differences in the 20 DPA fiber samples between Sea Island cotton (Xinhai 32) and upland cotton (17–24 or 62–33) were mainly caused by genes involved in the metabolic pathways of carbohydrates, fatty acids, and amino acids. However, differences in the plant hormone pathways were common to both comparisons involving L-20, and differences in secondary metabolite production were observed between the two cultivars of upland cotton.