Two different cotton species, Gh (accession TM-1) and Gb (accession 3-79), were grown synchronously under controlled conditions in a growth chamber (Conviron E-15, Controlled Environments Inc. N.D., USA) as previously described (Grover et al., 2025; Lee et al., 2025a, 2025b). Plants were grown individually in a two-gallon pot containing a potting mix (4:2:2:1 ratio of soil:perlite:bark:chicken grit). Growth chambers were set for 16 h days at 28°C and 500 µmol m^-² s^-¹ of light. Flowers were self-pollinated before noon every day and tagged. For 20 different time points (from 6 to 25 DPA), bolls were collected at mid-day and stored at -80° C. At each time point, three biological replicate bolls were used for glycome, proteome and transcriptome analysis.

From the cotton fiber, the CW and the buffer-soluble and alkali-soluble polysaccharide fractionswere extracted as per the established protocol (Avci et al.,2013; Swaminathan et al., 2024). In brief, using a razor blade and tweezers, fibers were isolated from the developing seeds, and care was taken to avoid the seed coat while harvesting the fibers. Harvested cotton fibers from a single boll were ground into powder using liquid nitrogen, and the CW was then extracted from the powdered fiber using solvents (Swaminathan et al., 2024). Polysaccharide fractions were sequentially extracted from the CW by first using a 50 mM CDTA:50 mM ammonium oxalate (1:1) buffer, followed by 4 M KOH to extract buffer-soluble and alkali-soluble polysaccharides, respectively (Avci et al., 2013; Swaminathan et al., 2024). The buffer-soluble and alkali-soluble extracts were adjusted to pH 7.0, and dialyzed (SnakeSkin Dialysis Tubing, MWCO: 3.5 kDa; ThermoFisher Scientific) for 4 days using only sterile double-distilled water to remove the salts (Avci et al., 2013). The dialyzed extracts were frozen, lyophilized to dryness, and weighed. The final pellet, which contained a mixture of amorphous and crystalline celluloses, was weighed. The crystalline cellulose content present in the final pellet was measured by using the Updegraff reagent (acetic acid: nitric acid: water, 8:1:2 v/v) (Updegraff, 1969) (Supplementary Data Sheet 1).

Glycome profiling was carried out using a standard protocol (Pattathil et al., 2010; Pattathil et al.,2012). Initially, the sugar content in each of the buffer-soluble and alkali-soluble polysaccharide fractions from each sample was estimated by the phenol-sulfuric acid method (Pattathil et al., 2012). Later, the polysaccharide sample was dissolved in water. An equal amount of 3 µg (50 µl/well from a 60 µg/µl solution) was added to each well in 96-well Enzyme-Linked Immunosorbent Assay (ELISA) plates (Costar 3595). The polysaccharides were allowed to bind to the bottom of the well by drying in an oven at 37°C. Glycome profiling was conducted using an ELISA method as described earlier (Swaminathan et al., 2024) with 71 polysaccharide epitope-specific antibodies selected based on the literature (Avci et al., 2013; Thorne et al., 2023). In brief, the ELISA procedure includes the following steps: the polysaccharide coated wells were blocked with nonfat dry milk solution, incubated for 1 h, washed once, incubated with glycome antibodies for 1 h, washed three times, incubated with peroxide-conjugated secondary antibody for 1 h, washed five times, incubated with TMB-peroxidase substrate solution for 40 min, the reaction was stopped by adding 1.0 N sulfuric acid and the color development was detected at 450 nm wavelength in an ELISA plate reader and the optical density (OD) reading from each sample were used for further analysis. These selected 71 antibodies detect the backbone and decorated side chains of the pectins, homogalacturonan (HG) and rhamnogalacturonan-I (RG-I), as well as the hemicelluloses, XGs and xylans (Xyls), and callose. The specific details of each antibody used in this study, along with the polysaccharide epitopes detected, are reported in the supplementary information (Supplementary Data Sheet 2; Supplementary Figure 1). Data from the ELISA for each epitope were obtained from either buffer-soluble oralkali-soluble polysaccharide fractions of the same sample (Supplementary Data Sheet 3). To distinguish between the two types of polysaccharide fractions, for the naming of epitopes, we used a naming convention where a suffix of “-C” refers to the buffer-soluble fraction (50 mM CDTA:50 mM ammonium oxalate buffer extracts) and “-K” refers to the alkali-soluble fraction (4 M KOH extracts).

The SOM is a robust unsupervised machine-learning method that has built-in Cosine similarity(statistical distance metrics) used for clustering analysis/pattern discovery (Kohonen, 1990). An SOM algorithm, based on an inner product distance metric, was applied to the ELISA-based glycome profiling results (as described in previously established protocols; Lee et al., 2025a; Yang et al., 2025) to cluster the abundance profiles of glycome epitopes from two separate interpolated datasets: buffer-soluble and alkali-soluble polysaccharide fractions from two different species of cotton (Wehrens and Kruisselbrink, 2018). We optimized the number of clusters by balancing between over−segmentation (leading to empty or single-member clusters) and under−segmentation (losing biologically meaningful groups), ultimately selecting a 12-cluster SOM. Thus, SOM grid structures consisted of 3 rows and 4 columns for both buffer-soluble and alkali-soluble polysaccharide fractions datasets. The resulting clusters represent groups of glycome epitopes with similar abundance trajectories. In the SOM analysis, the data from both species for each of the fractions were combined, and two different colors were used to maintain the identity of polysaccharide epitopes. The SOM was performed in R Studio (Version 2024.12.0.467, R Studio; R Core Team, 2022) using the supersom function in the kohonen package (Version 3.0.12, Wehrens and Kruisselbrink, 2018). The p-value threshold determining significance was 0.05. Based on the significant differences, the polysaccharide epitopes were grouped into mainly three different categories (presented in figures and supplementary information).

In parallel to glycome profiling, RNAseq was carried out for the cotton fibers collected simultaneously for the same days (6 to 25 DPA) (Grover et al., 2025). Briefly, a modified version of the Sigma Plant Spectrum Total RNA kit (Sigma-Aldrich) was used to extract RNA from the fibers. Library construction (NEBNext Ultra II RNA Library Prep Kit) and sequencing (as PE150; Illumina NovaSeq 6000) were conducted by the Iowa State University DNA facility. Data were processed using Trimmomatic version 0.39 (trimmomatic/0.39-da5npsr) (Bolger et al., 2014) from Spack (Gamblin et al., 2015) for read and quality trimming. The transcripts of each gene were annotated based on the published reference Gh genome (Chen et al., 2020) and Gb genome (Chen et al., 2020). The transcripts per million (TPM) output by Kallisto for each sample were imported into R/4.2.2 (R Core Team, 2022), and RNA-seq quality was assessed by consistency of the number of genes expressed over time and among replicates. The number of expressed genes per sample (TPM > 0) was plotted across developmental time using ggplot2, and visual outliers were discarded. RNAseq reads of cotton fibers collected from 6 to 25 DPA were deposited into the NCBI-SRA under PRJNA1099209 (G. hirsutum) and PRJNA1222456 (G. barbadense). Differential gene expression analysis was performed using DESeq2, and p-values were adjusted using the False Discovery Rate (FDR) Benjamini-Hochberg method (Benjamini and Hochberg, 1995). Differential gene expression was determined, where p-values were less than 0.05.

Glycosyltransferases and enzymes that decorate polysaccharides play key roles in the synthesis offiber CW polysaccharides. To compare the transcript levels of these enzymes between the two species,we first compiled a list of known or putative Arabidopsis thaliana glycosyltransferase genes based on our previous study and other studies (Swaminathan et al., 2024). Subsequently, the gene sequences of Arabidopsis proteins were used to identify the homologous genes of Gh and Gb from the Phytozome 13 (https://phytozome-next.jgi.doe.gov) and the Cotton Functional Genomics Database (https://cottonfgd.net) (Supplementary Data Sheet 4). Based on this information, transcript profiles of predicted glycosyltransferases and other polysaccharide-decorating enzymes in both Gh and Gb were generated from the transcriptome data (Grover et al., 2025). Comparisons of profiles of transcripts of these enzymes between Gh and Gb and their sub-genomes (At and Dt) were conducted.

Since the glycosyltransferases that synthesize polysaccharides are membrane proteins, themicrosomal protein fraction (P200) was isolated from intact Gb fiber from 6 to 25DPA, as described previously (Lee et al., 2025a). Briefly, apoplastic proteins and extracellular vesicles were removed from the intact ovules by dipping each intact ovule into 5 mL of pre-chilled microsome isolation buffer (MIB) for 10 min. Subsequently, the fiber tissues were isolated, homogenized in cold MIB using a Polytron homogenizer, and filtered through cheesecloth. From the filtrate, microsomes were enriched by ultracentrifugation (200,000 x g for 20 min at 4°C). The final pellet was mixed with 8 M Urea to denature membrane proteins, digested with trypsin, and desalted using C18 columns (Lee et al., 2025a, 2025b). The peptides were analyzed using Bruker’s TIMS-TOF HT (Bruker Daltonics GmbH) mass spectrometer (MS), which was coupled with a reverse-phase liquid chromatography system, nanoElute2 (Bruker Daltonics GmbH). The MS raw data were processed using the directDIA™ approach in Spectronaut software (Version 19.4, Biognosys) according to the vendor’s recommendations. The spectral library was constructed using a cotton FASTA file containing 55,237 entries and DPA-pooled DDA runs. Cysteine Carbamidomethylation was set as a fixed modification, while protein N-terminal acetylation and methionine oxidation were set as variable modifications. Trypsin was specified as the digestion enzyme, allowing a maximum of two missed cleavages. The FDR was set to 1% at both the peptide and protein levels. Quantification was based on MS2-level fragment ion extracted ion chromatograms (XICs), with protein abundances calculated using Spectronaut’s MaxLFQ algorithm. The glycosyl transferase proteomic data of Gb fibers presented in this study (Supplementary Data Sheet 5) is part of a larger whole proteome dataset. The whole Gb fiber proteomedataset consisting of 13,357 proteins including glycosyl transferases is presented in Supplementary Data Sheet 5.

The experimental data, including Pearson correlation coefficient (PCC) analysis, were subjected to statistical analysis using R software (Version 2024.12.0.467, R Studio; R Core Team, 2022). PCC value with >0.7 was considered the best correlation profile pattern. Statistical significance was set at p < 0.05, and p-values were adjusted using the FDR. The t-test was carried out to compare the polysaccharide content at different DPA to find out the statistical significance in R Studio (Version 2024.12.0.467, R Studio; R Core Team, 2022). Exact p-values are presented in the Supplementary Table, and p-values were adjusted using the FDR. The p-value threshold determining significance was 0.05.

The buffer-soluble fractions (50 mM CDTA:50 mM ammonium oxalate buffer extracts), alkali-soluble fractions (4 M KOH extracts), and cellulose were extracted from cotton fiber CW (6 to 25 DPA), lyophilized, weighed and averaged (Figure 1; Supplementary Data Sheet 1). As expected, the amount of cellulose and non-cellulosic polysaccharide fractions per boll increased gradually from 6 to 25 DPA in both species.

At 6 DPA, total cellulose content (consisting of both amorphous and crystalline cellulose) was 2.3 mg/boll and 2.1 mg/boll, for Gh and Gb, respectively (Figure 1; Supplementary Data Sheet 1). Cellulose content gradually increased in both species after this time; however, the rate of increase was slower in Gb. By 16 DPA, cellulose content reached 25.9 mg in Gh and 24.3 mg in Gb. After 16 DPA, the cellulose content of Gh started to increase rapidly (53.3 mg at 17 DPA; 124.2 mg at 19 DPA), coinciding with the SCW thickening stage (Li et al., 2013), and reached 461.6 mg at the end of 25 DPA. Interestingly, unlike Gh, the cellulose content of Gb increased gradually and slowly after 16 DPA and reached only about 45.9 mg at the end of 19 DPA, and after that, it started to increase rapidly (103.6 mg at 20 DPA; 162.6 mg at 21 DPA), coinciding with the SCW thickening stage, reaching 442.3 mg at the end of 25 DPA, a comparable level to that of Gh. The amount of estimated crystalline cellulose content of both species showed profile patterns similar to that of total cellulose content. Crystalline cellulose content was about 0.2 mg per boll at 6 DPA and about 9.0 mg at 16 DPA in both species (Figure 1; Supplementary Data Sheet 1). From 17 DPA onwards, crystalline cellulose content increased rapidly in Gh, but more slowly in Gb. At 20 DPA, crystalline cellulose content was 78.2 mg in Gh and 54.7 mg in Gb. After 20 DPA, this amount increased rapidly (at 21 DPA, 97.8 mg in Gh and 92.8 mg in Gb). At 25 DPA, the crystalline cellulose content was estimated to be 293.6 mg in Gh and 303.5 mg in Gb. Overall, in Gb there was a three-day delay in cellulose content accumulation relative to Gh.

The content of buffer-soluble polysaccharides at 6 DPA was lower in Gb (1.3 mg) in comparison with Gh (2.9 mg). The relative content of the buffer-soluble polysaccharide fraction gradually increased after 6 DPA, and the increment was faster after 17 DPA for Gh and 20 DPA for Gb, which coincided with their respective increase in cellulose content. The relative amount of buffer-soluble polysaccharide fractions reached 37.2 mg in Gh and 36.6 mg in Gb at the end of 25 DPA. The content of alkali-soluble polysaccharides was 2.1 mg for Gh and 2.0 mg for Gb at 6 DPA and 15.4 mg and 17.4 mg, respectively, at the end of 16 DPA. Similar to the buffer-soluble polysaccharide fraction, the relative content of alkali-soluble polysaccharides started to increase at a faster rate after 16 DPA, coinciding with an increase in the cellulose content, and reached 41.2 mg and 46.8 mg in Gh and Gb, respectively, at the end of 25 DPA.

We conducted comprehensive glycome, transcriptome, and proteome profiling of only Gb fibers (collected daily from 6 to 25 DPA). As described in our previous study (Swaminathan et al., 2024), glycome profiling of cotton fiber CW polysaccharide epitopes was performed using 71 monoclonal antibodies targeting diverse CW polysaccharide epitopes (Supplementary Figure 1; Supplementary Data Sheet 2), and the glycome profiling data of Gb fiber are presented in Supplementary Data Sheet 3. Then we curated a list of candidate glycosyltransferases (Supplementary Data Sheet 4). Glycosyltransferases are membrane proteins isolated in the microsomes during proteinpurification. So, in order to identify the glycosyltransferases responsible for synthesizing thepolysaccharides corresponding to the observed epitope patterns (Supplementary Data Sheet 3), we integrated transcriptome and proteome data from the microsomal (P200) fraction isolatedat each developmental stage (6 to 25 DPA) (Supplementary Data Sheet 5). PCC analysis was used to integrate glycome profiles with gene and protein accumulation data.

In this study, we focused on polysaccharide-synthesizing glycosyltransferases and associatedpolysaccharide-decorating enzymes, such as methyltransferases and acetyltransferases, that areresponsible for the biosynthesis and modification of CW polysaccharides analyzed in the glycome profiles. Among the 13,357 proteins detected in the P200 proteome of Gb fibers, we identified only 93 glycosyltransferases (Supplementary Data Sheet 5) that belong to the 550 putative glycosyltransferases identified by homology fromArabidopsis (Supplementary Data Sheet 4). About 38% of mRNA/protein pairs in P200 had PCC values > 0.50 (~15% > 0.7 PCC),suggesting that their protein abundance could be explained by transcriptional control to some extent(Supplementary Data Sheet 5). Also, the PCC analysis revealed significant correlations between a few of thepolysaccharide epitopes and the profile patterns of corresponding glycosyltransferases at bothtranscript and protein levels (Supplementary Data Sheet 6), with representative correlated profiles shown in Figure 2. Overall, from PCC analysis, higher correlations were observed between XGs and theirassociated glycosyl transferases, followed by HGs and Xyls (Supplementary Data Sheets 5, 6). A minimal correlation was observed between RG-I and callose, as well as between RG-I and the associated glycosyltransferases.

For XG biosynthesis or modification, four XG xylosyltransferases (XXTs) (GbXXT2-1Dt, GbXXT2-2At, GbXXT3-1At, and GbXXT3-2At/Dt) were detected in the proteome. These XXTs exhibited protein accumulation profiles that closely mirrored the abundance patterns of xylosylated XG epitopes (GbXG-XX and GbXG-XXXG), with strong correlations (PCC ≥ 0.7) and with profiles gradually decreasing from 6 to 25 DPA of fiber development (Figure 2; Supplementary Data Sheet 6). Notably, two out of the four XXTs (GbXXT2-1Dt and GbXXT3-1At) also displayed similar transcript profiles, with PCC values greater than 0.75. Similarly, the transcript levels of two of four XG galactosyltransferases (GbMUR3-2At and GbMUR3-2Dt) showed high correlations with PCC values of 0.82 and 0.92, respectively. Both MURs also exhibited strong correlations (PCC = 0.7 - 0.95) with galactosylated XG epitopes (GbXG-XLX and GbXG-L), as shown in Figure 2 and detailed in Supplementary Data Sheet 6.

The analysis revealed that many of Xyl-related glycosyltransferases, identified from both the proteome and transcriptome datasets, showed significant correlations with the profiles of xylan backbone (Xyl-BB) and branched Xyl epitopes (GbXyl-GlcA and GbXyl-MeGlcA) (Figures 2C, 2; Supplementary Data Sheet 6). Notably, glycosyltransferases such as GbIRX9-3At/Dt, GbIRX9-4At, GbIRX10-1At/Dt,GbIRX15-2At/Dt, GbFRA8-1At/Dt, GbGUX1-1At/Dt, GbGUX2-2Dt, and GbGXMT3-1At/Dt exhibited strongcorrelation (PCC > 0.7) with GbXyl-MeGlcA epitopes (Supplementary Data Sheet 6). Among them, xylosyltransferase GbIRX9-3At/Dt also showed a strong correlation with the corresponding transcript profiles (PCC = 0.8) (Figure 2; Supplementary Data Sheet 6). Also, a few methyltransferases (GbRWA1-1At/Dt, GbRWA1-4At/Dt) and two acetyltransferases(GbTBL32-1Dt, GbTBL32-2Dt) correlated with GbXyl-MeGlcA epitopes (Supplementary Data Sheet 6). For HG biosynthesis, GbQUA2-2At/Dt and GbTBR-5At/Dt showed strong correlations between their transcript and protein accumulation profiles (PCC > 0.7), and their protein levels were highly correlated with the abundance patterns of homogalacturonan (GbHG) epitopes (PCC > 0.85) (Figure 2; Supplementary Data Sheet 6). For RG-I biosynthesis, the profile of GbGal4-BB epitopes showed a strong correlation with the protein accumulation profile of GbGALS1-2Dt (PCC = 0.95) (Figure 2; Supplementary Data Sheet 6).

To find out the differences in polysaccharide epitope compositions and differences in temporal dynamics of polysaccharides between Gh and Gb fibers, ELISA absorbance data from glycome profiling of Gh and Gb fibers at 6–25 DPA (Supplementary Figure 1; Supplementary Data Sheet 3) were visualized as heat maps and grouped into SOM clusters. ELISA data for each epitopewere obtained from two different fractions of the same sample: the buffer-soluble and alkali-solublefractions (Supplementary Data Sheet 3). ELISA data from Gh and Gb were analyzed together, but the identities of all 71 epitopes were maintained separately for both species.

Heat maps were generated separately for buffer-soluble polysaccharide and alkali-soluble polysaccharide fractions from both species (Figure 3). In these heat maps, epitopes were matched one-to-one between Gh and Gb to highlight their differential distribution during fiber development (Figure 3). The heat maps of buffer-soluble polysaccharides showed that many of the highly branched RG-I, arabinogalactan (AG), XG, and some Xyl epitopes were present in lower amounts in Gb. Also, their quantity dropped at earlier stages of fiber development in Gb in comparison with Gh (Figure 3A). The heat maps of alkali-soluble fractions showed that many of the XG and Xyl epitopes were present at high levels and in comparable amounts between both species. In addition, the heat maps clearly showed that the profile patterns and content of some of the XG (XG-XX-1-K, XG-XX-2-K, XG-XLX-1-K, XG-XLX-2-K, XG-XLX-3-K, XG-XnLG-K), Xyl (Xyl-2Ar-1-K, Xyl-2Ar-2-K, Xyl-3Ar, Xyl-MeGlcA-2-K, Xyl-3-1-K), and of the RG-I (RG-I-1-4-K, RG-Ia-K) epitopes were prominently different between the two species (Figure 3B).

Next, the glycome profiles from Gh and Gb were combined and clustered using an unbiased SOM approach to investigate the dynamics of epitope compositions across fiber development. In total, 12 SOM cluster groups of glycome profile epitope patterns were generated. SOM analysis was mainly carried out to identify similarities and differences of polysaccharide epitope profiles between the two species (cluster G1-C to G12-C for buffer-soluble polysaccharide fractions and cluster G1-K to G12-K for alkali-soluble polysaccharide fractions) (Figure 4; Supplementary Data Sheets 7, 8). Further, the SOM grouping analysis helped us easily categorize epitopes at a finer levelbased on the profile patterns and the number of specific epitopes, as well as at the species level.Members present in each SOM cluster are provided in Supplementary Data Sheet 7. The SOM profiles for buffer-soluble polysaccharides showed that most polysaccharide epitopes had higher contents at 6 DPA, then gradually or rapidly decreased to lower levels at 25 DPA (Figure 4A). However, the total amount of buffer-extractable components continued to increase during fiber development (Figure 1A). This indicates that some polysaccharide epitopes remained undetected in buffer-soluble fractions of fiber CW by the antibodies used in our study. Epitopes in clusters G8-C, G11-C, and G12-C remained low throughout all DPAs. SOM clustering of alkali-soluble polysaccharides showed greater variability in epitope patterns, with some increasing or decreasing gradually or rapidly, while others peaked between 12 and 16 DPA during fiber development (Figure 4C).

Both the SOM and heat maps provided complementary details of the epitope patterns and their differential distribution between the fibers of the two species during development (6–25 DPA). These data allowed us to group epitopes into three major “categories” (Figures 5–7; Supplementary Data Sheet 8), namely 1) polysaccharide epitopes that were similar in content and have no temporal variability between the two species, 2) polysaccharide epitopes that had similar profiles through all 20 DPAs, but varied in their content between the two species, and 3) polysaccharide epitopes that had more noticeable differences in their content and profiles between two species.

The first category included 44 polysaccharide epitopes with comparable content and similar profiles in the two species, indicating no significant difference in some of these epitopes between the fibers (Figure 5; Supplementary Data Sheet 8A). This group contained epitopes present in HGs (both methyl esterified and de-esterified), RG-I-BB (RG-I backbone epitopes), Gal4-BB (secondary galactan backbone of RG-I), highly branched RG-I (RG-I-Gal2, RG-I-Gal3s), XG (XG-(XXXG)2, XG-(XXXG)3, XG-L, XG-FF), and Xyl (Xyl-BBs, Xyl-GlcA, Xyl-MeGlcA-1). Most of the epitopes in this category, from both Gh and Gb, were found in the same SOM groups. For instance, buffer-soluble HG epitopes (HG-Bbde-1-C/-2-C, HG-BBMe-C) and RG-I epitopes (RG-I-BB-1-C/-2-C/-3-C, RG-I-Gal3-3-C) of both Gh and Gb were present in the same SOM group G9-C (Figure 5; Supplementary Data Sheet 8A). Similarly, some alkali-soluble HG epitopes (HG-Bbde-1-K, HG-Bbde-2-K), RG-I epitopes (RG-I-BB-1-K, RG-I-Gal2-1-K, RG-Ic-2-K), XG epitopes (XG-(XXXG)2-2-K, XG-FF-2-K, XG-Fuc-K, XG-L-1-K, XG-L-2-K), and Xyl epitopes (Xyl-BB3-K, Xyl-GlcA-K) from both species were assigned in the same G9-K SOM cluster. Many other polysaccharide epitopes from both species were also grouped in the same SOM groups, e.g., RG-I-BB-3-K, Gal4-BB3-C, Gal4-BB3-K, Gal4-BB6-K, RG-Ic-1-C, RG-I-Gal2-1-K, RG-I-Gal3-5-C, RG-I-Gal3-7-K, Xyl-MeGlcA-1-K, and Cal-K. SOM analysis revealed subtle differences in overall patterns and abundance of epitopes between the species. For example, many of the Xyl backbone epitopes (Xyl-BBs) from Gh were found in the G10-K SOM group, whereas the same epitopes from Gb were found in the neighboring G9-K SOM group, although profiles in these two SOM groups were not statistically different (Figure 5; Supplementary Data Sheet 8A).

The other 48 epitopes fell into the second category, which had the same profiles but showed significant differences in content between Gh and Gb (Figure 6; Supplementary Data Sheet 8B). Many of the highly branched RG-I, AGs, XG, and some Xyl epitopes were included in this group. Despite having similar profiles, these epitopes showed significant quantitative differences between the two species, and the SOM analysis effectively distinguished them based on their quantitative differences. For example, in Gh, six XG epitopes (XG-XX-1-C, XG-XX-2-C, XG-(XXXG)2-2-C, XG-(XXXG)3-2-C, XG-FF-2-C, XG-XLX-1-C) and three highly branched RG-I epitopes (RG-I-Gal3-1-C, RG-I-Gal3-2-C, RG-I-Gal3-7-C) were assigned to the G1-C group (Figure 6; Supplementary Data Sheet 8B). In contrast, the same above epitopes from Gb were distributed to the G2-Cand G4-C groups. Similarly, the G5-C group contained only polysaccharide epitopes fromGh, such as branched RG-I and XG epitopes, whereas the same epitopes from Gb were predominantly distributed to the G2-C and G10-C groups. Interestingly, only four epitopes out of 48 were present in higher amounts in Gb in comparison with Gh, whereas the remaining 44 epitopes showed higher presence in Gh relative to Gb (Supplementary Data Sheet 8B). Overall, glycome profiling showed that Gb had significantly lower levels of highly branched RG-I, AG, and Xyl epitopes compared with Gh (Figure 5; Supplementary Data Sheet 8B).

Comparative heat maps and SOM analysis clearly differentiated 14 polysaccharide epitopes with distinct profile patterns between the two species, forming a third category (Figure 7; Supplementary Data Sheet 8C). Interestingly, three epitopes, RG-Ia-K (present in the group G4-K for Gh, G12-K for Gb), XG-XnLG-K (present in G3-K for Gh, G12-K for Gb), and RG-I-Gal1-4-K (present in G4-K for Gh, G11-K for Gb), were not abundant and their profiles were flat in Gb through 6 to 25 DPA. On the contrary, in Gh, the profiles of these polysaccharide epitopes peaked in the middle of the fiber development. The content of five XG epitopes, XG-XX-1 (present in G6-K for Gh, G2-K for Gb), XG-XX-2 (present in G5-K for Gh, G2-K for Gb), XG-XLX-1 (present in G6-K for Gh, G2-K for Gb), XG-XLX-2 (in G2-K for Gh, G11-K for Gb), and XG-XLX-3 (in G1-K for Gh, G6-K were Gb), and two Xyl epitopes, Xyl-2Ar-1 (in G4-K for Gh, G8-K for Gb), and Xyl-2Ar-2 (in G3-K for Gh, G8-K for Gb) were lower in content in Gb at the earlier DPAs but increased gradually to match the amount of Gh at later stages (Figure 7; Supplementary Data Sheet 8C). A reverse profile was observed for the Xyl epitope Xyl-MeGlcA-2 (present in G2-K for Gh and G10-K for Gb), with lower content in Gh that gradually increased to equal the amount present in Gb at later DPAs. In Gb, the amount of two of the Xyl epitopes, Xyl-3Ar (present in G11-K for Gh, G8-K for Gb) and Xyl-3-1 (in G2-K for Gh, G3-K for Gb), was comparable to Gh at earlier stages but drastically reduced after 15 DPA to reach a significantly lower levels in Gb by 25 DPA.

About 34 epitopes were detected at very low levels and showed noisy profiles and therefore not analyzed further; these epitopes were mainly present in the SOM groups, G12-C and G12-K (Figure 4; Supplementary Data Sheet 8D).

Here, we observed a three-day delay in the rapid accumulation of cellulose content in Gb fibers compared to Gh during the transition period prior to SCW thickening (Figure 1). This observation agrees with earlier studies reporting an extended fiber elongation phase at Gb compared to Gh (Tuttle et al., 2015; Hu et al., 2019; Pettolino et al., 2022). Such a delay in cellulose accumulation has been reported previously, even across different Gb accessions and growing conditions (Lee et al., 2015; Pettolino et al., 2022). Coincidentally, a similar delay in transcript accumulation was observed for ten SCW CESAs in Gb compared to Gh (Figure 8). In addition, transcript levels of two PCW CESAs in Gb were higher and persisted for a longer period than in Gh. Previous transcriptomic studies utilizing only three time points (10, 15, 20 DPA) revealed that only one of the PCW CESA (CESA6) had a higher transcript level in Gb relative to Gh (Hernandez-Gomez et al., 2015). Several other studies that examined fewer time points of fiber development also reported lower transcript levels and a few days’ delay in the accumulation of transcripts in some SCW CESAs in Gb relative to Gh (Lacape et al., 2012; Li et al., 2013; Hernandez-Gomez et al., 2015).

It was proposed that the delayed timing of cellulose accumulation might be one of the possible reasons for the better fiber quality of Gb relative to Gh (Avci et al., 2013; Lacape et al., 2012; Li et al., 2013; Rajasundaram et al., 2014; Pettolino et al., 2022). Additionally, it was well established that higher SCW CESA transcript accumulation leads to increased cellulose accumulation at the SCW developmental stage, resulting in the thickening/rigidification of the CW, which in turn halts CW expansion (Haigler et al., 2012). Our findings, combined with those of earlier reports, suggest that slower accumulation of cellulose (due to extended/higher PCW CESAs transcript accumulation combined with delayed lower levels of transcript accumulation of SCW CESAs for a few days) might have resulted in delayed rigidification of CW, which potentially resulted in longer period of elongation, and longer and thinner fibers of Gb. Even though cellulose accumulation was delayed by about three days and had a lower content in Gb at around 19 DPA compared to Gh (Figure 1), cellulose accumulation in Gb increased rapidly later. By the end of 25 DPA, both species had reached the same level of cellulose. Earlier, Li et al. (2013) reported that the CESA8 was the major player for SCW cellulose accumulation in both Gh and Gb fiber and significantly contributed to rapid cellulose accumulation at later stages in Gb during the SCW thickening stage. Interestingly, our transcriptome data also showed that the transcript level of only CESA8-B (Figure 8F) was highest and equal in both species in comparison to all other SCW CESAs (Figure 8; Supplementary Data Sheet 4). However, the amount of cellulose accumulated per unit fiber length might be lower in Gb compared to Gh, as reported earlier (Avci et al., 2013), and it might have resulted in longer fiber elongation time and longer and thinner fiber in Gb. In our future study, we will investigate fiber length, growth rate, cellulose accumulation and microfibril orientation to relate them to fiber quality traits.

Some unanswered questions remain in order to understand the whole dynamics and contribution of CESAs to fiber quality traits, such as 1) Would the temporal patterns of transcripts of CESAs that we observed directly correlate with protein/enzyme level and dynamics? 2) How much do the differentially accumulated PCW and SCW CESAs contribute to cellulose crystallinity? 3) What is the role of PCW and SCW crystalline cellulose in temporal plasticity, rigidity, elongation, and strength of fiber in both species? 4) How does the cellulose in PCW and SCW interact with other CW components and thus contribute to fiber quality?

Comparative analysis of the polysaccharide composition and related glycosyltransferases between Gh and Gb fiber revealed many interesting details that might be the cause for the fiber quality differences between these cotton species. First, it is worth pointing out the differences in the total amount of buffer-soluble and alkali-soluble polysaccharide fractions. In Gh, the content of buffer-soluble polysaccharides was higher than that of alkali-soluble polysaccharides at 6–7 DPAs, whereas in Gb the content of these two fractions were more comparable (Figure 1). However, the level of buffer-soluble polysaccharides reduced more rapidly in Gh than in Gb during development, and at later stages however, the relative proportion of buffer-soluble polysaccharides remained at the same level in the two species. These differences reflected in the differences in polysaccharide epitope distribution, discussed below.

Quantitative heat maps and SOM analyses of polysaccharide epitope profiles revealed significant differences between the two species, which guided a comparative analysis of the glycome and glycosyltransferase transcriptome data. The most interesting aspect of our findings was that some of the same polysaccharide epitopes from Gh and Gb fell into two different SOM groups, which suggests that polysaccharide distribution differs in fiber CW between the two species. Although some of these particular differences appeared insubstantial and statistically insignificant, even subtle differences, taken together, could potentially impact fiber growth and development. In particular, the polysaccharide epitopes from the second and third categories (Figures 6, 7) may be promising targets for further detailed investigations.

Our large-scale glycome profiling showed that most of the buffer and alkali-extracted fucosylated and non-fucosylated XG epitopes were very low in Gb in comparison with Gh (Figures 6, 7), which agrees with similar earlier comparative studies (Avci et al., 2013; Hernandez-Gomez et al., 2017; Guo et al., 2019). These earlier studies reported that the presence of XG was lower in cotton fiber middle lamellae (CFML, the special adhesive outer layer of PCW) and PCW from Gb than in Gh. The authors proposed that the role of CFML is to adhere to the adjacent fibers, which facilitates the formation of organized tissue-like fiber bundles around each seed. This orderly fiber packing facilitates the elongation of fibers to a maximum within the confined space inside each locule of the cotton boll (Singh et al., 2009). However, Avci et al. (2013) reported that CFML was not required for fiber elongation in Gb, since it continued to elongate rapidly even after lysis of CFML. Additional studies will be required to confirm the correlation between the amount of XG and the dynamics of CFML/PCW development in promoting fiber length in different species of cotton.

Interestingly, in our transcriptome profiling study, we found that the transcript levels of all glycosyltransferases related to XG biosynthesis (CSLCs/XXTs/MUR3s/XLT) (Figure 9; Supplementary Data Sheet 4) were lower in Gb than in Gh, which might have resulted in lower XG enzymatic activity and consequently less XG epitope content in Gb than in Gh. A previous study by Guo et al. (2019) using fiber at 10, 15, and 20 DPA showed that the transcript levels of a CSLC4 and a MUR3 were lower in Gb than in Gh. In our transcriptome data, we found that the transcripts of a significantly larger set of various XG-synthesizing glycosyltransferases from both At and Dt genomes accumulated less in Gb compared to Gh. It will be interesting in future studies to examine the contribution of all these enzymes to the fiber quality in Gb.

It was presumed previously that XG interacts with cellulose microfibrils through non-covalent hydrogen bonds that could potentially contribute to rigidity to the CW and restrain cell expansion (Cosgrove, 2022; McCann et al., 1990; Somerville et al., 2004; Park and Cosgrove, 2015; Pękala et al., 2023). Interestingly, some of the earlier studies showed that the short fiber mutants (Ligon lintless-1 & 2; Li₁/Li₂) had a higher accumulation of XG content during elongation stages due to higher expression of XG-synthesizing glycosyltransferases, which presumably resulted in an extreme reduction in fiber length (Shao et al., 2011; Naoumkina et al., 2017). Overexpression of a MUR3 (GhMUR3-2Dt; Gohir.D10G136900) protein in Gh cotton resulted in shorter and thicker fibers compared to wild type (Wu et al., 2024). All data obtained so far suggests that a higher amount of XG negatively correlates with fiber length, most likely impacting the “xyloglucan-cellulose interaction”, which may lead to an increase in the rigidity of fiber CWs. Thus, the results suggest that a lower amount of XG seems to loosen this interaction, which might have helped extend the elongation time and promote longer fibers in Gb in comparison with Gh.

Another hemicellulosic polysaccharide, xylans (Xyl), is the major component in SCWs, and due to its interaction with cellulose and lignin, Xyl is essential for SCW architecture and strength (Gille and Pauly, 2012). In our previous study, we demonstrated that some of the Xyl epitopes (Xyl-3Ar, Xyl-MeGlcA) and transcripts of associated glycosyltransferases rapidly peak at the transition stage in Gh fiber before the rapid SCW cellulose accumulation/CW thickening and the changes in microfibril orientation (Swaminathan et al., 2024). Here, we found that the amount of loosely bound buffer-extractable Xyls, mostly represented by the Xyl backbone epitopes (Xyl-BB) and a glucuronoxylan epitope (Xyl-GlcA), was lower in Gb than in Gh (Figure 6) and might contain undetectable Xyls with a high degree of substitution. The alkali-extractable arabinoxylans (Xyl-2Ar, Xyl-3Ar) and a methylated glucuronoxylan (Xyl-MeGlcA-2) showed different profiles between Gh and Gb (Figure 7), and Xyl-MeGlcA-2 was present in higher amounts in Gb than in Gh. Somewhat similar differences in Xyl-related epitope profiles were reported earlier (Avci et al., 2013; Hernandez-Gomez et al., 2015).

We compared the transcript levels of Xyl-synthesizing glycosyltransferases between the two species. Transcript levels of most of the IRXs/FRA8/PARVUSs and ESKs/TBLs/RWAs genes, which produce enzymes that are involved in Xyl backbone synthesis and decorating Xyl backbone with acetyl/methyl groups, respectively, were higher in Gh than in Gb. This correlates well with the presence of higher amounts of Xyl backbone epitopes in Gh than in Gb (Figure 10). On the other hand, the transcript levels of some of the GUXs and GXMTs, which are responsible for the synthesis of methylated glucuronoxylans, were higher in Gb than in Gh, which correlates with the higher content of Xyl-MeGlcA-2 epitopes detected in Gb (Figure 7).

Recent reverse-genetic studies in cotton with some of the Xyl-synthesizing glycosyltransferases (FRA8/IRX9/IRX10/IRX14/IRX15) from Gh showed that Xyls play an important role in cellulose microfibril orientation, SCW cellulose deposition, fiber length, and thickness (Li et al., 2014; Chen et al., 2020; Guo et al., 2023; Li et al., 2024). In our previous large-scale profiling study (Swaminathan et al., 2024), we also observed that three heteroxylan epitopes were highly correlated with the cellulose content, the transcript levels of SCW CESAs, the cellulosic microfibril orientation and the CW thickness phenotype of Gh fiber. The function of heteroxylans (glucuronoxylans and arabinoxylans) is usually related to the strengthening of CW by acting as a guiding scaffold for cellulose microfibril orientation/arrangement (Grantham et al., 2017; Smith et al., 2017; Crowe et al., 2021; Pfaff et al., 2024). The strength of the cotton fiber CW is an important quality essential for the textile industry processing. The differential profiles and contents of methylated glucuronoxylan (Xyl-MeGlcA) and arabinoxylans (Xyl-2Ar and Xyl-3Ar) observed in Gb could potentially contribute to its stronger fiber. In our proteome data, we also noticed that several IRXs, GUXs, GXMTs, and XATs were abundantly present and matched their corresponding transcript profiles (Figure 2; Supplementary Data Sheet 5). Also in our study, the transcript levels of GUXs and GXMTs, were found to be higher in Gb than in Gh. Additional research is still required to understand the temporal expression dynamics of IRXs, GUXs, GXMTs, XATs and heteroxylans accumulation in fiber and their impact on fiber quality. Interestingly, an earlier study of Arabidopsis mutants showing various levels of glucuronoxylan deficiency demonstrated that the levels of glucuronoxylans critical for the precise cellulose network formation and CW integrity of SCWs (Crowe et al., 2021). This might suggest that temporal synchrony in the presence of appropriate heteroxylans quality and quantity in combination with the other polysaccharides, affects the CW architecture and, thus, the final quality of the fiber.

Pectins are the major components of the PCW and CFML (up to 35% dry weight) of the cotton fiber (Kaczmarska et al., 2022). There are four main structural components of pectins: HGs, RG-I, rhamnogalacturonan-II, and xylogalacturonans. The cellulose and hemicellulose networks are embedded in a matrix of pectins and proteins in the PCW. HGs and RGs play a central role in regulating the viscosity/extensibility/gelling property of the CW matrix, thus controlling polysaccharide interactions, CW elongation, and cell growth (Hwang and Kokini, 1992; Yapo, 2011; Sousa et al., 2015; Zheng et al., 2020). Our results from glycome profiling showed that a methyl-esterified HG epitope (HG-BBMe-K) and most of the highly branched RG-I epitopes were less abundant in Gb relative to Gh. In contrast, the de-esterified HGs (HG-BBde), and RG-I backbone epitopes (RG-I-BB, Gal4-BB) were present in equal amounts in both species. Similar observations were reported earlier by Avci et al. (2013). Another comparative “omics” study (Liu et al., 2013) reported a lower amount of methyl-esterified HG in Gb compared to Gh during elongation time and proposed that pectin content might be an important factor influencing fiber elongation and quality. Our results suggest that the presence of a lower content of HGs and highly branched RG-Is, in addition to the changes in other polysaccharide epitopes during the PCW developmental/fiber elongation stage in Gb, could support its longer fiber phenotype.

Our comparative transcript analysis showed that transcript levels of most of the glycosyltransferases (that synthesize the HG epitopes), including methyl and acetyltransferases (GAUTs, GATLs, CGRs/QUAs, TBLs/TBRs) were present in higher amounts in Gh than in Gb. A similar trend was observed for RG-I-synthesizing glycosyltransferases, most likely resulting in a higher amount of pectin epitopes in Gh than in Gb. The cotton AG: GALT (Gohir.A04G041100/GhAG: GALT1-2At) protein was found to play a functional role in galactosylation of AG chains of RG-I/AGPs molecules (Qin et al., 2017). The AG: GALT1-2At RNAi-silenced cotton plants, had less galactose content in the RG-I and produced longer fibers as a result of CW loosening. Conversely, AG: GALT1-2At protein overexpressing cotton plants had the opposite effect. Overexpression of an HG- degrading polygalacturonase in Arabidopsis resulted in reduced HG content, which resulted in a longer hypocotyl and this might be because of increased CW loosening (Xiao et al., 2014). Pectate lyases (PEL) are another type of hydrolase involved in the degradation of HGs. Silencing of a cotton GhPEL resulted in accumulation of higher amounts of HGs in fiber CW during the elongation stage, which subsequently resulted in a shorter fiber phenotype in transgenic plants (Wang et al., 2010).

All these results indicate that higher amounts of pectin epitopes during the elongation stage lead to possible broader hydrogen bonding interactions (as reported in in vitro studies; Zheng et al., 2020; Chandel et al., 2022; Said et al., 2023) or calcium-mediated “egg-box” formations (Cosgrove, 2022) between the pectin backbone molecules, both can lead to increased firmness of CW, which could negatively impact the elongation time and length of cotton fibers. A more detailed comparative investigation of glycosyl hydrolases contributing to cotton fiber development will be conducted in our future study.

Expansins are CW-loosening proteins that act in a non-enzymatic manner to weaken the noncovalent bonds in hemicellulose-cellulose networks, resulting in slippage of cellulose microfibrils, and ultimately leading to cell expansion (Sampedro and Cosgrove, 2005) and one of the major factors responsible for cotton fiber elongation (Li et al., 2016; Lv et al., 2020). Therefore, here, in addition to glycosyltransferases, we examined the transcript levels of known expansins in both Gh and Gb fibers. Interestingly, the transcript levels of most of the expansins was higher in Gb relative to Gh (Figure 13; Supplementary Data Sheet 9). Additionally, the transcript levels of some of the expansins was significantly higher in Gb through the later stages (after 17 DPA), which coincides with the extended elongation stage in this cotton genotype. In contrast, in Gh, the transcript levels of most expansins decreased at 17 DPA, coinciding with the transition stage before SCW thickening and complete cessation of fiber elongation in this genotype. Our results suggest that the temporally controlled higher level of expansin expression plays an important role in extended fiber elongation and the development of superior fiber quality in Gb.

Previous reverse-genetics studies of EXPA8-2 (expansin with the highest transcript accumulation levels from both Gh and Gb) by Xu et al. (2013) and Li et al. (2016) showed that the cotton fiber length is increased in GbEXPA8–2 overexpressing lines and reduced in the silenced lines. Overexpression of GbEXPA8–2 protein also altered the expression of SCW-associated genes and delayed the accumulation of SCW cellulose, which resulted in extended elongation time of the PCW stage and long fibers in Gb. A separate transgenic experiment in Gh by overexpressing GhEXPA8–2 protein along with a BURP-domain containing protein resulted in the increased fiber length (Xu et al., 2013).