Fiber cell initiation in GhTM-1, Gb3-79, and GhMD17
Our study involved two cultivated cottons, G. hirsutum (GhTM-1) and G. barbadense (Gb3-79), and the fiber-less mutant GhMD17 (G. hirsutum background). These three cotton lines were chosen on the basis of their distinct fiber characteristics: GhTM-1 has both lint and fuzz fibers, Gb3-79 has almost exclusively lint fibers, and GhMD17 is fiberless (Figure 1A). Both GhTM-1 and Gb3-79 have similar lint fiber properties in mature seeds, except for the lack of fuzz fiber on the Gb3-79 seed. During fiber initiation stage (0 and 2 DPA), both GhTM-1 and Gb3-79 had fibers of ~10 mm in length at 0 DPA and ~50 mm at 2 DPA, with no significant difference between GhTM-1 and Gb3-79, although fiber growth was slightly faster in Gb3-79 than in GhTM-1 (Figure 1C). GhMD17 lacked any visible fiber cell initials during those stages (Figure 1B). Thus, GhTM-1 and Gb3-79 have a similar fiber initiation process, which is blocked in GhMD17.
RNA-seq analysis of fiber and epidermal cells separated by LCM
In wild-type ovules at 0 DPA and fiber-less mutants at both 0 and 2 DPA, LCM was used to cut the top layer of epidermal cells. For wild-type ovules at 2 DPA, LCM was used to cut along the epidermal surface to collect the protruding fibers, followed by a second cut collecting the underlaying epidermal cells (Figure 1D and Supplemental Figure 1A). LCM materials in each sample with two biological replicates were used for total RNA extraction (Supplemental Figure 1B), library construction, and mRNA sequencing. Sequencing reads were mapped onto the G. hirsutism genome v2.0 (GhTM-1 and GhMD17) and G. barbadense genome v1.0 (Gb3-79)  (Supplemental Table 1), and read counts were normalized to fragments per kilobase per million (FPKM) for gene expression analysis (see Methods).
Principal component analysis (PCA) showed a clear separation of gene expression variance by tissue types (PC1: 37% of variance), followed by species (PC2: 23% of variance) (Figure 2A). Epidermal samples between the 0 and 2 DPA stages were not separated. We also determined the Pearson’s correlation coefficient of all 8 samples (averaged values of two biological replicates) and performed hierarchical clustering analysis among these samples (Supplemental Figure 2A). This revealed three major groups: fiber samples (GhTM-1 and Gb3-79 at 2 DPA), mutant epidermal samples (GhMD17 at 0 and 2 DPA), and Gb3-79 epidermal samples (0 and 2 DPA). GhTM-1 epidermal sample at 0 DPA was grouped with Gb3-79 epidermal samples, while GhTM-1 epidermal sample at 2 DPA and mutant epidermal samples were in the same group. These results indicate that the gene expression differences in the ovules at early fiber stages between the tissue types and between species were suitable for this study.
At 0 DPA, the fiber cell initials develop in the wild type but not in the mutant. To identify differentially expressed genes in fibers at 0 DPA, we compared expression differences in the epidermal samples between the wild type (GhTM-1 and Gb3-79) and mutant (GhMD17). We identified 1,295 and 1,500 upregulated genes in GhTM-1 and Gb3-79 fiber cell initials, respectively, 981 of which (>60%) were shared in both, while 314 and 519 were specific to GhTM-1 and Gb3-79, respectively (Figure 2B, Supplemental Table 2 and 4). In addition, 2,101 and 2,582 genes were down-regulated in GhTM-1 and Gb3-79 fiber cell initials, respectively, 1,528 of which (>60%) were shared in both (Supplemental Figure 2B, Supplemental Table 2).
At 2 DPA, we compared gene expression of LCM fiber cell with epidermal cell tissues within the same genotype (GhTM-1 or Gb3-79) and between the fiber cell (GhTM-1 and Gb3-79) and epidermal cell tissues in the fiber-less mutant (GhMD17). We identified 1,861 and 2,046 up-regulated genes in GhTM-1 and Gb3-79 fiber cells, respectively, 743 of which (>35%) were shared in both, while 265 and 450 were unique to GhTM-1 and Gb3-79, respectively (Figure 2C, Supplemental Table 3 and 4). We also identified 2,978 and 4,269 down-regulated genes in GhTM-1 and Gb3-79 fiber cells, respectively, and 1,626 of which (>35%) were shared in both (Supplemental Figure 2C, Supplemental Table 3). These differentially expressed genes in the fiber cell tissues were used for further analysis.
Dramatic gene expression changes in early stages of cotton fiber cell development
Gene Ontology (GO) analysis of these differentially expressed genes showed overrepresentation of translation (ribosome activity), gene expression, and peptide, lipid, and fatty acid metabolic processes in the upregulated genes in fibers (Figure 3). These results may indicate that ribosome biosynthesis, required for protein synthesis as well as lipid and fatty acid metabolism, were active during fiber cell initiation. Moreover, the upregulated genes between the species were overrepresented in ATP and purine metabolic processes in GhTM-1, and ADP and pyridine (similar structure of pyrimidine) metabolic processes were overrepresented in Gb3-79. ATP is a critical energy resource for fiber growth. In cultured cotton ovules, fibers release ATP as they grow, and when the ectoapyrase activity is blocked by the addition of polyclonal anti-apyrase antibodies or small molecule inhibitors, the medium ATP level rises and fiber growth is suppressed; low concentrations of hydrolyzable nucleotides ATPgammaS/ADPbetaS stimulate fiber growth . Thus, the difference of ATP and ADP metabolic processes between GhTM-1 and Gb3-79 may contribute to longer fibers in Gb3-79.
Among down-regulated genes in fibers, they were overrepresented in response to hormones and glucose metabolic process at 0 DPA and DNA replication and the cell cycle at 2 DPA. These data may suggest that cell cycle and DNA replication have attenuated in fiber cells, which require active ribosome biosynthesis for translational regulation in fiber cells. Down-regulation of phytohormone-related genes in the fiber cells may suggest upregulation of these genes in ovules to support fiber elongation.
Allotetraploid cottons consist of A and D subgenomes (A and D) with 52,514 orthogroups representing 15,220 pairs of A and D homeologs in G. hirsutum and G. barbadense (Supplemental Table 5). Among them, ~20% of homeolog pairs were differentially expressed (Supplemental Figure 3A), with 10% more D homoeologs (1,441-1,550) expressed at higher levels (D>A) than the A homoeologs (1,250-1,402, A>D) in both fiber and epidermal cells (Supplemental Figure 3B, 3C). These data are consistent with the notion of dramatic gene expression divergence between the homoeologs, despite gene numbers and synteny are conserved between the A and D subgenomes in the allotetraploid cottons . This trend of biased expression is shared by 251 and 81 homeolog pairs that encode cell cycle regulation genes and ribosomal protein subunit genes, respectively, with expression bias for ~20% of homeolog pairs involved in cell cycle regulation and for ~15% of homoeologs encoding ribosomal protein subunits (Supplemental Figure 3D). Interestingly, for cell cycle regulation genes, 10% more D than A homoelogs had biased expression, whereas this number of homoeologs with expression bias was increased to 40-50% more in the D than in the A homoeologs for the ribosomal protein subunit genes (Supplemental Figure 3D). This high number of homoeologs with expression bias in the early stages of fiber cells may indicate that deep subfunctionalization of homoeologous genes during active developmental process of fiber cell initiation and elongation.
Translational regulation in early cotton fiber cell development
During fiber cell development, the cell cycle is predicted to attenuate in 20-30% of the epidermal cells, and those cells become fiber cell initials. Cell cycle arrest in a cell often leads to active ribosome biosynthesis . This is consistent with developmental process of fiber cell initiation, when the arrest of the cell cycle in the epidermal cells induces active ribosome biosynthesis for fiber cell initiation and elongation.
Arabidopsis has 231 annotated genes for cell cycle regulation, and cotton has a total of 845 putative cell-cycle regulation genes, including 366 involved in G1, 283 in G2, 122 in S, and 74 in M stages (Supplemental Table 6). Approximately 85% of all putative cell-cycle regulation genes were down-regulated in fiber cells (Figure 4A and Supplemental Figure 4A). There is a total of annotated genes in cotton for 681 ribosome protein subunits, consisting of large subunit, RPL and small subunit RPS, respectively (Supplemental Table 6). These multiple putative ribosomal subunit genes (618 genes) are grouped into 121 putative ribosomal proteins. In contrast to down-regulation of cell cycle related genes, approximately 75% of putative ribosomal protein genes were up-regulated, while 8% were down-regulated (Supplemental Figure 4B). Plants have more ribosomal protein genes than animals or yeasts , and many of those paralogues are differentially expressed during development and in UV-B responses . For example, in ArabidopsisRPL16A and RPL16B are expressed in developing pollen and early phase of lateral root initiation , a process similar to fiber cell initiation. The AtRPS5A gene is strongly expressed in dividing cells, whereas AtRPS5B expression correlates with cell differentiation rather than cell division . These results suggest that during fiber cell initiation, the cell cycle was attenuated, and ribosomal protein biosynthesis is activated and may regulate tissue-specific expression in the fiber and epidermal cells.
The cell cycle is regulated by Target of Rapamycin Complex1 (TORC1), which consists of three proteins TARGET OF RAPAMYCIN (TOR), RAPTOR, and LETHAL WITH SEC THIRTEEN PROTEIN 8 (LST8) [44, 45]. ErbB-3 EPIDERMAL GROWTH FACTOR BINDING PROTEIN1 (EBP1) is a TOR target. TORC1 is also a key regulator of many other biological processes such as DNA and RNA synthesis and autophagy and is conserved among plants and animals . In cotton, TOR, RAPATOR, LST8, and EBP1 have multiple members, two homoeolog pairs each for TOR, RAPTOR, and EBP1, and one homoeolog pair for LST8. The expression levels among multiple members of each gene were similar and averaged for the analysis (Supplemental Table 7). In fiber cells at 2 DPA, TOR was down-regulated in both GhTM-1 and Gb3-79 (Figure 4B), while RAPTOR and LST8 expression levels were similar to epidermal cells, and EBP1 was slightly up-regulated. In Arabidopsis, TOR is expressed in the primary meristem, embryo, and endosperm, but not in differentiated cells , while the mutant is embryonically lethal. Notably, the putative TOR target EBP1 ribosomal protein biosynthesis in humans and may regulated cell expansion and size in plants . Thus, in fiber cells down-regulation of TOR may inhibit the cell cycle, and strong expression of EBP1 could promote ribosome biosynthesis for fiber cell expansion.
To further test roles of ribosome biosynthesis and cell cycle regulation in fiber cell development, we used pharmacological assays in cultured GhTM-1 and GhMD17 ovules. Ribozinoindoles (Rbin-1) are potent and reversible triazinoindole-based inhibitors of eukaryotic ribosome biogenesis , and Roscovitine is a cyclin-dependent kinase inhibitor for cell cycles . Cotton ovules were isolated from fertilized flowers of GhTM-1 and GhMD17 at 0 DPA and incubated in BT medium for 21 days following a published protocol , with or without inhibitors. Inhibiting cell cycle by Roscovitine (2 mg/L) increased total fiber area and length at statistically significant levels, whereas impeding ribosome biosynthesis by Rbin-1 (40 mg/L) reduced the fiber area (Figure 4C, 4D), and the ovule size was not affected by either treatment (Figure 4D). Indeed, chemical inhibitors had no effects on ovule development in the GhMD17 that is fiberless (Figure 4C left). The altered fiber area could result from increased fiber length or number or both. To discern this, we further measured the number or density of fiber cell initials after two days of each treatment (Supplemental Figure 4). The result showed a similar trend of increasing and reducing fiber cell density by Roscovitine and Rbin-1 treatments, respectively. These results support the hypothesis that cell cycle attenuation and active ribosome biosynthesis promote fiber cell initiation and elongation.
Dynamic roles of phytohormones and MYB transcription factors in early fiber cells
Phytohormones play important roles in fiber cell development [2, 24], but the temporal regulation of gene expression in early fiber developmental stages is poorly understood. Using LCM, we unexpectedly uncovered genome-wide down-regulation of 43 phytohormone related genes with only a few up-regulated genes in early fiber cells (Figure 5A and Supplemental Table 8). This is probably because most genes are upregulated in the epidermal layer cells that support fiber cells. Abscisic acid (ABA) is known to inhibit fiber cell initiation and elongation . ABA in ovules at 0 DPA accumulates at much higher levels in the Ligon-lintless 1 mutant than in the wild type . This is consistent with the down-regulation of key regulators in the ABA pathway (Figure 4A), including phosphate type-2C HAB1 (a homeolog of ABA INSENSITIVE 1, ABI1; and ABI2)  and ABA receptor gene PYRABACTIN RESISTANCE-LIKE PROTEIN 9 (PYL9) (known to interact with ABI1 and ABI2) . Down-regulation of the ABA regulator and receptor genes may promote fiber cell initiation and elongation.
Gibberellins promote fiber cell elongation and increase ABA and auxin levels during fiber development to increase fiber quality (strength) . Over-expression of gibberellin biosynthesis gene GIBBERELLIN 20-OXIDASE1 (GA20OX1) leads to more and longer fibers , while their homeologs GA20OX2 and GA20OX3 are largely expressed in the ovule. Consistently, GA20OX2 was upregulated in the ovule and down-regulated in fibers (Figure 4A).
Ethylene regulates fruit ripening and inhibits vegetative growth [57, 58], and is reported to promote fiber elongation in cotton, which is mediated by very-long-chain fatty acid signaling pathway . However, the mechanism for ethylene in fiber initiation remains unclear. We found down-regulating of the ethylene signaling pathway genes like ETHYLENE INSENSITIVE 2 (EIN2) and EIN3 in fibers, which may suggest a negative role for ethylene in fiber cell initiation (Figure 4A).
Auxin promotes fiber initiation and accumulates in ovules prior to fertilization . Auxin is transported from ovule to fiber cells during fiber cell development . Consistent with this notion, most genes encoding AUXIN RESPONSE FACTORs (ARFs) and PIN-FORMED6 (PIN6), a key auxin transporter, were upregulated in the ovules and down-regulated in early fiber cells (Figure 4A). This is reminiscent of the finding that ARF2, ARF18 and PIN6 were expressed in ovules in the fiber initiation stage . Interestingly, although the majority of among ARF gene family members are upregulated in the ovules, a few such as ARF19a and ARF19b homoeolog pair, showed opposite expression patterns in the fiber cells (Supplemental Figure 6), suggesting functional divergence between homoeologs for fiber cell development in polyploid cotton.
Brassinosteroids promote both fiber cell initiation and elongation [27, 59]. A key regulator of brassinosteroid signaling, BRASSINAZOLE-RESISTANT 1 (BZR1), induces cotton fiber development . BSU1 is a dominant suppressor of BRI1 and is expressed in elongating cells in Arabidopsis . In cotton fiber cells, BSU1 was upregulated, whereas BZR1 and its semidominant suppressor gene (BES1)  were down-regulated (Figure 4A), suggesting a role for brassinosteroids in fiber cell initiation.
In addition to phytohormones, MYB transcription factors regulate fiber cell development in seeds and trichome development in leaves [36, 37, 63]. Among 18 MYB transcription factor genes examined, nine were upregulated, and the other nine were down-regulated in early fiber cells (Figure 4B and Supplemental Table 8). The upregulated genes included GhMYB25, GhMYB109, and cotton homologs of Arabidopsis MYB4, MYB66, MYB96, and MYB106 (Supplemental Table 8). Notably, cotton homologs of MYB106 (GhMYB25 in cotton) were up-regulated, while cotton MYB16 homologs were down-regulated in fibers (Figure 4B), while Arabidopsis MYB16 and MYB106 are in the same clade (Supplemental Figure 7). Similarly, cotton MYB4 was upregulated, and MYB6 was down-regulated in early fiber cells (Figure 4B and Supplemental Figure 4), but they are closed related MYB genes. These results may suggest functional divergence between closely related homologous genes during fiber cell development.