Deep Transcriptome Analysis Reveals Reactive Oxygen Species (ROS) Network Evolution, Response to Abiotic Stress, and Regulation of Fiber Development in Cotton

Reactive oxygen species (ROS) are important molecules in the plant, which are involved in many biological processes, including fiber development and adaptation to abiotic stress in cotton. We carried out transcription analysis to determine the evolution of the ROS genes and analyzed their expression levels in various tissues of cotton plant under abiotic stress conditions. There were 515, 260, and 261 genes of ROS network that were identified in Gossypium hirsutum (AD1 genome), G. arboreum (A genome), and G. raimondii (D genome), respectively. The ROS network genes were found to be distributed in all the cotton chromosomes, but with a tendency of aggregating on either the lower or upper arms of the chromosomes. Moreover, all the cotton ROS network genes were grouped into 17 families as per the phylogenetic tress analysis. A total of 243 gene pairs were orthologous in G. arboreum and G. raimondii. There were 240 gene pairs that were orthologous in G. arboreum, G. raimondii, and G. hirsutum. The synonymous substitution value (Ks) peaks of orthologous gene pairs between the At subgenome and the A progenitor genome (G. arboreum), D subgenome and D progenitor genome (G. raimondii) were 0.004 and 0.015, respectively. The Ks peaks of ROS network orthologous gene pairs between the two progenitor genomes (A and D genomes) and two subgenomes (At and Dt subgenome) were 0.045. The majority of Ka/Ks value of orthologous gene pairs between the A, D genomes and two subgenomes of TM-1 were lower than 1.0. RNA seq. analysis and RT-qPCR validation, showed that, CSD1,2,3,5,6; FSD1,2; MSD1,2; APX3,11; FRO5.6; and RBOH6 played a major role in fiber development while CSD1, APX1, APX2, MDAR1, GPX4-6-7, FER2, RBOH6, RBOH11, and FRO5 were integral for enhancing salt stress in cotton. ROS network-mediated signal pathway enhances the mechanism of fiber development and regulation of abiotic stress in Gossypium. This study will enhance the understanding of ROS network and form the basic foundation in exploring the mechanism of ROS network-involving the fiber development and regulation of abiotic stress in cotton.


Chromosome Mapping of the ROS Network Gene Family
The ROS genes were mapped and found to be widely distributed across the entire cotton chromosomes. In tetraploid cotton, G. hirsutum, the ROS genes were mapped in all the 26 chromosomes, in which the highest gene loci density were observed in chromosome At05 and its homeolog chromosome with 46 (18.11%) and 38 (14.28%) ROS genes, respectively, while the lowest gene loci density in chromosome At04 with only 10 ROS genes accounting for 1.94% of all the ROS genes identified in G. hirsutum ( Figure 1A). The gene distribution among the diploid cotton species, G. arboreum and G. raimondii were found to be present in each of the 13 chromosomes, in which the highest gene loci were detected in chromosome A10 (35,12.36%) and D09 (40,13.89%) in G. arboreum and G. raimondii, respectively ( Figure 1B,C). Moreover, the genes mapping in the individual chromosomes showed the tendencies of the genes aggregating on either upper or the lower arms of the chromosomes, similar observation have been previously observed in the distribution of the LEA genes in cotton, with high gene loci densities being observed in chromosomes A10 and D09 [27]. The detection of the high number of the ROS genes within the chromosomes showed that the three sets of chromosomes harbor significant genes with diverse roles in enhancing abiotic stress tolerance in plants. This result suggests that chromosome D09, chromosome Dt05 and its homeolog chromosome At05 play a core role in the evolution and dynamic balance in the ROS genes balance in cotton.

Phylogenetic Tree, Gene Duplication, and Synteny Analysis of the ROS Network Gene Family
In evaluating the evolution of the ROS network genes, the entire ROS network genes were aligned and phylogenetic tree constructed in order to understand the pattern and nature of their evolution. All the genes were classified into seventeen (17) different clades and each clade was made up of members of a particular family. The 17 clades were superoxide dismutase (SOD), ascorbate peroxidase (APX), monodehydroascorbate reductase (MDAR), dehydroascorbate reductase (DHAR), glutathione reductase (GR), catalase (CAT), glutathione peroxidase (GPX), ferritin (FER), respiratory burst oxidase homologue (RBOH), NADPH-like oxidase (FRO), Ubiquinol oxidase (AOX), peroxiredoxin (PrxR), thioredoxin reductase (Trx), Glutaredoxin (GLR), and peroxidase (PER) (Figure 2A). A total of 683 ROS network genes were found to have undergone duplications, in which 631 genes through segmental duplication while only 52 genes were tandemly duplicated (Table S1). The results were in agreement to previous findings in which a number stress responsive genes evolution were majorly governed by segmental type of gene duplication, for instance the LEA genes [27]. Gene duplication is considered as the main processes which do lead to the expansion of various gene families in organisms, and it is a prevailing feature in plant genomes [28,29]. Generally, a high number of gene copies in a gene family are maintained through large-scale segmental duplication or small-scale tandem duplication during evolution.

Orthologous ROS Gene Analysis and Evolution Prediction through Synonymous (Ks) and Non-Synonymous (Ka) Rate of Substituion
Allopolyploid cotton evolved following trans-oceanic dispersal of an A genome diploid to the Americas, where the immigrant underwent hybridization, as female, with a native D genome diploid similar to modern G. raimondii [30]. The formation of ROS network might be earlier than that of Synteny analysis, showed a total of 161, 166, 99, and 105 orthologous genes synteny blocks between At and A, Dt and A, At and D, and Dt and D, respectively ( Figure 2B-D). It is interesting to note that more ROS genes were observed between the syntenic blocks among the three cotton species, for instance, at most five genes were found within the syntenic regions between At and A, and the same number of genes were also observed between Dt and A. However, at most 20 genes have been found within the syntenic regions between Dt and D ( Table 1 and Table S2). The high number of ROS genes detected within the same syntenic regions showed that this syntenic regions have been highly conserved, and this indicates their integral role they play within the plants. This result supports the principle of asymmetric selection of genes within the genomes [8].  2  19  10  19  13  3  3  11  4  8  4  2  7  2  3  5  1  7  1  4  6  0  3  0  3  7  0  1  0  2  8  0  2  0  2  9  0  0  0  2  10  0  1  0  0  11  0  0  0  2  15  0  0  0  0  16  0  0  0  1  20  0  0  0  1 NGEB: number of genes per synteny block; * for example, there are 136 synteny blocks which include one ROS gene between At subgenome and A genome.

Orthologous ROS Gene Analysis and Evolution Prediction through Synonymous (K s ) and Non-Synonymous (K a ) Rate of Substituion
Allopolyploid cotton evolved following trans-oceanic dispersal of an A genome diploid to the Americas, where the immigrant underwent hybridization, as female, with a native D genome diploid similar to modern G. raimondii [30]. The formation of ROS network might be earlier than that of Gossypium diverged from its closest relatives and Allopolyploid cottons formation. Phylogenetic analysis (Table S3) revealed that, a total of 241 paired genes, were orthologous among the three cotton species, G. raimondii, G. arboreum and G. hirsutum, found to be either directly or indirectly involved in the ROS network. In the analysis of the evolutionary history of the tetraploid cotton, 243 genes were orthologous between G. arboreum and G. raimondii. G. arboreum contained 14 genes, but not in G. raimondii, while only five specific genes were found to be orthologs between G. arboreum and At subgenome. In the D genome, G. raimondii had 17 genes, but not in G. arboreum and only five specific genes were found to be ortholog between G. raimondii and Dt subgenome ( Figure 3A). This result suggested that 243 genes were formed evolved much earlier before the separation of the A and D cotton genomes, which are the G. arboreum and G. raimondii, respectively. Thirty-one (31) genes were found to have evolved before the separations of the two diploid cotton species, and even the emergence of the allotetraploid cotton. three genes were found to have originated from D, and introgressed into the Dt subgenome, of the allotetraploid cotton, while some element of gene loss was evident, this could have been occasioned by either loose of function or chromosomal rearrangement during the evolution of allotetraploid from the two diploid cotton species, A and D cotton genomes ( Figure 3B). The synonymous substitution value (K s ) peaks of orthologous gene pairs between At/A and Dt/D was estimated at 0.004 and 0.015, respectively. The distribution of 241 putative ortholog pairs between the three cotton species/genomes (mode K s = 0.0040-0.015) clearly indicated that the secondary peak for cotton, lies to the left of the ortholog peak, which represents a burst of gene duplications that occurred in cotton after the separation between the two lineages. Similar observation has been made in the analysis of the genome evolution of soybean, which showed that soybean genome also underwent whole genome duplication [31]. From our findings, the ROS genes network seemed to have occurred in tandem with the evolution of tetraploid cotton, approximately 0.8-2.8 million years ago (MYA) [2]. Both synonymous site (K s ) peaks of the ROS orthologous gene pairs between the two cotton progenitors, A and D genomes, and two subgenomes (At and Dt) was 0.045 ( Figure 3C,D). Moreover, we estimated the ROS genes network divergence time between the A and D genomes to have occurred approximately 8.6 MYA, in the results was in agreement to the divergence time of the two diploid cotton, G. arboreum and G. raimondii which occurred 5-10 MYA [32]. Additionally, the synonymous (K s ) and non-synonymous (K a ) substitution values ratio (K a /K s ) were less than one (1) that suggested that the ROS gene network evolution, was primarily governed by purifying selection, similar results were also obtained in the analysis of the 156 paralogous pairs of the LEA genes in cotton, in which over 99% of the orthologous gene pairs had the K a /K s ratio values of less than 1 [27].
introgressed into the Dt subgenome, of the allotetraploid cotton, while some element of gene loss was evident, this could have been occasioned by either loose of function or chromosomal rearrangement during the evolution of allotetraploid from the two diploid cotton species, A and D cotton genomes ( Figure 3B). The synonymous substitution value (Ks) peaks of orthologous gene pairs between At/A and Dt/D was estimated at 0.004 and 0.015, respectively. The distribution of 241 putative ortholog pairs between the three cotton species/genomes (mode Ks = 0.0040-0.015) clearly indicated that the secondary peak for cotton, lies to the left of the ortholog peak, which represents a burst of gene duplications that occurred in cotton after the separation between the two lineages. Similar observation has been made in the analysis of the genome evolution of soybean, which showed that soybean genome also underwent whole genome duplication [31]. From our findings, the ROS genes network seemed to have occurred in tandem with the evolution of tetraploid cotton, approximately 0.8-2.8 million years ago (MYA) [2]. Both synonymous site (Ks) peaks of the ROS orthologous gene pairs between the two cotton progenitors, A and D genomes, and two subgenomes (At and Dt) was 0.045 ( Figure 3C,D). Moreover, we estimated the ROS genes network divergence time between the A and D genomes to have occurred approximately 8.6 MYA, in the results was in agreement to the divergence time of the two diploid cotton, G. arboreum and G. raimondii which occurred 5-10 MYA [32]. Additionally, the synonymous (Ks) and non-synonymous (Ka) substitution values ratio (Ka/Ks) were less than one (1) that suggested that the ROS gene network evolution, was primarily governed by purifying selection, similar results were also obtained in the analysis of the 156 paralogous pairs of the LEA genes in cotton, in which over 99% of the orthologous gene pairs had the Ka/Ks ratio values of less than 1 [27].    (Tables S3 and S4). Moreover, we analyzed the homeolog and ortholog gene pair's expression levels at 10 and 20 DPA in the three cotton species, G. hirsutum, G. raimondii, and G. arboreum. It is interesting that the expression levels of genes encoding the SOD proteins, such as the CSD, FSD, and MSD proteins exhibited lower expression levels at 10 and 20 DPA in G. raimondii but were significantly higher in G. arboreum and G. hirsutum. The expression levels of most of the SOD genes were higher at 10 DPA than at 20 DPA ( Figure 4A,B). The result showed that the SOD genes are integral at 10 DPA, being a critical stage of fiber initiation and elongation, moreover, at the stage of 10 DPA, ABA inhibits fiber cell initiation and elongation [33]. Moreover, previous studies have shown that cotton annexin proteins; AnxGh1:AAR13288, AnxGh2:AAB67993, AnxGhFx:FJ415173, and AnxGhF:AAC33305 are present in higher amount in fibers of 10 DPA wild-type plants as compared to the fuzzless-lintless mutant [34]. This result indicated that SOD genes were highly expressed in taxa with long fiber (G. arboreum and G. hirsutum), while were either down regulated or partially expressed in cotton species with shorter fibers, such G. raimondii of the D genome. Additionally, the expression level of SOD genes showed up regulation at fiber elongation phase (10 DPA) compared to the secondary cell wall biosynthesis/thickening phase (20 DPA). The results obtained for the expression pattern of the ROS genes, correlated positively to previous findings on the expression patterns of some of the genes known to be highly involved in fiber development in cotton, such as APX, MDAR, GR, GPX, FER, AOX, PrxR, Trx, GLR, and PER [21,22,35]. The SOD gene family and other ROS scavenging genes are believed to be involved in cotton fiber development through the change in their expression level in the course of fiber development period. On the other hand, most of the genes with functional role in the regulation of ROS, including RBOH and FRO genes were up regulated at the secondary cell wall biosynthesis phase (20 DPA) relative to fiber elongation phase (10 DPA) in G. hirsutum, G. arboreum, or G. raimondii. In summary, ROS level might be low at fiber elongation phase, but high at secondary cell wall biosynthesis phase in G. hirsutum and G. arboreum. Compared with G. hirsutum and G. arboreum, although the expression model of ROS production genes was similar in G. raimondii, the expression level of ROS removing genes was lower. Fine regulation of steady-state levels of ROS are essential for proper cell elongation [36], and long and short fiber properties were regulated by ROS removing genes, such as SOD gene family. Moreover, the Ascorbate peroxidase 1 (APX1), has been found to be integral in the process of cotton fiber development, and that APX gene profiling in upland cotton, G. hirsutum has provided valuable information on the redox homeostasis in the process of fiber development [37]. Furthermore, when the cytosolic H 2 O 2 -scavenging enzyme, APX1, the whole chloroplastic H 2 O 2 -scavenging process stops in A. thaliana and the level H 2 O 2 increases resulting into protein oxidation [38]. Ferric Chelate Reductase-2 (FRO2), Ferric Chelate Reductase-5 (FRO5), and respiratory burst oxidase homologs (RBOH6) orthologous gene pairs, modulation ROS level in fiber development, were identified by GWAS analysis and RNA-seq analysis [8]. Furthermore, the RBOH do modulate the ROS and in turn promote root growth in Arabidopsis [39].

ROS Network Genes Response to Abiotic Stress
We downloaded data of ROS-gene expression levels from pubic RNA-Seq data of G. hirsutum (TM-1). After analyzing gene expression under cold, heat, dehydration, and salt stress treatment, a total of 135 ROS genes were founded to be upregulated by abiotic stresses (Table S5). Compared with the controls, 92 (68%), 79 (59%), 74 (55%), and 85 (62.9%) genes were differently expressed under cold, heat, dehydration, and salt stress, respectively. Twenty-nine genes were differentially expressed under all four different stress conditions, which included CSD, APX, MDAR, CAT, GPX, FER, RBOH, FRO, AOX, GLR, and PER. The majority of the differentially expressed genes (DEGs) showed more

ROS Network Genes Response to Abiotic Stress
We downloaded data of ROS-gene expression levels from pubic RNA-Seq data of G. hirsutum (TM-1). After analyzing gene expression under cold, heat, dehydration, and salt stress treatment, a total of 135 ROS genes were founded to be upregulated by abiotic stresses (Table S5). Compared with the controls, 92 (68%), 79 (59%), 74 (55%), and 85 (62.9%) genes were differently expressed under cold, heat, dehydration, and salt stress, respectively. Twenty-nine genes were differentially expressed under all four different stress conditions, which included CSD, APX, MDAR, CAT, GPX, FER, RBOH, FRO, AOX, GLR, and PER. The majority of the differentially expressed genes (DEGs) showed more expression changes under cold stress than under other abiotic stresses, suggesting that ROS genes were more sensitive to cold stress. ROS-producing genes were up regulated after cold, heat, dehydration, and salt treatments. FRO5 gene shown higher expression level than other ROS-producing DEGs and was upregulated after 3 and 6 h of post stress exposure. APX1 and MDAR1 genes, showed higher expression level than other ROS-scavenging DEGs, and were upregulated after 3 and 6 h ( Figure 5). These results proved that the ROS gene do respond positively towards abiotic stress during the early developmental stages. RBOH6 and APX1 genes were highly expressed in G. hirsutum D subgenome than At subgenome which suggested that some gene pairs of G. hirsutum D subgenome and A subgenome might be playing different roles in cotton growth and development [8]. expression changes under cold stress than under other abiotic stresses, suggesting that ROS genes were more sensitive to cold stress. ROS-producing genes were up regulated after cold, heat, dehydration, and salt treatments. FRO5 gene shown higher expression level than other ROSproducing DEGs and was upregulated after 3 and 6 h of post stress exposure. APX1 and MDAR1 genes, showed higher expression level than other ROS-scavenging DEGs, and were upregulated after 3 and 6 h ( Figure 5). These results proved that the ROS gene do respond positively towards abiotic stress during the early developmental stages. RBOH6 and APX1 genes were highly expressed in G. hirsutum D subgenome than At subgenome which suggested that some gene pairs of G. hirsutum D subgenome and A subgenome might be playing different roles in cotton growth and development [8].

Analysis of Morphological and Physiological Changes Under Salt Stress
Four semi-wild (Latifolium 130, Latifolium 32, Latifolium 40, and Marie-galante 85) and two cultivated (CRI12 and CRI16) Gossypium hirsutum accessions were evaluated for phenotypic traits under salt treatment. Among them, Marie-galante 85 (M85) was highly tolerant to salinity stress, as evidenced by significantly higher relative stem length, relative root water content, relative stem water content, and lower rate of area of damaged leaves (Table 2; Figure 6). Compared with the controls, proline content was significantly enhanced under salt stress in M85, but was significantly lower in CRI12. These results suggested that M85 is highly tolerant compared to CRI12, to salt stress. The content of MDA was not significantly changed between normal condition and salt stress condition in both CRI12 and M85. However, the concentration level of MDA was relatively lower in the leaf tissues of M85 compared to the level in CRI12. Proline plays an important role in maintaining osmotic balance and protecting some important enzymes in plant cells, thus mainly acts as osmoprotectant when plants are exposed to water deficit condition, the detection of high proline content in the leaf

Analysis of Morphological and Physiological Changes under Salt Stress
Four semi-wild (Latifolium 130, Latifolium 32, Latifolium 40, and Marie-galante 85) and two cultivated (CRI12 and CRI16) Gossypium hirsutum accessions were evaluated for phenotypic traits under salt treatment. Among them, Marie-galante 85 (M85) was highly tolerant to salinity stress, as evidenced by significantly higher relative stem length, relative root water content, relative stem water content, and lower rate of area of damaged leaves (Table 2; Figure 6). Compared with the controls, proline content was significantly enhanced under salt stress in M85, but was significantly lower in CRI12. These results suggested that M85 is highly tolerant compared to CRI12, to salt stress. The content of MDA was not significantly changed between normal condition and salt stress condition in both CRI12 and M85. However, the concentration level of MDA was relatively lower in the leaf tissues of M85 compared to the level in CRI12. Proline plays an important role in maintaining osmotic balance and protecting some important enzymes in plant cells, thus mainly acts as osmoprotectant when plants are exposed to water deficit condition, the detection of high proline content in the leaf tissues of M85, showed that the genotypes has higher adaptability to salinity compared to CRI12. Therefore, the same with salt (NaCl) stress, cotton is exposed to complex alkali-salt stress could result in osmotic stress. Malondialdehyde (MDA) is the final product of lipid peroxidation, and its content can reflect stress tolerance level of the plant [36]. However, in our study, the content of MDA was not significantly changed after salinity treatment. It suggested that alkali-salt stress has not resulted in lipid peroxidation of M85 and CRI12. Taken together, cotton regulation of complex alkali-salt stress was different from salt (NaCl) stress, and M85 had a higher tolerance to osmotic stress than CRI12. Reactive oxygen species (ROS) balance was also observed. The content of hydrogen peroxide (H 2 O 2 ) and the activity of superoxide dismutase (SOD), peroxidase (POD), catalase (CAT), ascorbate peroxidase (APX), and glutathione reductase (GR) were measured. The content of hydrogen peroxide was significantly enhanced in root and leaf of CRI12 under alkali-salt stress. Hydrogen peroxide was also increased in leaf of M85, but significantly decreased in root of M85 under alkali-salt stress. Furthermore, we found that the activity of POD, CAT, APX, and GR was not significantly changed in roots of M85, suggesting that ROS production and removal was balanced in roots of M85. The content of hydrogen peroxide was higher in leaf than that in root. In addition, the activity of POD was significantly increased in leaf of M85, but POD, CAT, APX, GR, and SOD was not significantly affected. This result suggests that the leaf of M85 was better to remove excessive accumulation of ROS. Additionally, the activity of SOD was significantly increased under alkali-salt stress in root of M85, and the activity of POD and CAT was significantly enhanced in root of CRI12, suggesting that redox enzymes was act on different conditions, flexibly [40].

RNA-Seq Analysis and ROS Network Response to Alkali-Salt Stress
After Illumina Hiseq2500 sequencing, a total of 39.62-57.16 M clean reads were generated for each sample, in which approximately 85.65%-89.61% of reads were mapped to the reference genome, and 77.79%-82.27% of clean reads were mapped to the unique locations of the reference genome (Table S6). Unique mapped reads were assembled and annotated with Cufflinks and were used for analyzing mechanism of plant response to alkali-salt stress. Based on a cutoff q-value ≤ 0.01 and fold change ≥ 2, different expressed genes (DEGs) were identified. The number of DEGs was higher at 12 h post alkali-salt stress treatment than that of 3 h in both roots and leaves in CRI12 and M85 moreover; the ROS network genes analysis in the same region by principal component analysis showed genes were aggregated at a particular region ( Figure 7A,B). The ROS balance system was found to highly enriched, the finding which was consistent with the previous studies [41]. Moreover, the majority of the DEGs were found to be enriched in response to oxidative stress (GO0006979), oxidation-reduction process (GO0055114), oxidoreductase activity

RNA-Seq Analysis and ROS Network Response to Alkali-Salt Stress
After Illumina Hiseq2500 sequencing, a total of 39.62-57.16 M clean reads were generated for each sample, in which approximately 85.65%-89.61% of reads were mapped to the reference genome, and 77.79%-82.27% of clean reads were mapped to the unique locations of the reference genome (Table S6). Unique mapped reads were assembled and annotated with Cufflinks and were used for analyzing mechanism of plant response to alkali-salt stress. Based on a cutoff q-value ≤ 0.01 and fold change ≥ 2, different expressed genes (DEGs) were identified. The number of DEGs was higher at 12 h post alkali-salt stress treatment than that of 3 h in both roots and leaves in CRI12 and M85 moreover; the ROS network genes analysis in the same region by principal component analysis showed genes were aggregated at a particular region ( Figure 7A,B). The ROS balance system was found to highly enriched, the finding which was consistent with the previous studies [41]. Moreover, the majority of the DEGs were found to be enriched in response to oxidative stress (GO0006979), oxidation-reduction process (GO0055114), oxidoreductase activity (GO0016717, GO0016491, and GO0046857), and peroxidase activities (GO0004601) ( Figure S1). We also observed H 2 O 2 accumulation in leaves and roots under alkali-saline stress. The result shows the importance of ROS networks in acclimation for alkali-saline stress in Gossypium hirsutum. DEGs were also enriched in protein serine/threonine kinase activity (GO0004674). This result supports that the serine/threonine kinases, act as one of signal transduction molecular, is an essential kinase for oxidative stress responses [42]. Genes of ROS network that control ROS levels, encoded most oxidation-reduction proteins. The correlation analysis (R 2 > 0.7) between different cotton accessions through orthologous genes expression levels analysis were significantly higher ( Figure 7C). Based on the analysis of the various ROS network genes, we further predicted the interlink between fiber and abiotic stress responsive ROS genes (Figure 8). We found that CSD1-6, FSD1-2, MSD 1-2, APX3, APX11, FRO5, FRO6, and RBOH6 have a functional role in enhancing fiber development, more specifically at 10 and 20 DPA, while CSD1, APX1, APX2, MDAR1, GPX4-6-7, FER2, RBOH6, RBOH11, and FRO5 are specifically involved in enhancing abiotic stress tolerance, more salt stress. Finally, ten genes were selected from three different libraries and validated through RT-qPCR, the expression levels were same as per the RNA sequencing data ( Figure 7D). (GO0016717, GO0016491, and GO0046857), and peroxidase activities (GO0004601) ( Figure S1). We also observed H2O2 accumulation in leaves and roots under alkali-saline stress. The result shows the importance of ROS networks in acclimation for alkali-saline stress in Gossypium hirsutum. DEGs were also enriched in protein serine/threonine kinase activity (GO0004674). This result supports that the serine/threonine kinases, act as one of signal transduction molecular, is an essential kinase for oxidative stress responses [42]. Genes of ROS network that control ROS levels, encoded most oxidation-reduction proteins. The correlation analysis (R 2 > 0.7) between different cotton accessions through orthologous genes expression levels analysis were significantly higher ( Figure 7C). Based on the analysis of the various ROS network genes, we further predicted the interlink between fiber and abiotic stress responsive ROS genes ( Figure 8). We found that CSD1-6, FSD1-2, MSD 1-2, APX3, APX11, FRO5, FRO6, and RBOH6 have a functional role in enhancing fiber development, more specifically at 10 and 20 DPA, while CSD1, APX1, APX2, MDAR1, GPX4-6-7, FER2, RBOH6, RBOH11, and FRO5 are specifically involved in enhancing abiotic stress tolerance, more salt stress. Finally, ten genes were selected from three different libraries and validated through RT-qPCR, the expression levels were same as per the RNA sequencing data ( Figure 7D).

Discussion
In plants, many different biology processes are controlled by ROS [43]. For example, fiber development and biotic and abiotic stresses have previously been found to be partially regulated by ROS in cotton [44,45]. During 5-10 MYA, Gossypium split into three major diploid lineages: the New World clade (D genome), the African-Asian (A, B, E, and F genomes) clade, and the Australian clade (C, G, and K genomes) [1]. Allotetraploid Upland cotton was formed approximately 1-2 MYA [46] and domesticated for at least 4000-5000 years ago [47]. During human domestication history, the central focus was the fiber quality and quantity, but in the recent past, researchers have found that fiber and fiber qualities are not only affected by the intrinsic nature of the cotton plant, but greatly influenced by the environmental conditions, and therefore in elucidating the mechanisms of fiber development and improvement of its qualities, multidimensional approach has so far been preferred by plant breeders. The change in environmental condition, occasioned by continued human interference, rainfall has become scarce and very erratic in nature, in addition to the emergence of other forms of abiotic stresses, such as salinity, heat, and cold, among others [48]. The sessile nature of plants, makes them to be much prawn to full effects of abiotic stresses, and thus in the course of plant evolution, various survival strategies have been adopted by plants, more so the evolution of various stress responsive plants transcription factors, in which ROS genes are among one of them. In our study, we found that ROS network was conserved in cotton. The Evolutionary history of ROS genes network in cotton was similar to Gossypium evolution. Thus, ROS network was positively selected during cotton evolution and the result of Ka/Ks analysis also supported this view. Gene order and colinearity of allotetraploid cotton are largely conserved between the A and D subgenomes and the extant D progenitor genome (G. raimondii) [47], and the genes colinearity of D09 (G. raimondii)

Discussion
In plants, many different biology processes are controlled by ROS [43]. For example, fiber development and biotic and abiotic stresses have previously been found to be partially regulated by ROS in cotton [44,45]. During 5-10 MYA, Gossypium split into three major diploid lineages: the New World clade (D genome), the African-Asian (A, B, E, and F genomes) clade, and the Australian clade (C, G, and K genomes) [1]. Allotetraploid Upland cotton was formed approximately 1-2 MYA [46] and domesticated for at least 4000-5000 years ago [47]. During human domestication history, the central focus was the fiber quality and quantity, but in the recent past, researchers have found that fiber and fiber qualities are not only affected by the intrinsic nature of the cotton plant, but greatly influenced by the environmental conditions, and therefore in elucidating the mechanisms of fiber development and improvement of its qualities, multidimensional approach has so far been preferred by plant breeders. The change in environmental condition, occasioned by continued human interference, rainfall has become scarce and very erratic in nature, in addition to the emergence of other forms of abiotic stresses, such as salinity, heat, and cold, among others [48]. The sessile nature of plants, makes them to be much prawn to full effects of abiotic stresses, and thus in the course of plant evolution, various survival strategies have been adopted by plants, more so the evolution of various stress responsive plants transcription factors, in which ROS genes are among one of them. In our study, we found that ROS network was conserved in cotton. The Evolutionary history of ROS genes network in cotton was similar to Gossypium evolution. Thus, ROS network was positively selected during cotton evolution and the result of K a /K s analysis also supported this view. Gene order and colinearity of allotetraploid cotton are largely conserved between the A and D subgenomes and the extant D progenitor genome (G. raimondii) [47], and the genes colinearity of D09 (G. raimondii) was conserved with Dt05 and At05 (G. hirsutum). Interestingly, more potential ROS genes were mapped to A10 (35, 12.36%), D09 (40,13.89%), At05 (46, 18.11%), and Dt05 (38, 14.28%) in G. arboreum, G. raimondii, G. hirsutum A subgenome, and G. hirsutum D subgenome, respectively. It suggested that the location of ROS genes was conserved during the domestication process and D09, Dt05, and At05 play a core role in regulated ROS dynamic balance. Above all, ROS plays an important role in cotton growth and development, because of ROS network genes were positive selection and conserved during evolution and domestication of cotton.
At least 15 gene families (SOD, APX, MDAR, DHAR, GR, CAT, GPX, FER, RBOH, FRO, AOX, PrxR, Trx, GLR, and PER) of the ROS network were identified in Gossypium. SOD, APX, and GPX gene families were identified and analyzed. Our results are consistent with the results of SOD gene analyses in previous studies [37] and SOD genes cover the chromosome at the same locations. However, two genes were annotated as copper chaperone for superoxide dismutase in G. arboreum and G. raimondii. Compared with a previous study, the member of GPX genes family was different because of different genomes (NBI_Gossypium hirsutum vs. BGI_Gossypium hirsutum) used for identifying GPX genes. In previous studies, ROS were found to regulate development, differentiation, redox levels, stress signaling, interactions with other organisms, systemic responses, and cell death in higher plants [49][50][51]. Wild cottons allocate greater resources to stress response pathways, while domestication led to reprogrammed resource allocation toward increased fiber growth, possibly through modulating stress-response network [52]. Modulating reactive oxygen species (ROS) production and scavenging could be regulating cotton fiber elongation [53], abiotic stress [41].
In previous studies, different developmental or environmental signals feed into the ROS signaling network and perturb ROS homeostasis in a compartment-specific or even cell-specific manner [12]. Then, the fluctuant ROS is perceived by receptor proteins (NPR, HSF), redox-sensitive transcription factors, or direct inhibition of phosphatases [54]. Serine/threonine protein kinase (OXI1) is involved in ROS sensing and the activation of mitogen-activated-protein kinases (MAPKs) [55] that plays a central role in ROS signal network [56][57][58]. ROS-producing or -scavenging pathway is activated to regulate ROS levels. Different developmental, metabolic and defense pathways are activated to regulate growth and development. In cotton, several families of transcription factors, such as MYB, WRKY, AP2/ERF, and NAC have been found to be involved in regulating plant response to abiotic stress [59], and some of them are also involved in ROS signal network [60].
Cotton fiber development is a complicated biological process. Many researches demonstrate that fiber development was regulated by ROS-mediated pathway. In the previous study, 119 association signals, 71 SNP associated with yield-related traits, 45 association with fiber qualities, and three associated with resistance to Verticillium wilt, were identified through genome-wide association study, and a total of 43 ROS genes were found [61]. Among the 43 genes, seven genes were within the associated SNP signals region of 200 kb (±100 kb), and one of the 43 genes, the following genes were found to have a direct involvement in the regulation of fiber development, for instance Gh_A13G1654 (GLR34) was at the downstream, and was highly associated to signals A13:75463332 (associated with fiber elongation), A13:75463913 (associated with fiber strength), and A13:75473149 (associated with fiber length); Gh_D05G2413 (FRO2) was at the downstream, at the boll weight association signal D05: 24,174,044 (associated with boll weight). This finding may partially support the speculation that ROS network associated with some different plant traits. ROS is involved in regulating plant cell expansion [62], which suggests ROS may also be involved in regulating cotton fiber development [17]. High concentration of ROS in wild-cotton fiber cells inhibits fiber elongation by regulating secondary-cell-wall synthesis [8]. Moreover, the H 2 O 2 , involved in cell-wall relaxation, is necessary molecule for cell elongation, however, higher ROS levels may halt cell elongation [63][64][65]. Compared with fiber elongation (10 DPA), ROS level was regulated by ROS network, was upregulated at secondary cell wall biosynthesis (20 DPA) [21]. It is interesting that the expression levels of the SOD genes and some other ROS-scavenging genes were higher G. hirsutum and G. arboreum, and they are known to have relatively longer fiber compared G. raimondii. Moreover, G. hirsutum and G. arboreum, exhibited lower ROS-scavenging enzymes, which implied that they had the higher ROS-producing ability than G. raimondii. The results therefore showed that perturbed ROS levels as a result of the SOD genes and some other genes could perhaps have a net effect on the regulation of various stages of fiber development in G. hirsutum, G. arboreum, and G. raimondii. In addition, fluctuating levels in ROS levels, caused by different expression level of ROS-scavenging genes and ROS-producing genes at fiber elongation and secondary cell wall biosynthesis phases, showed the integral role played by the ROS genes in fiber development.
Almost all the abiotic stress factors do lead to the development of secondary stress condition such as osmotic and oxidative stresses [66]. RNA-seq data, physiological, and phenotypic index analysis in the two upland cotton genotypes, M85 and CRI12 suggested that the ROS gene network is partially responsible to salinity stress tolerance. The majority of the ROS genes were expressed under abiotic stress, which demonstrated that ROS network was also interconnected with abiotic stress tolerance in cotton [67][68][69]. According to the expression values (FPKM) of the ROS network genes in G. hirsutum TM-1, the expression levels of the homologous genes in under various abiotic stress conditions were: 92 (68%), 79 (59%), 74 (55%), and 85 (62.9%) under cold, heat, dehydration, and salt stress conditions, respectively. A total of 29 genes were differentially expressed under four different stress conditions, in which out of the six genes, three genes were the ROS-producing genes, while the other three were the ROS-scavenging genes.

Phylogenetic Tree, Orthologous Gene Pairs and Values of K a , K s and K a /K s
Full-length sequences of ROS-network proteins were first aligned with the ClustalX program (http://www.clustal.org/clustal2/) [73], then phylogenetic trees were constructed by using neighbor-joining (NJ) method with 1,000 bootstrap replicas [74] and the Poisson model by using MEGA 7.0 software (http://www.megasoftware.net). Meanwhile, the orthologous gene pairs of the ROS network in A, D genomes, At and Dt subgenomes were searched by InParanoid software (http://inparanoid.sbc.su.se/cgi-bin/index.cgi). Additionally, the evolutionary rates K a , K s , and K a /K s ratio were estimated by K a K s _Calculator package (https://kakscalculator.herokuapp.com/). On the basis of the synonymous substitutions per year (λ) of 2.6 × 10 −9 for cotton, the divergent time of ROS orthologous gene pairs were estimated by T = K s /2λ × 10 −6 MYA [47].

Chromosomal Locations, Conserved Synteny Blocks and Gene Duplication
Chromosomal positions of the ROS genes were obtained from the annotation information of cotton genomes (G. arboreum, G. raimondii and G. hirsutum). The chromosomal distributions of genes were drafted on the cotton chromosomes according to the gene positions by Mapchart software [75]. The conserved synteny blocks between A genome and At subgenome, D genome and Dt genome in cotton were inferred using the MCSCANX program with the default parameters (http://chibba.pgml.uga.edu/mcscan2/) [76]. A local BLASTP program was run on different genomes with parameters (e −10 ). The syntenic relationships were illustrated with the CIRCOS program (http://circos.ca/) [77]. The duplicate gene classifier program of MCSCANX software was used to identify each ROS network gene duplication with the default parameters (OVERLAP WINDOW: 5). All the ROS network genes were classified into various types of duplications.

Plant Materials and Salt Stress Treatment
G. hirsutum var. marie-galante 85 (M85) and G. hirsutum cv. CRI12 were selected for the salinity treatment. Semi-wild cotton marie-galante was originally collected from Yucatan, Mexico and cultivated in the National Wild Cotton Nursery (Sanya, Hainan, China). G. hirsutum cv. CRI12 was widely planted in china. The seeds of M85 and CRI12 were first germinated at 28 • C in a 16 h light/8 h dark cycle. Then, seedlings were planted in the normal solutions for three weeks. The similar growing seedlings were selected for salt-alkali stress. The salt-alkali treatment was imposed by simulating the actual saline condition of Xinjiang province, a region commonly known for salt stress [78]. We examined the soil salt composition and simulated a similar environment to impose salt-alkali stress, the composition was made of 2.664 mg:g −1 of calcium chloride (CaCl 2 ), 0.179 mg:g −1 of sodium hydrogen carbonate (NaHCO 3 ), 9.94 mg:g −1 of sodium sulphate (Na 2 SO 4 ), 0.848 mg:g −1 of potassium sulphate (K 2 SO 4 ) and finally 3.587 mg:g −1 of hydrated magnesium sulphate (MgSO 4 ·7H 2 O) [79].

Measurement of Morphological, Physiological and Biochemical Traits
Samples of roots and leaves that collected at 48 h after salt-alkali stress treatment were used for measuring physiological traits. And, seedlings were also used for measuring phenotype at 48 h after salt-alkali stress treatment. Malondialdehyde (MDA) and Proline (PRO), the most important protective macromolecules in response to salt stress, were measured using the kits (Nanjing Jiancheng Bioengineering Institute). Hydrogen peroxide (H 2 O 2 ), superoxide dismutase (SOD), peroxidase (POD), ascorbate peroxidase (APX), catalase (CAT) and glutathione reductase (GR), important ROS and ROS-scavenging molecules, were measured using the corresponding assay kits (Nanjing Jiancheng Bioengineering Institute). The root length (RL), stem length (SL), roots fresh weight (RWC) and stem fresh weight (SWC) were measured after 48 h salt stress. The value of SPAD, one of photosynthesis indices, was measured using chlorophyll meter model SPAD-502.

RNA Extraction, cDNA Library Construction, and RNA-Seq
Roots and leaves were collected at 0, 3, 12 and 48 h time points after salt-alkali stress treatment. Then, collected samples were used for transcriptome sequencing. Total RNA was extracted from each cotton sample using TRlzol Reagent (Life technologies, CA, USA) according to the instruction manual. RNA integrity and concentration were checked using an Agilent 2100 Bioanalyzer (Agilent Technologies, Inc., Santa Clara, CA, USA). The mRNAs were isolated by NEBNext Poly (A) mRNA Magnetic Isolation Module (NEB, E7490) by New England BioLabs (Ipswich, MA, USA). The cDNA libraries were constructed by following the manufacturer's instructions of NEBNext Ultra RNA Library Prep Kit for Illumina (NEB, E7530) and NEBNext Multiplex Oligos for Illumina (NEB, E7500) (Illumina, San Diego, CA, USA). Briefly, the enriched mRNA was fragmented into RNAs with approximately 200 nt, which were used to synthesize the first-strand cDNA and then the second cDNA. The double-stranded cDNAs were performed end-repair/dA-tail and adaptor ligation. The suitable fragments were isolated by Agencourt AMPure XP beads (Beckman Coulter, Inc., Indianapolis, IN, USA) and enriched by PCR amplification. Finally, the constructed cDNA libraries were sequenced on a flow cell using an Illumina HiSeq™ 2500 sequencing platform.

Expression Analysis
Whole-transcriptome sequencing data for G. hirsutum were obtained from the NCBI Sequence Read Archive (SRA) (http://www.ncbi.nlm.nih.gov/sra) to analyze the tissue/organ-specific, stage-specific and stress-induced expression patterns of cotton ROS genes [47]. The details of the SRA are shown in (Tables S7 and S8). The row data of RNA-seq of M85 and CRI12 were separately analyzed and the clean reads were obtained by removing reads containing adapter, reads containing ploy-N and lower quality reads from raw data. At the same time, Q 30 [80], GC-content and sequence duplication level of the clean data were calculated. Raw sequences were transformed into clean reads after data processing. These clean reads were then mapped to the reference genome sequence. Only reads with a perfect match or one mismatch were further analyzed and annotated based on the reference genome. Tophat2 [81] software was used to map with reference genome. Base on the reference genome, using Cufflinks software [82], mapped reads had been assembled. Quantification of gene expression levels was estimated by fragments per kilobase of transcript per million fragments mapped (FPKM). Base on mapped reads, using FPKM as the index, each gene was estimated by Cuffquant and Cufform. EBSeq software was used to identify the differential expression genes by Fold change ≥ 2 and FDR ≤ 0.01 (False Discovery Rate). FDR was corrected using Benjamini-Hochberg method by p-value. Finally based on the RNA sequence expression data, ten (10) genes were selected and their expression validated through RT-qPCR in relation to salt-alkali stress. In each, three libraries were chosen and the gene specific primers were designed by using primer premier 5 (Table S9). The RT-qPCR analysis was carried out as described by Magwanga et al [83], with GhActin as the internal control gene.

Statistical Analysis
The experimental data on phenotype and physiology were statistically analyzed using the SAS version 8.2. Graphic presentations were performed using OrginPro 8.0 program and R packages.

Conclusions
Due to the ever-changing environmental conditions, the available land for cultivation is declining as a result of various environmental stress factors such as drought and salinity, among others [84]. The combined effect of drought and salinity has been estimated to cause huge losses in agricultural crops as over 6% of arable lands are saline [85]. Thus non edible crops such as cotton are currently being grown in harsher conditions, in soils with high levels of salt content, but due to the importance of the cotton plant, focus is aimed of developing a more robust and salt tolerant cotton cultivars. When plants are exposed to any form stress, the ROS equilibrium shifts leading to high concentrations which results in toxicity and eventually leads to total plant death [86]. In this research work, we carried out deep transcriptome analysis of the ROS network genes, evaluated their evolution pattern and their expression levels under abiotic stress and fiber development stages. A total of 515, 260, and 261 ROS candidate proteins were identified in G. hirsutum, G. arboreum, and G. raimondii, respectively. The high numbers of the ROS network genes were found to have evolved through segmental gene duplication compared to tandem, a result which was in agreement to evolution pattern of various stress responsive genes such the LEA genes which evolved mainly through segmental type of gene duplication [27]. Moreover, RNA sequencing and RT-qPCR validation revealed that the ROS genes had a putative role in enhancing abiotic stress tolerance in cotton and in turn promote fiber development in cotton. In relation to fiber development, the following genes were found to be highly upregulated, CSD1,2,3,5,6; FSD1,2; MSD1,2; APX3,11; FRO5.6; and RBOH6 while for salt stress tolerance, CSD1, APX1, APX2, MDAR1, GPX4-6-7, FER2, RBOH6, RBOH11; and FRO5 were found to have a significant role. This research lays a significant foundation for future exploration of the ROS genes in developing a more stress tolerant cotton genotypes and with superior fiber quality.
Supplementary Materials: Supplementary Materials can be found at http://www.mdpi.com/1422-0067/20/8/1863/ s1. Table S1. Characteristics of genes of ROS network in Gossypium hirsutum, G. arboreum, and G. raimondii; Table  S2. Synteny blocks of G.hirsutum and diploid cotton (G. arboreum and G. raimondii); Table S3. Orthologous gene pairs of G. hirsutum, G. arboreum, and G. raimondii; Table S4. Expression levels of ROS network genes of G. hirsutum (TM-1), G. arboreum, and G. raimondii fiber; Table S5. Expression levels of ROS network genes of G. hirsutum (TM-1) leaves under abiotic stress; Table S6. Summary of the results of RNA-seq and map to the G. hirsutum genome; Table S7. Gene accession number and samples information of transcriptome data used in our research; Table S8. Genes of ROS network within 2 Mb (±1 Mb) of association signals; Figure S1. GO enrichment analysis of M85 and CRI12 under alkali-salt stress.