Transcriptome analysis reveals gene expression associated with fuzz fiber initiation regulated by high-temperature in Gossypium barbadense

Background: Gossypium barbadense L. is the most important renewable source of textile fiber. Cotton fiber cell initiation and elongation are often affected by various environmental stimulus, such as high temperature. However, little is known about the underlying mechanisms of temperature regulating the fuzz fiber initiation. Results: In the present study, phenotypic observation revealed that high temperature (HT) accelerated the fiber development, improved fiber quality and induced fuzz fiber initiation. It has been proved that the fuzz fiber initiation was inhibited by low temperature (LT), and 4 days post-anthesis (DPA) was the key stage for fuzz fiber initiation. Based on comparative transcriptome analysis, a total of 43,826 differentially expressed genes (DEGs) were identified, of which 9,667 were involved in both fiber development and temperature response with 901 transcription factor genes and 189 genes related to plant hormone signal transduction. Further analysis of gene expression patterns revealed that 240 genes were involved in fuzz fiber initiation. Functional annotation revealed that the candidate genes related to fuzz initiation were significantly involved in asparagine biosynthetic process, cell wall biosynthesis and stress response. Furthermore, the expression trends of sixteen selected genes from the RNA-seq data were almost consistent with the results of qRT-PCR results. Conclusions: Our study revealed several potential candidate genes and pathways that related to fuzz fiber initiation induced by high-temperature and provided a new view of temperature-induced tissue and organ development in Gossypium barbadense. In this study, we explored the effect of temperature on G. barbadense fuzz density and believed that n2 may also be temperature-regulated.

(AD 2 ), whose annual yield is more than 90% of the cotton fiber production. G. hirsutum is widely cultivated all over the world due to its high yield, while G. barbadense is not high in yield but prized because of its excellent fiber quality [1]. Mature cotton seeds are covered with two different types of fibers, lint and fuzz, both of which are mainly composed of cellulose. The development of cotton fibers can be divided into four overlapping stages: fiber initiation, elongation, thickening of secondary cell walls, and dehydration [2]. Lint fiber initiation occurs from − 3 DPA to 3 DPA [3]. With the help of scanning electron microscope (SEM) and cotton fiber mutants, the fuzz fiber initiation stage was determined to be 4 DPA [4,5]. The first three stages of fiber development are hypersensitive to environmental stress [6]. Extreme temperature and other abiotic stresses during fiber development can significantly reduce cotton fiber yield and quality, and may even lead to falling bolls and squares [7][8][9].
In plants, studies of trichomes and root hairs in Arabidopsis thaliana has greatly expanded our understanding of cotton fiber initiation and elongation. Numerous genes have been found to be involved in trichome initiation in this model plant. In Arabidopsis, it has been widely believed that the MBW transcriptional complex (GL3/EGL3-GL1-TTG1) determines the fate of trichome cells by inducing the expression of GL2 [10]. In addition, GIS2, ZFP8, ZFP6 and ZFP5 were found to promote the expression of MBW complex members in Arabidopsis [11][12][13]. There are many similarities between cotton fibers and Arabidopsis trichomes in the regulatory mechanism of initiation and elongation development, which has been proven by some cotton homologs. GhMYB25-like and GhMYB25 are two R2R3 MYB genes that regulate fiber initiation and elongation. GhMYB25-like plays a major role in regulating fiber initiation and development. Silencing GhMYB25-like will inhibit cotton lint and fuzz fiber initiation development [14]. Overexpression of GhMYB25 in cotton increases the number of trichomes on leaves and the fibers number on ovule epidermis at anthesis [15]. As another R2R3 MYB transcription factor, GhMYB109 also specifically expressed during lint initiation and elongation development [16]. As a homeodomain-leucine zipper (HD-ZIP) transcription factor, GhHOX3 is revealed to play a major role in regulating fiber elongation. Overexpression of GhHOX3 increases fiber length, while silencing the gene reduces fiber length [17]. Overexpression of another HD-ZIP transcription factor, GhHD1, increased the number of initiating fibres, but had not other adverse effects on leaf hairs [18]. Cotton PROTODERMAL FACTOR1 gene (GbPDF1) is mainly expressed during fiber initiation and early elongation stages, and silencing of GbPDF1 will retard fiber initiation [19].
Cotton fiber initiation and elongation are not only dependent on their genetic specificity, but also regulated by environmental factors and endogenous hormones. Auxin, gibberellin and brassinolide have long been known to play important roles in plant cell expansion and elongation. As another phytohormone, ethylene is widely studied in fruit ripening, dormancy release, and stress response. In addition, ethylene also has been found to play an important role in regulating root hair and hypocotyl development [20,21]. Auxin, gibberellin, brassinolide and ethylene have positive effects on the fiber development, while abscisic acid and cytokinin inhibit the fiber cell development [22][23][24]. Jasmonic acid increases the number of leaf trichomes in Arabidopsis by promoting the expression of GL3 [25,26]. In cotton, JA-related metabolism promotes fiber initiation [27]. In addition to phytohormones, Ga 2+ and H 2 O 2 also regulate the fiber cells development [28].
Temperature is a major environmental factor affecting plant growth and development. Plants of different species require different growth temperatures, and different developmental stages or developing organs of the same plant also require favorable temperatures. High temperature stress is an abiotic stress common to all crops from subtropical origin. High temperature stress during flowering period can cause male sterility in some crops, such as rice and cotton [29,30]. Some molecular mechanisms of plant response to high temperature have been revealed, including the response of transcription factors and heat shock proteins, calcium ion and reactive oxygen species signal transduction, and response of endogenous plant hormones such as auxin and gibberellin [31][32][33]. These studies provide a model for temperature-regulated development of specific tissues or organs, such as the development of cotton fibers in high temperatures.
Most G. barbadense varieties have a low density of fuzz fibers, which is different from most G. hirsutum varieties whose seeds surface are completely covered by fuzz fibers. It is often observed that the fuzz density of the same G. barbadense cultivar varies with different environments, and the fuzz density in Xinjiang and Hainan Province is lower than that in Anyang City [34]. In this study, 121 G. barbadense cultivars were used for fuzz density evaluation. The temperature of full-bloom stage in three main growing sites in China was investigated. Accession L7009 was considered sensitive to environmental stimulus. We also verified the influence of temperature on fuzz fiber initiation of L7009 and when the fuzz fiber initiation occurred, and examined the effect of different temperatures on fiber quality. In order to understand the transcriptional regulation mechanism of temperature on the early development of lint fibers and fuzz fiber initiation, we conducted comparative transcriptome on ovules of L7009 in different temperatures. Then we performed functional enrichment analysis and coexpression analysis of differentially expressed genes from different profiles. These findings in this study extend our understanding of how temperature affects the development of cotton fibers, and also provide valuable data for studying the mechanism of temperature-regulated fuzz fiber development.

Temperature regulates fuzz fiber initiation and affects lint fiber development
To reveal the environmental influence on fuzz and lint fibers development, we investigated fuzz density (FD) and lint fiber quality of 121 G. barbadense cultivars in three cities (Korla, Anyang and Sanya). For convenience, we visually graded the fuzz density to eight levels (Additional file 1: Figure   S1). In Sanya, 23.97% of the cultivars showed low fuzz density (FD=0.5) and 20.66% showed high fuzz density (FD=3 or 3.5). However, in Korla and Anyang, 15.70% and 2.48% of the cultivars showed low fuzz density, while 29.75% and 60.34% showed high fuzz density, respectively (Fig. 1a). This suggests that the environment of Sanya is more likely to reduce fuzz density, and the environment of Anyang is more likely to increase the fuzz density. We also investigated the relationship between lint quality and fuzz density of 98 cultivars planted in Sanya. Three of the five fiber quality indexes showed an significant positive correlation between fuzz density and lint fiber quality (FL, r=0.59; FS, r=0.61; FE, r=0.47) (Additional file 2: Figure S2). To test whether temperature affects the fuzz density, we investigated the temperature of three environments during the blooming period of three years. The days with minimum temperature above 24℃ accounted for 64.6% of the total survey days (n=96) in Anyang, while 71.9% of the days (n=93) were below 22℃ in Sanya (Fig. 1b). Therefore, we infered that temperature may be the key environmental factor affecting fuzz density of G.
barbadense. We assumed that high temperature may increase the fuzz density, while low temperature may have the opposite effect on it.
To test our hypothesis, we simulated the ambient temperature of the natural environments (Fig. 1c).
higher than that in low-temperature environment (LT, 20~25℃) (Fig. 1d). This indicates that temperature is indeed the decisive environmental factor regulating fuzz initiation, and that high temperature increases the fuzz density. Based on the fact that fuzz density of L7009 was significantly affected by temperature, we designed an experiment to confirm the time when the fuzz initiation occurs. First, by transferring cotton plants from HT to LT (HTL), we found that the fuzz density at 4 DPA (or DT) was significantly lower than that at 5 DT or higher han that at 2 DT, while there was no significant difference when compared to 3 DT, indicating that 3~4 DPA was an important period for fuzz fiber initiation development (Additional file 3: Table S1; Fig. 2a-b). To make the results more credible, we transported another batch of cotton plants from LT to HT (LTH), and found that the fuzz density at 4 DT was significantly lower than that at 3 DT or higher han that at 6 DT, while there was no significant difference when compared to 5 DT, indicating that 4~5 DPA was also important for fuzz fiber initiation development (Additional file 3: Table S1; Fig. 2c-d). By comparing the fuzz density of HTL and LTH at the same stages, it was found that there was no significant difference between HTL and LTH at 4 DPA, while significant differences between any other stages were found, indicating that 4 DPA was the key stage for fuzz initiation, which was consistent with previous findings in G. hirsutum (Additional file 3: Table S1) [5]. For further study, the early and later stages of fuzz fiber differentiation were assigned to 1 and 7 DPA, respectively. With phenotypic observation, we found that high temperature could accelerate fibers initiation and development and increase the fuzz density on seed surface ( Fig. 3a-b). In addition, it was found that SFC_n, FS, HGW, FP, SS, Length_w, UQL, MR, FNC, SCNC, FL, SFC_w, Fineness, SCNMS and Micronaire in HT were significantly different from those in LT (p<0.05) (Additional file 4: Figure S3; Fig. 3c). Among them, six fiber traits in HT were significantly better than those under LT, and the differences in Length_n and IFC were not significant.
Transcriptome sequencing and comparative analysis of differentially expressed genes between different stages and environments development of L7009, we constructed cDNA libraries from 12 samples of three stages under two environments. The cDNA libraries were then sequenced using an Illumina HiSeq 4000 sequencing platform based on paired-end sequencing. All RNA-seq raw datasets were deposited in the NCBI database with a SRA accession number SUB6753336. We totally obtained 73.43 Gb clean data after mRNA sequencing for 12 samples with at least 6.00 Gb clean data for each sample. In each sample, more than 93.30% of bases score Q30 and above ( Table 1). The clean data were then mapped to the reference genome, with the mapping ratio varying from 90.62% to 95.17%. Based on the mapped results, alternative splicing prediction, genetic structural optimization and novel genes discovery were analyzed. Finally, 5,994 novel genes were identified, 5,180 of them are functionally annotated (Additional file 5: Table S2). Based on the alignment results, gene expression analysis was performed.
We utilized StringTie to estimate its expression level based on a maximum flow algorithm, and used  Table S3). The Pearson correlation coefficient between two biologically repeated samples of the same developmental stage in different environments fluctuated between 0.98 and 0.99, and the sample clustering also showed a good correlation between the two biological replicates (Fig. 4a). The results of principal component analysis (PCA) of the 12 samples were also consistent with the above analysis results (Fig. 4b).
Intra-environment and inter-environment comparisons of all expressed genes were performed to reveal the dynamic transcriptional changes during fuzz fiber initiation. In total, 26048 DEGs were identified under the thresholds of absolute log 2 ratio ≥ DEGs that were likely to be involved in both temperature response and fiber development (Additional file 7: Figure S4).

Gene expression trend analysis and PPI network construction
The 9667 DEGs were subjected to hierarchical clustering using Mfuzz package of R sortware and divided into 9 clusters (Additional file 8: Figure S5), in which 4756 (49.2%) genes with membership ≥ 0.7 were illustrated with a hierarchical clustering heatmap (Fig. 5a). The expression levels of 512 genes in cluster 5 were continuously up-regulated from H1 to H7 and L4 to L7 (Fig. 5b). Genes in clusters 6 (823) and 7 (240) were predominantly expressed at L1 and H4, respectively. According to our previous verification, fuzz fiber initiation occured at 4 DPA under HT (H4), but didn't occur in low temperature condition. Combined with phenotypic observation, we speculated that the genes in cluster 7 were mainly involved in fuzz fiber initiation induced by high temperature. In contrast, genes in cluster 5 and 6 were involved in lint fiber elongation and initiation under temperature regulation, respectively. However, most of the genes were involved in lint fiber initiation development, and only a small number of genes were related to the fuzz fiber initiation, which indicated that the regulation of lint fiber development is relatively complicated and preferentially.
The 240 genes with membership ≥ 0.7 in cluster 7 were mapped to the STRING database (https://string-db.org/) and the PPI network was visualized using the Cytoscape software to investigate the possible interaction between the differentially expressed genes (Fig. 6a). In the PPI network, the node genes with larger degree are defined as hub genes. In the cluster 7, P5CS1 (GB_D01G2569), HSPRO2 (GB_A10G2743) were defined as hub genes, which were considered as candidate genes involved in regulating fuzz initiation in high temperature.

Gene ontology and KEGG pathway analysis of DEGs
There is a partial overlap in the developmental stages of lint fiber elongation and fuzz fiber initiation.
Not only are their underlying gene expression patterns different, but the functions of related genes may also differ significantly. To further understand the functions of the DEGs in the three clusters, GO enrichment for each cluster was performed (Additional file 9: Table S4). To reduce the functional redundancy among GO-terms, the presence of overrepresented GO-terms were obtained using REVIGO program (http://revigo.irb.hr/) and visualized using treemap package of R software [35].
The DEGs in cluster 5 were significantly (FDR < 0.05) enriched into 30 GO-terms of biological process, of which the terms involved "fatty acid biosynthetic process", "monocarboxylic acid biosynthetic process", "lipid biosynthetic process", "microtubule-based process" and "intracellular protein transport" were the most significant (Additional file 9: Table S4). Forty significant terms for molecular functions were enriched, the most significant of which were redundant and involved "structural constituent of cytoskeleton", "GTPase binding", "UDP-glycosyltransferase activity" and "actin binding" (Additional file 9: Table S4). For cell component 16 terms were also enriched and the most significant involved "microtubule", "cytoskeleton" and "membrane" (Additional file 9: Table S4). The DEGs in cluster 6 were significantly enriched into 104 GO-terms of biological process, of which the supercluster terms involved "glycerolipid metabolism", "dephosphorylation" and "protein deacetylation" (Fig. 7). Fifty-one significant terms for molecular functions were enriched, the superclusters of which involved "phosphoric ester hydrolase activity" and "1-phosphatidylinositol-4phosphate 5-kinase activity" (Additional file 10: Figure S6). For cell component 10 terms were significantly enriched and "cytosol" was the main supercluster (Additional file 10: Figure S6). The DEGs in cluster 7 were significantly (p < 0.05) enriched into 63 GO-terms of biological process, of which the supercluster terms involved "asparagine biosynthesis", "response to stress" and "hemicellulose metabolism" (Fig. 7). Forty-one significant terms for molecular functions were enriched, the superclusters of which involved "xyloglucan:xyloglucosyl transferase activity" and "peroxidase activity" (Additional file 10: Figure S6). Nine significant terms for cell component were enriched, and the supercluster terms of treemap involved "extracellular region" and "external encapsulating structure" (Additional file 10: Figure S6). Genes in the PPI network of cluster 7 were mainly enriched in the GO-terms supercluster "response to stress" and "response to stimulus" (Fig.   6b).

Differentially expressed transcription factors involved in fuzz fiber initiation
Numerous studies have shown that the development of Arabidopsis trichome or cotton fiber was strictly regulated by several transcription factors [36][37][38][39][40][41][42]. To identify the differentially expressed TFs associated with fuzz fiber initiation induced by high-temperature, twenty-five differentially expressed TF genes assigned to thirteen families were detected in cluster 7.
In this study, totally 20 and 105 TF genes were identified in cluster 5 and 6, respectively (Additional file 12: Table S6). Of these, there were more genes involved in lint fiber initiation, and relatively few  Table S7). Five unigens were were from cluster 7, and these genes encode the components of auxin (GH3) and brassinosteroid (TCH4) signal transduction, respectively (Additional file 13: Table S7). There were more genes in cluster 6 than that in cluster 7. For example, no genes related to abscisic acid, cytokinine, and ethylene signal transduction were found in cluster 7, but several genes were identified in cluster 6. Although genes for the auxin and brassinosteroid signaling pathways were found in both c6 and c7, the components encoded by these genes were different, indicating that GH3 and TCH4 may play a key role in regulating the fuzz fiber initiation induced by high-temperature.
DEGs related to cotton fiber development According to previous reports, it is believed that the homologous genes related to trichome initiation of Arabidopsis in cotton may also be involved in the fiber initiation. We collected more than 100 genes associated with trichome initiation of Arabidopsis and identified 695 homologous genes in the G. barbadense genome. These 695 genes were compared with 9667 environment-related DEGs, and 56 genes were shared in the two gene sets. Based on the expression profile, we found that some of these genes were affected by temperature (Fig. 8). Among these temperature-responsive genes, the expression of the gene homologious to ARP3 (GB_D01G1926 and GB_A01G1791), NAC029 (GB_A12G1914), and MYB82 (GB_A08G2828) in L7009 showed an expression trend positively correlated with fuzz fiber initiation development. Hence these genes were identified as candidate genes that might be involved in fuzz fiber initiation induced by high-temperature.

Changes in antioxidant enzyme activity and ROS content during fuzz initiation
According to the GO analysis, some genes, such as RCI3 (GB_A03G1340), DOX2 (GB_A09G2023), PRX52 (GB_A03G0328), GPX4 (GB_D08G0779) and ERD10 (GB_A01G1053), potentially related to fuzz initiation were involved in stress response. Therefore, we measured the activities of four antioxidant enzymes and the contents of two ROS (Fig. 9

Verification of RNA sequencing data by qRT-PCR
To validate the RNA-seq data, quantitative real-time PCR (qRT-PCR) was performed for sixteen genes randomly selected from different expression profiles. The sixteen genes and their primer sequences were listed (Additional file 14: Table S8). In general, the relative expression levels based on qRT-PCR were consistent with the results of RNA-seq and the correlation between them was significant (R 2 =0.78) (Additional file 15: Figure S7; Fig. 10b). However, when we clustered the sixteen selected genes according to the RNA-seq data and qRT-PCR results, we found that GB_A01G1791 was not clustered with the other three genes of cluster 7, but with the genes of cluster 2 and 9 (Fig. 10a). We then calculated the correlation between the nine clusters and found that the Pearson correlation of cluster 2 and cluster 9 reached 0.94, indicating that the gene expression patterns from these two clusters are similar (Additional file 16: Table S9). These results indicated that the RNA-seq data was basically credible and accurate.

Discussion
Temperature is the key external factor that determines the fate of epidermal cells of G. barbadense ovules and the quality of lint fiber As an important fiber source plant, G. barbadense has been favored by scientific researchers and production industries for its excellent fiber quality traits. The quality of cotton fibers is sensitively affected by environmental stimulus [43,44]. The difference in fiber quality mainly comes from meteorological factors (e.g. temperature and rainfall) and soil factors [45]. Of these environmental stimulus, temperature is the most important factor affecting fiber quality. Low temperature stress, especially the mean daily minimum temperature (< 20 ℃), is the major abiotic stress limiting the formation of cotton fiber quality such as fiber length [46]. In this study, it led to significant differences in fiber quality when the ambient temperature was set to 20~25 ℃ and 28~36 ℃, respectively.
Interestingly, high temperature significantly increased the fuzz density (FD) when compared to a low temperature. This indicates that the change of FD is regulated by the ambient temperature.
Moreover, as the FD increased, FL, FS, and FE also increased accordingly. Since FD had a similar change trend with the quality traits of lint fibers, we speculated that the lint fiber quality must be regulated by the ambient temperature, that is, higher temperature will increase FL, FS and FE.
However, the underlying mechanism remains unknown and the relevancy between them is noticeable. It is reported that the elevated temperature can indeed accelerate fiber elongation process. Nevertheless, if the temperature is too high, it will shorten the fiber rapid elongation duration However, the n2 gene has still not been discovered, and the environmental stimulus affecting the function of the n2 gene remain unclear. In this study, we explored the effect of temperature on G.
barbadense fuzz density and believed that n2 may also be temperature-regulated.  [56]. AtMYC1, another bHLH gene, functions upstream of GL2 by interacting with MYB and TTG1 [57]. GhFP1 (a bHLH gene) as a positive regulator participates in controlling fiber elongation by activating BR biosynthesis and signaling transduction [58]. GhbHLH18 is negatively associated with fiber strength and length by activating peroxidasemediated lignin metabolism [59]. In summary, members of the bHLH family not only regulate trichome or fiber initiation, but also control lint fiber elongation and quality. Our results also indicated that bHLH family members were widely involved in the development of lint and fuzz fibers, which further indicates that this TF family members can not only positively regulate lint but also positively regulate fuzz development.

Major transcription factors involved in fiber development
In addition to bHLH and NAC, members of MYB and C2H2 also play an important role in regulating cotton fiber development. As mentioned earlier, the MYB genes play an important role in regulating the initiation and elongation of cotton fibers, such as GaMYB2, GhMYB109, GhMYB25, and GhMYB25like [14,15]. In this study, two MYB genes homologous with MYB116 and MYB305 were expressed in fuzz fiber initiation stages, indicating that they play a role in regulating fuzz development. C2H2 genes have also been reported to play an important role in determining the epidermal cell fate. As as a steroid 5α-reductase, is considered to be the rate-limiting enzyme in BR biosynthesis. The loss of this enzyme activity will lead to a reduction in BRs levels [61]. The GhDET2 is positively regulated during the process of fiber initiation and extreme elongation, and inhibition of this gene will restrain fiber initiation and elongation [62]. GhPAG1, as a homologous gene of CYP734A1, controls cotton fibers development by regulating the endogenous BR level through ethylene signal-mediated VLCFA synthesis [63]. GhBZR1, as a member of the BR signal transduction cascad, can be regulated by Gh14-3-3, and bind to the promoter of GhXTH1 to regulate fiber initiation and elongation development [64]. As a xyloglucan endotransglycosidase gene, the expression of TCH4 is induced by auxin, BR, and environmental stimulus and is mainly expressed in trichomes, lateral root primordia, and elongating hypocotyls [65]. In this study, we identified three homologous genes of TCH4, which were mainly present in cluster 7 expression pattern. These results indicate that temperature stimulation may control the fuzz initiation development by regulating the expression of genes in BR signal transduction pathways.
In addition, two GH3 genes related to auxin signal transduction pathway were also identified during fuzz fiber initiation (cluster 7). Studies have found that auxin plays an important role in the initiation and elongation development of cotton fibers [28]. GhARF2 and GhARF18 are two homologous genes of ARF in the auxin signal transduction pathway, which are mainly expressed during the fiber initiation process [66]. Overexpression of these two genes in Arabidopsis promotes the initiation of trichomes.
The homologous genes of IAA, GhAux8 and GhIAA16, are mainly expressed during the fiber development stage [67]. Of these, GhAux8 is preferentially expressed in elongating fibers, and GhIAA16 is expressed in fiber initiation stage. This shows that the IAA signal may play an important role in lint and fuzz fiber initiation stages. This suggests that these two hormones play different roles in the development of the two type fibers, and that the hormones involved in the lint fiber initiation may also be involved in the fuzz fiber initiation.

Maintenance of ROS homeostasis is sufficient to temperature-regulated fiber initiation
Reactive oxygen species (ROS) mainly include superoxide radical, hydrogen peroxide and hydroxyl radical, which can regulate plant cell proliferation and expansion, cell differentiation (such as root hair and lateral root), dedifferentiation and tissue regeneration [68]. Ascorbate peroxidase (APX), a hydrogen peroxide scavenging enzyme, is involved in regulating ROS homeostasis [69]. The expression of GhAPX1 at 5 DPA in the fl mutant ovule was much lower than that of the wild type, and application of H 2 O 2 induced GhAPX1 expression to promote fiber elongation. Enzyme activity and H 2 O 2 content assays showed that APX members except GhAPX1 were mainly involved in the late development of fibers (30 DPA) [70]. GhPOX1 encodes a class III peroxidase that functions during fiber elongation by mediating ROS production [71]. Studies have shown that homeostasis the H 2 O 2 and redox level is the core mechanism leading to inhibition of fiber initiation or elongation [72].
Increased hydrogen peroxide content was detected in 5 DPA ovules of cotton, and H 2 O 2 could promote fuzz fiber initiation in vivo [73]. In the present study, we found that a portion of the fuzz fiber initiation-related candidate genes were significantly enriched in the oxidative stress response pathway, in which two hydrogen peroxide scavenging enzyme genes, GbPrx52 and GbRCI3, were found. In Arabidopsis, AtPrx52 is involved in lignin synthesis as a hydrogen peroxide scavenging enzyme, while AtRCI3 is primarily involved in tolerance regulation under various abiotic stresses [74,75]. In this study, the content of H 2 O 2 and antioxidant enzymes activity at 4 DPA in the HT environment were lower than that in LT environment, indicating that the ROS level in the HT environment was maintained at a suitable state for the fuzz fiber initiation development, and the ROS homeostasis of LT environment may be disrupted to a certain extent, which may lead to the failure of the fuzz fiber initiation development. In summary, it is believed that the ROS homeostasis regulation must play a key role in fuzz fibers initiation stage.

Conclusions
This study emphasized that thermal stimulus has a significant effect on lint and fuzz fiber development, and validated the regulatory model that temperature determines the fuzz fiber initiation of G. barbadense. Low temperature could inhibite the fuzz fiber initiation and delay the lint fiber development, while high temperature could relieve the negative effects of low temperature on these phenotypes. Through comparative transcriptome analysis, a large number of DEGs (9667) involved in both fiber development and temperature response were identified, and three expression patterns closely related to the different development of the fibers were obtained, of which the genes in SC6 were differentially expressed between the two temperature environments during the fuzz fiber initiation. Functional enrichment analysis indicated that DEGs involved in asparagine biosynthetic process, cell wall synthesis and stress response were related to the fuzz fiber initiation induced by high-temperature. Transcription factors, plant hormone signals and ROS homeostasis maintenance may also be potentially involved in fuzz fiber initiation. These findings can better help us understand the morphology and molecular mechanism of fuzz fiber initiation regulated by high temperature.

Plant materials and temperature treatments
A total of 121 G. barbadense cultivars (provided by the Cotton Research Institute of the Chinese Academy of Agricultural Sciences) were used to investigate the effects of different environment on fuzz density and lint fiber quality. A cultivar named L7009 (provided by Zhejiang Agriculture and Forestry University) was used to evaluate the effects of temperature on fuzz density and fiber quality.
We simulated the high temperature environments with 25°C in the day and 20°C at night, and the low temperature environments with 36°C in the day and 28°C at night. Both environments had the same photoperiod (15 h light and 9 h dark) and light intensity. First, we transplanted the three-leaf seedlings into flowerpots with one plant per pot and then transferred them to the greenhouse. When it came to full-bloom stage, 120 cotton plants were equally divided into two groups and transferred into high and low temperature environment, respectively. Dated cards were hanged on the flowers, and the flowers in the HT environment were artificially pollinated everyday. When transferred into the cotton growing rooms on the 15th day, ovules of 0 to 7 DPA in the two environments were separately collected and stored in liquid nitrogen. Each biological replicate at each stage came from a mixed sample of three bolls. The seeds were harvested when the cotton bolls were split and dried. The seeds and lint fibers were manually separated, the fuzz density was investigated, and the fiber quality was measured. To reveal the time when fuzz fiber initiation occurred, we planted another batch plants of the L7009 cultivar, and the management was the same as the previous batch. To accurately track the temperature changes sensed by each seed during fiber development, we numbered sixty plants in each environment and dated cards were hanged on the flowers every day. After ten days cotton plants from two environments were exchanged. The flowers that bloomed on the day of changing the growing environment were recorded as 0 days after transfer (DT). The flower that bloomed on the day before changing the growing environment was recorded as 1 DT, and flower that bloomed on the first day after changing the growth environment was recorded as -1 DT. Flowers that bloomed every day were tagged with dated cards.

Fuzz phenotype determination and fibre quality measurement
The fuzz density of seeds was classified into eight grades by visual inspection. The seeds with a score of "0" were fuzzless, and the seeds with a score of "3.5" were covered with dense fuzz fibers. The fuzz density score was used to quantitatively characterize the fuzz phenotype of the seeds, and it was measured in the same way as previously reported [76]. After the fiber and seed were separated     Validation the time of fuzz fiber initiation. After being moved from HT to LT, the dynamic variation of fuzz density with time (a), and the significant difference of fuzz density between different treatment depths (b). After being moved from LT to HT, the dynamic variation of fuzz density with time (c), and the significant difference (p-value) matrix of fuzz density between different treatment depths (d).     Treemap of BP GO terms found in the three clusters. Each rectangle represent a significant GO-term.
Related GO terms are clustered together in superclusters of the same color and supercluster names are labelled in white. The sizes of rectangles are adjusted to reflect the relative corrected p-value.