A spontaneous mutation that changed fiber color from white to brown on one branch of the same plant was observed in a cotton (G. hirsutum L.) line of unknown background in Starkville, MS (Hinchliffe et al., 2016). The seeds from white and brown cotton bolls from this single plant were collected and designated as MC-WL and MC-BL. The lines differ in Lc1 locus only; therefore, they are near-isogenic lines (NILs). Both lines continued to maintain their corresponding fiber color through the following generations. The lines were further advanced by self-pollination single seed descent for three generations in New Orleans, LA. Pure lines of MC-WL and MC-BL were planted in the field in New Orleans, LA, in 2017 and grown using standard agricultural practices. Seeds for each cotton line were planted in rows measuring 20 m x 1 m with approximately 50 plants per row. Three rows were planted for each line, with each row used as a biological replicate for sample harvest. Developing cotton fibers were collected from the MC-WL and MC-BL cotton lines at 8, 12, 16, 20, 24, 32, 40 DPA, and mature fibers from open bolls. Mature seed cotton was hand-harvested and ginned using a laboratory roller gin.

Developing (8, 12, 16, 20, 24, 32, and 40 DPA) and mature (~50 DPA) cotton fibers in three biological replicates from MC-WL and MC-BL were used for metabolomic analysis. Fiber samples (15 mg each) were extracted using 8:1:1 acetonitrile: methanol: acetone ratio. Southeast Center for Integrated Metabolomics University of Florida performed sample extraction and metabolome analysis. Global metabolomics profiling was performed on a Thermo Q-Exactive Orbitrap mass spectrometer with Dionex UHPLC and autosampler. Catechins standards were run along with the samples to confirm their presence in the samples. All samples were analyzed in positive and negative heated electrospray ionization with a mass resolution of 35,000 at m/z 200 as separate injections. Separation was achieved on an ACE 18-pfp 100 x 2.1 mm, 2 μM column with mobile phase A as 0.1% formic acid in water and mobile phase B as acetonitrile. The flow rate was 350 μL/min with a column temperature of 25°C. Aliquot of 4 μL was injected for negative ions and 2 μL for positive ions.

Long chain alcohol standards, 1-Octacosanol and 1-Triacontanol, and Supelco 37-component FAME mix (C4 – C24) were obtained from Sigma-Aldrich (St. Louis, MO); MSTFA+1% TMCS was purchased from Thermo Scientific (Bellefonte, PA). Cotton waxes were extracted from 1 g of mature fibers with dichloromethane using an Dionex ASE 350 Accelerated Solvent Extractor (Thermo Fisher Scientific, (Waltham, MA). Extracts were taken to dryness under a nitrogen stream.

Esterification of fatty acids was performed as follows: 1 mg of dry extract was dissolved in 1 ml of diethyl ether and mixed with 20 µl of methyl acetate and 40 µl of 0.5M sodium methoxide in anhydrous methanol and incubated for 5 min; 30 µl of saturated oxalic acid in ether was added to precipitate glycerides; samples were centrifuged for 7 min and supernatants were transferred to new glass vials, dried in SpeedVac and diluted in 200 µl hexane; injection volume was 1 µl.

Derivatization of long chain alcohols was performed as follows: 10 µl of MSTFA + 1% TMCS were added to 1 mg of dry extracts; incubated at 50°C for 1 hour; 200 ul of hexane was added, mixed, and transferred into a vial for GC automatic sampling; injection volume was 100 split of 1 µl.

The fatty acids and long chain alcohols were analyzed by GC-MS using an Agilent 6890 GC, 5973 MS, and Gerstel MACH (Agilent LTM) fast GC adaptation, using an HP-88 LTM column (30 m × 0.18 mm, 0.18 μM; Agilent). The column heat gradient was 50°C (1 min hold) to 140°C at 25°C/min, followed by a 4°C/min gradient to 177°C and a 2°C/min gradient to 210°C (1.65 min hold). A pressure ramp program was set up as 29.2 psi (1 min hold) to 39.3 at 2.81 psi/min, followed by a 0.44 psi/min gradient to 43.4 psi and a 0.22 psi/min gradient to 47.1 psi (1.65 min hold), which kept a 1 mL/min helium flow rate. The inlet was set to a 50:1 split at 220°C, using an insert with glass wool (Agilent 5062-3587). The GC oven and MS transfer were set at 250°C. A 5 μL syringe was used in a 7683 autosampler. A mass range of m/z 35−500 was acquired.

Concentrations of wax components were determined using standard curves of commercially available FAMEs and long chain alcohols; the peak areas of target ions were used for concentration calculation; the peak areas were normalized by the area of the internal standard C19:0 for a specific sample. Analysis was performed in five technical replicates from two independent isolations of waxes from mature cotton fibers.

Total RNA was extracted using the Sigma Spectrum Plant Total RNA Kit (Sigma-Aldrich, St Louis, MO, USA) with on-column DNase I digestion following the manufacturer’s protocol. RNA quantity and integrity were evaluated as previously described (Hinchliffe et al., 2011).

Total RNA from 8 and 20 DPA developing fibers of MC-BL and MC-WL in two biological replicates were used for library preparation and RNAseq (LC Sciences, Houston, TX, USA). Sample preparation and library constructions were performed using TruSeq Stranded mRNA Library Prep Kit (Illumina Inc., San Diego, CA, USA) following the manufacturer’s protocols. Samples were sequenced using a HiSeq 2000 (Illumina Inc.) with 100 bp paired-end reads. Raw sequence reads were filtered for quality and trimmed by SICKLE (Joshi and Fass, 2011) and aligned to the draft G. hirsutum TM-1 reference genome (Zhang et al., 2015) with GSNAP software (Wu and Nacu, 2010). Reads mapping to annotated genes were counted using BEDTools software (Quinlan and Hall, 2010).

EDS and elemental mapping analyses of fiber surfaces in scanning electron micrographs (SEM) were conducted with an accelerating voltage of 15 kV using a Phenom G6 ProX SEM (Nanoscience Instruments, Phoenix, AZ, USA) equipped with an EDS detector. The sample was mounted on a stub using double-sided carbon tape, and then a seven nm-thick gold coating was placed onto the sample using a LUXOR gold sputter coater (Aptco Technologies, Nazareth, Belgium).

Statistical analyses of UHPLC-MS and RNAseq data of cotton fiber samples were performed using JMP Genomics 10 software (SAS Institute Inc. Cary, NC, USA). An ANOVA was conducted to find differentially accumulated metabolites or expressed genes between MC-BL and MC-WL as previously described (Naoumkina et al., 2013, 2014). The metabolites or genes that accumulated or expressed differently between two lines by at least 2-fold with a false discovery rate < 0.05 (Benjamini and Yekutieli, 2001) were considered significant. A two-way ANOVA of fatty acids was conducted with Prism 9 (GraphPad Software, Inc.).

Figure 1A shows images of open bolls of MC-WL and MC-BL near-isogenic lines (NILs). The mature fiber length of MC-BL was about 14 mm, whereas MC-WL fiber was about 27 mm. Developing (8, 12, 16, 20, 24, 32, and 40 DPA) and mature (50 DPA) fibers of white (MC-WL) and brown (MC-BL) cotton were used for metabolite analysis. Each sample (in three biological replicates) was analyzed in positive and negative heated electrospray ionization on a Thermo Q-Exactive Orbitrap mass spectrometer. A total number of detected peaks, 9140 (positive) and 5895 (negative), were used for principal component analysis (PCA) to explore the relationship between brown and white fiber samples in metabolite pools ( Figure 1B ). According to PCA, samples from developing fibers (8, 12, 16, 20, 24, 32, and 40 DPA) separated between white and brown cotton, while samples from mature fibers formed a distinct cluster ( Figure 1B ). Therefore, a substantial difference in metabolite composition (in the detectable range of molecular weights by the UHPLC-MS) was observed between MC-WL and MC-BL in developing fibers.

A two-way ANOVA determined that 5836 (ESI MS positive mode) and 4541 (ESI MS negative mode) metabolites were significantly different in accumulation between brown and white fibers ( Supplementary Data 1 ). Figure 2A shows the distribution of up-regulated and down-regulated compounds in MC-BL vs. MC-WL developing and mature fiber samples. The highest number of peaks was detected by both methods at 40 DPA and the lowest at 24 DPA of fiber development. There were significantly more up-regulated than down-regulated metabolites at each developmental time point.

Among detected peaks, 142 were non-redundant known metabolites, 70 of which were significantly (FDR < 0.05) up-regulated and 48 down-regulated ( Supplementary Tables 1, 2 ). Amino acids and their derivatives were the main category among down-regulated metabolites, whereas phenylpropanoids, organic acids, and sugars were the major categories among up-regulated metabolites ( Figure 2B ).

We analyzed RNAseq data between the NBC and white NILs from developing fibers at 8 and 20 DPA. Two strategically essential time points, including the peak of elongation (8 DPA) and secondary cell wall (SCW) biosynthesis (20 DPA), were selected for RNAseq due to cost efficiency, with two biological replicates minimally required for statistical analysis. A two-way ANOVA determined that 1691 (8 DPA) and 5073 (20 DPA) genes were significantly (FDR < 0.05) differentially expressed (DEGs) between NILs ( Supplementary Data 2 ). The count of up-regulated and down-regulated DEGs was 986 and 705, correspondingly, at 8 DPA ( Figure 3A ). However, at 20 DPA, almost all of the DEGs were down-regulated, 4993 genes ( Figure 3A ).

We used an integrated web-based Gene Ontology (GO) analysis toolkit agriGO to analyze the GO enrichment of DEGs (Tian et al., 2017). The GO analysis showed that the up-regulated genes at 8 DPA in the NBC fibers were enriched in secondary metabolism, phenylpropanoid, flavonoid, nucleotide, oxidoreduction, organic acid, carbohydrate, etc. biological processes, cell wall, plasmodesmata, vacuole, apoplast, cytosol, etc. cellular components, and oxidoreductase activity, heme and tetrapyrrole binding molecular functions, whereas down-regulated genes were enriched in plastid ribosome, plastoglobule, photosystem cellular components and chlorophyll binding molecular function ( Figure 3B ). At 20 DPA, the up-regulated DEGs in NBC fibers were enriched in secondary metabolic and oxidation biological processes, oxidoreductase activity, and transmembrane transported activity molecular function. In contrast, down-regulated genes were enriched in nucleic acid, chromosome, cell cycle, etc. biological processes, nucleus, cytoskeleton, nucleoplasm, etc. cellular components, and nucleic acid binding molecular function ( Figure 3B ).

Organic acids are synthesized as intermediates of metabolic pathways, representing fixed carbon’s transitory or stored forms. The citric acid cycle (TCA cycle) connects most metabolic pathways and is essential for energy production. Figure 4 shows changes in the MC-BL fibers in detected metabolites and transcripts of genes encoding enzymes involved in glycolysis and TCA pathways. Most genes from both pathways were up-regulated in MC-BL fibers at 8 DPA and down-regulated at 20 DPA. Genes encoding three enzymes from TCA pathway were significantly differentially expressed between fibers of MC-BL and MC-WL NILs, including aconitase (ACO), isocitrate dehydrogenase (IDH), and malate dehydrogenase (MDH); among them, the MDH was shown to control flux of TCA cycle (Zhang and Fernie, 2018). Glycerol-3-P and detected organic acids from TCA, including citrate, 2-ketoglutarate, succinate, and malate, were significantly (FDR < 0.05) increased in MC-BL compared to MC-WL fibers. In contrast, in MC-BL fibers, shikimate and most amino acids, including phenylalanine, were significantly decreased ( Figure 4 ). Thus, the data indicate that significant changes occurred in the primary metabolism of the brown compared to its white NIL.

Substantial carbon flux from primary metabolism was directed into the PA pathway. Both isoforms of flavan-3-ols, (+)-catechin and (-)-epicatechin, were detected in NBC and white developing and mature fibers. Figure 5A shows the structures of catechins and epicatechins with gallate modifications detected in fibers of both NILs. A representative total ion chromatogram of extract from fibers at 20 DPA reveals that detected catechins are present in brown and white fibers; the relative amount of catechin is multiple times greater than epicatechin and their corresponding gallate modifications ( Figure 5B ).

Figure 6 presents relative (based on peak areas) quantifications of catechins, epicatechins, and their gallate modifications in developing and mature fibers of NILs. Catechin and gallocatechin were the most abundant compounds among detected PA subunits. Gallocatechin-3-gallate and Epcatechin-3-gallate accumulated in small amounts, mainly in 8-12 DPA fibers, gradually decreasing during fiber maturation. Catechin, gallocatechin, epicatechin, and epigallocatechin accumulated the most during 16-24 DPA of fiber development in the NBC fibers and substantially less in white fibers. Concentrations of all detected PA subunits were significantly reduced in mature fibers, indicating that they converted into condensed tannins.

The composition of waxes extracted from mature fibers was compared between NBC and white NILs. Figure 7 shows a total ion chromatogram and quantification of ester derivatives of fatty acids extracted from NBC fibers. Palmitic acid (C16:0) and long-chain alcohols, 1-Octcosanol and 1-Triacontanol, were the main components of cotton fiber waxes. Palmitic acid is the building block for many other fatty acids and phospholipids and is widely present in different classes of plant lipids. Palmitic acid and 1-Triacontanol were significantly (FDR < 0.05) higher accumulated in NBC than white fiber waxes ( Figure 7B ). Expression levels of genes involved in biosynthesis fatty acids, including acetyl-CoA carboxylase, acyl-carrier protein, acyl CoA ligase, and 3-ketoacyl-CoA synthase, were significantly up-regulated in NBC developing fibers at 8 DPA but down-regulated at 20 DPA ( Supplementary Table 4 ).

It has been previously shown that the presence of various inorganics in raw cotton fibers can alter the thermal properties of cellulose (Nam et al., 2014; Hinchliffe et al., 2015; Nam et al., 2016, 2017). Elemental analysis of the surfaces of NBC and white mature fibers in SEM was conducted using EDS. Figure 8 shows the EDS results including color-coded elemental maps, EDS spectra, and elemental compositions in atomic percent. Among inorganic elements, only potassium was detected by EDS. The EDS spectra showed that potassium measured in NBC was three times greater than in white fibers ( Figures 8B, C ). To explain the increased potassium level in NBC fibers, we searched transcript data to identify DEGs involved in potassium transport. Two G. hirsutum genes homologous to Arabidopsis, the two-pore potassium (K+) channel (TPK3), were significantly (FDR < 0.05) up-regulated in NBC developing fiber at 8 DPA ( Figure 8D ). In Arabidopsis, the TPK3 regulates the proton gradient in the thylakoid membrane and, therefore, photosynthetic light utilization (Carraretto et al., 2013).