Zhngqing Song,Yu Chen, Chunyun Zhng, Jingxi Zhng, Xuehn Huo,,Yng Go, Ao Pn,c, Zhohi Du, Jun Zhou, Ynxiu Zho, Zhi Liu,c,Furong Wng,,*, Jun Zhng,,*

aKey Laboratory of Cotton Breeding and Cultivation in Huang-Huai-Hai Plain, Ministry of Agriculture and Rural Affairs, Cotton Research Center,Shandong Academy of Agricultural Sciences,Jinan 250100,Shandong,China

bCollege of Life Sciences,Shandong Normal University, Jinan 250014, Shandong,China

cCollege of Bioscience &Biotechnology, Hunan Agricultural University, Changsha 410128,Hunan,China

A B S T R A C T

1. Introduction

Cotton is an important crop used in the production of natural fibers as industrial raw materials.The cotton genus Gossypium contains ＞45 diploid(2n=26)and seven allotetraploid(2n=52)species, of which four are grown to produce fiber, including the two allotetraploids G. hirsutum L. and G. barbadense L. and the two diploids G. herbaceum L. and G. arboreum L. [1].Gossypium hirsutum, commonly referred to as upland cotton,is responsible for over 95% of annual cotton fiber output worldwide because of its high yield and wide adaptability. In contrast,G.barbadense,or sea island cotton,produces higherquality fibers with lower yield[2].Because these fiber qualities may have a genetic origin,using genetic resources to develop a high-yield upland cotton with excellent fiber qualities is a promising breeding strategy [3]. Several germplasm accessions with improved fiber quality that contain introgressed components from other Gossypium species [4-7], especially cultivated G. barbadense varieties [8-11], have been developed using wide crosses (interspecific hybridization).

Cotton fibers are single-celled trichomes that develop from ovule epidermal cells.They are an excellent model system for studies on plant cell elongation, as well as cell-wall and cellulose biosynthesis [12]. The process of cotton fiber development can be divided into four overlapping stages:initiation, elongation, secondary cell wall (SCW) synthesis,and maturation[12,13].During the elongation stage[3-23 days post-anthesis(DPA)],fiber cells rapidly elongate and gradually reach their maximum length [14]. From 15 to 40 DPA, large amounts of cellulose are synthesized in the SCW synthesis stage,making up ＞90%of the dry weight of mature fibers[15].Other noncellulose polysaccharides contributing lower content also play important roles in the synthesis of SCWs.These include hemicelluloses, which are indispensable for the normal assembly and mechanical strength of SCWs[16,17].

Cellulose is the largest structural component of mature fibers and is synthesized on the plasma membrane by the cellulose synthase complex, which includes 18 cellulose synthase (CESA) proteins [18]. In the Arabidopsis genome,there are 10 CESA genes,and CESA6 is involved in primary cell wall (PCW) cellulose biosynthesis [19,20]. The Arabidopsis genome encodes 30 cellulose synthase-like (CSL) genes,which are divided into six subgroups (CSLA, CSLB, CSLC,CSLD, CSLE, and CSLG) [20,21]. In Gossypium species, 228 CES/CSL genes have been identified and grouped into 11 subfamilies [22]. Although their functions await further elucidation,CSL genes that encode glycosyltransferase family 2 proteins are postulated to participate in the synthesis of non-cellulose polysaccharides, including the polysaccharide backbones of hemicelluloses[23].

Fiber initiation and development are controlled by plant hormones that play crucial roles during growth. Auxin and gibberellins(GAs)promote the growth of cotton fibers in vitro[24]. Auxins also play important roles in cell wall elongation and SCW deposition by loosening the cell wall [25]. Applied together, auxin and cytokinin promoted the expression of NCA transcription factors VNDs [26]. Abscisic acid (ABA)induced the expression of NST3/SND1, VND4, VND5, and VND2 [27]. The promoter of GbEXPA2, a cotton fiberpreferential gene that is strongly expressed during fiber development, was regulated by exogenous GA and ABA in Arabidopsis [28]. The transcription factor bri1-EMSSUPPRESSOR 1 (BES1) was activated by brassinosteroid receptor-like kinase BRI1 to promote the expression of the cellulose synthase gene CesAs [29]. In addition to secondary cell wall biosynthesis,ethylene responses may also indirectly influence fiber strength by promoting fiber-to-fiber interactions[30].

The genomes of G. raimondii[31,32],G. arboretum[33,34],G.hirsutum [11,35,36], and G. barbadense [11,37] have recently been sequenced, enabling studies of the molecular mechanism of fiber development. Recent research has identified several genes involved in the regulation of cotton fiber initiation and development. For example, GhHD1 regulates lint fiber initiation, and its deletion results in a delay in fiber initiation [38,39]. GhHD1 also interacts with GhHOX3 to promote fiber elongation by relieving the competitive binding of GA repressor DELLA to GhHOX3 [40]. In G. hirsutum,overexpression of GbEXPATR from G. barbadense results in longer, stronger but thinner fibers [41]. However, the molecular mechanisms that regulate fiber strength and secondary wall deposition remain unclear.Better understanding of these mechanisms would facilitate the genetic improvement of fiber quality.

The genetic and molecular basis of high-quality fiber traits can be investigated using genome-wide transcriptome profiles of relevant genes and pathways[14,42,43].Transcriptome analyses have shed light on the dynamic regulation of fiber development [30]. However, few studies have focused on germplasm with superior fiber quality carrying G. barbadense chromosomal fragments in a G.hirsutum background[14,42].

Development of chromosome segment substitution lines(CSSLs) is an effective strategy for improving cotton fiber quality by use of an exogenous gene pool. In CSSLs, the genome is partially substituted with one or few chromosome segments from the donor parent following distant hybridization [44]. In preliminary work in our laboratory, Luyuan343(LY343), a line with excellent fiber quality conferred by G.barbadense genomic introgressions [8,9] was used as a donor parent and backcrossed to recurrent parent Lumianyan 22(L22) to generate a single chromosomal segment substitution line (SCSSL). The single introgressed segment contains a cluster of major quantitative trait loci (QTL) increasing fiber length strength, and fineness and lying between markers DPL0757 and DC40182 on chromosome A07, and was named SL7.

The present study was undertaken to (1) confirm the association of the introgressed G. barbadense chromosomal segment with improved fiber quality and (2) investigate whether transcriptional regulation of morphogenesis and major components of the cell wall is associated with alterations in fiber quality,particularly fiber strength.

2.Materials and methods

2.1. Plant materials

All cotton materials were planted at the Linqing Experimental Station of Shandong Cotton Research Center (LES/SCRC) and on Hainan Island in winter.Lumianyan 22(L22)is a transgenic insect-resistant upland cotton cultivar developed by the Cotton Research Center of Shandong Academy of Agricultural Sciences. Luyuan343 (L343) is an introgressed line with excellent fiber quality produced by a natural cross between G.barbadense and an upland cotton[8,9].SL7 was obtained by backcrossing the donor parent L343 to the recurrent parent L22 for five generations, followed by selfing for two generations. To characterize the genetic background, 307 SSR markers evenly distributed across the genome were used to identify introgressed chromosomal segments(Fig.S1).

For evaluating fiber composition, fibers were harvested from ＞10 plants of L22 or SL7, mixed, and divided into three biological replicate samples of 3 g each. Cellulose, hemicellulose,and lignin content were determined by Aobaike Biotechnology Co.,Ltd.(Qingdao,China)using respectively anthronesulfuric acid colorimetry, hydrochloric acid hydrolysis, and titrimetry, following standard procedures. The results were represented as content per unit weight fiber(%).

2.2. RNA extraction, library construction, and transcriptome sequencing

Ovule samples at 0 and 5 DPA and fiber samples at 10,15, 20,and 25 DPA were harvested at the Linqing Experimental Station and frozen in liquid nitrogen. High-quality RNA extraction was performed using an RNAprep Pure Plant Kit(Huayueyang, Beijing, China). RNA degradation and contamination were monitored on 1% agarose gels. RNA purity and concentration were checked with a spectrophotometer(NanoPhotometer,Westlake Village,CA,USA).Approximately3 μg RNA per sample were used to generate sequencing libraries using a NEBNext Ultra RNA Library Prep Kit (NEB,Ipswich, MA, USA). Briefly, mRNA was purified using poly-T oligo-attached magnetic beads. Reverse transcription was performed with M-MuLV Reverse Transcriptase (RNase H-)(NEB, Ipswich, MA, USA), DNA polymerase I and RNase H.After adenylation of the 3′ends,DNA fragments with hairpin loop structure were ligated. cDNA fragments of 150-200 bp were selected with an AMPure XP system (Beckman Coulter,Brea, Calif., USA). PCR was performed and the products were purified with the AMPure XP system. Library quality was detected with the Agilent Bioanalyzer 2100 system (Agilent Technologies Inc., Santa Clara, Calif., USA). Clustering of index-coded samples was performed with a TruSeq PE Cluster Kit v3-cBot-HS (Illumina, San Diego, Calif., USA). The library was then sequenced on an Illumina HiSeq platform and 125-150 bp reads were generated.

Table 1-Fiber quality trait values in L22 and SL7.

2.3.Data quality control and mapping of reads to the reference genome

The original reads were first converted to sequencing reads,including the base quality of reads through base calling.Clean reads were obtained by removal of all low-quality reads,including those containing only adapters or poly-Ns,using the Illumina pipeline.The Q30 and GC contents of the clean reads were calculated. All subsequent analyses were based on the remaining clean reads.

The reads were aligned with the G. hirsutum reference genome sequence obtained from the CottonGen database(https://www.cottongen.org/)[36].Reference genome indexing and read alignments were performed using Bowtie 2.0.6 and TopHat 2.0.9(https://magic.novogene.com/),respectively.

2.4. Differentially expressed gene (DEGs) analysis

Genes differentially expressed between two fiber samples were identified with the DESeq R package (1.18.0) (https://magic.novogene.com/). Differential expression in the digital gene expression data was determined using statistical routines based on the negative binomial distribution.Genes with P-value ＜0.05 were considered to be differentially expressed.

2.5. GO and KEGG enrichment analysis of DEGs

The GOseq R package (https://magic.novogene.com/) was used to perform GO enrichment analysis of DEGs, GO terms with corrected P-value ＜0.05 were considered significantly enriched [45]. KEGG statistical enrichment of DEGs was evaluated with KOBAS software (https://magic.novogene.com/). DEGs were gathered into the KEGG pathways by matching in the database.

2.6. qRT-PCR analysis

Fig.1- Cellulose(A),lignin(B), and hemicellulose(C)content of SL7 and L22 mature fibers.* P ＜0.05. Error bars represent standard deviation.

qRT-PCR analysis of 12 DEGs was performed to test the reliability of the transcriptome sequencing data. Reverse transcription was performed using PrimeScript RT reagent Kit (TaKaRa, Kusatsu, Japan). qRT-PCR was performed on a Light Cycler 480 II (Roche, Basel, Switzerland) using SYBG Premix Ex Taq II (TaKaRa, Kusatsu, Japan). The primer sequences of 12 DEGs are shown in Table S1. The actin gene was used as internal reference (using the primers 5′-GGTGGTGTGAAGAAGCCTCAT-3′ forward and 5′-AATTTCACGAACAAGCCTCTGGAA-3′ reverse). Relative gene expression levels were calculated according to 2-ΔΔCT[46].

3. Results

3.1. Phenotypic evaluation

Fig.2-Heat map of DEGs from L22 and SL7 samples.Numbers 0,5,10,15,20,and 25 represent respectively 0,5,10,15,20,and 25 DPA.DEG:differentially expressed gene;DPA:days post-anthesis.

The fiber quality traits of SL7 and L22 were evaluated and compared to verify that the fiber quality of SL7 was superior to that of its recurrent parent L22. The fiber strength (FS) of SL7 was 11.8% higher than that of L22. SL7 also showed higher fiber length (FL) and fiber fineness (micronaire) (FM) than L22 (Table 1),indicating that SL7 inherited the excellent fiber properties of L343.

The cellulose contents in both L22 and SL7 reached 90%(Fig. 1-A). Differences between L22 and SL7 in cellulose and lignin contents were not significant, although these contents were slightly higher in SL7 (Fig. 1-A, B). However, the hemicellulose content was significantly higher in SL7 than in L22,probably resulting in an increase in SL7 fiber strength(Fig.1-C).

3.2. Transcriptome sequencing and alignment to the reference genome

RNA-seq libraries of ovule and fiber samples of L22 and SL7 at 0 and 5 DPA and fiber samples at 10, 15, 20, and 25 DPA were sequenced. Quality control yielded 1391 million clean reads,accounting for 95.5% of the raw reads (Table S2). Proportions of reads whose base calling error rate was ＜0.001 ranged from 87.1% to 93.4%, and GC content was ＞42.7%, indicating highquality RNA sequence(Table S2).

Fig.3- DEGs in L22 and SL7 and GO functional classification of DEGs.(A)Multiple comparisons of gene expressions among different fiber development stages of SL7 and L22.Red numbers represent upregulated genes and blue numbers represent downregulated genes.(B) Venn diagram of significant DEGs between L22 and SL7.Numbers 0,5,10,15,20, and 25 represent DPA.(C)GO functional classification of DEGs between L22 and SL7 at 10 DPA.DEG:differentially expressed gene;DPA:days post-anthesis.

Table 2-DEGs involved in cellulose and hemicellulose synthesis.

Table 3-DECs involved in cell wall synthesis.

Table 4-Numbers of transcription factor genes differentially expressed between L22 and SL7.

Of the 90.7% of clean reads mapped to the G. hirsutum reference genome,84.5%were uniquely mapped(Table S2).In RNA-seq analysis, FPKM (expected number of fragments per kilobase of transcript sequence per million base pairs sequenced)is used to show gene expression levels.About half of the mapped genes were obviously expressed in our samples(FPKM ＞1) at every fiber development stage (Table S3),indicating that the biological processes in fiber cells were still relatively active even in the later stage of secondary wall thickening. Pearson's correlation coefficient (PCC) analysis indicated that gene expression at the same fiber developmental stage between L22 and SL7 showed ＞90% correlation degree with a minimum of 0.915 at 25 DPA,indicating that the expressions of genes in L22 and SL7 at 25 DPA were more different than that at other stages(Fig.S2).Over the course of 0-25 DPA in the same line,the PCCs between 0 and 5 DPA were similar but showed sharp decreases at 10,15,20,and 25 DPA,possibly owing to the alterations in gene expression profiles in fiber and ovule cells.

3.3. DEG analysis and GO enrichment of DEGs

The read counts of all detected genes in each sample were used to identify DEGs (|log2fold| ＞1, P ＜0.05), and cluster analysis was performed to examine the expression patterns of DEGs at different developmental stages. The expression patterns of genes at 0 and 5 DPA were dramatically different from those at 10, 15, 20, and 25 DPA (Fig. 2),supporting the PCC findings. Interestingly, genes highly expressed in early ovule cells were generally suppressed in fiber cells undergoing rapid elongation and secondary wall thickening (Fig. 2). DEG expression profiles at 10 and 15 DPA also showed clear differences, indicating that the expression profiles of genes were altered during fiber development. There were also significant differences in DEG expression levels at the same stage between L22 and SL7,which might account for the regulatory mechanism of fiber quality improvement in SL7.

Multiple comparisons between fiber developmental stages were performed to identify changes in expression of DEGs.As many as 1437, 353, and 590 genes were induced in SL7 compared with L22 at 10, 20, and 25 DPA, respectively (Fig. 3-A), significantly greater than the numbers of downregulated genes. After 15 DPA, changes in gene expression levels were higher in SL7 than L22 (Fig. 3-A). There were many DEGs specific to SL7 in the secondary wall thickening period(5712 in 20 vs.15 DPA and 4955 in 25 vs.20 DPA),although many DEGs were common to L22 and SL7(Fig.3-B).

GO enrichment analysis was performed to identify genes involved in fiber development among the numerous DEGs.DEGs were classified by biological process, cellular component, or molecular function. We focused on DEGs at the fiber elongation and secondary wall thickening stages. In the biological process category, most DEGs were assigned to terms describing cellular process and metabolic process,especially at 15, 20, and 25 DPA (Fig. 3-C, Fig. S3-A-C). A few DEGs were involved in responses to stimuli and signaling at 10 and 15 DPA, indicating that the development of cotton fibers was associated with environmental changes through various pathways (Fig. 3-C, Fig. S3-A). For the category of cellular components, the majority of DEGs were located in the membrane or membrane-bound organelles. In the molecular function category, most of the proteins encoded by DEGs had binding and catalytic activity, suggesting that many were enzymes involved in various metabolic processes (Fig. 3-C,Fig.S3-A-C).

3.4. CSLs and other genes involved in cell wall syntheses are upregulated in SL7

We further examined DEGs involved in cellulose and cell wall syntheses. The expressions of one cellulose synthase A[Gh_D05G2313 (CESA6)] and nine cellulose synthase-like protein (CSL) were significantly upregulated compared with L22 in SL7 fiber cells at 10,15,20,and 25 DPA(Table 2). Several of these CLS genes were highly expressed in SL7 at more than two fiber developmental stages. For example, the expression of Gh_D02G0415(CSLG2)was upregulated in SL7 at 10,15,and 25 DPA, and Gh_D09G2165 (CSLE1) and Gh_D09G2166 (CSLG3)were also highly expressed in SL7 at 15 and 25 DPA compared with L22(Table 2).Because CSLs are involved in the synthesis of non-cellulose polysaccharides, these findings suggest that during fiber development in SL7, the biosynthesis of noncellulose polysaccharides such as hemicelluloses was more active than that in L22.Thus,the induced expression of these CSLs and CESA genes might be responsible for the observed superior fiber quality of SL7.

Other genes involved in cell wall synthesis were also induced in SL7 at 10,15,20,and 25 DPA,including members of xyloglucan endotransglucosylase/hydrolase (XTH), pectinesterase (PME), and endochitinase (Table 3). Eleven XTH genes were upregulated in SL7 during fiber development (Table 3),and five PME and four endochitinase genes were also highly expressed in SL7. Most of the DEGs involved in cellulose and cell wall synthesis were upregulated in SL7 compared with L22,suggesting that the high fiber quality of SL7 is attributable to higher metabolic activity.

In the molecular function category, a GO item contained 24-118 transcription factor DEGs(118 in 10 DPA,63 in 15 DPA,24 in 20 DPA,and 36 in 25 DPA),and more than a third of these transcription factors belonged to three families: WRKY, MYB,and NAC. Most of these transcription factor genes were upregulated in SL7 (Table 4). Forty-three MYB and 31 NAC genes were differentially expressed in SL7 compared with L22 at 10-25 DPAs, most of which were upregulated in SL7,especially at 25 DPA (Fig. 4-A, B). These findings suggest that these transcription factors are involved in the regulation of secondary wall synthesis. Thirty-five WRKY genes were differentially expressed between SL7 and L22 (Fig. 4-C), with most induced in SL7 at 10,15,and 20 DPA.The upregulation of WRKY, MYB, and NAC TFs suggests that secondary cell wall synthesis is more active in SL7 than in L22 and that there is some relationship between cotton fiber development and environmental response.

3.5. Hormone signal transduction pathways are involved in high-strength fiber formation in SL7

KEGG analysis was used to identify pathways in which the DEGs between L22 and SL7 participated.The 20 pathways with the most significant enrichment at every measured DPA time point were plotted (Fig. 5). The metabolic pathways were enriched the most by DEGs at every point but 10 DPA, when plant hormone signal transduction was the most significant pathway. At 15 DPA, the plant hormone signal transduction pathway still contained 33 DEGs(Fig.S4).Fiber development is controlled by a variety of plant hormones, and our results shed light on the regulatory mechanisms underlying the generation of excellent fiber quality.

Plant hormone signal transduction was investigated at 10 DPA during fiber elongation, and up to 70 DEGs were enriched in the plant hormone signal transduction (Fig. 5).These were divided into a series of hormone signaling pathways, and most genes were upregulated in SL7 (Fig. 6).Of the 70 DEGs, 46 were involved in the auxin signaling pathway, suggesting that auxin promotes fiber elongation.For example, the genes encoding AUX1, TIR1, AUX/IAA, and GH3 were possibly expressed in SL7 to promote cell enlargement (Fig. 6). In the brassinosteroid signaling pathway, the genes encoding BAK1, BRI1 and BKI1 were upregulated in SL7 (Fig. 6). At 15 DPA during secondary wall thickening, the genes involved in the auxin signaling pathway were downregulated in SL7 (Fig. S5), suggesting that auxin does not regulate the synthesis of secondary cell walls. The genes involved in the ethylene and ABA signaling pathways were also upregulated in SL7 as well as at 10 DPA (Fig. 6), suggesting that ethylene and ABA also promote secondary wall thickening. These results suggest that the exogenous chromosome segments in SL7 introgressed from G. barbadense might improve fiber quality via hormone signal transduction pathways.

3.6. DEG verification by qRT-PCR

To verify the DEGs identified by RNA-seq, several DEGs described in Tables 2-4 were selected and their expressions were measured by qRT-PCR (Fig. 7). The results showed that most of these DEGs were upregulated at the secondary wall thickening stage, including Gh_A08G0195 (XTH22),

Fig.5-KEGG pathway analysis of DEGs between L22 and SL7 at 10 DPA.Sizes of dots correspond to numbers of genes,and their colors correspond to q-values of pathway enrichment.DEG:differentially expressed gene;DPA:days post-anthesis.

Gh_A03G0887 (NAC029), Gh_A08G2417 (WRKY41), and Gh_A08G1691 (NAC022) (Fig. 7), which exhibited high expression levels at 20 and 25 DPA in both L22 and SL7,suggesting that the proteins encoded by these genes play important roles in secondary wall thickening. Moreover, in agreement with the RNA-seq results, these 12 DEGs were induced at certain DPAs in SL7, in contrast to L22,confirming the reliability of our RNA-seq results.

4. Discussion

Cotton fibers are differentiated from single cells in the ovule epidermis, and mature fibers are composed mainly of cellulose and small amounts of noncellulosic polysaccharides such as hemicellulose [15]. Fiber quality is determined by three main factors:length,strength,and fineness.Among the four cotton cultivars, G. barbadense shows superior fiber quality but lower yield, suggesting that its genome contains genetic resources for superior fiber quality. In this study, it was indicated that SL7 inherited the excellent fiber properties of L343. SL7 could be used as an ideal material for the identification of functional genes related to fiber development and further research on its underlying mechanisms.

Fig.6-Regulation of plant hormone signal transduction pathways in SL7 compared with L22 at 10 DPA.Red box indicates that gene expressions in the pathway were upregulated in SL7.Yellow box indicates that the genes of this node were both upregulated and downregulated in SL7.DPA:days post-anthesis.

RNA-seq is an effective method for detecting changes in the expression of relevant genes and pathways at the genomic level. Thousands of DEGs between L22 and SL7 during fiber development were identified (Fig. 2). Genes involved in cell wall syntheses were upregulated in SL7 (Table 2), including several CSLGs, members of the CSL family involved in the synthesis of non-cellulose polysaccharides including hemicellulose, which are synthesized in the Golgi apparatus and transported to the cell wall[22].The hemicellulose content in SL7 was significantly higher than that in L22, possibly accounting for the greater fiber strength of SL7. Cellulose content, the largest component of mature fiber, was not significantly different between SL7 and L22, in keeping with the similarity of expression of cellulose synthase genes.These results indicate that the arrangement and crosslinking of cellulose, which may be associated with hemicellulose content, are significant for fiber quality, especially fiber strength. A series of XTH and PME genes were also differentially expressed in SL7 and L22.In cotton,23 XTH cDNAs have been identified, and XTH regulates fiber elongation by catalyzing the cleavage and reconnection of cell wall xyloglucans [47-49]. PME and PME inhibitors (PMEIs) are thought to regulate biomechanical properties of the cell wall synergistically by controlling the methyl esterification of homogalacturonan in pectin [50]. The differential expression of these genes may also be an important reason for the fiber quality improvement in SL7.

Fig.7-qRT-PCR validation of 12 DEGs.Columns indicate the results of qRT-PCR;colored line segments represent the results of transcriptome sequencing.DPA,days post-anthesis.For qRT-PCR,*,P ＜0.05;**,P ＜0.01;Error bars represent standard deviation.DEG:differentially expressed gen.

The process of SCW synthesis is regulated by a variety of transcription factors,among which two families,the NAC TFs and MYB TFs, constitute the core of the regulatory network[51]. NAC and MYB TFs play vital roles in the regulation of SCW deposition. Five NAC TFs, namely NST1, SND1, NST2,VND6,and VND7,promote the synthesis of SCW by regulating the expression of a series of MYB TF genes, which directly regulate the expression of cellulose,hemicellulose,and lignin biosynthetic genes [28,51-53]. WRKY is a large plant-specific TF family with 74 members in Arabidopsis and participates in the regulation of various developmental processes and stress resistance [54]. RNA-seq analysis revealed that several transcription factor genes belonging mainly to the WRKY, MYB,and NAC families (Table 4) were differentially expressed in L22 and SL7 (Fig. 3-C, Fig. S3-A-C). These genes might influence fiber development by regulating the expression of genes encoding enzymes such as CSLs and XTHs,which were also differentially expressed between L22 and SL7.

Cotton fiber development is regulated by several plant hormones[13].At 10 DPA,70 DEGs between L22 and SL7 were enriched in signal transduction pathways of various plant hormones, including auxin, cytokinin, ethylene,brassinosteroid, and ABA (Figs. 5, 6). However, at 15 DPA,only 33 DEGs were involved in plant hormone signal transduction pathways (Fig. S4), indicating that plant hormones play key roles in the transition from the fiber elongation stage to the secondary wall thickening stage.A model(Fig.8)can be proposed to explain the possible mechanism of fiber strength increase in SL7. In SL7, exogenous chromosome segments introgressed from G. barbadense alter the expression profile,including the upregulation of plant hormone signal transduction pathways, which regulate the expression of TF genes such as WRKYs,MYBs,and NACs via a signal cascade.Then the TFs induce the expression of CSLs,resulting in an increase in hemicellulose content in SL7 compared with L22. Because these TFs are also involved in the response of plants to environmental stimuli,and GO analysis suggested that some DEGs were related to responses to stimuli, we propose that the improvement of fiber quality in SL7 is the result of the combination of genetic and environmental factors. This study suggests a molecular mechanism of cotton fiber development involving plant hormone signals, TFs, and transcriptional regulation of enzymes that participate in fiber synthesis.

Declaration of competing interest

Authors declare that there are no conflicts of interest.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (31601345 and 31671742),the Earmarked Fund for China Agriculture Research System(CARS-15-05),the Taishan Scholar Project of Shandong Province (ts201511070),and the Innovation Project in Shandong Academy of Agricultural Sciences(CXGC2016A01).

Fig.8- Model of the mechanism by which exogenous chromosome fragments in SL7 promote fiber development.The red segment represents exogenous chromosome segments introgressed from Gossypium barbadense in SL7.Arrows indicate regulatory relationships and dotted lines represent possible regulation.

Appendix A. Supplementary data

Supplementary data for this article can be found online at https://doi.org/10.1016/j.cj.2019.11.003.