Transcriptome Analysis of Ten Days Post Anthesis Elongating Fiber in the Upland Cotton ( Gossypium hirsutum ) Chromosome Substitution Line CSB 25

A chromosome substitution line, CS-B25, was developed by the substitution of chromosome pair 25 of Gossypium hirsutum TM-1 with the homologous pair of chromosome 25 from G. barbadense, a double haploid Pima 3-79 line. CS-B25 has improved fiber traits compared to its parent TM-1. To explore the molecule mechanisms underlying improved fiber traits, deep sequencing of total RNA was used to compare gene expression in fibers of CS-B25 and TM-1 at 10 days post anthesis (10-DPA). A total of 1872 differentially expressed genes (DEGs) were detected between the two lines, with 1175 up-regulated and 697 down-regulated in CS-B25. Gene Ontology (GO) enrichment analysis of the expression data by Generally Applicable Gene-set Enrichment (GAGE) and ReviGO indicated that the most prevalent Biological Process GO terms associated with DEGs included DNA-templated transcription, response to oxidative stress, and cellulose biosynthesis. Enriched Molecular Function GO terms included structural constituents of cytoskeleton, peroxidase activity, cellulose synthase (UDP-forming) activity, and transcription regulatory region sequence-specific DNA binding factors. GAGE was also used to find enriched KEGG pathways, and the highly represented pathways were Biosynthesis of Amino Acids, Starch and Sucrose Metabolism, Phenylpropanoid Biosynthesis, Protein Processing in Endoplasmic Reticulum, and Plant Hormone Signal Transduction. Many of the identified DEGs are involved in cyHow to cite this paper: Hsu, C.-Y., Arick II, M.A., Miao, Q., Saha, S., Jenkins, J.N., Ayubov, M.S., Abdurakhmonov, I.Y., Peterson, D.G. and Ma, D.-P. (2018) Transcriptome Analysis of Ten Days Post Anthesis Elongating Fiber in the Upland Cotton (Gossypium hirsutum) Chromosome Substitution Line CS-B25. American Journal of Plant Sciences, 9, 1334-1361. https://doi.org/10.4236/ajps.2018.96098 Received: March 27, 2018 Accepted: May 27, 2018 Published: May 30, 2018 Copyright © 2018 by authors and Scientific Research Publishing Inc. This work is licensed under the Creative Commons Attribution International License (CC BY 4.0). http://creativecommons.org/licenses/by/4.0/ Open Access


Introduction
Cotton fibers are seed trichomes and differentiated from the epidermal cells of a developing seed.The development of cotton fiber consists of four partially overlapping stages: initiation, elongation (primary wall synthesis), secondary wall thickening, and maturation [1].Fiber cells initiate from 3 days before anthesis to 3 days post anthesis (DPA).After initiation, fiber cells immediately enter into the elongation stage.Fiber elongation reaches peaks at 6 -12 DPA, but can last up to 22 -26 DPA.The secondary wall synthesis begins at 16 -19 DPA, reaches a maximal rate at 25 DPA, and lasts up to 40 DPA.Fiber maturation, the last stage of fiber development, occurs at 40 -50 DPA [1].Fiber cell initiation and elongation are the two most important developmental stages with regard to fiber quality [2].
Cotton provides the largest renewable natural fiber for the textile industry and is also an important source of food, feed, fuel and other products [3].Over 90% of the world cotton production is from Upland cotton, Gossypium hirsutum L. [4].Upland cotton is widely adapted and has a high percentage of lint fiber and a relatively high yield.Pima cotton (G.barbadense L.), the second most widely-grown cotton species, generally has a lower lint percentage and lower yield, but its fiber qualities are superior to Upland cotton in terms of length, strength, and fineness.The Pima fiber offers superior spinning and manufacturing performance and better textile products, and thus normally commands a 30% to 50% price advantage over fiber from high-yielding Upland cottons.Pima cotton, however, requires long growing seasons to produce profitable yields of high quality fiber.One of the major goals of cotton breeders is to develop Upland cotton with Pima-like elite fiber quality combined with high yield and adaptability.Combining the adaptability and yield of Upland cotton with the fiber quality of G. barbadense into an Upland cotton cultivar would maximize quality and productivity in an increasingly competitive market.However, crosses between G. hirsutum and G. barbadense at the whole genome level have not been highly successful due to technical and biological challenges associated with the conventional methods of interspecific introgression.An alternative approach to interspecific introgression is to use alien chromosome substitution lines where individual chromosomes or chromosome arms of G. barbadense are substituted into G. hirsutum.Seventeen interspecific chromosome substitution lines (CS-B lines) of Upland cotton (TM-1 background) containing whole chromosomes or chromosome arms of G. barbadense (Pima 3-79) chromosomes have been developed and released to the public [5].These lines are genetically identical except that each differs by the replacement of a specific homologous pair of chromosomes from Pima 3-79 (G.barbadense) into the Upland background.CS-B25, a line in which the chromosome 25 pair of TM-1 has been replaced with the homologous chromosome 25 pair from G. barbadense, exhibits several improved fiber quality traits including increased fiber length, increased fiber strength, and lower micronaire [6] [7] [8] [9] [10].Results from crosses of CS-B25 with selected USA improved cultivars also demonstrated that CS-B25 had several beneficial alleles with additive gene effects for improving fiber strength and decreasing micronaire in Upland cultivars [9] [10].Additionally, CS-B25 has some beneficial alleles for improving agronomic traits with dominance genetic effects for lint percentage, seed cotton, and lint yield [11].
One of the serious impediments in cotton genetic improvement is the paucity of information of genes and their association with important fiber traits.The CS-B lines, being isogenic lines with the recurrent parent TM-1 except for the substituted chromosome pair from the alien species of G. barbadense, provide excellent analytical tools to dissect genetically complex fiber traits at the molecular level.The CS-B22sh line (chromosome 22 short arm substitution from Pima 3 -79 into TM-1) was first used as a tool to identify differentially expressed transcripts in cotton [12].The authors reported that 36 genes were differentially expressed in 15-DPA fiber in CS-B22sh, including genes that encoded tubulins, actin, cellulose synthase, glycosyl hydrolase, glycoside hydrolase, glycoside pyruvate decarboxylase and genes involved in signal transduction and nucleic acid, protein and lipid metabolism.Another report showed that CS-B22sh had positive additive effect on decreasing micronaire and a positive homozygous dominance effect on fiber strength [9].
In this study a genome wide transcriptome analysis was used to identify differentially expressed genes (DEGs) in 10-DPA fibers in CS-B25 as compared to TM-1.These DEGs and KEGG pathways analysis have provided insight into the mechanisms of fiber development during the fiber elongation stage and offered novel candidate genes that can be utilized in the genetic improvement of fiber quality.

Extraction and Purification of Fiber RNA
Total RNA was isolated from 10-DPA fibers of TM-1 and CS-B25 lines using a modified hot borate method [13].Genomic DNA contamination in the RNA samples was removed by DNase I (Promega, Madison, Wisconsin, USA).After DNase treatment, other potential contaminants were removed from the RNA sample with a Qiagen RNasey Mini Kit (Qiagen, Valencia, California, USA).The quality of the total RNA sample was pre-checked on a 1% -1.2% (w/v) agarose gel by electrophoresis, and the RNA concentration was estimated with a Nanodrop 1000 Spectrophotometer (NanoDrop Technologies, Wilmington, Delaware, USA) and a Qubit® Fluorometer by using the Qubit RNA BR Assay Kit (Invitrogen, Carlsbad, California, USA).The quality of the total RNA was finally determined using both agarose gel electrophoresis and a Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, California, USA).

RNA Library Preparation, Illumina Sequencing and Data Collection
The isolated total RNA samples from 10-DPA fibers (1.5 µg from each sample) were used to construct strand-specific cDNA libraries using the TruSeq Stranded mRNA Sample Preparation Kit (Illumina, San Diego, California, USA) according to the manufacturer's instructions.Three biological replicate samples of the two cotton lines CS-B25 and TM-1 were used to generate 6 RNA-Seq libraries which were pooled and sequenced on one lane of an Illumina HiSeq 2000 sequencer (pair end sequencing, 50 bp per end) (Illumina, San Diego, California, USA).The raw reads were mapped to the G. hirsutum (AD1) NAU-NBI genome [14] using Tophat 2 [15].Next HTSeq [16] was used to estimate gene expression from the Tophat alignments.The mRNA sequencing reads were deposited into NCBI Sequence Read Archive (SRA) with accession numbers SRR6015224-6 and SRR6015221-3 for TM-1 and CS-B25, respectively.

Identification and Functional Analysis of Differentially Expressed Genes
EdgeR [17], a Bioconductor (http://www.bioconductor.org/)software package for differential gene expression analysis, was used to identify DEGs in CS-B25 10-DPA fibers via comparison with TM-1 10-DPA fibers.Generalized linear models were used to test for differential gene expression between CS-B25 and TM-1.Genes with an FDR adjusted p-value ≤ 0.05 were considered differentially expressed.GAGE [18], a Bioconductor software package for differential gene set analysis, was used to identify differentially expressed GO terms [19] and KEGG pathways [20].GO gene sets were compiled from the functional annotation available from CottonGen [21].REViGO [22] was used to summarize and visualize the enriched GO terms.Gene expression data was mapped to Arabidopsis KEGG pathways using gene homologies reported on CottonGen.rRNA was used as a reference gene to normalize the expression levels of RT-PCR products [23].All reactions were performed with three replicates.The 2 −ΔΔCt method [ΔCt = Ct (differentially expressed gene) − Ct (18S rRNA), ΔΔCt = ΔCt (CS-B25) − ΔCt (TM-1), 2 −ΔΔCt = Relative Expression] was used to calculate relative expression of differentially expressed genes [24].All primers were designed by using Primer3 plus and the sequences of the primers are listed in Table S1.

Analysis of SSR DNA Markers Associated with Chromosome 25
Seven chromosome 25-specific nuclear SSR markers associated with improved fiber traits of fineness (micronaire), length, and strength were selected from the previous QTL mapping of fiber quality loci (see Table S2).These QTLs identify genetic relationship among the selected improved G. hirsutum and G. barbadense accessions including the recurrent parent TM-1 and the donor parent Pi- ma 3-79 of the CS-B25 line.The seven chromosome-specific SSR markers were used to construct a dendrogram based on the average coefficient of similarity using JMP Genomics program (SAS TM ).The method of Egamberdiev et al. [25] was used to verify SSRs via PCR amplification and capillary electrophoresis.

Transcriptome Analysis of 10 DPA Fiber from CS-B25 and TM-1
Approximately 300 million total sequence reads of 50 bp were generated for this study.Analysis of sequencing data of CS-B25 and TM-1 from the six RNA libraries revealed that more than 80% of reads were mapped to the TM-1 reference genome.Table 1 shows the total reads statistics of the six libraries.Of the 70,478 genes annotated in the G. hirsutum reference, 57,366 had at least one sequence read alignment.A total of 1870 genes (Table S3) were found to be differentially expressed in 10-DPA fiber when comparing CS-B25 to the TM-1 control at the false discovery rate (FDR) cutoff of 5% (adjusted p ≤ 0.05) and an absolute value of log 2 ratio ≥ 1.Among them, 1173 and 697 genes were up-regulated and down-regulated in CS-B25, respectively.

Gene Function Annotation by GO Enrichment Analysis
GO enrichment analysis by GAGE and ReviGO indicated that Biological Process GO terms related to DNA-templated transcription, response to oxidative stress, and cellulose biosynthetic processes were enriched (Figure 1(a)).In addition, the Molecular Function GO terms including structural constituent of cytoskeleton, peroxidase activity, cellulose synthase (UDP-forming) activity, and transcription regulatory region sequence-specific DNA binding were predominant (Figure 2(a)).The 57,366 G. hirsutum genes with expression data mapped to 14,610 homologous Arabidopsis genes were used to find enriched KEGG pathways.Most enriched pathways were found to be related to biosynthesis of secondary metabolites, biosynthesis of amino acids, starch and sucrose metabolism, amino sugar and nucleotide sugar metabolism, and plant hormone signal transduction (Table 2).The genes differentially expressed in CS-B25 10-DPA fibers were functionally annotated, and representatives of these differentially expressed genes (DEGs) are listed in Table 3 and Table S3.

DEGs Involved in Hormone Synthesis and Signaling
Three phytohormones of gibberellin, ethylene, and auxin have been shown to play important roles in cotton fiber development [26].Several genes participating in gibberellin synthesis and in response to gibberellin stimulus were up-regulated in CS-B25 10-DPA fiber (Table 3 and Table S3).These included genes for gibberellin 20 oxidase (Gh_A04G0224, 11 fold upregulation), gibberellin 3-beta-dioxygenase 4 (Gh_A04G0228, 10.41 fold upregulation) and a gibberellin-regulated family protein (Gh_A06G0045, 4.97 fold upregulation).Genes involved in ethylene synthesis and signaling were also found to be up-regulated.

DEGs Involved in Cell Wall Synthesis
Many genes involved in primary cell wall (PCW) biosynthesis were up-regulated in CS-B25 10-DPA fiber cells (Table 3 and Table S3).These include genes coding for fasciclin-like arabinogalactan proteins 7 (Gh_A08G0457, 3338.37 fold and Gh_D08G0544, 211 fold) [34], non-specific lipid transfer proteins (LTPs) American Journal of Plant Sciences  [36].Other genes related to the cytoskeleton, including tubulins beta 8 (Gh_Sca012883G01) was also up-regulated 30.5 fold.It had been suggested that specific tubulins might affect microtubule organization and microfibril deposition in elongating fiber cells [37].An ABIL-1 gene encoding a subunit of the WAVE complex was highly expressed (1,262 fold upregulation).The WAVE complex is involved in regulation of actin and microtubule organization [38].Genes for LRR receptor-like serine/threonine-protein kinases (Gh_D10G2231, 294 fold; Gh_D12G0597, 21.11 fold; Gh_D02G0280, 23.75 fold) were highly or moderately up-regulated, and they are known to participate in signal transduction pathways [39], cellulose deposition [40], and maintenance of cell wall integrity [41].

Transcription Factors Involved in Fiber Development
Many MYB-type transcription factors, found to regulate plant cell wall biosynthesis [42] [43], were also up-regulated in CS-B25 (Table 3 and Table S3).55 fold), a negative regulator of ABA signaling with function in transition from active cell division to cell expansion, were also up-regulated.NAC domain-containing protein 67-like (Gh_A08G1086, 0.0032), however, was down-regulated.NAC transcription factors and NAC domain-containing proteins are a very large family in plants [45], and the SND protein classified as NAC domain protein 1 was reported in regulating secondary wall synthesis by targeting the MYB46 transcription factor [42]. WRKY transcription factors 21 (Gh_A11G1172, 16.78 fold) and 29 (Gh_A02G1042, 3.27 fold) were also up-regulated in CS-B25.

Dendrogram Based on the Chromosome 25-Specific Simple Sequence Repeat (SSR) Markers
Results from the dendrogram analysis divided all of the accessions into two broad groups of G. hirsutum and G. barbadense based on the average coefficient of similarity (IBD value indicating identical alleles), suggesting that the chromosome 25 is distinctly different at the molecular level and fiber traits between the two species accessions including the recurrent parent TM-1 and the donor parent Pima 3-79 of the CS-B25 line (Figure 4).Pima 3-79, the donor parent, was clustered with other improved G. barbadense accessions in the dendrogram based on the chromosome 25-specific nuclear SSR markers associated with improved fiber traits from previous reports.This observation suggests that the substituted alien species chromosome 25 pair in CS-B25 would have the potential to carry genes of improved fiber traits including those associated with micronaire, fiber length and fiber strength.The results from field experiments demonstrated that the substituted chromosome 25 pair from G. barbadense in CS-B25 had significant effects on improved micronaire (Figure 5), which is one of important fiber quality traits used in determining fiber price.

Discussion
Results from GO enrichment analysis by GAGE and ReviGO indicated that many genes were up-regulated in CS-B25 10-DPA fiber with the highest fold change (FC) in expression of the gene coding for a fasciclin-like arabinogalactan protein 7 (AGP, Gh_A08G0457, 3338.37 fold upregulation) in CS-B25 compared to TM-1 (Table 3 and Table S3).Using over-expression and RNAi approaches,  S2).The results suggested that chromosome 25 is distinctly different between G. hirsutum and G. barbadense in fiber traits and molecular markers.The scale is presented as the branch length distance.Huang et al. [34] reported that AGPs are involved in fiber initiation and elongation.They reported that overexpression of a fasciclin-like arabinogalactan protein (GhFLA1) in G. hirsutum promoted fiber elongation and caused an increase in fiber length.In contrast, the suppression of GhFLA1 expression slowed down fiber initiation and elongation and produced shorter mature fiber.In addition, expression levels of GhFLAs and the genes related to primary cell wall biosynthesis were notably enhanced in the GhFLA1 overexpression transgenic fibers, whereas the transcripts of these genes were dramatically reduced in the fibers of GhFLA1 RNAi cotton plants.Using immunostaining methods, Huang et al. [34] further showed that both AGP composition and primary cell wall composition were changed in the transgenic fibers.The CS-B25 had several different differentially expressed AGP-like transcripts upregulated (Table S3), suggesting a major role of these genes at the fiber elongation stage.Previous studies revealed that the CS-B25 had several improved fiber quality traits including fiber length, strength, and micronaire [6]  A large scale analysis of fiber gene expression in the elongation step had been previously conducted via cDNA array [48] [49] and deep RNA sequencing [50] [51].Ji et al. [48] used subtractive PCR and cDNA macroarrays to identify 172 genes that were up-regulated in 10-DPA elongating fibers.They found that genes encoding a V-ATPase catalytic subunit, arabinogalactan proteins (AGPs), a kinesin-like calmodulin binding protein, cytochrome P450-like protein I, and proline-rich protein 5 were up-regulated.They also found that the mRNAs for proteins involved in signal transduction, lipid metabolism (such as LTP and ketoacyl CoA synthase), and cell wall-loosening enzymes (xyloglucan endotransglycosylase XET and expansins) were abundant.Using microarray analysis, Shi et al. [49] found that the ethylene biosynthesis pathway was up-regulated during fiber elongation, indicating that ethylene plays a role in promoting fiber elongation.Via high-throughput RNA sequencing, Wang et al. [50] conducted a global analysis of transcriptomes of two cotton lines, a wild type and a fuzzless/lintless mutant, at the early fiber development.They identified many genes differentially expressed in 2 -8 DPA ovules, which included genes for a long-chain fatty acid elongation enzyme, 3-ketoacyl-CoA synthase, an ACC oxidase, NAC domain protein NAC2, P450 monooxygenase, WRKY transcription factor 45, MYB19, pectin esterase, chitinase, lipid transfer protein, and tubulin β-1 chain.Transcriptome analysis of short fiber mutant ligon lintless-1 had shown that DEGs involved in very-long-chain fatty acid biosynthesis and cell wall metabolism were differentially expressed in fiber bearing ovules at 3 and 8 DPA [51].In this study, the identified DEGs in 10-DPA CS-B25 fiber were very similar to those found by Ji et al. [48], Wang et al. [50], and Liang et al. [51].We found that many genes involved in primary cell wall synthesis were up-regulated in CS-B25 during the fiber elongation stage, which included fasciclin-like arabinogalactan proteins [34], non-specific lipid transfer proteins, cellulose synthase, UDP-glycosyltransferase, pectin methylesterases, xyloglucan endotransglucosylase (XET), glycine-rich cell wall structural proteins, cytosolic pyruvate kinase [46], and ascorbate peroxidase [47].Pectin methylesterases (PME), a multigene family, catalyze the demethylesterification of galacturonic acid units in pectins, which are the major components of plant cell wall.PME have been shown to play an important role in plant cell wall extension during stem elongation, pollen germination and pollen tube growth, and root development [52].XETs are generally considered to be important enzymes in cell wall loosing, but they have been also shown to function in secondary cell wall genesis during the expansion process [53] [54].Both PME and XET were found to be up-regulated in a cotton MD52ne line which had high fiber strength [55].Several cotton non-specific lipid transfer protein genes, including pFS6, LTP3, LTP12, were highly expressed in fiber cells [56] [57] [58].LTP proteins have been proposed to participate in lipid transport for primary cell wall synthesis in fiber cells [59].Recently, the functional role of cytosolic pyruvate kinase and ascorbate peroxidase in fiber elongation was elucidated.Zhang and Liu [46] showed that a cotton GhPK6 gene encoding a cytosolic pyruvate kinase was preferentially expressed in fibers, and its transcript level increased during the fiber elongation stage from 5 to 20 DPA and reached the maximum at the late elongation stage (20 DPA).Pyruvate kinase catalyzes the conversion of PEP to pyruvate in glycolysis, and pyruvate then undergoes oxidative phosphorylation in TCA cycle to generate the by-product, H 2 O 2 , a form of reactive oxygen species (ROS) which was also found to have elevated levels during fiber elongation [46].The GhPK6 activity was negatively correlated with the rate of fiber elongation [46].An ascorbate peroxidase, a reactive oxygen species (ROS) scavenging enzyme encoded by the GhAPX1 gene, was shown to be up-regulated in 10-DPA fibers [47].Both the cytosolic pyruvate kinase and ascorbate peroxidase were involved in fiber elongation regulation in a ROS-mediated manner by modulating ROS homeostasis.A tubulin protein of tubulin beta 8 (Gh_Sca012883G01) and an ABI-1 like protein (Gh_A08G1131) in the assembly of cytoskeleton were also up-regulated (Table 3 and Table S3).It had been suggested that specific tubulins might affect microtubule organization and microfibril deposition in elongating fiber cells [37].Three LRR receptor-like serine/threonine-protein kinases (Gh_D10G2231, Gh_D12G0597, and Gh_D02G0280) were highly up-regulated in CS-B25 (Table 3 and Table S3).It had been reported that these membrane-associated kinases were involved in signal transduction pathways [39], participating in cellulose deposition [40] and maintenance of cell wall integrity [41].An LRR receptor-like kinase GhRLK1 was reported to be involved in secondary cell wall formation during fiber development [60].Seventy-nine genes encoding Catharanthus roseus receptor-like kinases (CrRLK1Ls) had been identified in G. hirsutum TM-1 line.Some of them were found to be up-regulated in cotton fiber of chromosomal segment introgression lines (CSILs) with excellent fiber traits during the elongation stages from 10 to 20 DPA, but down-regulated at the early (5-DPA) and later (25-DPA) stages of fiber elongation.These observations revealed that RLKs positively affected fiber development [61].
The results of transcriptome analysis of 10 DPA fibers of the CS-B25 and TM-1 lines indicated that the synthesis and signaling pathways of GA, ethylene and auxin increased in CS-B25 when compared with TM-1.Many publications have indicated that ethylene modulates the elongation of cotton fiber cells.In this study, genes encoding the enzyme for ethylene synthesis, 1-aminocyclopropane-1-carboxylate oxidase (ACCO) and ethylene responsive transcription factors are up-regulated.Genes coding for ethylene overproduction 1-like protein (in chromatin remodeling and transcription regulator) and NAC transcription factors were also up-regulated.The gene encoding a RING E3 ligase XBAT32 (Gh_A08G0780), which targets the ACS proteins, was down-regulated in the CS-B25 fiber.This down regulation resulted in a slight increase of the ACS level.Auxin has been shown to promote fiber initiation and elongation and play an important role in fiber development [26].The transcript levels of auxin-responsive GH3 protein, auxin-responsive IAA14 protein, and auxin influx carrier family protein in CS-B25 10 DPA fiber were up-regulated, suggesting that they may have important functions in the fiber elongation.
Several encoding MYB transcription factors were found to be highly up-regulated in the CS-B25 elongating fibers.A total of 524 non-redundant MYB genes had been identified in the upland cotton genome [62].Among them, the G. hirsutum MYB genes MYB2, MYB3, MYB5, MYB7, MYB20, MYB25, MYB42, MYB43, MYB46, MYB52, MYB54, MYB69, MYB85, MYB103, and MYB109 were shown to be highly expressed in fiber [42] [62], suggesting that they play important roles in fiber development.GhMYB109 was reported to participate in fiber development, and the suppression of its expression resulted in the reduction of fiber length [63].GhMYB25 was shown to regulate trichome formation and act as a key factor in early cotton fiber development [64].We had previously reported that a cotton MYB gene GhMYB7 was differentially expressed in developing fibers, and GhMYB7 was involved in transcriptional regulation of a LTP3 gene encoding lipid transfer protein [65].The GhMYB7 was later shown to activate genes involved in secondary cell wall thickening and biosynthesis [66].Besides the MYB transcription factors, homeobox leucine zipper (HD-Zip) transcription factors are also up-regulated in CS-B25.HD-Zip proteins were shown to play important roles in cotton fiber differentiation [67] [68].Many NAC transcription factors were also highly expressed in CS-B25 10-DPA fibers.A total of 143 and 145 NAC genes have been identified in the genomes of the two diploid cottons G. arboretum and G. raimondii, respectively [69].
Among them, more than 30 had been shown to have higher expression levels in 10 -20 DPA fiber.NAC genes had been shown to control secondary cell wall thickening during fiber development [69].Our transcript profile observations during fiber development are similar to those of Al-Ghazi et al. [70].Similar results had also been observed in some of Upland cotton chromosome introgression lines (CSILs) carrying different G. barbadense chromosomal segments in the recurrent parent TM-1 [71] [72].These CSIL lines with longer and stronger fiber had increased expression in genes involved in cell wall biosynthesis, hormone signal transduction, and in genes encoding transcriptions factors during the fiber elongation stage.Both CS-B25 and chromosome segment introgression lines (CSILs) were useful resources for transcriptome profiling to study molecular mechanism of positive fiber traits.
Genome wide analysis has identified 116 and 102 WRKY genes in G. ramondii and G. hirsutum, respectively [73].We have only found that two WRKY transcription factors, WRKY 21 and WRKY 29, were moderately up-regulated in CS-B25.This observation is consistent with the results of Dou et al. [73], in that many GhWRKY transcription factors had high expression levels at 0, 20, and 25 DPA fiber and they are involved in fiber initiation and secondary wall synthesis.Using the homologous alignment method, 143 genes for NAC-transcription factors in the G. arboretum genome were identified [69].Many of them were highly expressed in 15 DPA fiber cells and might be involved in the fiber secondary wall thickening.
The genetic effects in the CS-B25 line could be due to any of the following reasons: 1) the genes on the substituted chromosome, 2) the genes on the re-maining chromosomes of its recurrent parent, and 3) an interaction between the gene(s) on the substituted chromosome and the remaining chromosomes of its recurrent parent [8].A comparative analysis of CS-B25 with TM-1 is expected to show the differential transcripts from the direct effect of alien species substituted chromosome or an interaction effect of the substituted chromosome and the remaining chromosomes of CS-B25.Results showed that the majority genes upand down-regulated in 10-DPA fiber in CS-B25 were not located on Pima 3-79 chromosome 25, which is mapped to Gh_D06G in the G. hirsutum (AD1) NAU-NBI genome used in this study (Table S3) [74].This observation suggested that chromosome interactions between the substituted G. barbadense chromosome 25 pair and other non-substituted chromosomes of TM-1 might be involved in potential gene regulation.This is important considering that molecular markers (including EST markers) are used to construct linkage map of important QTLs with the expectation that these markers are directly associated with the QTL of interest in the same chromosomal region and can be used in marker assisted selection using genomic DNA to expedite the breeding program.However, our results showed that some of the ESTs may not have direct association with the chromosome but can cause a chromosomal effect on the genes.For example, the substituted G. barbadense chromosome 25 may contain potential regulatory elements that interact with gene located on other chromosomes in CS-B25 to cause an effect.Our results showed that several transcription factor like elements were associated with differentially expressed transcripts in CS-B25 line.It would be interesting in future investigation to identify the roles of some inter-chromosomal interactions in the activation or repression of gene transcription.The role of long-range chromosomal interactions involved in gene regulation has been reviewed in details by different research groups [75] [76].The CS-B25 line can be used as a tool for the study of cotton chromosomal interactions and their role in transcriptional regulation because of the unique genetic background of CS-B25.The Chromosome Conformation Capture (3C), a high throughput technique developed by Dekker et al. [77], has been widely used to study chromosomal conformation/interaction.Using the 3C-based technology would allow us to study chromosomal interactions in CS-B25 and reveal how the interactions might affect up or down regulation of gene expression.
The results from the dendrogram and fiber QTL analyses suggested that chromosome 25 from G. barbadense is associated with potential valuable fiber traits (Figure 4).A previous report [8] also documented that CS-B25 had important association with several improved fiber traits, including decreasing micronaire, increasing fiber length and strength compared to TM-1 and some USA commercial cultivars, suggesting that CS-B25 would be an important genetic tool to discover many useful genes associated with fiber development.Cotton fiber, an elongated and thickened single cell of the seed epidermis in which the fiber cell wall polymers help to determine primarily the cotton fiber quality.Cotton fiber cell undergoes several phase of changes since its initiation: 1) extreme elongation (to >2.5 cm) via primary wall synthesis till about 10 days; 2) then wall thickening and primary wall remodeling; and 3) secondary wall thickening via deposition of nearly pure cellulose; and 4) final maturation and cell death processes before boll opening to reveal the fluffy fiber within the cotton bolls [78].Our results showed that CS-B25 had several up or down regulated differentially expressed transcripts in 10 DPA fiber associated with potential cytoskeleton or cell wall genes (Table 3).It is important to note that most of the improved cotton fiber quality traits are primarily controlled by the cytoskeleton or cell wall related genes during fiber development [77].These results suggested that future investigation with these differentially expressed transcripts, especially very high upregulated expression of fasciclin-like arabinogalactan protein 7 genes in CS-B25 10 DPA fiber cell, will have great potential to elucidate genes associated with fiber quality traits at the molecular level.

Conclusion
This study showed that the CS-B25 line is a valuable genetic resource to discern many novel differentially expressed genes due to the effect of the substituted alien species chromosome 25 or its interaction with other chromosomes in CS-B25.The CS-B25 line had several valuable fiber quality traits including micronaire and fiber length associated with the substituted chromosome 25 from G. barbadense.We discovered that several genes involved in hormone signal transduction, cell wall synthesis, and transcription factors in regulating cell wall synthesis are up-regulated in 10 DPA fiber of the CS-B25 line.These differentially expressed genes could be used as potential resources for fiber quality improvement.This research will facilitate our understanding of molecular mechanisms of fiber cell elongation and will allow us to develop new approaches for cotton fiber improvement.

2. 5 .
Quantitative Real-Time PCR Twelve DEGs including 8 up-regulated and 4 down-regulated in CS-B25 were selected for validation using qRT-PCR.One microgram of fiber 10-DPA total RNA was used as template for first-strand cDNA synthesis using the random hexamer (Promega, Madison, Wisconsin, USA) and M-MLV reverse transcriptase (Promega, Madison, Wisconsin, USA) according to the manufacturer's instructions.The synthesized cDNA was then used as template for qRT-PCR reactions using PowerUp SYBR Green Master Mix and a QuantStudio 5 Real-Time PCR System (Applied Biosystems, Waltham, Massachusetts, USA).Cotton 18S

Figure 1 .
Figure 1.A ReViGo scatter plot of GO analyses in biological process from (a) up-regulated and (b) down-regulated DEGs in the CS-B25 10-DPA fibers.Related GO terms are clustered and presented with bubbles with similar color shades.Bubble colors represent p-values and bubble sizes indicate the relative frequency of the GO terms.

Figure 2 .
Figure 2. A ReViGo scatter plot of GO analyses in molecular function from (a) up-regulated and (b) down-regulated DEGs in the CS-B25 10-DPA fibers.Related GO terms are clustered and presented with bubbles with similar color shades.Bubble colors represent p-values and bubble sizes indicate the relative frequency of the GO terms.

Figure 3 .
Figure 3. Quantitative RT-PCR analysis of differentially expressed genes in 10 DPA cotton fibers from CS-B25 and TM1 lines.Compared to TM1, eight selected up-regulated genes (a) and four selected down-regulated genes (b) identified in CS-B25 line from RNA-Seq analysis were validated with quantitative RT-PCR assays.At least three biological replicates were used in the assays, cotton 18S rRNA was used as a reference gene to normalize the expression level.The significant difference in the expression of each transcript between CS-B25 and TM1 was analyzed by t-test (*p ≤ 0.05; **p ≤ 0.005; and ***p ≤ 0.005).The error bars represent standard error of the mean.

Figure 4 .
Figure 4.The dendrogram analysis of cotton fiber traits-associated SSR markers.The dendrogram of a set of improved Gossypium hirsutum and G. barbadense cotton lines constructed using JMP Genomics program on the basis of association of fiber traits with seven chromosome 25 specific SSR markers (TableS2).The results suggested that chro-

Figure 5 .
Figure 5.Comparison of micronaire and 2.5% span length of CS-B25 fiber with TM-1 and Pima 3-79 lines.Some of the differentially expressed transcripts in CS-B25 might be associated with improved fiber quality traits.CS-B25 and TM-1 do not have significant differences in Fiber 2.5% Span length, although CS-B25 had average values higher than TM-1.The letters marked above bars indicate statistically significant differences (p ≤ 0.005 in the micronaire analysis and p ≤ 0.0005 in the analysis of 2.5% span length).The error bars represent standard deviation of the mean.
[7] [8][9] [10].These results were further confirmed from the top crosses with the selected USA improved cultivars that CS-B25 had several beneficial alleles with additive gene effects for improving fiber strength and micronaire in Upland cultivars[9] [10].Our results suggest that the improvement of fiber traits in CS-B 25 line was in part associated with very high upregulation of GhFLA genes in fiber cell due to the introgression of the substituted chromosome of Pima 3-79 into the CS-B25 line.The dendrogram results based on the SSR markers also separated all lines into two broad groups of G. hirsutum and G. barbadense.It will be interesting in future to investigate the relation of differentially expressed transcripts of fasciclin-like arabinogalactan proteins among the lines of these two species with reference to the improved fiber traits including elongation.This is important considering that fiber elongation is one of the most important economic traits for determining the fiber price in global textile market.C.-Y. Hsu et al.DOI: 10.4236/ajps.2018.960981350 Journal of Plant Sciences

Table 1 .
Total reads statistics of RNA-Seq data

Table 2 .
KEGG pathway analysis of KOs.
*FC represents fold change To test whether the gene expression results from RNA-seq were reliable, RNA C.-Y.Hsu et al.
DOI: 10.4236/ajps.2018.96098American Journal of Plant Sciences transcripts of 16 DEGs detected in CS-B25 were selected for qRT-PCR analysis.All the 16 transcripts were successfully detected by qRT-PCR and 12 transcripts