Comparative Transcriptome Analysis of Short Fiber Mutants Ligon-Lintless 1 And 2 Reveals Common Mechanisms Pertinent to Fiber Elongation in Cotton (Gossypium hirsutum L.)

Understanding the molecular processes affecting cotton (Gossypium hirsutum) fiber development is important for developing tools aimed at improving fiber quality. Short fiber cotton mutants Ligon-lintless 1 (Li1) and Ligon-lintless 2 (Li2) are naturally occurring, monogenic mutations residing on different chromosomes. Both mutations cause early cessation in fiber elongation. These two mutants serve as excellent model systems to elucidate molecular mechanisms relevant to fiber length development. Previous studies of these mutants using transcriptome analysis by our laboratory and others had been limited by the fact that very large numbers of genes showed altered expression patterns in the mutants, making a targeted analysis difficult or impossible. In this research, a comparative microarray analysis was conducted using these two short fiber mutants and their near isogenic wild type (WT) grown under both field and greenhouse environments in order to identify key genes or metabolic pathways common to fiber elongation. Analyses of three transcriptome profiles obtained from different growth conditions and mutant types showed that most differentially expressed genes (DEGs) were affected by growth conditions. Under field conditions, short fiber mutants commanded higher expression of genes related to energy production, manifested by the increasing of mitochondrial electron transport activity or responding to reactive oxygen species when compared to the WT. Eighty-eight DEGs were identified to have altered expression patterns common to both short fiber mutants regardless of growth conditions. Enrichment, pathway and expression analyses suggested that these 88 genes were likely involved in fiber elongation without being affected by growth conditions.


Introduction
Cotton fibers are single-celled trichomes that initiate from the ovule epidermal cells on or about the day of anthesis (DOA) [1]. Approximately 25% of the ovule epidermal cells differentiate into fiber cells during the initiation stage of cotton fiber development and subsequently undergo a period of rapid elongation known as the elongation stage [2,3]. The rate of fiber elongation peaks at approximately 6 to 12 days post-anthesis (DPA) and nears cessation around 22 DPA [4]. During peak elongation fiber cells can increase in length at rates of 2 mm/day or more depending on environmental factors and genotypes [5][6][7]. Beginning at 12-16 DPA and overlapping with the elongation phase is the secondary cell wall (SCW) biosynthesis stage. During this stage cellulose is synthesized and deposited between the primary cell wall and the plasmalemma [8,9]. Elongation and SCW biosynthesis continue until the fibers reach full length [25-35 mm in Upland cotton (Gossypium hirsutum L.) cultivars] [10], after which the cotton bolls open and the fibers desiccate under exposure to the environment. The environmental and genetic factors that influence the timing of these processes are shown to influence the development of desirable fiber traits such as lint yield and fiber quality [5,[11][12][13].
Understanding the molecular mechanisms of fiber development is essential for cotton researchers to devise strategies for developing cotton lines with superior fiber quality. Furthermore, it is important to identify key genes that could be genetically engineered to improve fiber properties. Toward these goals, scientists have been using fiber mutants to study the molecular and genetic mechanisms of fiber development [14][15][16][17]. Among them, two short fiber mutants Ligon-lintless 1 and Ligon-lintless 2 (Li 1 and Li 2 , respectively) were extensively studied by our group [18][19][20] and others [14,21,22] in order to develop a comprehensive understanding of the molecular and metabolomic mechanisms related to cotton fiber length development. In a near-isogenic state with the cotton cultivars Texas Marker-1 (TM-1) or DP5690, both the Li 1 and Li 2 mutants have seed fibers that are extremely short (,6 mm) as compared to wild type (WT) fibers that are typically longer than 25 mm in length [18,19,[23][24][25]. As a monogenic dominant trait, the short fiber phenotypes of Li 1 and Li 2 are similar (Fig. 1A). However, unlike the Li 2 mutant, which appears morphologically similar to the WT plants with the exception of short seed fibers, the Li 1 mutant exhibits pleiotropy in the form of severely stunted and deformed plants in both the homozygous dominant and heterozygous state (Fig. 1B) [18,23]. Although the Li 1 and Li 2 mutants have similar phenotype in cotton fiber, these two genes reside on different chromosomes with Li 1 on chromosome (Chr.) 22 and Li 2 on Chr.18 [18,19,26,27]. These two mutants, when taken in combination, provide an excellent experimental system to find both common and mutant locusspecific mechanisms related to fiber elongation.
Analyzing the microarray or RNA-seq data collected to date for the Li 1 or Li 2 mutant is limited in one aspect by the fact that very large number of genes showed altered expression patterns between a mutant and its WT near isogenic line (NIL). For example, previous microarray data obtained from fibers of Li 1 or Li 2 showed approximately 1,500 to 2,500 differentially expressed genes (DEGs) including many genes that may be affected or regulated by environments [18,19,21]. It is difficult to decipher which of the genes are truly vital and common to fiber elongation-related processes, and which are due to different environmental, genetic and physiological cues if using only one mutant in an experiment as reported in all the previous studies. In recognition of this issue, we conducted the present experiment with two mutant lines and two growth conditions in order to identify the genes that are differentially regulated in both mutants, with the goal of identifying the common molecular mechanisms involved in cotton fiber length development. First we took advantages of our unique NILs of the Li 1 and Li 2 mutants. As reported earlier [18,19], both Li 1 and Li 2 NILs were developed using the Upland cotton cultivar DP5690 as the recurrent parent. The Li 1, Li 2 and WT DP5690 are mutually near isogenic. Second, we conducted experiments in both field and greenhouse for Li 2 mutant, allowing the identification of genes impacted by environmental conditions and response to stress in this mutant. Third, we did a comparative analysis of transcriptome profiles between Li 1, Li 2 and WT to identify genes that had altered expression patterns in a short fiber mutant, regardless of the growth conditions (field or greenhouse) or the nature of the mutation (Li 1 or Li 2 ). Our major objective was to identify common genes (or molecular mechanisms) that were essential to the fiber elongation regardless of environment or a specific mutation. Herein we report our findings.

Plant Materials
The Li 1, Li 2 and WT DP5690 used in the present study were mutual NILs. Development of these NILs was described in our earlier reports [18,19]. For the greenhouse-grown Li 2 plants utilized in this study, growth and sample conditions were described in   [19], and the growing period was between October, 2009 and March, 2010. Each plant was grown in an 18.9 L pot. A commercial service provider periodically sprayed pesticides to control insects or diseases. Automatic drip irrigation was used throughout the growing season. For field grown plants, a total of 200 Li 1 , 100 Li 2 , and 100 WT DP5690 plants were grown in a field at the USDA-ARS Southern Regional Research Center, New Orleans, LA in the summer of 2012. The distance between two plants within a plot was 30 cm. The plot  distance was 45 cm. The soil type in this field was aquent dredged over alluvium in an elevated location to provide adequate drainage. Throughout the growing season, no pest control spray was applied. Supplemental sprinkle irrigation was provided when needed. Flowers were tagged and cotton boll sample collections were made before 10:00 a.m. and immediately placed on ice. To minimize environmental effects, boll samples were not collected from plants on the perimeter of the field and samples were only collected when 15-30 bolls were available for analysis. All samples of the same developmental stage were tagged and collected on the same day. Bolls were randomly separated into 3 replicates with 5-10 bolls per replicate. Bolls were then dissected, frozen in liquid nitrogen and stored at 280uC until further processing.

RNA Isolation from Cotton Fibers
RNA was isolated as previously described [19]. To separate the fibers from the ovules the samples were shaken vigorously enough to break fibers without damaging the ovules. Isolation of RNA was conducted using the Sigma Spectrum Plant Total RNA Kit (Sigma-Aldrich, St. Louis, MO) with on-column DNaseI digestion and used according to the manufacturer's instructions. RNA quantity was determined by using a Nanodrop 2000 spectrophotometer (NanoDrop Technologies Inc., Wilmington, DE). RNA integrity number (RIN) was determined for each sample using an Agilent Bioanalyzer 2100 and the RNA 6000 Nano Kit Chip (Agilent Technologies Inc., Santa Clara, CA). Only samples with RIN values of 7.0 or higher were used for expression analysis.

Reverse Transcription Quantitative Real-time PCR (RT-qPCR)
The experimental procedures and data analysis related to RT-qPCR were performed according to the Minimum Information for Publication of Quantitative Real-Time PCR Experiments (MIQE) guidelines [28]. We used RNA samples from 12 DPA fibers in two biological replicates for cDNA synthesis and in two technical replicates for qPCR. The cDNA synthesis reactions were performed using the iScript cDNA Synthesis Kit (Bio-Rad Laboratories, Hercules, CA) according to the manufacturer's instructions with 1 mg of total RNA per reaction used as template. The RT-qPCR reactions were performed with iTaq SYBR Green Supermix (Bio-Rad Laboratories) in a Bio-Rad CFX96 real time PCR detection system. The detail description of amplification parameters and calculation reported before [19]. Normalization of RT-qPCR data was performed by geometric averaging three internal control genes, including 18 S rRNA, ubiquitin-conjugating protein, and a-tubulin 4 [29]. The primer sequences of the 14 probe sets and the three internal control genes are listed in Table  S1.

Microarray Hybridizations and Data Analysis
The microarray technology used for this study was the commercially available Affymetrix GeneChip Cotton Genome Microarray (Affymetrix Inc., Santa Clara CA), comprised of 239,777 probes representing 21,854 cotton transcripts from a variety of expressed sequence tag (EST) databases. The source material for the EST data was derived from G. arboreum, G. barbadense, G. hirsutum, and G. raimondii. Labeling of the RNA was conducted using the Affymetrix GeneChip 39 IVT Express Kit and cotton genome hybridizations were conducted according to the manufacturer's protocols. Our earlier studies [18][19][20] indicated that a significant difference in both transcript profiles and fiber length measurement was observed at 12 DPA (peak elongation) between the Li 1 , or Li 2 and its WT NIL. Thus in this experiment, we used RNA samples from 12 DPA fibers (in two biological replicates) for microarray analysis. Probes sets demonstrating a two-fold or greater difference in expression levels between experimental samples were considered differentially regulated. Data normalization and the determination of statistically relevant deviations in expression patterns was performed as described [30].

Gene Annotation Analysis
Microarray data obtained from greenhouse-grown Li 2 plants, and field-grown Li 1 and Li 2 plants were first subjected to Venn analysis utilizing BioVenn [31] to determine which probes demonstrate consistent expression profiles between experimental sets. To assist in the identification of biological processes represented in the data, Gene-Ontology Enrichment Analysis (GOEA) was performed using the agriGO Singular Enrichment Analysis (SEA) [32] by comparing to the Gossypium raimondii reference genome sequence [33]. The statistical test method used was the Fisher's Exact test (significance level 0.05). Annotation of the probes was accomplished with Blast2Go [34]. For pathway analysis, MapMan software [35] was used to identify and illustrate pathways of interest using the January 12, 2013 Gossypium hirsutum     Table 1

Effects of Growing Conditions (Field and Greenhouse) on Gene Expression
The data obtained in these experiments allowed for an analysis of the environmental effects on transcriptome data, i.e.; a field vs greenhouse comparison. This was useful for our purpose as it allowed the identification of transcripts affected by variable environmental conditions, which could then be excluded from consideration of strictly fiber-related transcripts. Although it is known that many fiber-related genes are environmentally impacted [5], the exclusion of these genes permitted a more targeted analysis of the genes that are essential to fiber elongation. It also allowed us to investigate how the mutant line differed from its WT in responding to environmental stressors, providing insight into the interaction between stress response and fiber elongation. Microarray analysis was conducted on WT and Li 2 12 DPA bolls collected from cotton plants grown in both field and greenhouse conditions. This comparison identified 150 probes in the WT that were expressed higher in the greenhouse than in a field, and 754 probes that were expressed higher when grown in field conditions (numbers are sums shown in orange ovals of Fig. 2). Gene enrichment analysis of the probes higher in the field showed that genes related to nucleosome assembly (GO:0006334), flavonoid biosynthetic process (GO:0009813), and response to heat (GO:0009408) were enriched in the field conditions ( Fig. 2 and Table S2). Probes showing higher expression patterns in the greenhouse in WT were not enriched in any particular GO category to a statistically significant degree due to smaller number of probes, however did consist of ethylene and ap2 erf domaincontaining transcription factors, expansin, and probes related to NADH dehydrogenase. The results of a similar analysis conducted in the Li 2 mutant differed dramatically from what was observed in the WT, with 1,275 probes showing higher expression in the greenhouse and 1,136 probes showing higher expression in the field (numbers are sums shown in green ovals of Fig. 2). The probes showing higher expression in Li 2 in the greenhouse were enriched in lipid transport (GO:0006869), cellular nitrogen compound metabolic process (GO:0034641), and iron ion binding proteins (GO: 0005506), whereas, probes showing higher expression in the field were enriched in mitochondrial electron transport (GO:0006120) and response to reactive oxygen species (GO: 000302) ( Fig. 2 and Table S2). Only 124 probes showed the same expression profile in WT and Li 2 in both field and greenhouse conditions. This difference between WT and Li 2 in response to the environmental conditions is profound, as very few probes showed similar expression profiles between the two genotypes, (i.e., only 6 common probes were higher in the green house and only 118 were higher in the field).
Gene enrichment analysis for Li 2 specific greenhouse/field condition response identified two categories of probes that are of particular interest; probes related to mitochondrial electron transport (GO:0006120) and response to reactive oxygen species (ROS) (GO:0000302) were both enriched in Li 2 in the field relative to their expression in the greenhouse. Probes for genes in these categories were not enriched in WT. It has been suggested in previous studies that the inability of certain Gossypium species and G. hirsutum mutants to produce long fibers may be due to the inability to modulate ROS homeostasis [19,37,38]. The mitchondrial electron transport gene NADH dehydrogenase (GhiAffx.21609.1.S1_at, GhiAffx.53261.1.A1_at, GhiAffx.45916.1.A1_s_at) and NADH plastiquinone reductase (GhiAffx.4260.1.S1_at, GhiAffx.61308.1.S1_at, GhiAffx.9732.1.A1_at) were among those highly expressed in Li 2 in the Table 2. Common DEGs involved in vesicle (GO:0031982) in both Li 1 and Li 2 mutants regardless of growth conditions (P value, 1.7e-11; FDR, 3.1e-10).

DEGs Annotation
Ghi. 68  field condition. Because the field conditions presented both abiotic and biotic stresses that were absent or minimized in greenhouse conditions, it was likely that the field grown Li 2 plants needed to divert limited cellular resources to manage these additional stresses, leading to higher ROS accumulation and an even higher expression of ROS homestasis genes. Expression levels of mitochondrial-related genes were studied across all three microarray data sets. In field conditions Li 1 Fig. 3A). The 2 probes down-regulated (red squares in Fig. 3A) consisted of Gra.1550.1, which has sequence identity to a Choline transporter-related transcript (AT4G38640) and Ghi.648.1.a1_at, which has sequence similarity to an NADH dehydrogenase. The probe that was up-regulated in all 3 data sets in the Li mutants is homologous to the probe GhiAffx.45916.1, a subunit of NADP dehydrogenase (NAD2) (shown by a blue arrow in Fig. 3A). These results further demonstrated the potential relevance of ROS and stress response in fiber developing processes.
Chalcone synthase is the upstream enzyme in the flavanoids synthesis pathway responsible for production of secondary metabolites (flavonols, proanthocyanins and anthocyanin) that are often produced in response to stresses [39,40]. Here, the chalcone synthase-related probes in WT plants had significantly higher expression in field conditions than in the greenhouse as indicated by a significant enrichment in flavonoid biosynthetic processes (GO:0009813) (Table S2). Additionally, pathway analysis utilizing MapMan software indicated that under greenhouse condition, Li 2 had higher expression levels of chalcone synthase than WT (Fig. 3B). Thus, all of the plants in the field (WT, Li 1 and Li 2 ), and the Li 2 plants in the greenhouse had high levels of genes related to flavonol production relative to the WT plants in the greenhouse. This further supports the hypothesis that the Li 2 mutant, even in greenhouse conditions, were in a stressed state. Ghi.6103.1 that codes for chalcone synthase 3 (GhCHS3) (Fig. 3B, red squares) was the only identified flavonoid probe that decreased in Li 2 in both field and greenhouse conditions. The remaining probes that increased under field conditions demonstrate varying degrees of homology to Transparent Testa 4 (TT4) (8 probes) and TT5 (2 probes), both naringenin-chalcone synthases.
In brief, it was likely that field-grown plants were under more stressful conditions than greenhouse grown plants. Short fiber mutants (more specifically Li 2 mutant) commanded higher expression of genes related to energy production and transport such as mitochondrial electron transport or responding to ROS in order to fight against the stresses than the WT.

Identification and Annotation of DEGs Common to both Short Fiber Mutants Regardless of Growth Conditions
Although the Li 1 and Li 2 mutations are caused by different genes located on different chromosomes, both result in a similar short fiber phenotype. Thus it is possible that probes which are similarly altered in expression pattern between a mutant and WT in both field and greenhouse conditions are highly likely to be fiber elongation-related or specific genes. In all three data sets, 113 probes are commonly affected in the mutants in comparison with WT with 94 down-regulated and 19 up-regulated (Fig. 4). Of these 113 probes that showed altered regulation, 25 also showed to be differentially regulated between field and greenhouse, implying their altered regulation could be an environmentally controlled factor and were excluded from further consideration. The remaining 88 probes are likely the genes specific to fiber elongation. Among them, 66 genes were annotated by blastn or blastx and were used for further analyses (Table 1 for annotated  genes and Table S3 for full list).
Among the fiber elongation specific genes, actin depolymerizing factor (ADF5) (GhiAffx.33535.1) was decreased in all conditions analyzed. ADF family proteins have previously been shown to affect the 3-dimensional structure of actin filaments [41] and to alter the disassociation rates of actin subunits [42]. The downregulation of GhADF1 in transgenic cotton show altered fiber length and strength [43] and GhADF5 has been shown to localize to elongating cells of the root stem in A. thaliana [44]. As shown in Table 1, actin (Gra.1375.1.A1_at) and tubulin (Ghi.8448.1.S1_-x_at) were commonly down-regulated in short fiber mutants. It was proposed that actin and microtubule play important roles in fiber elongation [45]. Functional analysis showed that actin was indeed required for fiber elongation [46].
Two probes (Ghi.3235.1 and Ghi.9654.1.) showing altered regulation in all conditions are different regions of the same gene, UGT73C14, which codes for an UDP glycosyltransferase that glycosylates ABA in vivo and in vitro [47]. One probe (GhiAffx.53295.1.A1_at) also codes for UDP-glucuronosyl and UDPglucosyl transferase, but has not been characterized in the context of elongation and warrant further investigation. Ghi.6548.1.S1 codes for a MAP kinase-like protein that was previously identified as preferentially expressed from elongating fibers in G. hirsutum by subtractive PCR but remains otherwise uncharacterized [48]. Two probes, GhiAffx.43008.1 and GhiAffx.24789.1 code for SAUR (Small Auxin Up RNA) family genes, which comprise the largest family of auxin-responsive genes. However, some members have also reportedly been associated with cell expansion [49]. This is potentially evidence of a substantial and specific link between auxin regulation and fiber elongation.
Two probes that decreased in all conditions analyzed code for either plant invertase or pectin methylesterase (PME) inhibitors (GhiAffx.33585.1 and Ghi.5186.1). Invertases catalyze the hydrolysis of sucrose, and PMEs inhibit the enzyme that catalyzes the de-methylesterification of pectins, a process important in initiation and for the polysaccharides role in the primary cell wall of cotton fibers [50,51]. Li 1 and Li 2 have previously demonstrated a decrease in de-esterified pectin localized to the primary cell wall during elongation stages as indicated by antibody staining [52]. Previous studies have also demonstrated that there are at least 5 PMEs showing fluctuating expression profiles throughout fiber initiation, elongation, and transition to secondary cell wall synthesis, indicating that different PME regulated different fiber development processes [53]. It was further demonstrated that the timing of individual PME gene expression differed between the longer fibered Pima (G. barbadense) and Upland (G. hirsutum) species in such a manner that suggested the timing of PME activity affected the onset of the transition stage from elongation to secondary cell wall synthesis, in turn affecting length of the fiber [53]. Since it is suggested that PME activity generally corresponds to longer fibers [54][55][56], it is possible that the increased expression of PME inhibitor(s) identified in our study play a role in inhibiting mid or late elongation stage PMEs, thus low levels of inhibitor would induce early transition to secondary cell wall synthesis. This would effectively decrease the elongation period. In support of this, the PME inhibitor probed by Ghi.5186.1 also exhibited decrease expression ratios of 0.42 in Li 1 and 0.58 in Li 2 at 3 DPA (data not shown).

Gene Ontology Analysis of the DEGs Common to both Short Fiber Mutants
To identify the potential biological processes governing differential expressions of genes that affect fiber length development of the two short fiber mutants, we analyzed the identified 88 DEGs using agriGo SEA analyses. The GO enrichment analysis classified 9 DEGs as a group that is involved in vesicle transportation (GO:0031982) in cotton fiber cells (Fig. 5). The vesicle plays an important role of carrying membrane components to the growing site of elongating cells [57,58]. Among the 9 vesiclerelated DEGs (Table 2), fasciclin-like arabinogalactan protein (GhFLA1, Ghi.68.1.A1_s_at) was reported as an essential cotton fiber gene for fiber elongation and primary cell wall biosynthesis [59].

Co-expression Analysis of Annotated DEGs
Beginning with our list of genes showing altered regulation in both mutants in all three data sets and no variability between field and greenhouse (Table 1), we wanted to identify genes that were previously identified as being co-expressed, and further support their roles in a common pathway and their involvement in the phenotype of the mutant plants. Subjecting the aforementioned list to the ATTED-II Arabidopsis co-expression database revealed a putative gibberellic acid-regulated pathway (Fig. 6). Several GDSL esterase/lipase family genes have been found with giberellinresponsive elements (P-box and GARE motif) in their 59 upstream regulatory regions [60]. These genes are characterized by the presence of a conserved motif and have been shown to be involved in a wide range of functions, including stress response and development [61,62], however to our knowledge this is the first report of a possible link with fiber development. The only gene upregulated in the identified pathway, 4-coumarate-coenzyme A ligase, has been shown to play a role in lignin deposition in the plant cell wall of several species, which has exhibited a decrease in plant height when the gene was suppressed [63,64]. In cotton fibers, lignin is a structural component of the primary cell wall, and genes that are involved in lignin deposition are up regulated in parallel with fiber elongation [65]. Finally, an additional topic that is re-visited in this analysis is the role of ROS-related genes in fiber development identified by this and other studies. Experimental evidence has demonstrated that the membrane bound glyoxal oxidases, which decreases in the mutants examined here, are known to produce hydrogen peroxide [66,67], which is then used as a cofactor for lignin-synthesizing peroxidases [68]. Thus, the genes represented in Fig. 5 could help explain the apparent perturbation to ROS-related processes affected in the Li 1 and Li 2 mutants.

Cell-wall Precursors
Because our analysis was focused on cell wall development we utilized MapMan software to determine which steps in the synthesis of cell wall-related polysaccharides are affected in the mutants. Fig. 7 illustrates that the conversion of myo-inositol to Dglucaronic acid by myo-inositol oxygenase (MIOX), an early step in the synthesis UDP-D-glucaronic acid, was affected. Between the 3 experimental groups, there are 4 separate probes that were upregulated (Fig. 7 Likely related to this, MapMan software also revealed the alternative pathway for UDP-D-glucaronic acid synthesis contains probes for genes that were down-regulated (red dots). The conversion of UDP-D-glucose to UDP-D-glucaronic acid is mediated by UDP-glucose 6 dehydrogenase, which was downregulated in Li 2 /greenhouse (Gra.2095.2) and Li 1 /field (Ghi.8750.1) (red arrows). This gene was unaffected in Li 2 /field, however GhiAffx.64086.1, which codes for phosphofructokinase 3 and mediates the formation of D-Glucose-1-P, was downregulated, likely indicating an alternative site of regulation (red arrow). The regulation of UDP-D-glucaronic acid is a vital step, as it serves as the common precursor for arabinose, xylose, galacturonic acid, and apiose residues in the cell wall. It has been reported in null Arabidopsis mutants that the limited supply of myo-inositol generally prevents the effective induction of this pathway as an alternative to UDP-glucose formation [69]. Thus it remains a strong possibility that the UDP-D-glucose levels and other polysaccharides necessary for cell wall development were perturbed and insufficient in the Li 1 and Li 2 mutant fiber tissues .

Verification of Microarray Results by RT-qPCR Analysis
To test the reliability of microarray data, RT-qPCR analysis was performed for a subset of 14 genes selected from Table S3 showing altered regulation in Li 1 , Li 2 (field) and Li 2 (greenhouse) as compared to the wild type. Expression of selected genes was tested at 12 DPA of fiber development, the common time point between greenhouse and field experiments. Overall, the results of RT-qPCR analysis were consistent with results of microarray for 14 selected genes (Table S4).

Conclusions
The development of cotton fibers is a very complicated and poorly understood biological process, however understanding this process is vital for the targeting of genes to use in the creation of value-added crop. We have developed a genetic model system consisting of two short fiber cotton mutants, Li 1 and Li 2 , which when combined together with their near-isogenic WT line allows for the study of genes and processes specific to fiber elongation. Here we analyzed multiple transcriptome profiles obtained from Li 1, Li 2 short fiber mutants and their WT grown under different environmental conditions. We classified the differentially expressed genes into two groups: one was mainly affected by environmental conditions, and the other was largely regulated by Li 1 and Li 2 . Our results provide new insight to how environmental factors affect fiber elongation by transcriptional regulation.
Further, the short list of 88 genes required for fiber elongation without being affected by environmental conditions would warrant further investigation in hope to identify targets for improving cotton fiber property.