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). It was found that the fuzz density in high-temperature environment (HT, 28~35℃) was significantly 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
To reveal the molecular regulation mechanism of ambient temperature on the fuzz fiber initiation and 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 FPKM (Fragments Per Kilobase of transcript per Million fragments mapped) to measure the transcript or gene expression levels. A total of 63,113 transcripts were obtained with FPKM values > 0 in at least one sample, of which 44,864 were expressed with FPKM values ≥ 1 (Additional file 6: Table S3).
Table 1 Summary of sequencing data for different fuzz fiber developmental stages in high and low temperatures.
Sample
|
Raw Reads (×107)
|
Raw Bases (×109)
|
Clean Reads (×107)
|
Clean Bases (×109)
|
GC (%)
|
Q30 (%)
|
Mapped Reads (%)
|
H1-1
|
4.26
|
6.38
|
4.18
|
6.27
|
45.94
|
93.71
|
92.62
|
H1-2
|
4.15
|
6.22
|
4.07
|
6.10
|
45.63
|
93.83
|
94.86
|
H4-1
|
4.23
|
6.34
|
4.14
|
6.22
|
46.19
|
93.46
|
94.49
|
H4-2
|
4.19
|
6.28
|
4.11
|
6.16
|
45.81
|
93.61
|
95.17
|
H7-1
|
4.11
|
6.17
|
4.03
|
6.05
|
46.13
|
93.32
|
94.54
|
H7-2
|
4.25
|
6.37
|
4.17
|
6.25
|
45.80
|
93.56
|
94.89
|
L1-1
|
4.10
|
6.16
|
4.02
|
6.03
|
45.25
|
93.50
|
91.86
|
L1-2
|
4.10
|
6.15
|
4.01
|
6.01
|
45.54
|
93.62
|
90.62
|
L4-1
|
4.14
|
6.21
|
4.05
|
6.08
|
46.98
|
93.30
|
92.65
|
L4-2
|
4.17
|
6.25
|
4.09
|
6.13
|
45.76
|
93.45
|
94.36
|
L7-1
|
4.17
|
6.25
|
4.08
|
6.13
|
45.92
|
93.33
|
95.05
|
L7-2
|
4.08
|
6.11
|
4.00
|
6.00
|
45.84
|
93.56
|
92.36
|
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 log2 ratio ≥ 1.0 and p-value < 0.01. In HT, 14,490 DEGs were found at 1 DPA vs 4 DPA (3393 up-regulated, 11,097 down-regulated), 2160 DEGs at 4 DPA vs 7 DPA (1230 up-regulated, 930 down-regulated), and 12,716 DEGs at 1 DPA vs 7 DPA (3628 up-regulated, 9088 down-regulated) (Fig. 4c-d). In LT, 5694 DEGs were found at 1 DPA vs 4 DPA (1473 up-regulated, 4221 down-regulated), 4290 DEGs at 4 DPA vs 7 DPA (1667 up-regulated, 2623 down-regulated), and 17,748 DEGs at 1 DPA vs 7 DPA (3352 up-regulated, 14,396 down-regulated) (Fig. 4c-d). Compared to LT, the HT samples had 2480 DEGs at 1 DPA (951 up-regulated, 1529 down-regulated), 7355 DEGs at 4 DPA (3072 up-regulated, 4283 down-regulated), and 4164 DEGs at 7 DPA (2876 up-regulated, 1288 down-regulated) (Fig. 4c-d). This indicated that the fuzz fiber initiation development is more susceptible to temperature stress at 4 DPA. Of these DEGs, the genes differentially expressed between LT and HT treatments at the same stages were identified as temperature responsive genes, and those identified at different stages under the same temperature environment were defined as possibly related to fiber development. Therefore, we identified 9667 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), KING1 (GB_A12G2167), STZ (GB_D05G2149), BCAT-2 (GB_A05G2168), ACX4 (GB_A13G1045), and 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-4-phosphate 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).
At the same time, we performed KEGG pathway enrichment analysis on the DEGs in the three clusters (Additional file 11: Table S5). DEGs in cluster 5, 6 and 7 were significantly (p < 0.05) mapped into 11, 7 and 14 KEGG pathways, respectively. This indicated that there might be more pathways involved in fuzz fiber initiation compared with lint fiber development in high temperature. The most significant pathways in cluster 5 were "Phagosome (14 DEGs)", "Fatty acid elongation (5 DEGs)", "Metabolic pathways (59 DEGs)" and "Amino sugar and nucleotide sugar metabolism (9 DEGs)" (Fig. 5c). The most significant pathways in cluster 6 were "Circadian rhythm (7 DEGs)", "Plant hormone signal transduction (18 DEGs)" and "Brassinosteroid biosynthesis (3 DEGs)" (Fig. 5c). And the top three pathways with the most DEGs in cluster 7 were "Biosynthesis of secondary metabolites (29 DEGs)", "Phenylpropanoid biosynthesis (8 DEGs)" and "Metabolic pathways (32 DEGs)" (Fig. 5c).
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-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 genes involved in fuzz fiber initiation. Different genes of four families, such as NAC, bHLH, MYB and C2H2, were found to be participated in lint fiber early development and fuzz fiber initiation in high temperature, respectively. However, transcription factors of four families, such as AP2/ERF(ERF), C2C2-LSD, HD-ZIP and PLATZ, were mainly participated in lint and fuzz fiber initiation. In addition, six genes distributed in five families (RAV1, BLH1, HB40, LBD1 and OFP11) were exclusively enriched in cluster 7. According to the PPI network, STZ related to Cys2/His2-type zinc-finger proteins was identified as the key transcription factors of cluster 7, which might play a pivotal role in fuzz fiber initiation induced by high-temperature.
DEGs involved in plant hormone signal transduction pathway
Previous studies have found that cotton fiber development is also regulated by some genes related to plant hormone signal transduction. We identified 189 genes related to hormone signal transduction from 9667 DEGs involved in both temperature response and fiber development. Eighteen unigenes were from cluster 6, including 7 auxin (AUX1, IAA and SAUR), 5 abscisic acid (PYR/PYL, PP2C, SnRK2 and ABF), 3 cytokinine (CRE1, AHP and B-ARR), 1 ethylene (CTR1) and 2 brassinosteroid (CYCD3) signal transduction component-encoding genes (Additional file 13: 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). The results showed that the H2O2 content at 4 DPA was significantly higher than 1 or 7 DPA under HT environment, and the content at 4 DPA under LT environment was significantly higher than that of HT environment. The content of OFR in the early development of fibers under HT environment was generally higher than that in LT environment, indicating that the two ROS, OFR and H2O2, might play different roles in fuzz fibers initiation development. The activity of the four enzymes showed a downward trend from 1 to 7 DPA under different temperature environments. Among them, the enzyme activities of SOD, CAT and APX at 4 DPA in HT environment were generally higher than those in LT environment. This indicates that the regulation level of ROS homeostasis is an important cause of the difference during fuzz fiber initiation in different temperatures.
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 (R2=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.