Comprehensive analysis of WOX genes uncovers that WOX13 is involved in phytohormone-mediated fiber development in cotton

Background The WOX (WUSCHEL-RELATED HOMEOBOX) gene family encodes a class of transcription factors that are unique to green plants, where they are involved in regulating the development of plant tissues and organs by determining cell fate. Although the importance of the WOX gene is well known, there are few studies describing their functions in cotton. Results In this study, 32 WOX genes were found in Gossypium hirsutum. Phylogenetic analysis showed that WOX proteins of cotton can be divided into three clades: the ancient, intermediate, and WUS clades. The number of WOX proteins in the WUS clade was greater than the sum of the proteins in the other two clades. Our analysis revealed that 20 GhWOX genes are distributed on 16 cotton chromosomes and that duplication events are likely to have contributed to the expansion of the GhWOX family. All GhWOX genes have introns, and each GhWOX protein contains multiple motifs. RNA-seq data and real-time PCR showed that GhWOX13 gene subfamily is specifically expressed at a high level in cotton fibers. We also identified putative GA, NAA, and BR response elements in the promoter regions of the GhWOX13 genes and GhWOX13 transcripts were significantly induced by GA, NAA, and BR. Conclusions Our data provides a useful resource for future studies on the functional roles of cotton WOX genes and shows that the GhWOX13 genes may influence cotton fiber development. Our results also provide an approach for identifying and characterizing WOX protein genes in other species. Electronic supplementary material The online version of this article (10.1186/s12870-019-1892-x) contains supplementary material, which is available to authorized users.

initiating vascular primordium of the cotyledons and emerging leaves, and over-expression of AtWOX1 caused small leaves, dwarf plant and low fertility [11,12]. Ectopically over-expression of TaSF (STENTOFOLIA), the ortholog of AtWOX1, results in transgenic plants with broadened leaves, accelerated flowering and increased chlorophyll content [13]. WOX3 has a highly conserved function during the recruitment of founder cells to form the lateral domains of vegetative and floral organs in Arabidopsis and maize. A loss-of-function AtWOX3 mutant exhibited a degeneration of the lateral stipules and sepals [14]. AtWOX4 is required for auxin-dependent procambium differentiation and/or maintenance, which shows that the WOX gene family influences the auxindependent regulation of lateral plant growth [15,16]. Rice WOX4 is responsible for meristem maintenance and reduction of OsWOX4 expression led to severe malformed leaf primordia [17]. AtWOX5, which is induced by auxin, is a negative factor to regulate IAA homeostasis to maintain stem cells in root apical meristem (RAM) [18]. In Populus, WOX5a was mainly expressed in the adventitious root (AR) tip and lateral root tip. Over-expression of PtoWOX5a led to increased AR number, decreased AR length and leaf number [19]. Also, over-expression of WOX11 and WOX12 significantly increased AR number and decreased ectopic roots in poplar [20]. A homolog of AtWOX11, OsWOX11, interacts with H3K27me3 demethylase JMJ705 to promote the expression of downstream gene during shoot development [21]. In addition, OsWOX11 is also involved in lateral root initiation, root hair formation, and responses to abiotic stresses [22,23].
Cotton is the important economic crop in the world, which can be used for textile. Cotton fiber is the single seed coat epidermal cell. Quality and yield of fiber depends mainly on three biological processes; fiber initiation, fiber elongation, and secondary cell wall deposition [24][25][26]. Many studies have focused on fiber cell growth in cotton. So far, genome sequences of three cotton species have been completed, including the allotetraploid G. hirsutum. (AADD, 2n = 52) and the two diploid progenitors G. arboreum (AA; 2n = 26) and G. raimondii (DD, 2n = 26) [27][28][29][30][31]. These genome sequences have enabled the identification of WOX genes at a genome-wide scale in cotton.
The roles of WOX proteins have been welldocumented in Arabidopsis, rice, and maize, however, the functions of WOX proteins in cotton, especially in fiber development, are largely unknown. In this study, we identified 32 cotton WOX genes and found that they cluster into the previously-known ancient, WUS, and intermediate phylogenetic clades. Chromosomal distribution analysis revealed that 20 GhWOX genes are distributed on 16 cotton chromosomes. Analysis of gene expression patterns showed that GhWOX13 genes are specifically expressed at high levels in the cotton fiber. In addition, we found that all of the GhWOX13 genes had putative GA, NAA, and BR response elements in their promoter regions and can be induced by these hormones. Taken together, the results of our study provide a useful resource for future studies on the functional roles of cotton WOX genes, and demonstrated that the GhWOX13 gene products may act to influence cotton fiber development. We also provide an approach for identifying and characterizing WOX proteins in other plant species.

Identification of WOX gene family proteins in Gossypium
To identify WOX family members in Gossypium, 15 AtWOX protein sequences were used as queries for blast searching of the cotton genome database with the E-value cutoff of 0.001. The candidate GhWOX proteins were manually inspected to ascertain that the candidates contained the conserved DNA-binding homeodomain using the InterPro protein sequence analysis & classification tool. As a result, we identified 32, 18, and 19 WOX protein coding genes in G. hirsutum, G. arboreum, and G. raimondii, respectively (Table 1). We found that the number of WOX genes identified in the tetraploid G. hirsutum was lower than the combined number from its diploid progenitors suggesting some loss after polyploid formation. Compared with Arabidopsis, multiple WOX genes, such as WOX3, WOX4, WOX13 and WUS, were expanded in Gossypium. On the contrary, WOX7 and WOX8 were missing in three cottons, whereas they were found in Arabidopsis, suggesting that WOX7 and WOX8 genes may have been lost in Gossypium after they diverged from the common ancestor of Gossypium and Arabidopsis.
Phylogenetic relationships of the WOX proteins from G. hirsutum, G. arboreum, G. raimondii, and A. thaliana To evaluate the evolutionary and phylogenetic relationships of the WOX gene families among species, a total of 89 WOX protein sequences, including 32 GhWOXs, 18 GaWOXs, 19 GrWOXs, and 15 AtWOXs sequences, were used to construct an unrooted phylogenetic tree based on the alignment using the Neighbor-Joining (NJ) method as implemented in MEGA 6.0 with 1000 bootstrap replicates (Fig. 1). According to the phylogenetic relationship of WOX family, the GhWOX family proteins were also divided into three clades: the WUS clade, the ancient clade and the intermediate clade, The WUS clade, which was the largest clade, contains 23 GhWOX proteins. The ancient clade has six GhWOX proteins. The remaining members belong to the intermediate clade.

Chromosomal locations of the GhWOX genes
We investigated the chromosomal location of the 32 GhWOX genes based on the annotation of the G. hirsutum genome. Genomic localization of the GhWOX genes showed that 20 GhWOX genes map to 16 chromosomes, while 12 GhWOX genes map to scaffold sequences that are not assigned to chromosomes (Fig. 2). Cotton chromosomes At_02, At_09, Dt_08, and Dt_09 harbor two GhWOX genes each, while all other chromosomes have only one. However, no GhWOX genes were found on chromosomes At_06, At_07, At_10, At_12, Dt_04, Dt_ 06, Dt_07, Dt_12, and Dt_13. We analyzed duplication events in the GhWOX genes mapped to chromosomes, and identified three types of gene duplication: tandem duplication, dispersed duplication, and segmental duplication (Additional file 2: Table S2). Except for GhWOX12_At and GhWUS1b_Dt which were derived from dispersed and tandem duplications, respectively, the other mapped GhWOX genes all resulted from segmental duplications This result indicated that segmental duplications were the main contributors to the expansion of GhWOX genes in G. hirsutum.

GhWOX gene structure and domain analysis
In order to better understand and explore the diversity of the WOX genes in G. hirsutum, we analyzed the conserved motifs and exon/intron structures of the GhWOX genes. The gene structures of GhWOX was generated by aligning the cDNA sequences and corresponding genomic DNA sequences. An unrooted phylogenetic tree was obtained with predicted GhWOX protein sequences (Fig. 3a). Exon/intron structure analysis showed that most GhWOX genes contain two introns, and the others contain no less than one intron, suggesting that GhWOX gene structure is conserved (Fig. 3b).
We also analyzed the conserved motifs of the GhWOX proteins using the InterPro protein sequence analysis server, and five distinct conserved motifs were identified. The composition patterns of the motifs tend to be consistent with the phylogenetic tree. Motifs 1 and 2 are present in all GhWOX protein subfamilies, motif 4 is present in 24 GhWOX proteins except for GhWOX9, GhWOX10, GhWOX13, and GhWOX14. Motif 3 is present in GhWOX10, GhWOX13, and GhWOX14, and motif 5 is found only in GhWOX1 and all members of the GhWOX4 subfamily (Fig. 3c). The GhWOX proteins from the ancient clade contain three identical motifs, which could indicate that the proteins in this clade have the same function, and that the functions are redundant. In addition, a new motif was found that is specific to the GhWOX1_At protein, although it shared three identical motifs with GhWOX6A. The two proteins are very close in the phylogenetic tree, suggesting that they may have a shared function. However, a new motif could contribute a novel function to the GhWOX1_At protein.

Expression patterns of the GhWOX genes in various tissues
To investigate the tissue-specific expression profiles of 32 GhWOX genes, their relative expression levels in different organs, including young roots, young stems, young leaves, flowers, and 0 dpa ovules were analyzed from RNA-seq data. As shown in Fig. 4, the data indicate that expression of the GhWOX genes varies in the different organs and tissues; GhWOX2, GhWOX4, GhWOX5, GhWOX9, GhWOX11, and GhWOX13b_Dt are highly expressed in young roots, GhWOX4, GhWOX5, GhWOX11, and GhWOX13b_At are highly expressed in young stems, and GhWOX9_At, GhWUSa_ At, and GhWUSb_Dt are highly expressed in young leaves, implying that the GhWOX genes have different roles in the growth and development of cotton. These results are consistent with those of previous reports on the function of individual WOX genes, which suggests that WOX gene functions are conserved in different species. Interestingly, all of the GhWOX genes have nearly identical expression levels in the ovules. We further investigated the expression profiles of GhWOX genes and found that none of them displayed a significant change during the fiber development stage except for GhWOX13a_Dt, GhWOX13b_At, and GhWOX13b_D., which were highly expressed in fibers at 10 and 15 dpa (Fig. 5a). Furthermore, we verified the RNA-seq data using quantitative real-time PCR (qRT-PCR) experiment. As a result, GhWOX13 were consistent with the RNA-seq data. The transcripts of GhWOX13a_Dt and GhWOX13b_At gradually increased and reached a peak in the fibers at 10 dpa (Fig. 5b). These results indicate that GhWOX13 genes may play an important role in cotton fiber development.
GhWOX13 gene expression is stimulated by GA, BR, and auxin In order to identify the putative regulatory elements that control expression of the WOX genes in cotton, we analyzed the promoter regions of the GhWOX genes. We found that most of the GhWOX gene promoters contained at least one putative GA-response element (i.e. GARE motif, P-box, and TATC-box) except for GhWOX1_At, GhWOX3a_At, GhWOX3b_At, GhWOX3b_Dt, GhWOX4b_Dt, GhWOX4c_ At, GhWOX12_Dt, GhWOX13a_Dt, and GhWOX14_Dt (Additional file 1: Figure S1). We next analyzed the putative BR-response element and found that all GhWOX genes had putative BR-response elements in their promoter regions (Additional file 1: Figure S2). We also searched for putative auxin-related cis-elements in the promoter regions of GhWOX genes, and found that GhWOX3a_ Dt, GhWOX4b_Dt, GhWOX4c_At, GhWOX5_Dt, GhWOX9_At, GhWOX13b_At and GhWOX3b_Dt contain a putative auxin response element (AuxRE) in their upstream promoter regions (Additional file 1: Figure S3). Since the GhWOX13 genes are specifically expressed in cotton fibers, we analyzed the promoter regions of the GhWOX13 genes and identified putative GA, NAA, and BR responsive cis-elements in the promoters (Fig. 6a, c and e). In order to determine whether GhWOX13 genes are induced by GA, NAA, and BR, the The results showed that the expression of GhWOX13a_Dt and GhWOX13b_Dt was up-regulated after treatment with GA for 48 h, the GhWOX13b_At mRNA level was significantly increased after treatment with GA for 24 h, but the transcription of GhWOX13a_At showed no change in response to GA treatment (Fig. 6b). GhWOX13a_Dt, GhWOX13b_At and GhWOX13b_Dt were up-regulated, with expression peaking after 24 h of BR treatment, while transcription of GhWOX13a_At showed no response to BR (Fig. 6d). Expression of GhWOX13b_At and GhWOX13b_Dt was induced by NAA, while the expression of GhWOX13a_At, GhWOX13a_Dt showed no effect (Fig. 6f). These results indicate that the GhWOX13 genes may play an important role in cotton fiber elongation mediated by plant hormones.

Discussion
Cotton is the one of the most significant economic crops cultivated worldwide, and it is an important source of fiber and oil. Whole genome sequencing of the G. arboreum, G. raimondii, and G. hirsutum genomes has been completed recently, and a number of genes associated with fiber development and leaf formation have been discovered by combining genome analysis with highthroughput transcriptome sequencing. However, the functions of the WOX genes in cotton, especially as they relate to fiber development, are largely unknown. In this study, we identified 32, 18, and 19 WOX genes in G. hirsutum, G. arboreum, and G. raimondii, respectively. Although the total number of WOX genes in G. hirsutum is much higher than in the other two species, the number of genes in the At and Dt genomes of G. hirsutum are less than the numbers in G. arboreum or G. raimondii. Gene copy numbers were largely reduced during the cotton-specific polyploidization event [32- b Expression profiling of the GhWOX13 subfamily genes during cotton fiber development. Gene expression data were obtained from qRT-PCR assays with three independent replicates. The G. hirsutum ubiquitin gene GhUBQ7 was used to normalize the expression of the GhWOX genes 34]. Our results are consistent with the finding that the rate of gene loss in allotetraploid cotton was higher than that in the two ancestral diploid species. Previous research has shown that WOX gene family can be divided into three major clades; the ancient clade which is present in green algae and all lineages of land plants, while intermediate and WUS clades are only found in ferns and seed plants [35]. In this study we identified 32 WOX genes representing all three major clades in the allotetraploid cotton species G. hirsutum. Yang et al. identified WOX gene family members in Upland cotton and also found that WOX genes could be naturally classified into three clades [36], which is consistent with our finding. Four of the GhWOXs belonged to the intermediate clade and six belonged to the ancient clade. The number of genes in the WUS clade was much higher than in either the intermediate or ancient clades, which is consistent with a previous report [37]. In A. thaliana, Gene structure analysis predicted that 15 GhWOX genes contain one intron, 15 GhWOX genes contain two introns, and two GhWOX genes (GhWOX1_At and GhWOX6_Dt) contain three introns. GhWOX1_At and GhWOX6_Dt have nearly identical gene structures within the WOX1/6 subgroup, although the intron phases and lengths differ. Similarly, the Arabidopsis genes AtWOX1 and AtWOX6 are in one phylogenetic clade and share a similar gene structure with three introns. Protein sequence motif analysis revealed that all GhWOX subfamilies contain motifs 1 and 2. Motif 3 is only found in GhWOX10, GhWOX13, and GhWOX14, and these proteins belong to the ancient clade ( Fig. 1), indicating that motif 3 is conserved in the most ancient group of WOX proteins. These results suggest that the specific functions of proteins in the different subfamilies may be due to the presence of specific sequence motifs. The patterns of introns and motifs, which correlate well with the phylogenetic clades, strongly support the close evolutionary relationships among the GhWOX genes within each of the subfamilies.
The WOX proteins have been shown to act a pivotal part in numerous developmental processes, such as organ formation, embryo patterning, stem cell maintenance [38,39]. Our data show that many GhWOX genes have different expression patterns with the homologous genes in Arabidopsis. For example, in contrast to WOX2 expression in the Arabidopsis embryo, we found that GhWOX2a_At/2a_Dt/2b_At/2b_Dt, are specifically expressed in the young root, and exhibit an expression profile similar to that of WOX2 in Populus tomentosa, which is mainly expressed in the roots [40]. We also noted that the duplicated genes show varying expression levels; GhWOX13a_Dt and GhWOX13b_Dt are primarily expressed in flowers, while GhWOX13b_At is expressed at a high level in young stems. These results suggest that some of WOX genes might acquire diverse functions during evolution.
GhWOX13a_Dt, GhWOX13b_At, and GhWOX13b_Dt are highly expressed in the cotton fiber, and their expressions gradually increase during fiber elongation. We analyzed the promoters of the GhWOX13 genes and identified putative GA, NAA, and BR response cis-elements in the promoters; we also found that the GhWOX13 genes are up-regulated after treatment with GA, BR, or NAA. AtWOX13 have high expression levels in inflorescences, floral buds, emerging lateral roots and root tip [41,42]. The wox13 mutant showed slightly wider fruits with a reduced number of lateral roots in Arabidopsis. Our data indicate that GhWOX13 genes may involved in phytohormone mediating fiber elongation in cotton.
Taken together, our work not only comprehensively identified the WOXs in G. hirsutum, but also found that GhWOX13 transcripts were expressed during cotton fiber development and can be significantly induced by GA, NAA, and BR. Our research provides a clue for further uncovering the precise roles of WOX genes in cotton fiber development.

Conclusions
Gossypium hirsutum contained 32 WOX proteins and can be divided into three clades: the ancient, intermediate, and WUS clades. The number of WOX proteins in the WUS clade was more than that in the remaining two clades. Multiple introns and motifs were found in GhWOX genes. GhWOX13 genes were specifically and highly expressed in cotton fiber cells and they transcripts can be induced by GA, NAA, and BR. Our results also provide an approach for identifying and characterizing WOX genes in other species.
Promoter analysis, multiple sequence alignment, and phylogenetic analysis Sequences of 2000 nucleotides prior to the start codon were extracted from the genomic sequences of the GhWOX genes, and the cis-acting regulatory elements were predicted by the program PlantCARE online (http://bioinformatics.psb.ugent.be/webtools/plantcare/ html/). The gene IDs of the WOX genes from the above three species are shown in Additional file 2: Table S1. The full predicted amino acid sequences of the WOX genes from the four plant species were aligned using ClustalW. A Neighbor-Joining (NJ) phylogenetic tree was constructed by MEGA 6.0 software with pairwise deletion option, poisson correction model and uniform rates. Bootstrap test were carried out with 1000 replicates for evaluating the statistical reliability of the phylogenetic tree [43].
Gene structure, protein motif analysis, and genomic location The Gene Structure Display Server (GSDS, http://gsds. cbi.pku.edu.cn/) was used to draw the exon-intron structures of GhWOX genes based on the full-length genomic sequences and the corresponding coding sequences obtained as described above [44]. InterPro, protein sequence analysis & classification (http://www. ebi.ac.uk/interpro/scan.html), was used to identify conserved domains in the GhWOX proteins [45]. Chromosomal position information of each GhWOX gene was acquired from the genome annotation, and the relative locations were shown on the respective chromosomes in a top-to-bottom orientation.

Plant materials and hormonal treatments
Gossypium hirsutum ('Xuzhou 142') seed was obtained from Institute of Cotton Research of Chinese Academy of Agricultural Sciences (Anyang, China). Cotton plants were grown in a climate-controlled greenhouse with a 16 h light/8 h dark photoperiod at 30°C. Fresh cotton seed fibers were harvested from bolls at various days post anthesis (DPA) at the indicated time points and then immediately frozen in liquid nitrogen as previously reported [46].
The plant hormone treatments were performed using previously-described methods [46]. Briefly, selfpollinated cotton flowers were collected at 1 DPA (days post anthesis). Subsequently, bracts, sepals and petals were removed and the ovaries were sterilized using 10% hypochlorous acid for 15 min, and then rinsed with sterile distilled water for seven times. The ovules were carefully dissected from the ovaries under sterile conditions and immediately floated on the liquid medium  [46]. Then the ovules were cultured at 30°C in darkness without agitation. For GA and NAA treatment, ovules incubated for 12 h, 24 h, and 48 h were collected together with fibers because of the difficulty of separating them from each other. For BR treatment, 6 h, 24 h, and 5 d ovules were sampled together with fibers. The samples were frozen in liquid nitrogen immediately after collection and stored at − 80°C prior to use in the experiments. For each treatment, three individual samples were collected and the analysis was performed on the three biological replicates. GA, BR, and NAA were purchased from Sigma-Aldrich.

RNA extraction and qRT-PCR analysis
Total RNA was extracted using the E.Z.N.A.® Plant RNA Kit (OMEGA) according to the manufacturer's instructions. The tissues were disrupted and homogenized. The column flow-through at the very last step was mixed with the membrane-binding solution and then loaded into the HiBind RNA Mini column. Finally, RNA was washed with RWC buffer and RNA wash buffer to remove proteins, polysaccharides, and salts. The total RNA was treated with DNaseI to remove gDNA contamination. First-strand cDNA was then synthesized using the PrimeScript RT Reagent Kit following the manufacturer's protocol (TaKaRa, Dalian, China).
The qRT-PCR assays were performed using a Bio-Rad Real Time PCR detection system (Bio-Rad CFX96 Touch) [47]. The SYBR Green qRT-PCR reactions contained 10 μl SYBR® Premix Ex Taq™ II (Takara), 0.5 μl of 10 μM primers, and 2 ng cDNA template. The final volume was adjusted to 20 μl with ddH 2 O. The PCR cycling conditions were 95°C for 30 s, followed by 40 cycles of 95°C for 5 s and 60°C for 30 s. Relative gene expression was calculated using the 2 −ΔΔCt method [47]. The experiment was performed with three biological replicates and each biological replicate was performed with three technical replicates. Primers used for qRT-PCR analysis are shown in Additional file 2: Table S3. The cotton ubiquitin gene GhUBQ7 was used as the internal control for normalization of gene expression in each qRT-PCR experiment.

Additional files
Additional file 1: Figure S1. Identification of putative GA-related cis-elements in the promoter regions of the GhWOX genes. Figure S2. Putative BR response elements identified in the GhWOX gene promoter regions. Figure S3. Detection of putative auxin-related cis-elements in the GhWOX gene promoter regions. (DOCX 3317 kb) Additional file 2: Table S1. Analysis of the G. hirsutum WOX gene family members and their orthologs in the AA and DD genome cotton species. Table S2. Analysis of duplication events in the G. hirsutum WOX genes mapped to chromosomes. Table S3. PCR primers used for analysis of duplication events in the G. hirsutum WOX13 gene family members