LCM and RNA-seq analyses revealed roles of cell cycle and translational regulation in early stages of cotton ber cell development

Cotton bers provide a powerful model for studying cell differentiation and elongation. Each cotton ber is a singular and elongated cell derived from epidermal-layer cells of a cotton seed. Efforts to understand this dramatic developmental shift have been impeded by the diculty of isolating ber cells from epidermal cells. Here we employed laser-capture microdissection (LCM) to separate these cell types. RNA-seq analysis revealed transitional differences between the ber and epidermal-layer cells at 0 or 2 days post anthesis. Specically, down-regulation of putative cell cycle genes was coupled with upregulation of ribosome biosynthesis and translation-related genes, which may suggest their respective roles in ber cell initiation and elongation. Indeed, the amount of bers in cultured ovules was increased by cell cycle progression inhibitor, Roscovitine, and decreased by ribosome biosynthesis inhibitor, Rbin-1. Moreover, many phytohormone-related genes were upregulated in the ovules and down-regulated in the bers, suggesting their spatial-temporal effects on ber cell development. Key cell cycle regulators were predicted to be epialleles, and MYB-transcription factor related genes displayed expression divergence between bers and ovules, implying their effects on ber traits. breeding


Background
Cotton is the largest renewable source of the textile ber in the world and is an important oil crop [1,2].
Most widely cultivated cottons are allotetraploid, consisting of two sets of chromosomes from different origins across the old and new world, with the "mother" (A-genome like species) from Africa and the "father" (D-genome like species) from the Americas [3]. The polyploidization took place in the new world 1-1.5 million years ago [4], and resulted in ve extant species that diversi ed over 300,000-600,000 years [4,5]. Over the last 8,000 years, two of them, Upland (Gossypium hirsutum L.) and Pima (G. barbadense L.) cottons were independently domesticated in NW S. America and the Yucatan Peninsula of Mexico, respectively, under strong human selection, leading to the modern annualized crops [6]. To date, Upland cotton dominates in ~ 94% of total cotton production, with the remainder (~5%) being primarily Pima (with a superior ber quality) [1], and other diploid cottons that are adapted to certain growth environments [7].
Cotton ber is seed hair, and each individual cotton ber initiates from a single epidermal-like cell from the ovular epidermis, which elongates dramatically to reach up to ~5 cm in length, one of the longest singular cells in the plant kingdom [2]. Each cotton seed has over 20,000 semi-synchronically developed ber cells, representing 20-30% of cells in the epidermal layer of the seed [8]. There are two types of ber, lint and fuzz bers [9,10], and only the lint bers are ginned and used in textiles. Lint ber cells are initiated before or on the day of anthesis (0 Day Post Anthesis, DPA), which induces production and transportation of the phytohormones auxin and gibberellin to promote ber cell elongation [11]. The fuzz ber is produced approximately 5 days after the development of lint bers and is much shorter than the lint ber. The ber cells undergo overlapping stages of elongation, cellulose biosynthesis, and maturation [12]. During the elongation phase, primary cell wall is synthesized, and the cell extends through some combination of anisotropic expansion and tip growth [10]. Once bers have reached full length, they transition to the stage of cellulose and secondary cell wall biosynthesis, which is accompanied by the maturation phase where bers, now comprising ~95% pure cellulose, become metabolically inactive, leading to cell death, which happens before the boll opens [10].
In developing ber cells, the cell cycle is arrested at G1 or S phase [13]. Cell differentiation and cell cycle regulation are linked in plants [14], which involves cyclin-dependent kinases (CDKs) [15], phytohormones (auxin and gibberellins), and key transcription factors. Cell cycle exit without mitotic division often leads to endoreduplication, which is common in eudicot plants [16]. Previous studies showed that DNA content was increased by endoreplication during early ber development [17,18], which is accompanied with small RNA and DNA methylation changes in early ber cells [19][20][21]. In human cancer cells, tumor suppressor p53 has been shown to induce cell-cycle arrest in response to impaired ribosome biogenesis [22], and vice versa, cell cycle arrest during tumor metastasis is fueled by upregulation of ribosome biogenesis [23].
Phytohormones and transcription factors also play crucial roles in ber cell development [2,[24][25][26]. Auxin, brassinosteroids, ethylene, and gibberellins promote or induce ber cell development [27][28][29], whereas cytokinin and abscisic acid inhibit ber cell development [30]. Previous studies have shown that auxin accumulates in ovules before fertilization, peaking around 2-3 DPA, with the level gradually decreasing through the elongation stage [11], while abscisic acid gradually accumulates in ber cells during the initiation and elongation stages, with an accumulation peak around 10 DPA, followed by a decrease in later stages [31].
MYB transcription factors are required for the development of leaf trichomes and seed hair or ber cells [32][33][34][35]. Several MYB transcription factor families are shown to regulate ber cell initiation and elongation [33][34][35][36], while some MYB transcription factor members do not have an obvious role in cotton ber cell development [37].
Regulatory networks involving phytohormones and transcription factors during the early stages of ber cell development remain poorly understood, partly because it is technically di cult to isolate ber and epidermal cells during ber cell initiation. Here we performed RNA-seq analysis using the ber and epidermal cells that were separated by Laser Capture Microdissection (LCM) in early stages (0 and 2 DPA). We compared transcriptomes of LCM samples between the wild-type and ber-less mutant (both in the Upland background) and between Upland and Pima cottons, which revealed expression differences in the genes associated with cell cycle regulation and ribosome activity. Using pharmacological inhibitors for these pathways, we showed that cell cycle regulation and ribosome biogenesis affect ber cell initiation and elongation in the ovules cultured in vitro. Moreover, phytohormone related genes were upregulated in the ovules and down-regulated in early bers, suggesting a spatial-temporal role in ber cell development. Key cell cycle regulators were predicted to be epialleles, and MYB-transcription factor related genes displayed expression divergence between bers and ovules. Together, these results support a new role for cell cycle and translational regulation in ber cell initiation and elongation, which is accompanied by spatial-temporal regulation of transcriptional and phytohormonal networks.
Thus, GhTM-1 and Gb3-79 have a similar ber initiation process, which is blocked in GhMD17.

RNA-seq analysis of ber and epidermal cells separated by LCM
In wild-type ovules at 0 DPA and ber-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 bers, 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) [4] (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 coe cient 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: ber 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 ber stages between the tissue types and between species were suitable for this study.
At 0 DPA, the ber cell initials develop in the wild type but not in the mutant. To identify differentially expressed genes in bers at 0 DPA, we compared expression differences in the epidermal samples between the wild type (GhTM-1 and Gb3-79) and mutant (GhMD17). We identi ed 1,295 and 1,500 upregulated genes in GhTM-1 and Gb3-79 ber cell initials, respectively, 981 of which (>60%) were shared in both, while 314 and 519 were speci c 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 ber 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 ber cell with epidermal cell tissues within the same genotype (GhTM-1 or Gb3-79) and between the ber cell (GhTM-1 and Gb3-79) and epidermal cell tissues in the ber-less mutant (GhMD17). We identi ed 1,861 and 2,046 up-regulated genes in GhTM-1 and Gb3-79 ber 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 identi ed 2,978 and 4,269 down-regulated genes in GhTM-1 and Gb3-79 ber cells, respectively, and 1,626 of which (>35%) were shared in both (Supplemental Figure 2C, Supplemental Table 3). These differentially expressed genes in the ber cell tissues were used for further analysis.
Dramatic gene expression changes in early stages of cotton ber 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 bers ( Figure 3). These results may indicate that ribosome biosynthesis, required for protein synthesis as well as lipid and fatty acid metabolism, were active during ber 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 ber growth. In cultured cotton ovules, bers 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 ber growth is suppressed; low concentrations of hydrolyzable nucleotides ATPgammaS/ADPbetaS stimulate ber growth [38]. Thus, the difference of ATP and ADP metabolic processes between GhTM-1 and Gb3-79 may contribute to longer bers in Gb3-79.
Among down-regulated genes in bers, 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 ber cells, which require active ribosome biosynthesis for translational regulation in ber cells. Down-regulation of phytohormone-related genes in the ber cells may suggest upregulation of these genes in ovules to support ber 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 ber 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 [4]. 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 ber cells may indicate that deep subfunctionalization of homoeologous genes during active developmental process of ber cell initiation and elongation.
Translational regulation in early cotton ber cell development During ber cell development, the cell cycle is predicted to attenuate in 20-30% of the epidermal cells, and those cells become ber cell initials. Cell cycle arrest in a cell often leads to active ribosome biosynthesis [39]. This is consistent with developmental process of ber cell initiation, when the arrest of the cell cycle in the epidermal cells induces active ribosome biosynthesis for ber cell initiation and elongation.
Arabidopsis has 231 annotated genes for cell cycle regulation, and cotton has a total of 845 putative cellcycle 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 downregulated in ber 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 [40], and many of those paralogues are differentially expressed during development and in UV-B responses [41]. For example, in ArabidopsisRPL16A and RPL16B are expressed in developing pollen and early phase of lateral root initiation [42], a process similar to ber cell initiation. The AtRPS5A gene is strongly expressed in dividing cells, whereas AtRPS5B expression correlates with cell differentiation rather than cell division [43]. These results suggest that during ber cell initiation, the cell cycle was attenuated, and ribosomal protein biosynthesis is activated and may regulate tissue-speci c expression in the ber 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 [46]. 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 ber 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 [47], 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 [48]. Thus, in ber cells down-regulation of TOR may inhibit the cell cycle, and strong expression of EBP1 could promote ribosome biosynthesis for ber cell expansion.
To further test roles of ribosome biosynthesis and cell cycle regulation in ber 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 [49], and Roscovitine is a cyclin-dependent kinase inhibitor for cell cycles [50]. Cotton ovules were isolated from fertilized owers of GhTM-1 and GhMD17 at 0 DPA and incubated in BT medium for 21 days following a published protocol [38], with or without inhibitors. Inhibiting cell cycle by Roscovitine (2 mg/L) increased total ber area and length at statistically signi cant levels, whereas impeding ribosome biosynthesis by Rbin-1 (40 mg/L) reduced the ber 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 berless ( Figure  4C left). The altered ber area could result from increased ber length or number or both. To discern this, we further measured the number or density of ber cell initials after two days of each treatment (Supplemental Figure 4). The result showed a similar trend of increasing and reducing ber cell density by Roscovitine and Rbin-1 treatments, respectively. These results support the hypothesis that cell cycle attenuation and active ribosome biosynthesis promote ber cell initiation and elongation.

Dynamic roles of phytohormones and MYB transcription factors in early ber cells
Phytohormones play important roles in ber cell development [2,24], but the temporal regulation of gene expression in early ber 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 ber cells ( Figure 5A and Supplemental Table 8). This is probably because most genes are upregulated in the epidermal layer cells that support ber cells. Abscisic acid (ABA) is known to inhibit ber cell initiation and elongation [51]. ABA in ovules at 0 DPA accumulates at much higher levels in the Ligon-lintless 1 mutant than in the wild type [52]. 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) [53] and ABA receptor gene PYRABACTIN RESISTANCE-LIKE PROTEIN 9 (PYL9) (known to interact with ABI1 and ABI2) [54]. Down-regulation of the ABA regulator and receptor genes may promote ber cell initiation and elongation.
Gibberellins promote ber cell elongation and increase ABA and auxin levels during ber development to increase ber quality (strength) [55]. Over-expression of gibberellin biosynthesis gene GIBBERELLIN 20-OXIDASE1 (GA20OX1) leads to more and longer bers [56], while their homeologs GA20OX2 and GA20OX3 are largely expressed in the ovule. Consistently, GA20OX2 was upregulated in the ovule and down-regulated in bers ( Figure 4A).
Ethylene regulates fruit ripening and inhibits vegetative growth [57,58], and is reported to promote ber elongation in cotton, which is mediated by very-long-chain fatty acid signaling pathway [28]. However, the mechanism for ethylene in ber initiation remains unclear. We found down-regulating of the ethylene signaling pathway genes like ETHYLENE INSENSITIVE 2 (EIN2) and EIN3 in bers, which may suggest a negative role for ethylene in ber cell initiation ( Figure 4A).
Auxin promotes ber initiation and accumulates in ovules prior to fertilization [9]. Auxin is transported from ovule to ber cells during ber cell development [25]. 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 ber cells ( Figure 4A). This is reminiscent of the nding that ARF2, ARF18 and PIN6 were expressed in ovules in the ber initiation stage [55]. 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 ber cells (Supplemental Figure 6), suggesting functional divergence between homoeologs for ber cell development in polyploid cotton.
In addition to phytohormones, MYB transcription factors regulate ber 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 ber 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 bers ( 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 ber 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 ber cell development.

Discussion
In this study, we used LCM to successfully dissect ber and epidermal cells from the cotton ovules during ber cell initiation and investigated gene expression changes between ber cell initials and epidermal layer cells. Our data support a new model of ber cell initiation and elongation, involving regulation of cell cycle and ribosome biosynthesis activities, as well as phytohormones and MYB transcription factors ( Figure 6).
Fiber cells are initiated in 20-30% of epidermal cells at or prior to 0 DPA and before fertilization. Fiber cell initials differentiate from epidermal layer cells and involve a rapid and dramatic transition from actively dividing ovular cells to elongation of new singular ber cells. At the cellular and molecular levels, cell cycle is predicted to arrest and coupled with active biosynthesis of ribosome for translational activities including protein synthesis. This rapid process of cell differentiation is similar to the origin of human tumor cells, where cell cycle arrest and ribosome biosynthesis are intrinsically correlated [22,23]. The tumor suppressor p53 has been shown to induce cell-cycle arrest in response to active ribosome biogenesis for rapid cell growth [22], and the transition associated with development and tumor metastasis is accelerated by upregulation of ribosome biogenesis during G1/S arrest [23]. By this analogy, although the cause and consequence are unknown, some unknown factor(s) may trigger cell cycle arrest in pre-ber cells to induce active ribosome biosynthesis to prepare for rapid cell initiation and elongation. Considering that each ber cell does not divide, abundant ribosomes may promote translational regulation and protein synthesis, leading to rapid ber cell elongation.
The TOR signaling pathway integrates nutrient, hormone, growth factor and environmental inputs to control cell cycle, proliferation, metabolism, and growth in diverse organisms, including plants [46]. Loss of TOR function results in cell cycle arrest in G1 phase in human and yeast [64], which is consistent with the observation that cotton bers cells are arrested at G1 stage [13,65]. We found in early stages of ber cells, TOR expression is downregulated (Figure 4), which may induce ribosome biosynthesis and translation rate through the TOR signaling pathway [44]. Moreover, EBP1 is upregulated in the cotton ber cells, which can positively regulate ribosome biosynthesis as part of the TOR signaling pathway and auxin regulation in plants [48]. It is likely that the TOR signaling pathway may regulate cell cycle and ribosome biosynthesis in the ovule epidermis for ber cell initiation as well as elongation in early ber cells.
Phytohormones are essential for ber cell initiation and development [2,24,51]. While some phytohormones such as ABA and cytokinin have a negative role in ber cell development, many others such as auxin, gibberellins, and brassinosteroids positively regulate ber cell initiation and development. It is unexpected to nd genome-wide down-regulation of phytohormone related genes such as ARFs, PIN6, BRZ1, EIN2, EIN3, and GA20OX2 in early ber cells. This result may suggest an important role for spatial and temporal action of phytohormones in ber cell development, which is revealed by the LCM approach in this study. Many phytohormone related genes such as ARFs and PIN6 are strongly expressed in epidermal layer cells [55] (Figure 5). Auxin may actively accumulate in the ovules and is transported into young ber cells, re ecting a spatial and temporal role of auxin in ber cell initiation. Indeed, expression of cotton ARF2:GUS in Arabidopsis is localized in the trichome base but not in the body or tip [66], consistent with its spatial and temporal role in leaf trichome development as in cotton ber cell initiation. Alternatively, there may be crosstalk between phytohormone pathways and other factors such as cell cycle regulation and ribosome biosynthesis to regulate ber cell development. For examples, inhibiting TOR signaling pathway decreases ABA synthesis [67], and auxin activates TOR signaling pathway and promotes translational regulation [68] and some other associated pathways [69].
MYB transcription factors are required for the development of leaf trichomes as well as seed hair or ber cells [32][33][34][35]. Cotton has two major groups of MYB transcription factor genes, one expressed in trichomes or leaves and another largely in ovules and ber cells [37]. Up and down-regulation of MYB transcription factor related genes may re ect different roles of these genes in trichome and ber cell development in cotton. In Arabidopsis, MYB16 and MYB106 (GhMYB25 in cotton) are closely related genes that promote cuticle formation in the epidermal cells including trichomes [70] and promote both trichome branching and expansion [71], while the cotton homolog GhMyb25 regulates cotton ber cell development [18].
Thus, upregulation of MYB106 and down-regulation of MYB16 in early ber cells may re ect their functional divergence in trichome and ber cell development. This expression divergence between closely related homoeologs are common to many other genes including those in phytohormone, cell cycle, and ribosome biosynthesis pathways. The mechanism for this expression divergence may involve small RNAs including siRNAs and miRNAs, DNA methylation, and other epigenetic modi cations [19][20][21]33].
Coincidently, the cell cycle regulators TOR and RAPTOR have been identi ed as epialleles, which are demethylated in the cultivated cotton species but not in the wild relatives [72]. This result underscores important role for epigenetic regulation of gene expression changes in the selection and domestication of agronomic traits including ber cell development, possibly ber yield and quality. Our study has provided valuable genomic resources and new insight into the roles of cell cycle regulation, ribosome biosynthesis as well as spatial-temporal action of phytohormone and MYB transcription factors in early stages of ber cell development. Understanding this complex regulation during ber cell initiation and elongation may help us develop new breeding and biotechnological tools to improve cotton ber yield and quality.

Plant materials
Two cultivated species GhTM-1 (G. hirsutum) and Gb3-79 (G. barbadense), and the ber-less mutant GhMD17 (G. hirsutum background) were grown in soil in a greenhouse under natural sun light at The University of Texas at Austin. Flowers of cotton were tagged on the day of (post) anthesis (0 DPA) and collected at 1-3 pm to minimize diurnal effects. Immediately after ower collection, ovules were manually dissected from owers under a dissection microscope. Flowers were collected in March 2016 for LCM RNA-seq analysis and in August 2018 for ovule culture assays. Unless noted otherwise, three biological replicates were used for each analysis.
Laser-capture microdissection (LCM) and mRNA-seq library preparation Ovules were collected from owers at 0 and 2 DPA from GhTM-1, Gb3-79, and GhMD17, each with two biological replicates. Tissue xation was performed as described in [73]. Brie y, cotton ovaries were cut using a razor, and individual ovules were removed using forceps, xed in a solution containing 75% ethanol and 25% acetic acid under vacuum for 15 min, and incubated at 4 o C for overnight. The samples transferred into fresh 75% ethanol and 25% acetic acid for 15 min on ice and microwaved at 250 W for 15min using a PELCO BioWave (TED PELLA Inc., Redding, CA); the process was repeated three times each with a new solution. The samples were then dehydrated by ethanol series (70%, 80%, 90%, 100%, and 100%) with microwaving at 58 o C, 300 W, for 1.5 minutes at each step. Samples were then transferred to butanol in two steps (50%:50% ethanol:butanol followed by 100% butanol, each time microwaved at 58 o C, 300W, for 1.5min. Finally, samples were transferred into molten para n (50%:50% butanol:melted para n, then 100% molten para n, with microwaving at 58 o C, 250 W, for 10 min for each step). Molten para n was changed four times, with the sample vials immersed in a hot water bath, and microwaved at 58 o C, 250W, for 30 min after each change. Finally, the samples were embedded in para n tissue blocks. Para n-embedded tissues were sectioned by a microtome (at 7mm) and mounted on a pen-foil slide (ZEISS, Oberkochen, Germany). Fiber and epidermal cells were isolated and collected using Zeiss PALM Laser Microdissection with AdhesiveCap 500 (ZEISS, Oberkochen, Germany). The LCM ber and epidermal tissues were subject to RNA extraction by RNA Aqueous Micro Total RNA Isolation Kit (Thermo Fisher Scienti c, Waltham, MA), following the manufacturer's protocol. The total RNA was treated with DNase treatment and quali ed by a Bioanalyzer (Agilent, Santa Clara, CA) with RNA Integrity Number (RIN) > 6. An aliquot of total RNA (~100 ng) was used for RNA-seq library construction using the NEB Next Ultra RNA Library Prep Kit (New England Biolabs, Ipswich, MA). mRNA-seq libraries were sequenced by paired-end reads (150 bp) using an Illumina HiSeq 2500 (Illumina, San Diego, CA). mRNA-seq data analysis of LCM samples LCM mRNA-seq reads were mapped to G. hirsutum TM-1 genome v2.0 (GhTM-1 and GhMD17) and G. barbadense genome V1.0 (Gb3-79) whole genome [4] using Tophat 2.1.1 [74]. Uniquely mapped reads were extracted and analyzed by Cu inks 2.2.1 [75] to determine gene expression levels using Fragments Per Kilobase Million (FPKM). Orthologs of G. hirsutum and G. barbadense (52,514 genes) expressed in these samples were used for further analysis. Values from two biological replicates were averaged to determine differentially expressed genes (DEGs) using the criteria of >5 FPKM, ANOVA (P < 0.01), and fold changes >1.5 (wild type epidermal vs. mutant epidermal; ber vs. wild type epidermal) or >2 ( ber vs. mutant epidermal).
Principal component analysis (PCA) was performed using log2-transformed gene expression values with singular value decomposition via the prcomp function in R (https://www.r-project.org).
Averaged gene expression values from two biological replicates were used for correlation coe ciency analysis. Pearsons's correlation coe cients were determined using the cor function in R (https://www.rproject.org).

Ovule culture
Ovules were collected from fertilized owers at 0 DPA of GhTM-1 and GhMD17. Ovule culture was conducted using a modi ed protocol as previously described [38]. Brie y, ovaries were isolated from owers and surface sterilized in 100% ethanol for 15 min, and further dissected to expose the ovules. Individual ovules were collected and oated in BT medium [51] with or without 40mg/L Rbin-1 (Axon Medchem, Groningen, Netherlands) [49] or 2mg/L Roscovitine (#R7772-1MG, MilliporeSigma, Burlington, MA) [50] under dark conditions at 28 o C. For ber initialization density counts, after 2 days, cultured ovules were brie y dried on a laboratory wipe before being gently placed on a microscope slide covered a UV curable adhesive (All-Ways Adhesives, Bensenville, IL). These slides were exposed to a hand-held UV torch for 30 seconds to set the adhesive before removing the ovules, and the resulting impressions were imaged via stereomicroscope SMZ1500 (Nikon, Tokyo, Japan). Fiber density was calculated by de ning a region of interest in ImageJ, measuring the area and counting the ber cell impressions.
For area measurements, cultured ovules were incubated for 21 days and then transferred to a glass plate where a stream of distilled water from a lab wash bottle was used to gently splay bers. Excess water was then removed using laboratory wipes, and photos were taken. Total ber and ovule areas were calculated and analyzed using ImageJ (https://imagej.nih.gov/ij/).

Consent for publication
Not applicable.

Availability of data and material
The mRNA sequencing data have been deposited in the NCBI Gene Expression Omnibus (GEO) database, https://www.ncbi.nlm.nih.gov/geo (accession GSE152995).

Competing interests
The authors declare no con ict of interest in this work.