Genome-wide analysis of WOX genes in upland cotton and their expression pattern under different stresses

Background

WUSCHEL-related homeobox (WOX) family members play significant roles in plant growth and development, such as in embryo patterning, stem-cell maintenance, and lateral organ formation. The recently published cotton genome sequences allow us to perform comprehensive genome-wide analysis and characterization of WOX genes in cotton.

Results

In this study, we identified 21, 20, and 38 WOX genes inGossypium arboreum(2n = 26, A₂),G. raimondii(2n = 26, D₅), andG. hirsutum(2n = 4x = 52, (AD)_t), respectively. Sequence logos showed that homeobox domains were significantly conserved among the WOX genes in cotton,Arabidopsis, and rice. A total of 168 genes from three typical monocots and six dicots were naturally divided into three clades, which were further classified into nine sub-clades. A good collinearity was observed in the synteny analysis of the orthologs from At and Dt (t represents tetraploid) sub-genomes. Whole genome duplication (WGD) and segmental duplication within At and Dt sub-genomes played significant roles in the expansion of WOX genes, and segmental duplication mainly generated the WUS clade. Copia and Gypsy were the two major types of transposable elements distributed upstream or downstream of WOX genes. Furthermore, through comparison, we found that the exon/intron pattern was highly conserved betweenArabidopsisand cotton, and the homeobox domain loci were also conserved between them. In addition, the expression pattern in different tissues indicated that the duplicated genes in cotton might have acquired new functions as a result of sub-functionalization or neo-functionalization. The expression pattern of WOX genes under different stress treatments showed that the different genes were induced by different stresses.

Conclusion

在当前的工作,WOX基因,分为三种clades, were identified in the upland cotton genome. Whole genome and segmental duplication were determined to be the two major impetuses for the expansion of gene numbers during the evolution. Moreover, the expression patterns suggested that the duplicated genes might have experienced a functional divergence. Together, these results shed light on the evolution of the WOX gene family, and would be helpful in future research.

WUSCHEL-related homeobox (WOX), one of the sub-clades of homeodomain (HD) superfamily, is a plant-specific homeobox (HB) transcription factor, characterized by a short stretch (60–66 residues) of amino acids that form a DNA-binding domain named as homeodomain [1]．Previous reports have shown that WOXs assay a wide variety of important roles in development and growth process of plants, such as in embryonic patterning, stem cell maintenance, and organ formation [2,3]．InArabidopsis thaliana, which is a model plant, the evolution and function of WOXs have been well studied and characterized. As of date, 16 WOXs have been discovered by whole genome identification [4]; these have been clustered into three clades, namely, the ancient, WUS, and intermediate clades, by phylogenetic analysis [5]．WOXs of lower plants, such as those of green algae and non-vascular moss,Physcomitrella patens, belong only to the ancient clade, whereas those of higher plants, such asArabidopsis, sorghum, maize, and rice, are present in all the three clades [1,2]．AtWUSexpressed in shoot apical meristem, ovule, and anther has been proven to play important roles in stem cell maintenance [6]．Ectopic expression ofAtWUSin upland cotton was demonstrated to promote callus dedifferentiation resulting in the formation of embryo callus [7]．AtWOX3, as well as its orthologs in maize,NS1(narrow sheath1) andNS2(narrow sheath2), perform a highly conserved function during the recruitment of founder cells to form lateral domains of vegetative and floral organs [8]．AtWOX1was demonstrated to express in the initiating vascular primordium of the cotyledons during heart and torpedo stages, and its overexpression resulted in defects in the development of meristem and dwarf phenotype [9]．AtWOX4was determined to be an essential regulator in auxin-dependent cambium stimulation that regulates lateral plant growth [10]．AtWOX5, induced by turanose and auxin, was shown to play vital roles in a correct pattern of root-formation through mediation of auxin homeostasis and by maximizing the auxin content in a restricted area [11]．AtWOX6/PSF2was revealed to play important roles in ovule development via the regulation of cell proliferation of the maternal integuments and through differentiation of the megaspore mother cell [12]．WOX6/HOS9was also shown to be an essential component, mediating cold tolerance inArabidopsisthrough a CBF-independent pathway [13]．AtWOX2andAtWOX8were demonstrated to be essential for determination of the boundary between the cotyledons through the activation of threeCUPSHAPED COTYLEDON(CUC) genes [14]．STIMPY/AtWOX9was found to be crucial for the growth of vegetative shoot apical meristem through maintenance of cell division and prevention of premature differentiation [15]．Also,STIMPY/AtWOX9determined the fate of meristem by activation of the cytokinin signaling pathway in meristematic tissue [16]．WOX11, being a key regulator of shoot-borne crown development, was associated with the activation of crown root emergence and growth via directly repressingPR2to regulate cell proliferation during crown root development [17]．AtWOX11directly responds to an auxin maximum, induced by wounding in and around the procambium, and likeAtWOX12, it positively upregulatesLATERAL ORGAN BOUNDARIES DOMAIN 16and29genes, resulting in the initiation of a leaf procambium or the transition of its neighboring parenchyma cells to root founder cells [18]．AtWOX13mainly expresses in meristematic tissues including replum and promotes replum during the development of theArabidopsis thalianafruits [19]．AtWOX14andWOX4, regulatory elements downstream ofPHLOEM INTERCALATED WITH XYLEM(PXY), regulate vascular cell division instead of vascular organization, playing crucial roles in stem formation [20]．

Cotton is an important fiber crop, which provides the natural renew fiber for textile industry [21]．The roles of WOXs have been well-documented in embryogenesis inArabidopsis; however, the functions of WOXs in cotton, especially in somatic embryogenesis, are largely unknown, thus far. The completion of genome sequencing in cotton allows comprehensive identification and analysis of WOXs in cotton [22,23,24,25,26]．We, therefore, conducted a thorough investigation on WOXs in cotton; the study included the identification of gene families, phylogenetic tree analysis, as well as the analyses of segmental duplication, gene structure, chromosome location, and expression pattern.

Sequence identification

The complete genome sequence data of three cotton species,Gossypium arboreum(BJI, version 1.0),G. raimondii, (JGI, version 2.0), andG. hirsutum(NAU, version 1.1, BJI, version 1.0), available from COTTONGEN (http://www.cottongen.org) [27] were used. The rice (version 7.0), sorghum (version 2.1), cacao (version 1.1), poplar (version 1.1), and maize (version 1.1) genome sequence data were retrieved from JGI (https://phytozome.jgi.doe.gov/pz/portal.html). The amino acid sequences of WOXs fromArabidopsis thalianawere acquired from TAIR 10 (http://www.arabidopsis.org); these were used as query sequences to search theG. arboreumprotein database for candidate sequences employing blastp program. Thereafter, Interproscan 56.0 [28] was used to search for the HB domain (IPR001356) in the candidate sequences, and eventually the WOX sequences were identified. WOXs in rice, sorghum, cacao,G. hirsutum, andG. raimondiiwere identified using the same method as employed forG. arboreum。

Conserved sequence and phylogenetic analysis

Multiple sequence alignment was performed using Clustal X 2.0. For sequence logo analysis, the conserved HB domain sequences of WOXs from rice,Arabidopsis, and upland cotton were aligned, and the multiple alignment result was submitted to an online tool, WEBLOG [29], for generating the logos. For phylogenetic analysis, the full-length WOX sequences fromArabidopsis, rice, sorghum, cacao,G. arboreum,G. raimondii, andG. hirsutumwere aligned, and MEGA 7.0 [30] was used to construct a neighbor-joining (NJ) tree. Bootstrap method was used to test the tree with 1000 replicates. Substitution was evaluated by Poisson model using the default parameters. To validate the phylogenetic tree, constructed using the NJ method, the minimum-evolution method was also used. The bootstrap method was used to test the tree with 1000 replicates.

Chromosome location and collinearity analysis

The gene loci of WOXs were extracted from the annotation gff3-file. Mapchart was then used to obtain the chromosome location [31]．All the protein sequences of upland cotton were included in a local database using Basic Local Alignment Search Tool (BLAST). The entire protein sequences were used as queries to search the above-mentioned database with an e-value of 1e-5. The blastp result was analyzed by MCSCAN to produce the collinearity blocks across the whole genome. The collinearity pairs belonging to WOX family were extracted to draw a collinearity map within WOXs by CIRCOS software [32]．

Calculation of Ka/Ks values

The amino acid sequences from segmentally duplicated pairs and orthologous pairs were first aligned using Clustal X 2.0; thereafter, the aligned sequences were converted to the original cDNA sequences using the PAL2NAL program [33] (http://www.bork.embl.de/pal2nal/). The CODEML program of the PAML package [34] was used to estimate the synonymous (Ks) and nonsynonymous (Ka) substitution rates.

Annotation and analysis of transposable elements

Thede novoprediction and homolog search method based on Repbase [35] were used in the present study to identify the repeat content. For thede novoanalysis, PILER-DF, RepeatModeler, and LTR_FINDER [36,37] were used to predict the transposable elements (TEs) in the genome. For the analysis using the homology-based approach, the known TE library was used; the TEs were identified at the DNA level with RepeatMasker [38] using Repbase TE. To analyze the function of TEs in the expansion of the WOX family, we identified the TEs located 10,000 and 2,000 bp upstream and downstream of the gene and made statistics of the different types of TEs (mutator-like transposable element (MULE), hAT, CACTA, helitron, retrogenes, and retrotransposons) present.

Gene structure analysis

ArabidopsisandG. hirsutumsequences were aligned with Clustal ✕ 2.0, and MEGA 7.0 [30] was used to construct an NJ tree using the method and parameters as described above. The exon positions were acquired from the bed-file, and they were displayed by an online tool, GSDS 2.0 [39]．

Transcriptome data analysis and gene expression heatmap

The raw data of RNA-seq was downloaded from the NCBI Sequence Read Archive (SRA: PRJNA248163). Tophat and cufflinks were used to analyze the RNA-seq expression, and the gene expressions were uniformed in fragments per kilobase million (FPKM) [40]．The expression of WOXs was extracted from the total expression data. Heatmap was drawn by Genesis software [41]．

Real-time PCR

Cotton seeds of TM-1 were obtained from the Institute of Cotton Research of the Chinese Academy of Agricultural Sciences. The cotton (TM-1) seeds were germinated on a wet filter paper for 3 days at 28 °C, and then transferred to a liquid culture medium [42]．3 leaf阶段,幼苗治疗with 10% PEG 6000 and 300 NaCl; the true leaves were sampled at 0, 1, 3, 6, and 12 h of the treatment and were immediately frozen in liquid nitrogen and stored at -80 °C. The total RNA was extracted from the seedlings using RNAprep Pure Plant Kit (TIANGEN, Beijing, China), as per the manufacturer’s instructions. The first strand cDNA was synthesized using a PrimeScript® RT reagent kit (Takara, Dalian, China). SYBR Premix Ex TaqTM II (Takara) was used for PCR amplifications. The cotton histone 3 (GenBank accession no. AF024716) was used as an internal control.

Gene identification and homeobox domain retrieval

We used WOXs fromArabidopsis thalianaas queries for searching the rice, sorghum, poplar, maize, cacao,G. arboreum,G. raimondii, andG. hirsutumdatabases using blastp program and hits with e-values of 1e-5 were considered significant. In preliminary analysis, we recognized 14 candidates in rice, 14 in sorghum, 18 in poplar, 20 in maize, 14 in cacao, 26 inG. arboreum, 31 inG. raimondii, 50 inG. hirsutum(NAU), and 33 inG. hirsutum(BJI). Thereafter, PROSITE (http://prosite.expasy.org/) and InterProscan 56.0 (http://www.ebi.ac.uk/interpro/) were used to search for the HB domain in the obtained sequences and 14, 12, 12, 19, 11, 21, 20, 38, and 33 genes were confirmed as WOX family members in rice, sorghum, poplar, maize, cacao,G. arboreum,G. raimondii,G. hirsutum(NAU), andG. hirsutum(BJI), respectively (Additional file1:Table S1). According to newly sequenced A genome database (Unpublished) by PacBio RS II [43], we found the previous annotations of Cotton_A_11936 and Cotton_A_11937 were not accurate, and the total WOX genes inG. arboreumshould be 20 , so we usedGaWOX1to represent them for further study. Comparing the genes from the twoG. hirsutumgenomes (NAU and BJI), we found that they were highly similar, and the genes from NAU contained all the genes from BJI; therefore, we took genes from NAU for most of the analyses (Additional file2: Table S2). The total number of WOX genes identified in the two diploid species,G. arboreum(AA) andG. raimondii(DD), was higher than that in the tetraploidG. hirsutum(AADD), which is derived from hybridization of progenitors resemblingG. arboreumandG. raimondii。The number of WOXs in the two diploid cotton species (G. arboreumandG. raimondii) was found to be much higher than that in cacao [11], poplar [12],Arabidopsis[16], sorghum [12], and rice [14], indicating that the WOX family in cotton has undergone enlargement during the evolution.

Conserved amino residues within homeobox domains

The WOX gene family is a plant-specific clade of HD-containing superfamily, typically characterized by the presence of a conserved HB domain in the full-length sequence. To investigate the presence of homologous domain sequences and the degree of conservation of each residue in the HB domains, we performed multiple sequence alignment to generate sequence logos of the HB domain inG. hirsutum,Arabidopsis, and rice. The sequence logos revealed that the residue distribution in the HB domains was highly similar in these three plants (Fig.1a-c). Some amino acid residues in the HB domain, for instance, Q, L, and E in helix 1, P, I, and L in helix 2, and I, N, V, F, V, W, F, Q, N, R, and R in helix 3, were highly conserved. In contrast, the amino acid residues in the loop and turn were more variable; for example, the twentieth amino acid was blank due to the insertion of an extra amino acid in some atypical HB domains (Fig.1a-c). Therefore, our results demonstrate that the HB domain sequences from WOXs are highly conserved among the typical dicot and monocots species.

Phylogenetic analysis and nomenclature of WOX genes

To determine the evolutionary relationship of WOX genes among cotton (G. hirsutum,G. arboreum, andG. raimondii) and other species, we constructed a phylogenetic tree by MEGA 7.0 using the NJ method. The cotton WOX genes were named based on the phylogenetic analysis. Ga, Gr, Gh, and At were used as prefixes before the names of WOX genes fromG. arboreum,G. raimondii,G. hirsutum, andArabidopsis, respectively. Moreover, following rules were also considered for the nomenclature: 1) cotton WOXs were named after their orthologs inArabidopsis; 2) if there were no orthologous counterparts inArabidopsis, the cotton WOXs were named based on their homologs in the same clades; and 3) “a” and “b” were appended to the gene names to distinguish the relatively recent paralogs in a particular lineage. The phylogenetic tree showed that WOX genes could be naturally classified into three clades, namely, the ancient clade, WUS clade, and intermediate clade (Fig.2). To validate the phylogenetic tree constructed using the NJ method, we also used the minimum-evolution method to construct a tree. As show in Additional file3: Figure. S1, WOXs were divided into three clades as shown in Fig.2。Although, there were differences between the topologies of the two trees, the member within the subclades and the topology within the subclades were relative stable, which indicated that the NJ tree could be used for further analysis.

Compared to those of WUS and intermediate clades, the WOX genes of ancient clade were probably the earliest to diverge because the ancient clade separated the WUS and intermediate clades in the midpoint root. We found that the number of genes in the WUS clade (95) was greater than the sum of genes in ancient (28) and intermediate clades (45). To further investigate the evolutionary relationship and to predict the gene functions, we divided WOXs into nine sub-clades, named α through ι; six sub-clades were included in the WUS clade, two in the intermediate clade, and one in the ancient clade (Fig.2; Additional file3: Figure S1). Furthermore, we found that all the sub-clades comprised of dicot and monocot species. It is noteworthy that the genes within the sub-clades clustered with a dicot- or monocot-specific pattern. The number of WOXs in each species was variable within the sub-clades. For example, in the sub-clade ε, only one member was present inArabidopsis, sorghum, rice, and cacao, each, but each of the other four species contained two members.

Compared to other species, WOX genes in cotton showed a closer relationship with that in cacao because they always clustered closely to each other in the phylogenetic tree. However, their gene number were not similar within the sub-clades, and in most cases one cacao gene corresponded to two homologous genes ofG. arboreumandG. raimondii, whereas in some sub-clades, one cacao gene only had one corresponding gene inG. arboreumandG. raimondii。For example, in sub-clades θ,Thecc1EG029776had two orthologs inG. arboreumandG. raimondiieach, whereas in sub-clade γ,Thecc1EG012755andThecc1EG010888had only one ortholog inG. arboreumandG. raimondii, respectively.

As shown in Fig.2and Additional file3: Figure S1, almost all the orthologs from A genome and At sub-genome tended to form an orthologous pair at the branch end; same was the case with the orthologs from D genome and D_tsub-genome, indicating that the orthologs from At-A or Dt-D had a more closer relationship.

Gene enlargement and synteny analysis

G. hirsutum, the typical allotetraploid, is an ideal material for studying the effect of naturally occurring polyploidy [26]．To study the locus relationship of orthologs between the At and Dt genomes, we investigated the gene locus on chromosome and performed synteny analysis. The synteny analysis revealed that most of the WOX loci were highly conserved between the At and Dt sub-genomes (Fig.3). We also found that 38 WOX genes inG. hirsutum位于20染色体,distribu吗tion was uneven; for example, there were no WOX genes on chromosomes A04, D04, A06, D06, A09, and D09. Except for A02-D02, A03-D03, A11-D11, and A12-D12, the number of genes located on the chromosome in At sub-genome was the same as that on its homologous chromosome in Dt sub-genome, indicating that some genes might have been lost during the evolution or incomplete sequencing of genome might have resulted in the identification of less number of genes than were actually present; for example,GhWOX1_Atwas located on A12, but no corresponding ortholog was found in the Dt sub-genome. We further checked the genome sequence released by BJI [24]; unfortunately, we identified only one un-oriented orthologous gene (CotAD_21819), which was located on scaffold506.1 (Additional file2: Table S2). Therefore, we took eleven genes in the collinearity block aroundGhWOX1_At(fromGh_A12G2428toGh_A12G2439) to align with the cotton genome and coding sequencing database to determine its collinearity block in the Dt sub-genome. As shown in Additional file4: Figure S2a, we found that there was a corresponding collinearity block on D12 containing 11 genes, among whichGh_D12G2554corresponded toGhWOX1_At。Thereafter, the sequences from the promoter and genomic regions (from 58,423,052 to 58,426,551 bp) were extracted from the genome to align with the A12 chromosome. As shown in Additional file4:图开通中,我们发现有一个删除(998 bps) between the promoter and genomic regions on D12, which resulted inGh_D12G2554losing the typical HB domain of the WOX genes (Additional file4: Figure S2c). We also found thatGhWOX13b_Atwas located on A02, but its corresponding ortholog,GhWOX13b_Dt, was located on D03. To confirm this, we checkedGhWOX13bin another genome (BJI) and found two corresponding genes (CotAD_49496andCotAD_38760); however, both of them were located on the unanchored scaffold (Additional file2: Table S2). Therefore, we used the protein sequence ofGhWOX13bto align with the newly sequenced A genome protein database (Unpublished) by PacBio RS II [43),并发现WOX13b也位于Chr02,相应的是哪一个chromosome of A02 (Additional file5). We found that the collinearity around theGhWOX1 3b_Dtlocus was conserved between the Dt03 and D03 genomes ofG. raimondii; therefore, we presumed that a chromosomal translocation has been accomplished between Chr02 and Chr03 before cotton polyploidization forming an allotetraploid (Fig.3; Additional file6: Figure S3). MCSCAN was used to identify the duplicate gene types. Almost half of the WOX genes were determined to be singletons, whereas a total of 12 genes were observed to have undergone segmental duplication. No tandem and proximal duplications (in the nearby chromosomal region but not in the adjacent region) were found among WOXs. Eight WOX genes were present as dispersed genes (duplications other than segmental, tandem, and proximal) (Additional file7: Table S3).GhWOX4b_AtandGhWOX4b_Dtseem to have been more active during evolution because both of these formed two collinearity gene pairs within At and Dt sub-genome, respectively.

在漫长的进化历史,重复的genes might have experienced functional divergence, including nonfunctionalization (loss of original functions), neofunctionalization (acquisition of novel functions), or subfunctionalization (partition of original functions) [44,45]．To investigate whether Darwinian positive selection was associated with the WOX gene divergence after duplication, the non-synonymous divergence levels (Ka) versus synonymous divergence levels were calculated for 25 homologous pairs. Based on the Ka/Ks ratio, we can presume the selection pressure for the duplicated genes. It is generally accepted that a value of Ka/Ks = 1 indicates that the genes are pseudogenes with neutral selection, Ka/Ks < 1 indicates that the duplicated genes have a tendency to purify, and Ka/Ks > 1 shows an accelerated evolution with positive selection. In our study, we found that the Ka/Ks ratios from 15 gene pairs were smaller than 0.5 and those from seven gene pairs were between 0.5 and 1.0. Only three gene pairs (GhWOX13a_At-GhWOX13a,GhWOX13b_At-GhWOX13b_Dt, andGhWUS1b_At-GhWUS1b_Dt) had Ka/Ks larger than 1, and these gene pairs might have experienced relatively rapid evolution following duplication (Table1). Because most of the Ka/Ks values were smaller than 1.0, we presumed that the Cotton WOX gene family has undergone strong purifying selection pressure with limited functional divergence that occurred after segmental duplications and whole genome duplication (WGD).

The TEs are spread throughout the genome, and many of them are located in the vicinity of host genes [46]．Under abiotic or biotic stress, the TEs can be activated or repressed [47]．To explore whether TEs were involved in the WOX family expansion, we usedde novoprediction and homolog search methods to identify the TEs in the whole genome, and the TEs close to the WOX genes were taken out. As show in Table2, when checking the 2,000-bp region around the gene locus of WOX, only two retroelements (L1 and Caulimovirus) were found (Additional file8: Table S4). We then broadened the scanning region to 10,000 bp upstream and downstream of the genes, respectively, and thirty one TEs were identified including five DNA transposons and 26 retroelements. All the five DNA transposons belonged to CMC-EnSpm family, whereas the retroelements were made up of Caulimovirus (1), Copia (13), and Gypsy (12) (Additional file9: Table S5). Upon further inspection of the distribution of TEs within the 2,000-bp region, only one L1 was identified downstream ofGhWOX4_Dtand one LTR/Copia was found to be located downstream ofGhWOX2a_At。Within the 10,000-bp region around the gene locus, four CMC-EnSpm TEs were found to be located upstream ofGhWOX4a_Dt, one CMC-EnSpm was located downstream ofGhWOX8_At, three Copia elements were located downstream ofGhWOX13a_DtandGhWOX2a_Dt, two Copia TEs were located upstream ofGhWOX4_AtandGhWOX6_At, one Copia was located downstream ofGhWOX3a_At,GhWUS1bandGhWOX3a_Dt, five Gypsy TEs were located upstream and downstream ofGhWOX13a_Dt, three Gypsy TEs were located downstream ofGhWOX4_At, two Gypsy TEs were located downstream ofGhWOX4_At, and one Gypsy element was located downstream ofGhWOX13a_AtandGhWOX3b_At。We noticed that most of the TEs correlated with the presence of duplicated genes, which suggested that TEs, especially the retroelements, played important roles in the WOX family expansion. Compared to TEs, simple repeat sequences are more abundant, most of which locate downstream or upstream of the genes, and only 39 were located in the genomic region. The length of simple repeat sequence was very variable, which might play important roles in the divergence of gene function after duplication.

Analysis of gene structure and homeodomain location

Gene structure is closely related to its function and, together with phylogenetic analysis, it can reflect the phylogenetic relation among the WOX genes. To further study the phylogenetic relationship betweenArabidopsisandG. hirsutum, an NJ tree was generated with MEGA 7.0 usingArabidopsisand cotton WOX protein sequences (Fig.4a), and positions of exons and introns in their genes were determined (Fig.4b). Because the number of genes used for generating the phylogenetic tree described earlier (Fig.2) was different from the one shown in Fig.4a, the topologies of the two trees were different; however, the gene members within the sub-clades were mostly the same. As shown in Fig.4a, WOX genes from cotton were grouped into nine sub-clades;AtWOX15could not be divided into any of the sub-clades. As evident from Fig.4a, WOXs in α, ε, and ζ sub-clades of cotton might have undergone duplication because oneArabidopsisWOX in these sub-clades matched more than two orthologs from cotton (Fig.4a) and this speculation was confirmed through collinearity analysis within At and Dt sub-genomes (Fig.3). The exon number inArabidopsisranged from one to four, with the average being 2.9375, and in cotton, it ranged from two to four, with an average of 2.7297. In general, the gene structures among most of the sub-clades were conserved; for instance, in ζ, δ, and β sub-clades, the members in each sub-clade consisted of two exons. On comparing the exons in the orthologs from the same position inArabidopsisand cotton, we found that the intron lengths were more divergent between the two; for example, the intron length ofAtWOX5was much shorter than that of its ortholog in cotton. In the θ sub-clade, 6-bp-long exons were found at the fourth exon locus ofAtWOX8andAtWOX9but they were not found inGhWOX8_At/DtandGhWOX9_At/Dt, which might be due to an intron insertion in the third exon of the ancestor ofAtWOX8andAtWOX9or these exons might have been lost in cotton during the evolution. We noticed that the intron lengths ofGhWOX5_AtandGhWOX5_Dtwere different. Therefore,GhWOX5genomic sequences were extracted from the chromosome sequences to perform multiple alignments. We found many indels in the genomic sequence and statistically analyzed the indels more than 2 bps in length. The result showed that the introns ofGhWOX5in the At and Dt sub-genomes were less conserved than the exons. This was because of the fact that of the eight indels present in the genomic sequences, seven were located in the introns and only one was present in the exons (Additional file10: Figure S4). The typical conserved domains (HB domain) are marked in the exons with orange color. The lengths of HB domain ranged from 177 to 183 bp, and were considerably conserved. The HB domain locus in the same sub-clade was significantly conserved; for example, in α, γ, θ, and ι sub-clades, the HB domains were located in the second exon, and in δ, ε, ζ, β, and η sub-clades, they were located in the first exon (Fig.4b). OnlyGhWOX10_Dtwas used in our study, and its ortholog in At (Gh_A11G2876) was observed to have lost the third exon and part of the second exon, which resulted in the HB domains being incomplete. Further investigation revealed a point mutation at 1433 bp, turning ‘CAG’ to ‘TAG’, which lead to a premature termination ofGh_A11G2876protein (Additional file11: Figure S5). To validate the above results, we searched the BJI genome; however, no corresponding orthologs were identified. Therefore, we used the RNA-seq aligned data (Bam file) to verify the result. As shown in Additional file12: Figure S6, we found the clean reads mapped to the reference genome very well, and according to the mapping rate, we clearly found that there were two introns inGhWOX10_Dt。However, there was only one intron inGh_A11G2876。We enlarged the mapping results near the mutation site and clearly saw thatGh_A11G2876acquired a termination codon at 92,69,533 bp. The point mutant might have affected theGh_A11G2876gene structure, indicating thatGh_A11G2876might have a function different from that ofGhWOX10_Dtor might have lost its biological function.

Gene expression pattern in different tissues and under multiple stresses

因为与生物相关的基因表达logical function, we inspected the expression patterns of different WOX genes. RNA-seq data were downloaded from NCBI and analyzed. As shown in Fig.5a, we found that WOXs were widely expressed in the vegetative (root, stem, and leaf) and reproductive (torus, petal, stamen, pistil, calycle, and -3, -1, 0, 1, 3, 5, 10, 20, 25 and 35 days post-anthesis (DPA) ovule) tissues as well as in the fiber (5, 10, 20, and 25 DPA), indicating that WOXs have diverse biological functions and work in different tissues. We found that some WOXs did not express in the vegetative tissues and had a very low expression levels in the reproductive tissues. For instance, we could not detectGhWOX2_At/Dt,GhWOX3a_At/Dt, andGhWOX3b_At/Dt expression in root, stem, and leaf, and very low expression levels were detected in some of the reproductive tissues. On comparing the expression patterns of the orthologs between At and Dt, we found that the expression patterns and the levels of expression of the two were not always the same; for example,GhWOX8_Dtwas expressed in root and leaf butGhWOX8_Atwas not.GhWOX8_Dthad higher expression levels in 20, 25, and 35 DPA ovule and 10 DPA fiber compared to that inGhWOX8_At. GhWOX13a_At/DtandGhWOX13b_At/Dtwere not only close to each other in the phylogenetic tree (Fig.2and Fig.4a), but also had similar expression patterns (Fig.5), suggesting that they have a similar biological function.

Cotton faces multiple abiotic stresses during its growth and development. Therefore, a comprehensive analysis of the expression pattern of WOXs was performed in the present study. No obvious changes in the expression levels were observed for more than half of the WOXs under hot, cold, salt, and PEG 6000 conditions.GhWOX10_Dtresponded to multiple stress treatments. The expression ofGhWOX13a_At/DtandGhWOX13b_At/Dtwas strongly induced by multiple stresses, indicating that they might take part in response to stress and their expression might be regulated by stress.

Examination of the expression of WOX genes by qPCR

To verify the expression profiles of the WOX genes obtained by RNA-seq data, qPCR was performed using the leaves from plants treated with PEG6000 and NaCl. Eight genes, including three orthologous pairs, that were presumed to be highly expressed under PEG6000 or NaCl treatment based on the RNA-seq data, were selected for qPCR, (Fig.5b; Fig.6). It was difficult to distinguish between the orthologous pairs by qPCR because the sequences of the orthologous genes (GhWOX13a_At-GhWOX13a_Dt,GhWOX13b_At-GhWOX13b_Dt, andGhWOX4a_At-GhWOX4a_Dt) were highly similar. Therefore, we designed primers to amplify the orthologous pairs together (Additional file13: Table S6). The pattern of gene expression assessed by qPCR showed a similar tendency with that detected using the RNA-seq data, as shown in Fig.5band Fig.6。The expression ofGhWOX1_Atwas down-regulated under PEG 6000 and NaCl treatment. In contrast, the expression levels ofGhWOX10_Dt,GhWOX13b_At/DtandGhWOX13a_At/Dtwere increased under both PEG 6000 and NaCl treatment, indicating that they might play vital roles in the stress response and could be the candidate genes for further study on cotton stress biology.

Previous analysis of the WOX gene family has been performed in rice, sorghum, maize,Arabidopsis, and poplar [2]．However, there were no reports on the analysis of this family in cotton. In this study, we performed a comprehensive identification of WOX genes inG. hirsutum,G. arboreum, andG. raimondii, mainly focusing on the allotetraploid cottonG. hirsutum, with the aim of understanding the roles of this gene family in cotton, in future studies.

Cotton WOX family underwent enlargement during the evolution

Upland cotton is one of the most important cash crops worldwide, providing more than 90% of the natural renew fiber for textile industry. It is also an ideal material for studying the effects of natural polyploidy [48]．A-genome and D-genome diploid cotton were native to Africa and Mexico, respectively, and they diverged about 5–10 million years ago. About 1–2 million years ago, A-genome cotton resemblingG. arboreumandG. herbaceumhybridized with D-genome cotton resemblingG. raimondii, which was followed by chromosome doubling, and the eventual formation of nascent AtDt allopolyploid, including the upland cotton [22,48]．

A total of 38 WOX genes were identified in the upland cotton, which is more than the number of this gene reported in most of the other species whose genome has been sequenced [4]．One of the most important features was that the upland cotton experienced polyploidization, resulting in WGD. The At and Dt donors of the upland cotton are close relatives, with almost equal number of orthologs; therefore, the nascent duplication resulted in a number of WOX doublings [49]．Polyploidization is an important event during the evolution of flowering plants that might play important roles in the adaptation of plants to new environment [50]．Although the total number of upland cotton dramatically increased, gene loss also happened during the evolution of upland cotton, as evident from the comparison of the number of genes in At or Dt genomes with that in A- (G. arboreum) or D-genome (G. raimondii). Gene loss always happened during the rapid arrangement of genomic sequences after hybridization and following chromosome doubling during polyploidization [51]．Compared to the paleopolyploid maize andBrassica, cotton displays fewer changes in their genomic sequences [52,53]．

Although polyploidy was the main contributor in duplication, segmental duplication has also been responsible for the increase in the WOX gene number. Segmental duplication is one of the most important impetuses for evolution; it occurs most frequently in plants since most plants are diploidized polyploids and retain large amount of duplicated chromosomal blocks within their genomes [54]．The WRKY genes in soybean underwent large scale segmental duplication, with the duplicated genes exhibiting differential expression and functional divergence [55]．热休克蛋白在res扮演重要的角色ponse to drought stress in sesame; evolutionary analysis showed that segmental duplication was the primary force promoting the expansion of heat shock promoting genes [56]．In our study, 12 out of 38 genes were associated with segmental duplication, indicating that segmental duplication played significant roles in the expansion of WOX genes, and differential expression suggested that the duplicated genes might have experienced functional divergence.

Although we did not find tandem duplication in our study, it is also a basic contributor to gene expansion, which arises from unequal crossing over, chromosomal anomalies, transposon insertions, and other reverse-transcriptase-mediated processes [57]．After duplication, the new gene would be redundant with the previously existing one. The redundancy has been considered a driving force for evolutionary innovation [57]．一些模型被提出d for understanding the gene duplication; these include the neofunctionalization, DDC-subfunctionalization, escape from adaptive conflict (EAC), and dosage-balance models, which provide a theoretical framework for further studying the process [57]．

Cotton WOXs have been highly conserved during the evolution

WOX is a plant-specific sub-clade with few members of HD-containing superfamily, which consists of a leucine zipper (HD-ZIP), plant homeodomain associated to a finger domain (PHD finger), the distinctive Bell domain (Bell), zinc finger associated to a homeodomain (ZF-HD), knotted-related homeobox (KNOX), and WUSCHEL-related homeobox (WOX) [58]．Previous studies have shown that HDs (HB-domain) might have diverged before the separation of the branches forming plants, animals, and fungi [59]．The members of HD-containing superfamily were different not only in the sequence of HB domain, but also in their size, HB location, other related domains as well as in their structures [58]．WOX predominantly has HB domain, and can be distinguished by the phylogenetic relatedness of its homeodomains [1]．It has been reported that HB domain contains some conserved amino acids located in the helices; for example, it contains Q, L, and Y in helix 1 and I, V, W, F, N, K, R, and R in helix 3 [2]．Our results demonstrate that the above amino acids were also conserved in cotton. Previous studies have indicated that the amino acids in loops were less conserved [2], which is in consonance with our results. Previous studies have also shown that WOXs could be divided into three clades [1,2,60], as was observed in the present study for cotton (Fig.2). In a study by Zhang et al., WOX family was further classified into nine sub-clades, as was the case in the present study [2]．We observed thatAtWOX6was not divided into any sub-clades (Fig.2) because it had a bootstrap value lower than 50. We found that in the study by Zhang et al.,AtWOX6was clustered into the same clade asAtWOX1despite the bootstrap value not supporting this classification. However, as shown in Fig.4, we found thatAtWOX6andAtWOX1might have originated from the same ancestor, a fact supported by a high bootstrap value. Therefore, we presumed that the increasing number of sequences in Fig.2would make the total number of sequences more divergent, resulting inAtWOX6being not stable between Fig.2and Fig.4。Because the annotation is gradually improving,AtWOX15(AT5G46010.1) was identified as a novel member that was not reported in previous studies [1,2]．We found that it belonged to the intermediate clade, but it could not be divided into any sub-clade because as per our definition in the present study, a new sub-clade should have at least two members. The monocot- and dicot-specific cluster pattern suggested that the main functions of WOX genes had been determined before the monocots and dicots split. Previous studies have shown that in ζ (NS/WOX3) sub-clade, the genes from monocots formed two separate branches [2], which was consistent with our results. However, we found that there were two pairs of duplicated genes (GhWOX3a_At-GhWOX3b_AtandGhWOX3a_Dt-GhWOX3b_Dt) in cotton, with no orthologs in poplar, which indicated that poplar might have lost these genes during evolution or this might have resulted from an incomplete sequencing result. In the report of Zhang et al., it was thought that the subgroup B (corresponding to sub-clade α in our study) might be an ancient subgroup, because only one member from each species was present in this clade [2]．However, our data did not support this opinion because we found thatGhWOX4_At/Dt,GhWOX4a_At/Dt, andGhWOX4b_At/Dtwere duplicated genes (Fig.3). Cacao and cotton diverged from a common ancestor about 18–58 million years ago, and both of them had undergone the ancient duplication event together; cotton had subsequently experienced a nascent duplication event again [25], our data also indicated that cacao and cotton were close relative and probably derived from the same ancestor. One cacao WOX gene should theoretically correspond to two orthologs in the diploid cotton, but in fact some genes inG. arboreumandG. raimondiimight have been lost during the evolution [51]．It is generally believed thatG. hirsutum(2n = 4x = 52, (AD)_t) was reunited by hybridization of an A-genome species resemblingG. arboreum(2n = 26, A₂) with a D-genome material resemblingG. raimondii(2n = 26, D₅), followed by chromosome doubling [24]．Our data suggested that WOXs in A genome and At sub-genome had a common ancestor, and those in D genome and Dt sub-genome had a common ancestor, which was consistent with the above hypothesis, indicating that cotton WOXs were highly conserved during the evolution.

Expression of duplicated WOX genes

Previous studies on the expression and function of WOX genes have indicated that the WOX family members play crucial roles in key developmental processes of plants, such as in embryonic patterning, stem-cell maintenance, and organ formation [1]．Our data (Fig.5a) suggest that most of the WOX genes had very low expression levels or they did not express in selected tissues, the reason for which might be that WOX family members mainly function in the process of embryogenesis because of which their expression is restricted. Despite the expression levels and the number of expressed genes being different among the different tissues (Fig.5a), except forGhWOX8_At/Dt, the orthologous gene pairs did not show A- or D-ortholog bias, which suggests that they may have conserved functions. Segmental duplication was one of the most important impetuses for increasing the diversity at the molecular level. After duplication, the coding regions could have acquired new regulatory context through the acquisition/deletion of tissue-specific enhancers and repressors, causing the spatial and/or temporal change in the expression pattern of the duplicated gene, promoting diversification of gene functions, like subfunctionalization (acquisition of part of the function of a pre-existing gene) or neofunctionalization (acquisition of a new function) [61]．AtWOX4works in procambium development, which is associated with vascular patterning and leaf complexity [62]．In upland cotton, we observed the expansion ofWOX4through WGD and segmental duplication, and the duplicated genes showed different expression levels in root and stem, which indicated that these genes might have acquired subfunctionalization or a novel function [57]．

InArabidopsis thaliana, WOX8 protein not only positively regulates early embryonic growth, but also interacts withAtCLE8to promote seed growth and the overall seed size [63]．In the present study,GhWOX8_Dt, the ortholog ofAtWOX8in upland cotton, was observed to express abundantly during the seed development, indicating thatGhWOX8_Dtmight have similar function in seed development asAtWOX8has inArabidopsis。However,GhWOX8_Atmight have a divergent function in cotton seed development, because it had a different expression pattern in ovule compared to that ofGhWOX8_Dt, suggesting thatGhWOX8smight have experienced functional divergence after duplication andGhWOX8_Dtcould be a key candidate gene in cotton seed development that should be studied further.

AtWOX11interacts withWOX12; these two genes are involved in the first-step of cell fate transition during thede novoroot organogenesis [18]．GhWOX12was expressed specifically in root, which suggested thatGhWOX12might play important roles in root development; however, confirmation as to whether it works together withGhWOX11in regulating thede novoroot organogenesis would need future investigation.GhWOX11andGhWOX12had similar expression patterns asGhWOX8_Dthad in ovule; however, whether these genes have similar functions in seed development would need to be ascertained in further studies.

AtWOX13is a key component in the regulation of floral transition and root development, and has a high expression level during primary and lateral root initiation and development, in gynoecium, and during embryo development [64]．We found that its orthologs in cotton,GhWOX13a_At/DtandGhWOX13b_At/Dt, were dramatically expressed in the root, in reproductive organs, as well as during the embryo development, indicating thatWOX13 genes in cotton andArabidopsishave a conserved expression pattern, and might have a similar biological function. We also noticed that the other homologs ofWOX13in cotton from the same sub-clade also expressed in all the assessed tissues, which suggested that the genes in the ancient clades might be more active during cotton growth and development. Although most of the reports on WOXs determined their roles in embryogenesis, they might play important roles in stress response. The mutanthos9-1, which was generated by a T-DNA insertion into the fourth exon ofAtWOX6, was reported to be hypersensitive to freezing treatment [13]．Three WOXs from paper mulberry were induced by cold exposure [65]．A previous study showed thatOs08g14400was up-regulated under desiccation and salt stress,Os03g20910was up-regulated under salt stress, andOs01g63510was down-regulated under salt stress [61]．In our study, some of the WOX genes were also induced by different stress, indicating that they may meditate the stress response. We found that the duplicated genesGhWOX4_At/Dt,GhWOX4a_At/Dt, andGhWOX4b_At/Dt responded differently when exposed to different stresses, suggesting that they have different functions under different stresses.GhWOX13b_At/Dtwas induced by multiple stresses, indicating thatGhWOX13b_At/Dtmight be the node for multiple stress regulation. Although WOX genes were shown to have different expression levels under multiple stresses, there are no reports on the validation of the function of cotton WOXs in stress. Therefore, future studies should focus on determining the function of WOX genes in stress response.

Previous studies have illustrated that members of the WOX gene family play significant roles in the regulation of plant growth and development by determining the cell fate. The results of the present study indicate that WOX genes are highly conserved among cotton and other plant species. Furthermore, whole genome and segmental duplication have probably been the two major ways for gene amplification during the expansion of the WOX family in upland cotton. Moreover, the duplicated genes in cotton seem to have experienced functional divergence because the duplicated genes showed different expression patterns in different tissues and organs. In addition, some WOX gene members are likely to be involved in the mediation of stress response. Our results will not only deepen the understanding of the evolutionary processes of cotton WOX genes, but would also be helpful in formulating further functional genomic studies of WOX gene family in cotton.

BLAST:

Basic Local Alignment Search Tool

CUC:

Cupshaped cotyledon

DPA:

Days post-anthesis

EAC:

Escape from adaptive conflict

FPKM:

Fragments per kilobase million

G. arboreum:

Gossypium arboreum

G. hirsutum:

Gossypium hirsutum

G. raimondii:

Gossypium raimondii

HB:

Homeobox

HD:

Homeodomain

NJ:

Neighbor-joining

PXY:

Phloem intercalated with xylem

TE:

Transposable elements

WGD:

Whole genome duplication

WOX:

WUSCHEL-related homeobox

Not Applicable.

Funding

This work was supported by The Major Program of Joint Funds (Sinkiang) of the National Natural Science Foundation of China (Grant No. U1303282) and a grant from Henan Province (No. 162300410160). The funding body supported the funds for research and was not involved in the design, data collection, and analysis, decision to publish, or preparation of the manuscript.

Availability of data and materials

The phylogenetic trees were deposited in treebase (http://treebase.org) under following URL:http://purl.org/phylo/treebase/phylows/study/TB2:S21071。All another data generated or analyzed during this study are included in this published article and its Additional files.

Author notes

Affiliations

1. Xinjiang Research Base, State Key Laboratory of Cotton Biology, Xinjiang Agricultural University, Urumqi, 830052, China

   Zhaoen Yang, Wenqiang Qin & Fuguang Li

2. Institute of Cotton Research, Chinese Academy of Agricultural Sciences, Anyang, 455000, China

   Zhaoen Yang, Qian Gong, Wenqiang Qin, Zuoren Yang, Yuan Cheng, Lili Lu, Xiaoyang Ge, Chaojun Zhang, Zhixia Wu & Fuguang Li

Contributions

FL and Z-E.Y conceived and designed the study; XG, YC. LL and Z-R.Y carried out experiments; Z-E.Y, QG and WQ prepared Figures; YC, LL and Z-R.Y analyzed and interpreted the data; Z-E.Y, QG and WQ prepared the manuscript; XG, CZ and ZW participated in the design of the experiments, wrote part of the manuscript and critical review for intellectual content of it. All the authors have read, edited, and approved the current version of the manuscript.

Corresponding author

Correspondence toFuguang Li。

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and Permissions