RNA-seq analysis of short fiber mutants Ligon-lintless-1 (Li1) and – 2 (Li2) revealed important role of aquaporins in cotton (Gossypium hirsutum L.) fiber elongation

Cotton fiber length is a key determinant of fiber quality for the textile industry. Understanding the molecular basis of fiber elongation would provide a means for improvement of fiber length. Ligon lintless-1 (Li1) and Ligon lintless-2 (Li2) are monogenic and dominant mutations, that result in an extreme reduction in the length of lint fiber to approximately 6 mm on mature seeds. In a near-isogenic state with wild type (WT) cotton these two short fiber mutants provide an excellent model system to study mechanisms of fiber elongation. We used next generation sequencing (RNA-seq) to identify common fiber elongation related genes in developing fibers of Li1 and Li2 mutants growing in the field and a greenhouse. We found a large number of differentially expressed genes common to both mutants, including 531 up-regulated genes and 652 down-regulated genes. Major intrinsic proteins or aquaporins were one of the most significantly over-represented gene families among common down-regulated genes in Li1 and Li2 fibers. The members of three subfamilies of aquaporins, including plasma membrane intrinsic proteins, tonoplast intrinsic proteins and NOD26-like intrinsic proteins were down-regulated in short fiber mutants. The osmotic concentration and the concentrations of soluble sugars were lower in fiber cells of both short fiber mutants than in WT, whereas the concentrations of K+ and malic acid were significantly higher in mutants during rapid cell elongation. We found that the aquaporins were the most down-regulated gene family in both short fiber mutants. The osmolality and concentrations of soluble sugars were less in saps of Li1 – Li2, whereas the concentrations of malic acid, K+ and other detected ions were significantly higher in saps of mutants than in WT. These results suggest that higher accumulation of ions in fiber cells, reduced osmotic pressure and low expression of aquaporins, may contribute to the cessation of fiber elongation in Li1 and Li2 short-fiber mutants. The research presented here provides new insights into osmoregulation of short fiber mutants and the role of aquaporins in cotton fiber elongation.


Background
Cotton is the major source of natural fibers used in the textile industry. Apart from its economic importance, the cotton fiber provides a unique single-celled model system to study cell elongation and cell wall biogenesis in the absence of cell division [1]. Cotton fiber development consists of four distinct but overlapping stages, including fiber initiation, elongation, secondary cell wall biosynthesis, and maturation [1]. Each cotton fiber is a single cell that initiates from the epidermis of the outer integument of the ovules at or just prior to anthesis [2]. Fiber elongation starts on the day of anthesis and continues for about 3 weeks before the cells switch to intensive secondary cell wall cellulose synthesis. Lint fibers of the economically important Gossypium hirsutum generally grow about 30-40 mm in length. During peak elongation fiber cells can increase in length at rates of 2 mm per day or more depending on environment and genotype [1][2][3]. The fiber cells elongate up to 3000 fold during 3 weeks which makes them the fastest growing and longest single cell known in higher plants [4]. Understanding the molecular basis of fiber elongation would provide a means for cotton breeders and researchers to improve the fiber length while maintaining yield and other cotton characteristics.
Genetic mutants are useful tools for studying the molecular mechanisms of fiber development. Our laboratory uses two short fiber mutants, Ligon lintless-1(Li 1 ) and Ligon lintless-2 (Li 2 ) as a model system to study fiber elongation [5][6][7][8][9][10]. Both Li 1 and Li 2 are monogenic and dominant mutations, resulting in an extreme reduction in the length of lint fiber to approximately 6 mm on mature seeds [11,12]. Both mutations are located in the D T subgenome of G. hirsutum: the Li 1 gene is on chromosome 22 [8,13,14], whereas the Li 2 gene is on chromosome 18 [5,10,14,15]. Cytological studies of cotton ovules did not reveal much difference between mutants and their nearisogenic WT lines during initiation and early elongation up to 3 DPA [5,13]. In a fiber developmental study Kohel and co-authors observed that the elongation pattern is similar and restricted in both, Li 1 and Li 2 fibers [16]. However, unlike the normal morphological growth of the Li 2 plants, the Li 1 mutant exhibits pleiotropy in the form of severely stunted and deformed plants in both the homozygous dominant and heterozygous state [8,11,12]. The near-isogenic lines (NILs) of Li 1 and Li 2 with the elite Upland cotton variety DP5690 previously used in our research [5,8] provide an excellent model system to study mechanism of fiber elongation.
In our previous report we used a microarray approach to identify common genes related to fiber elongation, those with altered expression as a result of the Li 1 and Li 2 mutations, growing in the field and a greenhouse [7]. We found a relatively small number; 88 genes were differentially regulated in both short fiber mutants, which may be due to limitations of microarray technology. RNA-seq offers a larger dynamic range of quantification, reduced technical variability, and higher accuracy for distinguishing and quantifying expression levels of homeologous copies than microarray [17]. RNA-seq can provide a more comprehensive and accurate transcriptome analysis of cotton fiber development by using the reference genome sequence of Gossypium raimondii Ulbr. [18].
In this study we used a RNA-seq approach for the same goal: to determine fiber elongation related genes affected in both mutants growing in the field and a greenhouse. We found a larger number of differentially regulated genes common to both mutants, and from those the major intrinsic proteins were significantly over-represented among the down regulated genes. We measured the osmolality and concentrations of major osmotic solutes in sap of fiber cells. Although the osmolality and the concentrations of soluble sugars were less in saps of both short fiber mutants than in WT the concentrations of K + and malic acid were significantly higher in saps of mutants than in WT during rapid elongation time. The higher concentrations of malic acid and ions suggest limited uptake of water into fiber cells of mutants that can be result of down regulation of major intrinsic proteins.

Sources of variability in RNA-seq data
We examined genome-wide gene expression in elongating cotton fiber cells at 8 DPA in Li 1 , Li 2 mutants and WT under different growing conditions, in the field and greenhouse. The time point 8 DPA was selected because our earlier research revealed significant transcript and metabolite changes between the Li 2 and WT NILs during this time of fiber development [5,6]. Approximately 1.06 billion 100 bp reads from 13 libraries, including 9 libraries from field grown plants (this work) and 4 libraries from greenhouse grown plants (previously reported [9]), were trimmed with Sickle [19] and mapped to transcripts from the G. raimondii genome reference sequence [18,20]. The results of mapping reads are provided in Additional file 1.
Principal component analysis (PCA) was applied to explore relationships in gene expression among the samples. According to PCA, the samples from the near-isogenic lines and from the same lines growing in the field and a greenhouse are separated, indicating effects of the mutations and growth conditions on gene expression ( Figure 1A). To further investigate the proportion of variation in gene expression explained by each factor, a principal variance components analysis (PVCA) was run on the same data set. This approach first reduces data dimensionality with PCA, and then fits a mixed linear model to each principal component with variance components analysis (VCA). The largest source of variability in fiber transcriptome was the variance component L (the near-isogenic lines; weighted average proportion of 56.4%), whereas the variance component E (environmental factor) explained 13.8% of the total transcriptional variance ( Figure 1B).

Differential gene expression analysis
An ANOVA model for gene expression was specified in which the measured level of gene expression in Li 1 and Li 2 under different growth conditions was compared with gene expression in corresponding WT. The ANOVA analysis of transcript data is provided in Additional file 2. We found that 4,128 genes were significantly (FDR q-value < 0.05) up-regulated in field grown Li 1 fibers, whereas only 2,144 genes were up-regulated in field grown Li 2 fibers and 3,442 genes were up-regulated in greenhouse grown Li 2 fibers (Figure 2A). The largest amount of down-regulated genes 2,536 was detected in field grown Li 1 fibers, whereas 1,740 and 1,914 genes were down-regulated in field and greenhouse grown Li 2 fibers, consequently. Only small portions of these genes were common among up-regulated (531) and down-regulated (652) in all tested conditions by ANOVA model (Figure 2A). In the following gene set enrichment analysis we focused only on these common genes since our objective was to identify fiber elongation related genes common between short fiber mutants growing in the field and a greenhouse.
MapMan ontology was used for gene set enrichment analysis [21]. Two main categories (electron transport and transport) were overrepresented among up-regulated genes and five main categories (transport, enzyme families, cell wall, cell and development) were overrepresented among down-regulated genes in Li 1 -Li 2 developing fibers. Figure 2B shows only sub-categories from the above mentioned main categories which are significantly (Chi-square, p < 0.05) overrepresented in the Li 1 -Li 2 fiber transcriptomes. Particularly, NADH dehydrogenase, cytochrome c and alternative oxidase were significantly (p < 0.0001) overrepresented subcategories in electron transport, whereas ABC transporters and transport of amino acids were overrepresented subcategories Li 1 -Li 2 up-regulated genes. The most significantly (p < 0.0001) overrepresented sub-categories in Li 1 -Li 2 down-regulated genes were: major intrinsic proteins and transport of sulphate in transport category; and the plastocyanin -like enzyme family.
Genes categorized into transport functional category were overrepresented among up-regulated and downregulated pools of genes; however, proportions of gene family members of transporters were different among up-regulated or down-regulated genes. Significantly upregulated and down-regulated transporters in Li 1 -Li 2 mutants growing in the field and a greenhouse are shown in Tables 1 and 2. Major intrinsic proteins, sulphate and phosphate transporters were present only among pool of down-regulated genes, whereas proportions of amino acids and ABC transporters were significantly higher among pool of up-regulated genes. The sugars transporters were not significantly more abundant among up-regulated than down-regulated genes.

Major intrinsic proteins
Major intrinsic proteins or aquaporins were one of the most significantly (p < 0.0001) over-represented gene family among down-regulated genes in Li 1 -Li 2 fibers. Aquaporins facilitate the efficient transport of water and other small molecules across membranes in plants and other organisms [22]. Cotton aquaporins form a large family of proteins phylogenetically divided into five subfamilies including: plasma membrane intrinsic proteins (PIP), tonoplast intrinsic proteins (TIP), NOD26-like intrinsic proteins (NIP), small basic intrinsic proteins (SIP), and the recently identified X (or unrecognized) intrinsic proteins (XIP) [23]. To assess which subfamily members of aquaporins were affected by Li 1 -Li 2 mutations: first, we conducted phylogenetic analysis of G. raimondii genes annotated as aquaporins; and second, evaluated their expression level in Li 1 -Li 2 developing fibers. The analyzed G. raimondii aquaporins clustered into five main clades (marked by empty squares) representing the above mentioned subfamilies (Additional file 3). The members of subfamilies PIP (7 genes), TIP (4 genes) and NIP (2 genes) were down-regulated in Li 1 -Li 2 developing fibers (marked by black triangle in Additional file 3). The most highly induced aquaporins in WT fibers, for which transcript levels were dramatically reduced in Li 1 -Li 2 mutants, were tested by RT-qPCR. In most cases results of RT-qPCR analysis were consistent with results of RNA-seq analysis ( Figure 3). There were a number of aquaporins which showed increased transcript level only in greenhouse grown Li 2 (Additional file 4), indicating interactive response to Li 2 mutation and growth conditions. However, relative   expression level of those genes was considerably less compared with WT expressed aquaporins as shown in Figure 3 (1,500 reads in greenhouse Li 2 induced vs. 500,000 reads in WT expressed).

Osmotic concentrations and solutes in saps of Li 1 and Li 2 fiber cells
We measured the osmotic concentration and calculated osmotic pressure of the sap of cotton fiber cells. The sap solution represents the average osmotic concentration of the vacuole, the cytoplasm, and the apoplast (i.e. freespace solution) of the fiber cells. In fiber cells the vacuole occupies approximately 90% of the cell volume [4]; therefore the measured osmotic concentration values largely represent the solute concentration of the vacuoles. The calculated osmotic pressure in sap of WT fibers was steadily high during rapid fiber elongation, at 3 -16 DPA, and significantly dropped during the transition to the cell wall biosynthesis stage (Figure 4). The pattern of osmotic pressure in sap of Li 1 fibers was similar with pattern in WT; although the osmotic pressure was significantly lower (p < 0.05) at 3 -8 DPA. In sap of Li 2 fibers the osmotic pressure was significantly lower than in WT at 3 -5 DPA, but higher at 24 DPA. Soluble sugars, K + , and malate are major active solutes in elongating fibers, to which are often attributed 80% of the fiber sap osmolality [4,24,25]. To assess which osmotic solutes altered in the Li 1 and Li 2 developing fibers we measured the concentrations of sugars, malic acid, and ions in fiber sap solutions ( Figure 5). Concentrations of hexoses (D-glucose and D-fructose) were significantly less in sap of Li 1 and Li 2 fibers compared to WT during rapid fiber cell expansion (at 5 -16 DPA). The level of sucrose was low during elongation at 3 -16 DPA in sap of all near-isogenic lines; however, at 20 -24 DPA the concentration of sucrose significantly increased in Li 1 and Li 2 , but not in WT fiber. Surprisingly, the concentrations of malic acid and K + were significantly (p < 0.001) higher in sap of Li 1 and Li 2 fibers comparing to WT during elongation ( Figure 5). The concentrations of Na + were not significantly different in saps of Li 1 , Li 2 and WT. We also measured the concentrations of Ca +2 and phosphorus, which were significantly higher in saps of mutants compared to WT.

Discussion
Experimental design for identification of fiber elongation related genes In this study we compared the transcriptomes of developing fibers of two short fiber mutants and their WT NIL growing in the field and a greenhouse. The mutated genes of the Li 1 and the Li 2 are yet to be discovered. A defect in the Li 1 gene affected a number of traits (dwarf deformed plants and short fiber phenotype), while the defect in Li 2 gene affected only fiber length. Therefore, the Li 1 and Li 2 , most likely, are different types of genes; their alterations interrupt different parts of a complex biosynthetic process, but in both cases cause a short fiber phenotype. Both Li 1 and Li 2 mutations have an enormous effect on the fiber transcriptomes; the largest source of variability in the fiber transcriptome data was due to mutations (56.4%; Figure 1B). However, altered expression of many genes in Li 1 -Li 2 transcriptomes can be result of chain-reactions to adverse effects of the causative mutation, and is not necessary directly related to fiber elongation process. Also it is known that many fiber-related genes are environmentally regulated [26]; in our experiment the environmental factor contributed 13.8% to the data variability ( Figure 1B). Therefore, to reduce noise in the data we selected common regulated genes between Li 1 /wt and Li 2 /wt grown in the field and Li 2 /wt grown in a greenhouse. This approach allowed the identification of transcripts directly related to fiber elongation process regardless of far downstream effects of the mutations and environmental conditions.

Gene set enrichment analysis
We found a large number of differentially expressed genes common to both mutants (Figure 2A). To gain insight into biological processes altered by Li 1 -Li 2 mutations we used MapMan ontology for gene set enrichment analysis. Consistent with our previous microarray study, mitochondrial electron transport functional category was over-represented among up-regulated genes in short fiber mutants [7]. Enrichment of the cell wall functional category was expected among down-regulated genes and described for Li 1 and Li 2 in our previous reports [5][6][7][8]. However, strong downregulation of major intrinsic proteins in short fiber mutants was not noticed before in our microarray studies, probably due to limitations of microarray techniques. Here, we found  that the major intrinsic proteins were the most downregulated gene family in both short fiber mutants; their role in osmoregulation of Li 1 -Li 2 fibers is discussed below.

Osmoregulation in short fiber mutants
The rapid expansion of fiber cells requires high turgor pressure and cell wall relaxation [4,25,27]. The force of turgor pressure is related to the osmotic potential and to the transport coefficient for water uptake [28]. The maintenance of sufficient osmoticum to compensate for dilution effects resulting from the influx of water is an important component of sustainable cell expansion [27].
In the fiber sap of short fiber mutants we detected significantly lower osmotic pressure than in WT. The reduced osmotic pressure in Li 1 -Li 2 may not be sufficient to maintain rapid and sustainable cell expansion and may cause short fiber phenotype. Soluble sugars, K + and malic acid are considered as major active solutes in rapidly expanding fiber cells [4,24,25]. We detected lower concentrations of glucose and fructose in sap of short fiber mutants than in WT that correlate with lower osmotic pressure, suggesting sugars are the main solutes to positively impact turgor in fiber cells. Sucrose was almost undetectable in mutants and WT fibers during the rapid elongation phase (3)(4)(5)(6)(7)(8)(9)(10)(11)(12)(13)(14)(15)(16). In developing fiber cells, sucrose is degraded into hexoses by sucrose synthase in the cytoplasm and acid invertase in the vacuole [24,29,30].
We tested the expression levels of sugars transporters in mutants because their regulation may cause a reduced supply of sugars in developing fibers. However, the number of up-regulated sugars transporters in Li 1 -Li 2 was higher than down regulated: 4 versus 2 genes, correspondingly (Tables 1 and 2). Therefore, the transport of sugars is unlikely altered in short fiber mutants. In our previous report we observed significant reductions in the levels of detected free sugars, sugar alcohols, sugar acids, and sugar phosphates in the Li 2 metabolome; also biological processes associated with carbohydrate biosynthesis were significant down-regulated in the Li 2 transcriptome [6]. Consequently, detection of low amount of sugars in sap of Li 1 -Li 2 fibers might be the result of reduced de novo synthesis of sugars in mutants. The driving force for the transport and accumulation of ions into the protoplast and vacuole is provided by the plasma membrane and vacuolar H + -ATPases [27,31]. We did not detect the plasma membrane and vacuolar H + -ATPases among common Li 1 -Li 2 up-regulated or down-regulated pools of genes. Numbers of calcium, potassium and other metal transporters were not significantly different between pools of up-regulated and down-regulated genes in short fiber mutants; except for sulphate and phosphate transporters which were present among down-regulated genes only (Tables 1 and 2). Thus, ion transport in Li 1 -Li 2 is unlikely to be affected  by the mutations and proceeds normally as in wild type plants. The higher concentrations of malic acid, K + and other inorganic ions detected in sap of Li 1 -Li 2 can be explained by reduced influx of water into fiber cells of mutants ( Figure 5). Since malic acid and K + (major osmotic solutes) cannot restore the balance of water uptake into developing Li 1 -Li 2 fibers, there is another factor, which might be crucial for osmoregulation of cotton fibersthe major intrinsic proteins ( Figure 6). The major intrinsic proteins or aquaporins were the most overrepresented gene family among down-regulated genes in both short fiber mutants ( Table 2). The expression level of some members of PIPs and TIPs at 8 DPA of fiber development in WT was enormous, up to 500,000 reads ( Figure 3). It has been indicated in a number of studies that the osmotic water permeability (or hydraulic conductivity) is controlled by the activity of aquaporins. For instance, Javot and co-authors showed that Arabidopsis PIP2;2 is highly expressed in several root cell types, and that, by comparison to WT plants, the hydraulic conductivity of corresponding knock-out mutants (pip2;2) was reduced by 14% [32]. The hydraulic conductivity of pip1;2 mutants and pip2;1 and pip2;2 double mutants was decreased by 20% and 40% respectively, compared to that of WT [33,34]. A link between aquaporins and cell growth has also been shown in different species. Virus-induced silencing of rose PIP2;1 resulted in a reduction in size of cells and petal expansion [35]. Over-expression of a cauliflower TIP1-GFP fusion in tobacco suspension cells or of ginseng TIP in Arabidopsis leaves led to an increase in cell size [36,37]. Vacuole regeneration and cell expansion were accelerated in protoplast prepared from BY-2 cells overexpressing the NtTIP1;1 [38]. Knockdown of expression of GhPIP2 genes by RNA interference in G. hirsutum markedly inhibited fiber elongation [39]. Thus, the reduced expression of aquaporins in short fiber mutants may reduce the influx of water into fiber cells and slow down the elongation process ( Figure 6).

Conclusions
Here, we used an RNA-seq approach to determine common fiber elongation related genes in developing fibers of Li 1 and Li 2 mutants growing in the field and a greenhouse. We found that the aquaporins were the most down-regulated gene family in both short fiber mutants. The osmolality and concentrations of soluble sugars were less in saps of Li 1 -Li 2 , whereas the concentrations of malic acid, K + and other detected ions were significantly higher in saps of mutants than in WT. These results suggest that higher accumulation of ions in fiber cells, reduced osmotic pressure and low expression of aquaporins, may contribute to the cessation of fiber elongation in Li 1 and Li 2 short-fiber mutants.

Plant materials
Two mutant lines Li 1 and Li 2 in a near-isogenic state with the WT upland cotton line DP5690 were developed in a backcross program at Stoneville, MS as described before [5,8]. The growing period for the greenhouse grown Li 2 plants was between October, 2009 and March, 2010; planting and growth conditions were previously described [5]. For the field grown plants, a total of 150 Li 1 , 100 Li 2 , and 100 WT plants were grown in a field at the USDA-ARS Southern Regional Research Center, New Orleans, LA in the summer of 2013. All samples of the same developmental stage were tagged and collected on the same day. Cotton bolls were harvested at 3, 5, 8, 12, 16, 20, and 24 DPA. Bolls were randomly separated into 3 replicates with 15-30 bolls per replicate.

RNA isolation and reverse transcription quantitative polymerase chain reaction (RT-qPCR)
Total RNA was isolated from detached fibers [40] using the Sigma Spectrum Plant Total RNA Kit (Sigma-Aldrich, St. Louis, MO) with the optional on column DNase1 digestion according to the manufacturer's protocol. The concentration of each RNA sample was determined using a Nano-Drop 2000 spectrophotometer (NanoDrop Technologies Inc., Wilmington, DE). The RNA quality for each sample was determined by RNA integrity number (RIN) using an Agilent Bioanalyzer 2100 and the RNA 6000 Nano Kit Chip (Agilent Technologies Inc., Santa Clara, CA) with 250 ng of total RNA per sample. RNA from each of the above mentioned time-points was used for RT-qPCR analysis. A detailed description of reverse transcription, qPCR and expression analysis was previously reported [9]. Sequences of primers used for qPCR are listed in Additional file 5.  Deep sequencing and differential gene expression Library preparation and sequencing were performed by Data2Bio LLC (Roy J. Carver Co-Laboratory, Ames, Iowa). The libraries were sequenced using 101 cycles of chemistry and imaging, resulting in paired end (PE) sequencing reads with length of 2 × 101 bp. For the greenhouse grown Li 2 plants samples were sequenced in two biological replicates (sequencing and data were described elsewhere [9]). For the field grown Li 1 , Li 2 and WT plants RNA samples from cotton fiber at 8 DPA were sequenced in three biological replicates. The 8 DPA was chosen because it is the peak of fiber elongation phase according to our earlier studies [5,6,8]. RNA-seq expression analysis was conducted following the PolyCat pipeline [20]. Briefly, all reads were aligned to the JGI Gossypium raimondii reference genome [18], then the PolyCat software assigned each categorizable read to either the A or D subgenome based on an index of homeoSNPs. We followed two adjustments previously described [10], particularly: 1) we only counted exonic reads; 2) we used the ratio of Aassigned to D-assigned reads to proportionally divide the total number of mapped reads for each gene which ensures that unassigned reads contribute to the total expression of genes. The data normalization and ANOVA process were conducted as previously described [9]. Principal component analysis (PCA) was conducted using JMP Genomics 7 software (SAS Institute Inc., Cary, NC, USA). Transcript data for each sample were used as continuous variables.
Principal variance component analysis (PVCA) reduces the dimensionality of the data set with PCA, and then fits a mixed linear model to each principal component to partition variability with variance components analysis (VCA). A summary of variance components across all principal components is constructed as a weighted average of the individual estimates, using eigenvalues as weights [41].

Osmotic concentration measurement of cotton fiber saps
Saps were collected from cotton fiber cells from the each time point mentioned above. From short fiber ovules fibers were collected by shaking ovules frozen in liquid nitrogen. Fibers from long fiber ovules were pulled off by forceps. For each replicate, about 200 mg of frozen fibers were thawed on ice. Sap was separated by centrifugation at 5,000 g for 1 min [30]. The osmotic concentration of cotton fiber sap was measured using a vapour pressure osmometer, VAPRO-5600 (WESCOR INC., South Logan, UT, USA) in three biological replicates as previously described [30]. The osmotic pressure (MPa) was calculated from the osmotic concentration using the equivalence 2.44 MPa per 1 Osmol•kg −1 (25°C) [42].

Measurement of sugars and malic acid concentrations in cotton fiber saps
Concentrations of sucrose, D-fructose and D-glucose in cotton fiber saps were measured using enzyme assay kit