The genome sequences of J. curcas were retrieved from the Jatropha Genome Database (http://www.kazusa.or.jp/jatropha/; JAT_r4.5; Hirakawa et al., 2012). The genome sequences of H. brasiliensis (Rahman et al., 2013), M. esculenta (Prochnik et al., 2012), and R. communis (Chan et al., 2010) were obtained from the National Center for Biotechnology Information (NCBI, https://www.ncbi.nlm.nih.gov/).

Redundant sequences were removed from the results obtained using the following two methods described below to identify WOX genes in four Euphorbiaceae species. (1) The hidden Markov model (HMM) file (PF00046) associated with WOX protein family-related domains was downloaded from the Pfam database (http://pfam.xfam.org/) (Cao et al., 2017). Using HMMER 3.0 software to analyze WOX proteins, the HMM file (PF00046) of four species of Euphorbiaceae was screened with an E-value cut off of 0.001. (2) 15 WOX protein sequences from A. thaliana were used as queries for alignment with the total protein sequences of the four Euphorbiaceae species by BLASTP with an E-value of 0.0001 to confirm the accuracy. The candidate sequences of WOX genes were detected using SMART (http://smart.embl-heidelberg.de/) to verify the presence of homeodomains. Ultimately, each WOX gene was assigned a unique name by BLASTP.

The properties of WOX proteins, including their molecular weight (MW), isoelectric point (pI), instability index (II), aliphatic index (AI), were computed with ExPASy (https://web.expasy.org/protparam/; Artimo et al., 2012).

All sequences from the four Euphorbiaceae species were subsequently aligned by DNAMAN software (version 6.0) to visualize the results. Moreover, multiple sequence alignment of all WOX protein sequences was performed by ClustalX 2.0 (Thompson et al., 1997). An interspecific phylogenetic tree containing the species was generated using the Neighbor-Joining (NJ) method with MEGA X software, including 1,000 bootstrap replicates (Kumar et al., 2018). In addition, the exon and intron structures of all WOX genes were obtained from the online Gene Structure Display Server (GSDS; http://gsds.cbi.pku.edu.cn; Hu et al., 2015).

The conserved motifs of WOX proteins from the four Euphorbiaceae species were analyzed by MEME online server (http://meme-suite.org/). The maximum number of identified motifs was set to 25 (E-value = 0.0001; Li et al., 2018).

The codon usage bias of WOX genes in the four Euphorbiaceae species was analyzed using CodonW software (version 1.4.2). The relative synonymous codon usage (RSCU) of WOX genes in these species was calculated in such a way that RSCU value >1 represents a codon with positive bias; an RSCU value = 1 indicates no bias (Wang Z. et al., 2021). The relative frequencies of synonymous codons (RFSCs) in the J. curcas WOX genes (JcWOXs) were also calculated. When the RFSC value exceeds 60% or is 0.50 times higher than the average frequency of synonymous codons, the codon is a high-frequency codon (Zhou et al., 2007a). Furthermore, a comparative analysis of the codon usage frequency among WOX genes of the four Euphorbiaceae species and four representative species, A. thaliana, Nicotiana tabacum (N. tabacum), Populus trichocarpa (P. trichocarpa), and O. sativa, was performed. The ratio ranged from 0.50 to 2.00 reflects that the codon usage bias of WOX genes is highly similar to that of the representative plant species (Zhou et al., 2007b).

The genomic sequences located 2000 bp upstream from the initiation codon (ATG) were considered as the promoter fragments. The promoter sequences of JcWOXs were downloaded from the Jatropha Genome Database (Hirakawa et al., 2012). Meanwhile, PlantCARE online sever (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/) was used to analyze potential cis-acting elements (Lescot et al., 2002).

J. curcas RNA transcriptome data from three different tissues (leaves, roots, and seeds) were downloaded from the Sequence Read Archive (SRA) database to investigate the spatial expression characteristics of JcWOXs (Wu et al., 2015; Zou et al., 2016). The spatial expression profiles of each identified JcWOX in three different tissues were determined using three transcriptome datasets. The spatial characteristics of 10 JcWOXs were explored using the R library heatmap, and the trends were ultimately presented as color changes. Detailed information of the three transcriptome datasets was provided in Supplementary Table S1.

The RNA transcriptome data collected at seven developmental stages (14, 19, 25, 29, 35, 41, and 45 days after pollination-DAP) of seeds were downloaded following the spatial expression profiles of JcWOXs (Jiang et al., 2012). The temporal characteristics of the 10 JcWOXs were explored by generating the heatmap. The detailed data information of the seven transcriptomes is presented in Supplementary Table S1.

J. curcas calli were selected as experimental materials to test whether WOX gene has functional conservation in regulating callus proliferation. The third leaf from the top of Guangxi J. curcas was selected as the explant. Next, the explants were disinfected by washing with water for 0.50–1.00 min, followed by soaking in 70% ethanol for 15 s and 3% NaClO solution for 12–16 min. Callus induction was performed by inoculating the explant in dorsal contact medium supplemented with MS, BAP (+0.80 mg/L), TDZ (+0.60 mg/L) and NAA (+0.10 mg/L). Explants were cultured in the dark at 25°C. The whole cycle of callus culture was approximately 42 days, which was divided into three stages: 14 days (S1), 28 days (S2), and 42 days (S3), as illustrated in Supplementary Figure S1. With obvious differences in the growth states and callus characteristics were selected as samples. Three biological replicates of each sample were analyzed. All samples were collected and cryo-preserved in liquid nitrogen.

Total RNA was extracted from J. curcas callus was extracted using the FastPure plant total RNA isolation kit (RC401). The RNA concentration was determined using a Nanodrop-2000 spectrophotometer (Thermo, Inc.). The integrity of total RNA was detected by performing gel electrophoresis. First-strand cDNAs were synthesized using the reverse transcriptase method with the HiScriptR III 1st Strand cDNA Synthesis Kit (+ gDNA wiper) (R312-01/02). Additionally, the samples in triplicate were analyzed with qRT-PCR utilizing AceQ qPCR SYBR Green Master Mix (without ROX), as previously described by Wang D. et al. (2020). Four reference genes (JcGAPDH, JcEF1α, JcActin, and JcTUB8) were selected as candidate reference genes according to Zhang et al. (2013), and JcGAPDH and JcActin were determined as reference genes by using semi-quantitative RT-PCR (Liu et al., 2010). Based on the coding DNA sequences (CDSs) of WOX gene from J. curcas, primers were designed using Oligo7 and SnapGene (Supplementary Table S2).

Both HMMER and local BLAST searches were performed simultaneously with the HMM file (PF00046), and the sequences of 15 members of the WOX family in A. thaliana were used as templates to identify all possible WOX proteins in the four Euphorbiaceae species and confirm the accuracy of the identification. Ultimately, 59 WOX genes were identified in the genomes of the four Euphorbiaceae species such as, 20 in H. brasiliensis, 10 in J. curcas, 18 in M. esculenta, and 11 in R. communis. These WOX genes were named by BLASTP (Supplementary Table S3).

Next, the candidate WOX gene sequences were isolated, followed by analyzing the physicochemical properties of their corresponding proteins utilizing by ExPASy. As shown in Supplementary Table S3, their protein lengths and predicted MW varied little, and few differences in their pI were observed. For H. brasiliensis, the WOX proteins were 185–398 residues long, corresponding MW ranged from 21,346.99 to 44,058.10 Da, and their pI values ranged from 5.63 (HbWOX11b) to 9.51 (HbWOX4b). For J. curcas, the corresponding WOX proteins ranged in length from 190 (JcWOX7) to 392 (JcWOX1) aa, their MW ranged from 21,747.21 to 43,900.70 Da, and their pI values ranged from 5.15 (JcWOX14) to 9.51 (JcWOX4). WOX proteins in M. esculenta were 182 to 391 residues long, with the pI values ranging from 5.61 (MeWUSb) to 9.40 (MeWOX4a). For R. communis, the WOX proteins ranged in length from 192 to 401 residues long, and their pI values ranged from 5.35 (RcWOX14) to 9.42 (RcWOX1a). Therefore, the physicochemical properties of the WOX proteins in the four Euphorbiaceae species were similar. The hydropathicity index represents the hydrophilicity of a protein (Li et al., 2016). The AI ranged from 44.34 to 82.58, indicating that WOX proteins are thermally stable. As shown in Supplementary Table S3, the hydrophilicity index of 59 members of the WOX protein families was negative, suggesting that the members of the four WOX protein families of the Euphorbiaceae species were hydrophilic.

According to the visualized of sequence alignment (Supplementary Figure S2), the homeodomain of WOX protein of A. thaliana and the four Euphorbiaceae species presented high conserved structures. The number of amino acid residues in WOX homeodomain fluctuated slightly, from 63 to 66. These results proved that the WOX proteins of four Euphorbiaceae species were highly conserved, further indicating that the identified WOX genes were correctly identified.

In order to compare the evolutionary relationships of WOX genes of A. thaliana and the four Euphorbiaceae species, NJ phylogenetic trees of 74 WOX proteins from A. thaliana (15), H. brasiliensis (20), J. curcas (10), M. esculenta (18), and R. communis (11) were constructed with MEGA X software to compare the evolutionary relationships of WOX genes from A. thaliana and the four Euphorbiaceae species. It was observed from Figure 1 that WOX genes from the five species were classified into three major clades, which is consistent with the previously known clades. For instance, ancient clade (WOX10, WOX13, and WOX14 subfamilies), intermediate clade (WOX8, WOX9, and WOX11), and modern/WUS clade contains eight subfamilies (WUS and WOX1-7). The number of WOX genes in the modern/WUS clade (38 genes, 64.4%) was greater than that in the intermediate (12 genes, 20.3%) and the ancient clades (nine genes, 15.3%). Moreover, phylogenetic tree analysis reflected that the homologous genes lies in close proximity, suggesting that the WOX genes in Euphorbiaceae species were highly evolutionarily conserved. Furthermore, members of the WOX transcription factors in H. brasiliensis and M. esculenta were located close to each other in the evolutionary tree.

Structure analysis of 59 WOX genes was performed with GSDS. WOX gene is composed of exons, introns, and untranslated regions (UTRs), as illustrated in Figure 2. The number of exons in each WOX gene ranged from two to four. The WOX genes members in the same subclade usually presented similar exon–intron patterns. For example, the members of the closely related modern clade i.e., HbWOX7, JcWOX7, HbWOX5, MeWOX5, MeWOX7, and RcWOX7 exhibit similar gene lengths with two exons. In addition, members of the ancient clade had two intron insertion sites, while members of the intermediate and modern/WUS clades exhibited different intron insertion patterns, and the number of introns ranged from one to four. Moreover, five genes (HbWOX9, MeWOX9a, MeWOX9b, HbWOX13b, and MeWOX13) were more than 3 kb. Interestingly, most JcWOX and RcWOX members were located in adjacent regions and had similar gene structures, showing that the similarity between J. curcas and R. communis species is high and the kinship is close. The same phenomenon was also found in most of the HbWOX and MeWOX members. Collectively, these results indicated that similar gene structures also reflect the conservation of WOX gene members and their evolutionary relationships.

The WOX family typically contains additional conserved motifs that likely have involved in different functions. 25 conserved motifs were identified among the 59 WOX proteins in the four Euphorbiaceae species using the MEME online portal to obtain a comprehensive understanding of the structural features and relationships of the WOX proteins. The phylogenetic tree (Figure 3) was divided into ancient, intermediate, and modern/WUS clades. Among the 25 putative motifs, motif 1 and motif 2 were present in all WOX protein sequences, indicating that motifs 1 and 2 are characteristic domains of the proteins encoded by 59 WOX genes. Structural analysis demonstrated that proteins corresponding HbWOX3a and MeWOX3 exhibit only two motifs while HbWOX1d and MeWOX1c possess 10 motifs. With the exception of JcWOX3, motif 5 was present at a high frequency in nearly all WOX proteins in the modern/WUS clade, but this motif was not present in the intermediate clade or ancient clade, consistent with the distribution of the WUS box (TLXLFPXX). In particular, motif 4 existed only in the ancient clade. Furthermore, the conserved motifs of JcWOX and RcWOX were highly similar, as observed for HbWOX and MeWOX. The WOX members within the same clade had similar motif structures, suggesting that the homologous WOX genes in different plant species are closely related. These results are consistent with previously reported phylogenetic tree analysis and could further strengthens the classification of WOX subfamily members.

CodonW software was used to study the codon usage bias in the 59 WOX genes from the four Euphorbiaceae species. As shown in Table 1, the codon usage of WOX genes across the Euphorbiaceae species was quite conserved. According to the RSCU values, strong commonalities were detected in WOX genes from the four Euphorbiaceae species. 19, 23, and two identical codons displayed positive bias, negative bias and no bias, respectively; the remaining 20 codons differed in bias. The strongest positive bias in all species presents same codons (AGA), however, the codons with the strongest negative bias were different such as J. curcas (CGC) and M. esculenta (CGC), H. brasiliensis (GCG) and R. communis (ACG). An analysis of the RFSCs of the JcWOXs revealed four identical high-frequency codons from WOX genes of the four Euphorbiaceae species: GCA, GAU, CCA, and AGA. However, some differences in the high-frequency codons were observed. For instance, CUU, CAA, and ACU were high-frequency codons in M. esculenta but not in the other three species. Additionally, the codon usage frequency ratios between the four Euphorbiaceae species and four model organisms were analyzed (Table 2). The ratio of six codons (GCU, CGC, UCG, UAA, UAG, and UGA; GCU, GAU, CGC, UAA, UAG, and UGA) were greater than 2.00 and lower than 0.50 in the JcWOX genes than in the A. thaliana and P. trichocarpa genes, respectively. Comparative analysis reflected, five similar codons (GCC, CGG, UAA, UAG, and UGA) in RcWOX genes and P. trichocarpa genes, five codons (GCC, CGC, UAA, UAG, and UGA) showed differences in the comparison between the MeWOX genes and N. tabacum genes, and six codons (GCG, UUG, CGG, UAA, UAG, and UGA) showed differences in the comparison between the HbWOX genes and P. trichocarpa genes. Hence, we considered P. trichocarpa to be the optimal choice for comparing exogenous WOX genes in J. curcas, R. communis, and H. brasiliensis and N. tabacum the optimal choice for comparing exogenous WOX genes in M. esculenta.

As shown in Figure 4, the cis-acting elements identified in this study were divided into three types based on their functions such as stress responses, hormone responses, and growth and development. The largest group of cis-acting elements was related to stress at 104 among the JcWOXs. For most of JcWOXs, the number of identical cis-acting elements usually ranged from one to three. A large number of cis-acting elements were associated with growth and development—10. The most abundant cis-acting elements were MYCs (related to abscisic acid (ABA) induction), of which 26 existed in the promoters of nine JcWOX genes, excluding JcWOX6. Moreover, AREs (related to the anaerobic stress response) were the most widespread cis-acting elements in the JcWOXs (present within all 10 JcWOX gene members), suggesting that JcWOXs are widely involved in reactions related to anaerobic stress. In addition, JcWOX14 contained the most stress-related cis-acting elements, four MYB elements (related to the drought stress response) and four LTR elements (JcWOXs related to the low-temperature stress responses), which reflected that JcWOX14 plays an indispensable role in the drought stress response and low-temperature response. With respect to hormone-related promoter elements, JcWOX1 contains the largest number six types and 11 elements—suggesting that JcWOX1 is likely involved in hormone-related responses. ABREs (related to the ABA reaction element) were the most widespread cis-acting elements in the JcWOXs (present within eight JcWOX gene members, excluding JcWOX2 and JcWOX7), which demonstrates that JcWOXs are widely involved in reactions related to ABA-related reactions. In terms of cis-acting elements related to growth and development, JcWOX1 and JcWOX4 contained the most—seven types and 10 elements. In addition, the G-box (related to the light response element) was present within seven JcWOX gene members, excluding JcWOX2, JcWOX7, and JcWOX11. Interestingly, JcWOX13 contains four GT1 motifs that participate in the light response.

We investigated previously published transcriptome profiles of various tissue types to determine the critical role of WOX genes in the development of J. curcas. Based on a heatmap of these data, we found that four genes (JcWOX4, JcWOX7, JcWOX9, and JcWOX11) in the roots and five genes (JcWOX1, JcWOX2, JcWOX6, JcWOX13, and JcWOX14) in the seeds. Nevertheless, WOX genes in J. curcas were expressed at low levels in the leaves (Figure 5). Based on these results, the JcWOXs might play more significant roles in the growth and development of roots and seeds than in leaves. In summary, the expression patterns of JcWOXs in different tissues respectively are conducive to identifying functional genes expressed in J. curcas.

Based on the analysis of the expression of WOX genes and their variation trends across seven J. curcas seed development stages (Figure 6), all JcWOXs were expressed at relatively high level during the initial stage (19 DAP, 25 DAP, and 29 DAP) and at relatively low level during the fast oil accumulation stage (35 DAP, 41 DAP, and 45 DAP). The expression of four WOX genes (JcWOX3, JcWOX6, JcWOX11, and JcWOX13) markedly decreased between 29 and 35 DAP. Interestingly, the expression of two WOX genes (JcWOX6 and JcWOX13) decreased continuously during all seven stages. In general, the expression of WOX genes tended to decrease during the seven stages, hence, we deduced that JcWOX6 and JcWOX13 play a significant anti-regulatory role in seed development.

The expression levels of JcWOXs were analyzed in calli collected at three stages (S1, S2, and S3) via qRT-PCR to better understand their functions. We performed qRT-PCR on WOX genes of J. curcas; the findings for some genes with polar expression levels were discarded (JcWOX1, JcWOX2, JcWOX3, JcWOX6, and JcWOX9). The results of qRT-PCR experiments (Figure 7) showed that among the five WOX genes of J. curcas, the expression of four WOX genes (JcWOX4, JcWOX7, JcWOX13, and JcWOX14) did not display obvious variation in the three biological repetitions. The expression of the JcWOX11 gene was significantly upregulated at S1. Moreover, the gene variation in the expression of JcWOX11 was most obvious. These findings provide new insights and a comprehensive understanding of the characteristics of JcWOXs for functional validation in the future.