Comparative transcriptome analysis of cotton fiber development of Upland cotton (Gossypium hirsutum) and Chromosome Segment Substitution Lines from G. hirsutum × G. barbadense

Background How to develop new cotton varieties possessing high yield traits of Upland cotton and superior fiber quality traits of Sea Island cotton remains a key task for cotton breeders and researchers. While multiple attempts bring in little significant progresses, the development of Chromosome Segment Substitution Lines (CSSLs) from Gossypium barbadense in G. hirsutum background provided ideal materials for aforementioned breeding purposes in upland cotton improvement. Based on the excellent fiber performance and relatively clear chromosome substitution segments information identified by Simple Sequence Repeat (SSR) markers, two CSSLs, MBI9915 and MBI9749, together with the recurrent parent CCRI36 were chosen to conduct transcriptome sequencing during the development stages of fiber elongation and Secondary Cell Wall (SCW) synthesis (from 10DPA and 28DPA), aiming at revealing the mechanism of fiber development and the potential contribution of chromosome substitution segments from Sea Island cotton to fiber development of Upland cotton. Results In total, 15 RNA-seq libraries were constructed and sequenced separately, generating 705.433 million clean reads with mean GC content of 45.13% and average Q30 of 90.26%. Through multiple comparisons between libraries, 1801 differentially expressed genes (DEGs) were identified, of which the 902 up-regulated DEGs were mainly involved in cell wall organization and response to oxidative stress and auxin, while the 898 down-regulated ones participated in translation, regulation of transcription, DNA-templated and cytoplasmic translation based on GO annotation and KEGG enrichment analysis. Subsequently, STEM software was performed to explicate the temporal expression pattern of DEGs. Two peroxidases and four flavonoid pathway-related genes were identified in the “oxidation-reduction process”, which could play a role in fiber development and quality formation. Finally, the reliability of RNA-seq data was validated by quantitative real-time PCR of randomly selected 20 genes. Conclusions The present report focuses on the similarities and differences of transcriptome profiles between the two CSSLs and the recurrent parent CCRI36 and provides novel insights into the molecular mechanism of fiber development, and into further exploration of the feasible contribution of G. barbadense substitution segments to fiber quality formation, which will lay solid foundation for simultaneously improving fiber yield and quality of upland cotton through CSSLs. Electronic supplementary material The online version of this article (10.1186/s12864-017-4077-8) contains supplementary material, which is available to authorized users.


Background
On the global scale, cotton is not only one of the major crops, but also the most important plant producing natural fibers for the textile industry. Being composed of 46 diploid (2n = 2× = 26) and 5 allotetraploid (2n = 4× = 52) species [1], only 4 species of Gossypium were cultivated worldwide, namely G. herbaceum, G. arboreum, G. hirsutum and G. barbadense. The allotetraploid cotton species, Upland cotton (G. hirsutum) and Sea Island cotton (G. barbadense), derived from a natural hybridization event between A-genome and Dgenome 1-2 million years ago [2], contribute over 95% of cotton fiber yield. Upland cotton possesses the characteristics of high yield and moderate fiber quality, while Sea Island cotton is famous as its premium fiber quality performance and low fiber productivity. Facing the diminishing of arable land while increasing of human population, to develop novel varieties combining high productivity of G. hirsutum and excellent fiber quality of G. barbadense remains a preferential alternative in cotton breeding [3]. Chromosome segment Substitution lines (CSSLs), of which the genome was composed mostly of the recipient parent and minimal chromosome segments substituted from the donor parent, provided useful tools to take full advantages of both Upland and Sea Island cotton through marker assisted-selection (MAS) and conventional breeding procedures such as hybridization, backcross and selfing. To provide ideal materials for further genome research and crop improvement through MAS, CSSLs have been intensively applied to Quantitative Trait Locus (QTL) mapping for yield, quality, disease resistance and stress tolerance in a plenty of crops such as cotton [4][5][6][7][8][9][10][11][12][13], tomato (Lycopersicon esculentum), rice (Oryza sativa) and wheat (Triticum aestivum) [14][15][16][17].
Cotton fiber, a single-celled seed trichome developed from ovule epidermal cells, is an ideal model for investigating the mechanism of cell elongation, cell wall, and cellulose biosynthesis, which undergoes four overlapping developmental stages: initiation (−3 to +3 days post anthesis, DPA), elongation (3 to 23 DPA), secondary cell wall synthesis (16 to 40 DPA) and maturation (40 to 50 DPA) [18][19][20][21][22]. Fiber quality traits were collectively determined by performance of each of the aforementioned stages, so any disorders in these stages will significantly impact the final formation of the quality traits of fiber length, micronaire, strength, elongation and uniformity [23]. Along with the rapid development of highthroughput sequencing technologies, the genome sequencing of G. raimondii [24], G. arboretum [25], G. hirsutum [26,27] and G. barbadense [28,29] have been reported, which doubtlessly enhanced the cognition of polyploid properties of Gossypium species, including size variation of functional genome and significance of agronomic traits, and would further promote the research activities of structural genetics and functional genomics. As an effective tool to reveal the molecular mechanism of particular biological processes through concentrating on the overall level of gene expression and regulation, transcriptome analysis provided a novel high-throughput method for comprehending fiber growth system in combination with gene expression comparison between the different cotton species or between diverse developmental stages [30][31][32][33][34][35]. As a complicated biological process, fiber development was regulated by plant hormones such as auxin [36,37], gibberellins [38], brassinosteroid [39], ethylene [40,41], abscisic acid [42,43] and cytokinin [44,45]. Furthermore, carbohydrate, lipid, lignin, flavonoid and phenylalanine metabolisms have also been proved to play significant roles in fiber development [46][47][48][49][50] and some related genes are identified, such as GhKCS13 [51], GbEXPA2 [52], GbEXPATR [53] and so on. Despite all the above-mentioned results, the molecular mechanism of fiber elongation and Secondary Cell Wall (SCW) formation remains unclear, making it necessary for further investigations.
In the present study, two CSSLs, MBI9915 and MBI9749, having some substituted G. barbadense chromosome segments, were selected to conduct transcriptome analysis owing to their excellent fiber performance, and the recurrent G. hirsutum parent CCRI36 was also sequenced as the genome background and low fiber quality control. The CSSLs were developed from several generations of backcrossing and selfing of a cross between CCRI36, a G. hirsutum cultivar, and Hai 1, a G. barbadense cultivar, with G. hirsutum as the recurrent parent. Two of the most important fiber properties, fiber length and strength, were formed largely on the process of fiber elongation and SCW biosynthesis, respectively [54], therefore, the developing fibers were sampled at 10, 15, 20, 25 and 28 DPA for transcriptome sequencing and analysis. Aiming at exploring the differentially expressed genes (DEGs) relevant to fiber growth and inquiring into the possible contribution of chromosome substitution segments from Sea Island cotton to fiber quality formation, multiple comparisons among the RNA-seq data from different fiber developing stages of the three materials (two CSSLs and CCRI36) were performed. Plenty of DEGs were identified based on GFOLD analysis. The DEGs were then underwent GO functional enrichment and KEGG metabolic pathway analysis followed by accuracy validation using quantitative real-time PCR (qRT-PCR). The results indicated that the candidate genes bound up with fiber formation would contribute to the comprehension of the mechanism of cotton fiber elongation and SCW synthesis, and further make a significant difference in the cotton breeding and genomics research.

Results and Discussion
The phenotypic traits and substitution background of CSSLs In the study, the descriptive statistic characteristics of fiber yield and quality traits of MBI9915, MBI9749, the recurrent parent CCRI36, and the donor parent Hai 1 were shown in Table 1. With regard to the fiber quality performance, fiber length and strength of the two CSSLs were superior to those of CCRI36. In contrast, the boll weight and lint percentage of the two CSSLs were better than those of Hai 1, indicating that the two CSSLs showed a significant improvement in fiber yield as compared to the donor parent. The two CSSLs, which had significant fiber quality performance, would provide excellent materials for the further research of fiber development and functional gene identification.
In the high-density linkage map of CSSLs population derived from CCRI36 × Hai 1 which was described previously [12], 527 Simple Sequence Repeat (SSR) markers were selected at a marker interval of approximately one SSR marker per 10 cM genetic distance to identify chromosome substitution segments in the two CSSLs. The results showed that there were 19 substituted segments from G. barbadense and 97.90% genomic background recovery to G. hirsutum in MBI9915, while 19 substituted segments and 96.10% background recovery in MBI9749 ( Table 2). Heterozygocity of the substituted segments in the two CSSLs were analyzed (Table 2) and the results indicated that there was less than 4% of the total substituted segments of G. barbadense in the two CSSLs were identified, which indicated that the CSSLs were eligible for our intended studying of transcriptome analysis for functional dissection of the substituted segments on fiber quality formation.

Transcriptome sequencing and alignment to the reference genome
With the purpose of systematic identification of the critical genes affecting the cotton fiber development, 15 RNA-seq libraries of fiber samples were constructed and analyzed, which included the ones of fiber samples at 10, 15, 20, 25 and 28 DPA in MBI9915, MBI9749, and CCRI36. After removal of low-quality reads, a total of 750.433 million clean reads (approximately 88.18 Gb data) were obtained, of which the average reads per sample were 47.029 million (Table 3). In the alignment analysis, the proportion of clean data that were mapped to the genome of G. hirsutum reached to 91.93-94.30%, together with not less than 89.02% of the Q30 and over 44.79% of the GC content, implied a reliable quality of the RNA-seq result.
A total of 10,198 genes were identified, which were expressed in no less than half of the 15 libraries. So as to test the correlations between the experimental samples, Pearson Correlation Coefficient (PCC) analysis was adopted to analyze the expression of the genes obtained from these 15 libraries (Fig. 1). The results indicated that gene expressions in all the three lines at the same stage of fiber development showed more than 90% similarities except for the genes identified between MBI9749 and MBI9915/CCRI36 at 28DPA. Such gene expression similarities were higher at both the early (10 and 15DPA) and later (25 and 28DPA) stages than at the stage of 20 DPA, which hinted a significant gene expression alteration at the stage of 20DPA during fiber development. A gradual decreasing trend was also observed in similarities between 10DPA and each of the remained later stages in all the three lines, suggesting a remarkable alteration of gene expression profiles along with the stages of fiber development. All the findings of the gene For further confirming the relationship of the three lines at different stages, Principal Component Analysis (PCA) was performed on the above-mentioned expressed genes (Fig. 2). The results indicated that the expressed genes at the same stage in the three lines could be aggregated except for MBI9749 at 20DPA and at 28DPA. The similar gene expression patterns during the same development stage the PCA analysis indirectly proved, to some extent, the reliability of our RNA-seq data.

Analysis of differentially expressed genes (DEGs)
In order to identify DEGs in fiber development, the fragments per kb per million of the mapped reads (FPKM) value was employed to measure the gene expression quantity. The FPKM values were underwent GFOLD (Generalized Fold Change) algorithm to identify DEGs between MBI9915/MBI9749 and CCRI36 at same fiber development stage. Totally 1801 DEGs were obtained, including 898 down-regulated ones and 903 upregulated ones (Additional file 1). A cluster analysis of the FPKM values of these DEGs in the three lines was carried out, which resulted in a dynamic change of DEGs as shown in a heat map (Fig. 3). Interestingly, a high similarity of expression pattern of DEGs in CCRI36 and the two CSSLs at 10DPA was identified, which was in consistency with the above-mentioned PCC results. From an overall point of view, MBI9915 and MBI9749 had a higher similarity in gene expression profiles which was in contrast to that of CCRI36, indicating that the CSSLs presented a certain degree of difference with respect to their recurrent parent. The dynamically varied DEGs identified from the comparisons among the three lines might reveal a key gene expression regulating mechanism in fiber development and quality formation.
To comprehend the fiber development processes and identify the putative functional genes, the 1801 DEGs were classified into 37 GO terms based on biological process, cellular component and molecular function (Fig. 4). In the molecular function, most DEGs were sorted into catalytic activity (485 genes, 26.91% of the 1801 DEGs) and binding (478 genes, 26.53%). In the cellular component category, the most abundant subcategories were cell (223 genes, 12.38%) and organelle (186 genes, 10.32%). Under the biological process category, metabolic process (520 genes, 28.86%), cellular process    (340 genes, 18.87%) and single-organism process (307 genes, 17.04%) were the most abundant subcategories, while only one gene was categorized into growth process, namely XLOC_018364, which belongs to COBRA-like extracellular glycosyl-phosphatidyl inositolanchored protein family. According to previous studies, COBRA-Like (COBL) genes encoding a plant-specific glycosylphosphatidylinositol (GPI) anchoring protein have been found to play pivotal roles in the orientation of microfibrils and cellulose crystallinity status, and ultimately make a difference on fiber development and quality formation [55,56]. Referring to the expression quantity in the current RNA-seq data, XLOC_018364 was found to express throughout all the stages in this study (from 10 to 28DPA) in the three lines in spite of a lower expression level at early stages (10 and 15DPA).
While the FPKM value of XLOC_018364 in CCRI36 gradually decreased after reaching the peak at 20DPA, but the maximum FPKM value in MBI9915 and MBI9749 appeared at 28DPA and 25DPA, respectively. Considering the apparent differences between the CSSLs and their recurrent parent, the high activity of COBL gene in the two CSSLs at later developmental stages, might be due to the G. barbadense chromosomal segments harbored in them. For further uncovering the potential functions of these DEGs, GO enrichment and KEGG pathway analysis were respectively performed on the up-regulated and downregulated genes with Blast2Go software. The 903 upregulated genes were annotated and categorized into 41 GO terms based on biological process, cellular component and molecular function (Fig. 5). Under the biological processes, most up-regulated DEGs were concentrated on the response to oxidative stress (44 genes, 4.88% of the up-regulated 903 genes), response to auxin (39 genes, 4.32%), cell wall organization (37 genes, 4.10%) and lignin biosynthetic process (35 genes, 3.88%). In the cellular component category, mitochondrion (76 genes, 8.43%) was the most abundant subcategory followed by anchored component of membrane (27 genes, 2.99%) and anchored component of plasma membrane (24 genes, 2.66%). While in molecular function category, oxidoreductase activity and oxidizing metal ions (15 genes, 1.66%), oxygen binding (12 genes, 1.33%) and amino acid transmembrane transporter activity (10 genes, 1.11%) were the principal subcategories. The 898 down-regulated genes were annotated and categorized into 28 GO terms based on biological process, cellular component and molecular function (Fig. 6). The translation (58 genes, 6.46% of the 898 down-regulated genes), regulation of transcription and DNA-templated (34 genes, 3.79%) and cytoplasmic translation (33 genes, 3.67%) enriched most of the down-regulated DEGs under the biological process, while the dominant subcategories in cellular component and molecular function were cytosolic small ribosomal subunit (35 genes, 3.90%), followed by respiratory chain complexI (12 genes, 1.34%) and protein binding (118 genes, 13.14%), structural constituent of ribosome (98 genes, 10.91%), respectively.
To identify the DEGs between CCRI36 and two CSSLs, multiple comparisons among the different fiber developmental stages were performed as shown in Fig. 7. The results indicated that various sets of DEGs were identified between each of the CSSLs and CCRI36. Based on these results, the DEGs of the two CSSLs at each stage were compared and the results were shown in Fig. 8.

Analysis of gene temporal expression patterns
For the sake of exploring the temporal patterns of these DEGs in the two CSSLs and CCRI36, the Short Timeseries Expression Miner software (STEM) was employed. The results indicated that the expression of these DEGs in all the three lines was assorted into 8 profiles by STEM (Additional file 2, Fig. 9). Each profile represents a set of genes with a similar expression pattern. The largest expression profile of the DEGs in CCRI36 was profile 1 (598 genes, 41.50% of 1441 genes), followed by profile 7 (400 genes, 27.76%), profile 2 (208 genes, 14.43%) and profile 6 (89 genes, 6.18%). However, a similar expression trend was identified in the two CSSLs. In MBI9915, the highest number of DEGs was clustered into profile 1 (685 genes, 50.22% of 1364 genes), followed by profile 7 (282 genes, 20.67%), profile 6 (111 genes, 8.14%) and profile 4 (100 genes, 7.33%). While in MBI9749, the highest number of DEGs was clustered into profile 1 (692 genes, 49.96% of 1385 genes), followed by profile 6 (207 genes, 14.95%), profile 7 (190 genes, 13.72%) and profile 4 (71 genes, 5.13%).
The DEGs in profile 1 from the three lines were selected to perform a GO-term analysis to try identifying putative functional genes. All the DEGs in profile 1 were classified into 3 categories including biological process, cellular component and molecular function (Fig. 10). In the biological process, metabolic process (200 genes), cellular process (150 genes) and single-organism process (89 genes) were the dominant subcategories in CCRI36. While same trend was identified in the two CSSLs; in MBI9915, 232, 173, 97 DEGs were identified in metabolic process, cellular process and single-organism process respectively and in MBI9749, the same processes were identified to have 241, 174, 107 DEGs, respectively. In the molecular function category, the catalytic activity was the most abundant subcategory in the three lines, while cell was the most abundant subcategory in the cellular component. Based on the above results, the same category/subcategory was identified to have a varied number of DEGs implied that the different gene expression profiles might be responsible for the difference of fiber quality between CCRI36 and the two CSSLs.
Instead of above-mentioned processes, there had been increasing evidence that ROS played pivotal roles in fiber   [57,58]. As the second messenger during the biological processes, reactive oxygen species (ROS) has been reported to mediate cellular signal transduction and homeostasis regulation by forming the unique biological oxidation and reduction cycle [59], specifically participate in accommodating stress resistance and tip growth including pollen tube growth and root hair elongation in plants [60][61][62][63][64]. Therefore, a total of 29 DEGs in the subcategory "oxidation-reduction process" in biological process were selected to make another cluster analysis based on their FPKM values as shown a heat map (Fig. 11). An expression trend of gradual declining from 10DPA to 28DPA in CCRI36 were identified, while similar tendencies were also observed in MBI9915 and MBI9749 with a timing delay. Currently, 2 peroxidases (XLOC_036121 and XLOC_010767) were identified with high-expression quantity at 10DPA in all the three materials, which is consistent with a previous report that a large accumulation of ROS occurred at the cotton fiber initiation and elongation stage [65]. However, higher expression levels of the two genes were identified in two CSSLs than in CCRI36 at 10DPA, which indicated that fiber cell elongation obtained more enhancements in these genes in the CSSLs, which probably might result in the regulation of the quality formation in the end.
Furthermore, 2 leucoanthocyanidin dioxygenases (XLOC_019010 and XLOC_053890) and 2 flavanone 3hydroxylases (XLOC_029725 and XLOC_065014) associated with flavonoid pathway were identified in the aforesaid "oxidation-reduction process" subcategory. With strong antioxidant activities, flavonoids, which were relevant to the formation of pigment, participated in plant color alteration and eliminated the reactive oxygen species generated and accumulated under adverse stresses, mainly participated in plant resistance reactions [66] and in the natural coloration of cotton fiber [67].
Some studies reported that flavonoid pathway and its associated genes exhibited remarkable activities in the fiber development of wild and color cotton species [68][69][70][71]. Several flavonoid relevant genes in Upland and Sea Island cotton, such as F3H, CHS, ANS and LDOX were found with high level expression during the fiber growth, especially at the early stages of fiber development. Therefore, they were also believed to take part in the control of plant development [72,73]. Researches indicated that flavonoid pathway might mediate the cotton fiber development with a negative selection in the process of artificial domestication [74,75]. The four flavonoid genes identified in our RNA-seq results arised with a high-level expression at 10DPA, which is consistent with the above-mentioned research findings, indicating flavonoid pathway might participate in the fiber development, especially at the early stages.

DEG validation by qRT-PCR
After GO term analysis and KEGG annotation, twenty DEGs, including the ones involved in Starch and sucrose metabolism (XLOC_073805, XLOC_051736), Flavonoid biosynthesis (XLOC_070258, XLOC_053137) and Pentose and glucuronate interconversions (XLOC_055253), were randomly selected to confirm the validity of RNAseq data ( Table 4). The house-keeping gene β-Actin was selected as the reference gene during the validation analysis. The results indicated that expression patterns of 17 genes (85% of 20 genes) were highly consistent with the RNA-seq data, while the remained 3 genes (15%) were partially consistent (Table 4, Fig. 12). These results highly confirmed the reliability of the RNA-seq results.

Conclusion
Multiple comparisons among MBI9915, MBI9749, and CCRI36 during various fiber developmental stages The results showed that less than 4% of the total substituted segments of G. barbadense in the two CSSLs were heterozygous, which indicated the eligibility of the CSSLs for our intended studying of transcriptome analysis. A high proportion of 91.93-94.30% of clean RNAseq reads were mapped to the reference genome of G. hirsutum implied a reliable quality of the sequencing result. A total of 10,198 genes were identified and 1801 DEGs were identified through GFOLD analysis, including 898 down-regulated ones and 903 up-regulated ones. The up-regulated DEGs were mainly enriched in response to oxidative stress and auxin, cell wall organization, and lignin biosynthesis processes while the down-regulated ones in the processes of translation, regulation of transcription, DNA-templated, and cytoplasmic translation. Temporal expression patterns analysis of DEGs revealed that they were mainly enriched in profile 1, profile 7 and profile 6. GO-term analysis of DEGs in profile 1 indicated that these DEGs were mainly enriched in subcategories of metabolic process, cellular process, and single-organism process in biological process in all the three lines. Referring to the previous studies and bioinformatics analysis of the genes identified in the oxidation reduction process subcategory revealed that these genes also played important roles in fiber development in CSSLs in our current study. Our study not only offered novel insights into the developmental mechanism of fiber elongation and SCW biosynthesis, but also discussed the potential contribution of substitution segments to fiber quality performance, which undoubtedly is of great value for cotton breeding and genomics researches.

Plant material
Two Chromosome Segment Substitution Lines (CSSLs) and their G. hirsutum recurrent parent CCRI36 planted in the experimental farm of the Institute of Cotton Research, Chinese Academy of Agricultural Sciences (Anyang, Henan Province) were used as plant materials in the current study in 2014. The CSSLs were constructed from a cross between the paternal parent Hai 1 (G. barbadense) and the maternal parent CCRI36 (G. hirsutum) grown in Anyang in 2003. After five generations of backcrossing with CCRI36 as the recurrent parent and three generations of consecutive selfing, the BC5F3 CSSLs were obtained in 2009. Through plant-torow method, BC5F3:4 rows were planted in Anyang (Henan province) in 2010 for fiber quality evaluation. In 2011, BC5F3:5 lines were planted in Anyang (Henan province), Liaoyang (Liaoning province) and Shihezi (Xinjing province) for further evaluation. The details of CSSLs construction were shown in Fig. 13.
Based on the high-density linkage map constructed from the aforementioned population of CSSLs [12], 527 Simple Sequence Repeat (SSR) markers were selected to confirm the genetic compositions of CSSLs   at an approximate interval of one SSR maker per 10 cM [11]. Owing to the superior fiber performance and relatively clear substituted chromosomal segments composition, MBI9915 and MBI9749 were selected to conduct the current transcriptome study, with their recurrent parent CCRI36 was used as a reference. Due to concentrating on studying fiber length and strength, our emphasis was laid on the fiber elongation and SCW biosynthesis, which made us to sample the developing fibers at 10, 15, 20, 25 and 28DPA to conduct RNA-seq.
Flowers were tagged on the day of anthesis, which were carried out consecutively for three days as one biological replicate. Three-ten developing bolls from the three lines were harvested by 10:00 a.m. at 10, 15, 20, 25 and 28 day post anthesis (DPA) in each replicate, then immediately immersed in ice. Fiber samples were dissected from the developing bolls, frozen into liquid  nitrogen, and stored at −80°C, which were divided into two portions, one for RNA-seq and the other for qRT-PCR. Totally three replicates were tagged and sampled.

RNA extraction, library construction and transcriptome sequencing
High-quality RNA extraction was performed from frozen fiber tissue using the RNAprep Pure Plant Kit (Tiangen, Beijing, China). RNA degradation and contamination was confirmed through 1% agarose gel electrophoresis. After being assessed the RNA integrity by the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, CA, USA), RNA concentration was calibrated by Qubit® RNA Assay Kit in Qubit® 2.0 Flurometer (Life Technologies, CA, USA). For minimizing experimental errors, after calibration, RNA samples extracted from the three biological replicates of each stage were mixed with same quantity and volume to generate one sample and were calibrated again for transcriptome sequencing in this study. A total amount of 3 μg RNA sample was used for transcriptome library sequencing, which was constructed using NEBNext® Ultra™ RNA Library rep Kit for Illumina® (NEB, USA) following manufacturer's recommendations. Firstly, mRNA purification and fragmentation were respectively carried out by poly-T oligo-attached magnetic beads and divalent cations under elevated temperature (94°C for 5 min) in NEBNext First Strand Synthesis Reaction Buffer (5×). Next, the first strand cDNA synthesis was accomplished using random hexamer primer and M-MuLV Reverse Transcriptase (RNase H-) and the second strand cDNA was synthesized using DNA Polymerase I and RNase H. After purification with AMPure XP system (Beckman Coulter, Beverly, USA), the cDNA segments of 150-200 bp in length were utilized to construct the cDNA libraries through end-repair, adaptor ligation and PCR amplification. At last, the AMPure XP system and Agilent Bioanalyzer 2100 system were separately used to purify the PCR products and to evaluate the library quality. In total, 15 cDNA libraries of 3 cotton lines with 5 stages (10,15,20,25,28DPA) were sequenced on the Illumina HiSeq™ 2000 sequencing platform.
Data quality control and reads mapping to the reference genome Through Base Calling, the raw sequence data were firstly transformed to Sequenced Reads, including reads sequences and corresponding base quality. The raw reads in FASTQ format were subsequently processed by the Illumina pipeline to obtain the clean data by filtering all low-quality reads, such as the reads containing only adapter or poly-N whilst the Q30 and GC contents were calculated. All the downstream analyses were based on the clean data with high quality.
The clean reads were mapped to the G. hirsutum reference genome [27], together with gene model annotation files, which was downloaded from the CottonGen database (http://www.cottongen.org). Bowtie v2.0.6 and TopHat v2.0.9 was respectively applied to build the index of the reference genome and to make the reads alignment [76].

Differentially expressed genes analysis
In order to comprehend the status of differential gene expression in the developmental processes of cotton fiber, 15 separate libraries were established from 5 fiber development stages (10,15,20,25,28DPA) in the 3 lines, thereby producing comparison groups of the different stages of the same line and the different lines of the same stage. Gene expression levels were subject to transcript abundance, which was calculated as the reads mapped to reference genes or exons. The FPKM (fragments per kb per million of the mapped reads) parameter was employed to count the transcript data [77], whilst the algorithm GFOLD (Generalized Fold Change) was adopted to produce biological meaning and the rankings of DEGs [78]. As the recurrent parent of MBI9915 and MBI9749, CCRI36 was used to compare the two CSSLs to generate GFOLD value. The genes that their absolute GFOLD value >1 was defined as DEGs. Furthermore, Gene Ontology (GO) enrichment analysis of DEGs was performed by Blast2go software with corrected P-value ≤ 0.5 as the cutoff. KEGG pathways of DEGs was tested by KOBAS software [79,80].

Comparison of expression patterns of DEGs
The Short Time-series Expression Miner software (Carnegie Mellon University, USA) was adopted to identify the temporal expression patterns of genes along with the growth stages of cotton fiber. GO enrichment analysis was also employed to identify the potential functional genes or to characterize the expression profiles identified by STEM analysis.

Validation of RNA-seq by qRT-PCR
qRT-PCR was employed to verify the veracity of the transcriptome data with twenty randomly selected DEGs. The specific primers of these DEGs were designed by Primer-BLAST of online tool of NCBI. RNA samples were prepared as described in "RNA extraction, library construction and Transcriptome sequencing" section. The cDNAs synthesized using a TranScript All-in-One First-Strand cDNA Synthesis SuperMix for qPCR Kit (Transgen Biotech, Beijing, China). qRT-PCR was performed following the protocol of TransStart Top Green qPCR SuperMix kit (Transgen Biotech, Beijing, China), on the ABI 7500 fast Real-Time PCR System (Applied Biosystems, USA). The housekeeping gene of β-Actin gene was used as the reference to normalize the relative expression levels, with its primer sequences: F: 5′-ATCCTCCGTCTTGACCTTG-3′, and R: 5′-TGTCC GTCAGGCAACTCAT-3′. The qRT-PCR carried out in 20 μL system at the following conditions: one cycle of 94°C for 30s; 40 cycles of 94°C for 5 s, 60°C for 34 s, and one cycle of 60°C for 60s. Three biological and technical replicates were performed to validate the results in qRT-PCR tests. The relative gene expression level was quantified by the 2 -ΔΔCt method [81].