Genome-Wide Identification and Expression Profiling of the WOX Gene Family in Citrus sinensis and Functional Analysis of a CsWUS Member

WUSCHEL-related homeobox (WOX) transcription factors (TFs) are well known for their role in plant development but are rarely studied in citrus. In this study, we identified 11 putative genes from the sweet orange genome and divided the citrus WOX genes into three clades (modern/WUSCHEL(WUS), intermediate, and ancient). Subsequently, we performed syntenic relationship, intron-exon organization, motif composition, and cis-element analysis. Co-expression analysis based on RNA-seq and tissue-specific expression patterns revealed that CsWOX gene expression has multiple intrinsic functions. CsWUS homolog of AtWUS functions as a transcriptional activator and binds to specific DNA. Overexpression of CsWUS in tobacco revealed dramatic phenotypic changes, including malformed leaves and reduced gynoecia with no seed development. Silencing of CsWUS in lemon using the virus-induced gene silencing (VIGS) system implied the involvement of CsWUS in cells of the plant stem. In addition, CsWUS was found to interact with CsCYCD3, an ortholog in Arabidopsis (AtCYCD3,1). Yeast one-hybrid screening and dual luciferase activity revealed that two TFs (CsRAP2.12 and CsHB22) bind to the promoter of CsWUS and regulate its expression. Altogether, these results extend our knowledge of the WOX gene family along with CsWUS function and provide valuable findings for future study on development regulation and comprehensive data of WOX members in citrus.

Several previous studies indicated that the WUS gene is one of the most important genes in the WOX gene family and involved in numerous important developmental processes including size of shoot meristem, somatic embryo, as well as adventitious shoot and lateral leaf formation [5,38,39]. For example, AtWUS is crucial for shoot apical meristem maintenance to replace the function of WOX1 and PRESSED FLOWER (PRS) [5]. The ectopic overexpression of AtWUS in tobacco is involved in stem cell fate and lateral leaf formation [40]. In Medicago truncatula, the AtWUS homolog HEADLESS (HDL) is also involved in leaf development [31]. WUS gene functions downstream of the CLAVATA3 (CLV3) signaling pathway [41]. The WUS gene is expressed in the organizing center and enhances CLV3 expression in stem cells. Likewise, CLV3 negatively regulates meristem size by suppressing WUS expression [42,43]. WUS interacts with CYCLOIDEA 2 (CYC2) and regulates reproductive organ development (ovary, stigma, and style) in Chrysanthemum morifolium [44]. In the embryonic columella, WOX5 and CYCD3;3/CYCD1;1 facilitate cell proliferation, and CYCD3 plays a major role in normal cell division [45]. The WUS orthologs STERILE AND REDUCE TILLERING (MOC3/SRT1) and TILLERS ABSENT1/MONOCULM 3 are involved in bud formation and female fertility of rice [46]. The WUS gene regulates histone acetylation and interferes with HISTONE DEACETYLASE (HDAC) activity, which stimulates the auxin signaling pathway in stem cells [38]. In addition, WUS expression is indirectly repressed by AGMOUS (AG) and stimulates expression of the zinc finger TF C2H2 type (KNUCKLES), which in turn suppresses WUS expression directly or indirectly involved in the maintenance of floral meristem cells [5]. Thus far, functional characterization of WUS genes has been studied in other plants but is rarely studied in citrus.
The above-mentioned studies report that WOX TFs primarily affect plant development by regulating the expression of downstream genes. Notably, WOX protein cloning and functional analyses were predominantly focused on model plants such as Arabidopsis. However, there have been relatively few studies on the regulatory and genetic development role of the WOX family in citrus. The availability of the citrus genome database gives us a valuable genetic resource to study specific sweet orange genes [47,48]. A thorough screening of the citrus database allowed us to find the evolutionary, regulatory, and developmental role of WOX orthologs in sweet orange. In the current study, we identified 11 putative WOX members in Citrus sinensis. Tissue-specific expression patterns and co-expression profiles of CsWOXs under water deficit and floral inductive conditions were comprehensively studied. The subcellular localization, transactivation activity, and DNA-binding ability confirmed that CsWUS may be a TF. Moreover, CsWUS overexpression and the virus-induced gene silencing (VIGS) assay revealed new insights into floral organ development, stem cell activity, and leaf development in citrus. In addition, yeast two-hybrid assays and DNA-protein interactions confirmed the complex involvement of CsWUS in developmental regulatory networks. Our data provide an evolutionary, co-expression, and spatial expression analysis of the WOX gene family and new perspectives that contribute to the function of the WOX family in citrus.

Genome-Wide Identification and In Silico Subcellular Localization Prediction of CsWOX Gene Family
To identify citrus WOX genes, Arabidopsis and rice WOX proteins were used as queries and all resulting sequences were retrieved from the sweet orange database (http://citrus. hzau.edu.cn/cgi-bin/orange/blast) using BLASTP. After removing sequence redundancies of the same protein, a total of 11 potential WOX proteins were identified as being allied with CsWOX proteins (Table 1). We further confirmed that all of these CsWOX proteins contained the homeodomain (PF00046 and SM00389). The 11 putative WOX proteins corresponding to the gene were named according to their physical location (from top to bottom) on chromosomes 1-8 (Table 1). Notably, one gene (CsWOX10) is located on an unknown chromosome. The coding sequence (CDS) length of CsWOX genes varied from 582 bp (CsWOX7) to 1104 bp (CsWOX2), encoding polypeptides of 193-367 amino acids in length, with a predicted molecular weight range of 15,437.5-40,473.1 Da and a theoretical isoelectric point (pI) ranging from 5.4 to 11.5 (Table 2). In addition, the subcellular localization of CsWOX proteins was predicted ( Table 2). The 3D structure of proteins was also predicted (Supplementary Figure S1). The predicted locations of three CsWOX proteins (CsWUS, CsWOX3, and CsWOX5) were found to be nuclear localized. The remaining members of CsWOX were projected to be localized in chloroplast, mitochondria, cytoplasm, or plastid.   To explore the phylogenetic relationship of WOXs between citrus and other model plants, a phylogenetic tree was resolved using 39 WOX family members from Arabidopsis (15 genes), rice (13 genes), and sweet orange (11 genes). The phylogenetic distribution showed that all CsWOX members were grouped into three clades, modern/WUS, intermediate, and ancient, consistent with previous WOX family distribution schemes [2]. The modern/WUS clade was the largest in this phylogenetic tree, containing 18 members: four members from sweet orange, six from rice, and eight from Arabidopsis. The intermediate clade was the second largest and included two from sweet orange, six from rice, and four from Arabidopsis members. The ancient clade had five members from sweet orange, one from rice, and three from Arabidopsis. Additionally, we also explored the orthologous relationships among sweet orange, rice, and Arabidopsis WOX families. These included 13 orthologous genes, putative orthologs with sweet orange proposed based on the phylogenetic tree were as follows: CsWOX5/AtWOX10, CsWOX1/AtWOX13, CsWOX6/AtWOX11, CsWOX2/AtWOX9, CsWOX8/AtWOX1, CsWOX7/AtWOX5, CsWOX9/AtWOX2, CsWOX1/ OsWOX9B, CsWOX6/OsWOX11/12, CsWOX2/OsWOX9D, CsWOX8/OsNS1/2, CsWOX7/ OsWOX5, and CsWOX9/OsWOX3 ( Figure 1A).

Gene Structure and Synteny Analysis of CsWOX Genes
To further observe the structural diversity of the WOX genes in sweet orange, an exon-intron diagram of the CsWOX genes was created with reference to their genomic and coding sequences; the number of exons varied from two to four and the number of introns from one to three in CsWOXs ( Figure 1C). Furthermore, MCScan was used to identify duplicate gene types ( Figure 1B). Almost all CsWOX genes were singletons. In a syntenic block (Citrus sinensis and Arabidopsis), each member belonged to the same subfamily and phylogenetic group.  To gain further insight into the expression of sweet orange WOX genes, we carried out in silico analysis of potential cis-elements of each member of the CsWOX family conferring responsiveness to plant hormones (ABRE, TGA-element, P-box, CGTCA-motif, TGACG-motif, GARE-motif, TCA-element, and TATC-motif) within their promoters. Ciselements related to endosperm and meristem expression (CAT-box and GCN4-motif) and biotic/abiotic stress response (TC-rich repeats, LTR, and MBS) were also found in the promoters of most CsWOXs (Figure 2A). Cis-elements involved in light responsiveness and anaerobic induction (ARE) were found in the promoter region of almost all CsWOXs. Meanwhile, the cis-elements of whole CsWOX family were divided into different categories based on function prediction [50]. The result showed that the largest category of cis-elements was hormone-related, followed by stress and development-related cis-elements (Supplementary Figure S2). However, many motifs have not yet been functionally characterized and whether these motifs confer unique functional roles to CsWOXs remain to be further investigated.

Motif Analysis of CsWOX Family
In addition, the conserved motifs of CsWOX proteins were analyzed by the MEME Suite. A total of 15 conserved motifs were predicted in the CsWOX proteins ( Figure 2B). The size of the identified motifs ranged from 6 to 50 amino acids ( Figure 2B). The results show that groups classified by phylogenetic analysis shared similar conserved motif compositions. However, there were some differences. For example, motif 3 was present in all members of the CsWOX family; motif 15 was found in just two members, CsWOX5 and CsWOX9; motif 4 was found in CsWOX5 and CsWOX1; motif 9 was found in three members, CsWOX10, CsWOX3, and CsWOX4, of the ancient clade; and motifs 7 and 5 were distinctly found in the intermediate and modern/WUS clades, respectively ( Figure 2B). To some extent, these specific motifs may lead to the functional differences of WOX genes in sweet orange.

Expression of CsWOX Genes in Different Citrus sinensis Tissues and under Floral Inductive Water Deficit Conditions
To further determine the function of CsWOX genes, their tissue-specific expression was evaluated ( Figure 3A). The results revealed that the majority of CsWOX genes were expressed in multiple tissues, although some showed similar expression trends in different tissues. For example, CsWUS showed the highest expression in the flower and stem. CsWOX1 and CsWOX2 were constitutively expressed in all tissues, implying that they are involved in multiple developmental stages, while other CsWOX genes were differently expressed, suggesting that they have tissue specificity. CsWOX4 was highly expressed in stem and apical meristem but had no expression in fruit. CsWOX3, CsWOX8, and CsWOX6 showed similar expression patterns in apical meristem ( Figure 3A). However, CsWOX8 was expressed in the flower and CsWOX6 had relatively high expression only in apical meristem tissue, indicating that it may be involved in apical meristem cell maintenance. CsWOX9 was highly expressed in the flower and had low expression in the stem and leaf, indicating functional redundancy in flower development. CsWOX10 had maximum expression in the leaf and relatively less expression in the flower. CsWOX7 had high expression in the root. However, CsWOX5 had no expression in fruit ( Figure 3A). Collectively, this indicates that these genes may be involved in the maintenance of stem cells and in organ development and differentiation. To investigate the potential role of CsWOX genes under floral inductive water deficit conditions, we evaluated the expression patterns using previously reported RNA-seq data in lemon [51]. RNA-seq analysis was performed at stage 1, one week before the water deficit; stage 2, one week after the beginning of the water deficit; and stage 3, one week after the release of the water deficit. We categorized 10 CsWOX genes into three clusters and excluded the CsWOX7 gene as it has no expression in this specific condition ( Figure 3B). Cluster 1 consisted of four genes (CsWOX3, CsWOX4, CsWOX1, and CsWOX5).
These genes were subsequently induced at the beginning of the water deficit condition, while CsWOX3 had wide expression at stage 3 ( Figure 3B). Such genes may be crucial for citrus flower bud differentiation. Cluster 2 genes (CsWUS, CsWOX6, and CsWOX10) were suppressed at stage 2 and induced at stage 3. This cluster showed upregulation in the expression of genes involved in vegetative growth after water recovery. Cluster 3 genes (CsWOX8, CsWOX9, and CsWOX2) were suppressed throughout the water deficit and water recovery treatment. Four genes (CsWOX4, CsWOX1, CsWOX6, and CsWUS) were differentially expressed depending on probability ≥ 0.8 and absolute value of the log 2 ratio of ≥1 as a threshold. These findings show that CsWOX genes have narrow expression under water deficit conditions ( Figure 3B).

Co-Expression Analysis of CsWOX Genes under Water Deficit Floral Initiation
Co-expression analysis was done under floral inductive water deficit conditions to determine the probable role of CsWOX genes in citrus ( Figure 4A). Co-expression analysis generally clustered 11 CsWOX genes into two modules, together with 1638 differential expressed genes via 42,176 interactions. For instance, module 1 was the widest co-expression group, consisting of nine CsWOX genes (CsWUS, CsWOX8, CsWOX9, CsWOX4, CsWOX2, CsWOX5, CsWOX6, CsWOX10, and CsWOX1). Genes from module 2 were co-expressed uniquely with CsWOX3 ( Figure 4A). The co-expression network and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis of CsWOX genes found a potential role for them in water deficit conditions and growth. The genes evaluated in the co-expression network were enriched in terms of the secondary metabolite biosynthetic process, response to wounding, water deprivation, or response to water as well as active transmembrane transporter activity ( Figure 4B). The co-expression network was divided into two modules, as illustrate above, and it was important to check respective functions performed by different modules. In addition, the Gene Ontology (GO) analysis showed that these genes from various modules had diverse functions in various facets of plant development.

Overexpression Analysis of CsWUS in Tobacco and Gene Silencing Analysis in Lemon
CsWUS consists of functional domains with the dual function as a suppressor or activator similar to AtWUS, as previously reported [5,52]. In fact, CsWUS, as a member of the WUS clade, plays an important role during the growth and development of plants [24], yet its functional role is rarely studied in citrus. Therefore, we focused on functional characterization of CsWUS during development and growth of citrus. Sequence alignment and phylogenetic analysis showed that Cs1g25270 has high similarity with Arabidopsis compared with other citrus CsWOX genes, and thus was named CsWUS (Supplementary Figure S3). To investigate the functional role of CsWUS, it was over-expressed in tobacco driven by the 35S promoter. A total of eight independent transgenic lines (CsWUS-OE) were obtained ( Figure 5A). All CsWUS transgenic plants showed a similar phenotype ( Figure 5A). CsWUS-OE transgenic plants had leaf lamina curled up inside, reduced leaf area, and stunted plant growth ( Figure 5B). The gynoecium of the flower in CsWUS-OE lines was markedly smaller compared with the wild type (WT) ( Figure 5C), which suggests that overexpressed WUS may restrict the ability of female reproductive organ development. It is worth noting that CsWUS-OE plants had empty seed pods with no seed development compared to WT seed pods ( Figure 5D). Furthermore, the relative expression of the endogenous NtWUS ortholog gene was performed in CsWUS-OE transgenic plants, the results showed that the expression of NtWUS was significantly up-regulated compared with the wild-type (Supplementary Figure S4). It is known from overexpression analysis that CsWUS may be involved in aberrant cell division, stem cell fate, and organ development. To further investigate the functional role of CsWUS in citrus, the VIGS approach was used to knock down the CsWUS in lemon. Transcript abundance of CsWUS in the positive VIGS plants, designated as TRV (tobacco rattle virus)-CsWUS, were repressed compared to the control seedlings (designated as TRV) that were only infiltrated with empty vector ( Figure 5F). Two weeks following VIGS knockdown, conspicuous differences were noticed in plant morphology with thorn like meristem outgrowths present on the plant stems ( Figure 5E). These results suggested that CsWUS may be involved in thorn development activity in citrus.

Identification of Interacting CsWUS Proteins
CsWUS is generally associated with other proteins and TFs to form a transcription complex that regulates the expression of downstream genes directly or indirectly. To recognize the function of CsWUS protein, we evaluated its interaction partners via yeast two-hybrid screening. The conserved domain of CsWUS was cloned into the PGBKT7 vector. CsWUS-BD had no activation capability by the self-activation detection experiment; we used it as bait to perform yeast two-hybrid screening against the citrus mix cDNA expression libraries. Five interaction proteins were identified in at least two independent screens, indicating CsWUS interacting protein ( Figure 6A). These protein interactions are reported for the first time with CsWUS protein interaction: CsCYCD3 (Cs3g23120) ortholog of AtCYCD2;1 encodes cyclin-dependent kinase; zinc finger protein (Cs1g13130), which is an ortholog of Arabidopsis; miraculin-like protein (Cs7g08260) encodes soybean Kunitz super-family ortholog of AtMLP-like protein; Shaggy related protein kinase CsSK7 (Cs2g04660) ortholog of AtSK7 encodes catalytic domain of serine/threonine kinase and glycogen synthase kinase 3; and Cs14-3-3 protein (orange1.1t01991) encodes the 14-3-3 domain ubiquitous class of phosphoserine threonine binding protein. Generally, CsWUS may be related to various complexes composed of these interacting proteins.
Among these five interacting proteins, CsCYCD3 protein has a crucial role in enhancing cell division, while CsWUS protein is involved in stem cell fate. Sequence alignment and phylogenetic analysis showed that CsCYCD3 has high similarity with Populus trichocarpa CYCD protein compared with CYCD3 from other plants ( Figure 6B). Sequence alignment showed that CsCYCD3 has high similarity with Arabidopsis AtCYCD2;1 ( Figure 6C). Hence, it is hypothesized that CsWUS interacts with CsCYCD3 and regulates reproductive growth and stem cell fate. To further confirm the interaction between CsWUS and CsCYCD3, we performed BiFC analysis in tobacco ( Figure 6D). The findings show that the interaction between CsWUS and CsCYCD3 takes place in the nucleus. Therefore, we may speculate that the interaction between CsWUS and its co-partner CsCYCD3 is involved in tissue proliferation. These results provide a basis for further evaluation of these proteins regarding undifferentiated growth in citrus.

Sub-Cellular Localization, CsWUS Transcription Activation Analysis and the Identification of RAP2.12 and CsHB22 Transcription Factor in Citrus sinensis
To identify the subcellular localization of CsWUS, its full-length open reading frame (ORF) sequence was cloned into the pRI101 vector under 35S promoter, resulting in an in-frame fusion protein of CsWUS:GFP. The results show that the tobacco epidermal cell expressing GFPs showed cytoplasmic and nuclear staining, while CsWUS:GFP can also be detected in the whole cell, similar to a positive control ( Figure 7B). To determine the transcription activity of CsWUS, it was fused with GAL4 DNA-binding domain (GAL4BD) and tested in yeast AH109 ( Figure 7A). The result shows that CsWUS may act as a transcriptional activator. Consequently, the DNA binding ability of CsWUS gene confirms that CsWUS binds with (TAATTCA) motif and verifies the DNA binding ability of CsWUS through yeast one-hybrid ( Figure 7C).
To identify TFs that regulate CsWUS, a yeast one-hybrid assay was performed using ProCsWUS as bait. The 2kb CsWUS promoter fragment was cloned and 0.5 kb core sequence (from -579 to -1003) was inserted into pAbAi and used as bait. We obtained 20 positive clones; only two genes (Cs1g16690 and Cs3g22190) were found after putative re-streaking on high stringency medium supplemented with 100 mM Aureobasidin A (Papdi, #207) [53] ( Figure 7E). Cs1g16690, belonging to the AP2/ERF TF family, and Cs3g22190, belonging to the zinc finger homeodomain TF family, were identified from the citrus genome database. Meanwhile, a number of Cis-elements were found on ProCsWUS, including AP2/ERF (GGCGGCC) elements, which have been recognized by Cs1g16690. ZF-HD (TGATTAG) elements have been shown to be recognized by Cs3g22190 ( Figure 7D).
The results of comparison between cDNA and genomic DNA sequences reveled that Cs3g22190 is located on chromosome 3. Alignment and phylogenetic analysis showed that this gene has high similarity with Arabidopsis ZF-HD homeobox protein 22 (At4g24660), and thus was named CsHB22 (Supplementary Figure S5). It is composed of a 630 bp full-length ORF encoding a 209 amino acid putative protein. The CsHB22 consists of ZF-HD dimer protein and homeobox domain, consistent with previous reports on HB22 protein. A comparison was made between Cs1g16690 cDNA and genomic DNA located on chromosome 1. Alignment and phylogenetic analysis showed that this gene has similarity with Arabidopsis BRELATED TO AP2.12 (AT1G53910) and it was named CsRAP2.12 ( Supplementary Figure S6). It is composed of 1170 bp, ORF encoding 326 amino acids. CsRAP2.12 consists of AP2-ERF domain and ORC2 super-family. To further confirm CsRAP2.12 and CsHB22 binding to the CsWUS promoter, we investigated whether CsRAP2.12 and CsHB22 activated or suppressed ProCsWUS in vivo by performing dual luciferase assay on tobacco leaves. In this study, CsRAP2.12 and CsHB22 were used as effectors and two constructs consisting of ProCsWUS were used as reporters ( Figure 7F). The results show that co-transformation of effectors and reporters significantly elevated the promoter activity of CsWUS ( Figure 7F). Taken together, this suggests that these two TFs may be involved in phenotypic complementation of CsWUS and alter growth compared with the control by regulating the intrinsic hormonal and developmental pathways.

Expression of CsRAP2.12, CsHB22, and CsCYCD3 in Different Citrus Sinensis Tissues and under Floral Inductive Water Deficit Conditions
To further investigate the spatial expression pattern of CsRAP2.12, CsHB22, and CsCYCD3, their tissue specific expression was evaluated in different Citrus sinensis tissues including the leaf, flower, fruit, stem, apical meristem, and root ( Figure 8A-C). The results showed that CsRAP2.12 showed a higher level of expression in roots, leaves, and stems, and a lower level of expression in flowers, fruits, and apical meristem ( Figure 8A). CsHB22 was relatively highly expressed in the leaf, flower, and fruit ( Figure 8B). Moreover, we also performed tissue specific expression analysis of CsCYCD3 ( Figure 8C). CsCYCD3 was mainly expressed in the leaf, stem, and apical meristem ( Figure 8C). It is worth noting that all three genes showed high levels of expression in the leaves. Overall, these three genes presented a very broad expression pattern, implying that they may play multiple roles in the growth and development of sweet oranges. To explore the potential role of CsRAP2.12, CsHB22, and CsCYCD3 under floral inductive water deficit conditions, we also evaluated their expression patterns using previously reported RNA-seq data in lemon [51]. The result showed CsRAP2.12 was highly expressed after water recovery while CsHB22 is not induced under floral inductive water deficit conditions ( Figure 8D). In addition, CsCYCD3 was more highly expressed at stage 3 than stage 2 ( Figure 8D). These results suggest that they may also play an important role under floral inductive water deficit conditions.

Discussion
The WOX family plays a crucial role in shoot apical meristem and embryonic development, stem cell activity, and various other developmental processes in plants [5,22,44]. Due to the significance of their functions, research on plants has become more urgent and widespread. Previous genome-wide analysis of the WOX gene family was done in some important plant species [13,18,44,54]. However, there is no report of genome-wide analysis of the WOX family in citrus. In this study, we identified 11 putative CsWOX genes in the citrus reference genome. Based on their phylogenetic relation with Arabidopsis and rice WOX proteins, CsWOX proteins were split into three clades (modern/WUS, intermediate, and ancient), consistent with already reported classifications in different plant species [21,24,55]. However, synteny between CsWOX genes and their Arabidopsis homologs was less than expected.
Identifying the role of CsWOX genes by tissue-specific expression and co-expression analysis is an important and useful tool. For example, CsWOX1 a putative ortholog of AtWOX13 has highly expressed in the flower, root, and apical meristem. The WOX13 mutant exhibited slightly wider fruits with a reduced number of lateral roots in Arabidopsis [9]. However, GhWOX13 in cotton is involved in hormonal mediation of fiber elongation [56]. In this study, tissue-specific expression analysis suggests that CsWOX1 performs a similar function to AtWOX13 [9]. CsWOX7 is a putative ortholog to AtWOX5 in Arabidopsis, and CsWOX5 was shown to be involved in root meristem maintenance [12]. Further studies of AtWOX5 orthologs in maize (ZmWOX5), poplar (PtoWOX5a), and rice (OsWOX5) showed a similar conserved function and molecular mechanism in root development [57,58]. Tissue-specific expression analysis showed that CsWOX7 is highly expressed in root tissue, indicating that it may be involved in root development of sweet orange. CsWUS was found as a putative ortholog of AtWUS and OsWUS. WUS is a critical regulator and encodes a homeodomain protein that is mandatory for stem cell activity [59]. Nevertheless, AtWUS regulates maintenance of stem cells in floral meristems, and its expression level influences the number of flower organs that develop [60]. AtWUS is a bi-functional TF that acts to repress stem cell regulation and activate floral patterning [42]. Previous findings explained the conserved functioning of WUS genes in Gossypium hirsutum, Arabidopsis, Medicago truncatula, Glycine max, Triticum aestivum, Coffea canephora, Nicotiana tabacum, Oryza sativa, and Chrysanthemum morifolium [5,13,31,39,40,44,61,62]. High expression of CsWUS in the flower and stem suggests that it may be linked with maintenance of stem cell and flower organ development.
Interestingly, CsWUS-OE in tobacco exhibited a developmental role and induced ectopic growth. Overexpressed plants had malformed leaves, stunted growth with limited gynoecium development, and no seed development. Similar results were observed with AtWUS overexpression in tobacco [40]. CsWUS consists of functional domains considering its dual function as an activator similar to AtWUS, as previously reported [5,52]. Results show that the WUS gene has a conserved developmental role in Arabidopsis and sweet orange. However, gene silencing of CsWUS revealed conserved functioning in stem cell activity involved in thorn development in lemon.
To further investigate the regulation mechanism of CsWUS, we performed a yeast two-hybrid analysis using CsWUS as a bait, and five interacting protein partners of CsWUS were identified. Phenotypic functional complementation experiment with CsWUS confirmed involvement in aberrant cell division and stem cell activity. CsCYCD3 ortholog in Arabidopsis enhances cell division and plays a vital role in tissue proliferation, such as in the meristem and young leaves [9,63]. Cyclins are well conserved in functioning, and therefore have been comprehensively recognized in plants. Previous studies reported that CYCD3 defines distinct developmental zones and is locally regulated by CYCLOIDEA [64]. In addition, CYCLOIDEA (CYC2-like) and WUS regulate reproductive organ development in Chrysanthemum morifolium [44]. Furthermore, we found physical interaction between CsCYCD3 and CsWUS in the tobacco nucleus. Therefore, we may predict that CsCYCD3 and CsWUS complex might be involved in stem cell activity and promote ectopic growth in citrus. The conserved mechanism of CsCYCD3 and CsWUS in citrus needs to be studied further. Subsequently, other interacting proteins, such as 14-3-3 protein 6, shaggy related protein kinase, zinc finger homeodomain, and MLP-like protein are involved in various developmental, hormonal signaling, and stress responses in Arabidopsis [65][66][67][68]. These proteins' triggering mechanisms still need to be further elucidated.
In addition, we identified two TFs (CsRAP2.12 and CsHB22) binding to the CsWUS promoter by yeast one-hybrid library screening. These results indicate that the expression of CsRAP2.12 and CsHB22 activates the expression of CsWUS by binding to its promoter. CsRAP2.12 is an ortholog of AtRAP2.12 and belongs to the ERF-VII TF family. AtRAP2.12 and its orthologs shared a conserved AP2 domain that was mandatory for protein-DNA interaction. AtRAP2.12 involved in hypoxia tolerance and reduced growth in the presence of oxygen in Arabidopsis [53]. The ubiquitin-dependent N-end rule pathway for protein degradation, functions as an oxygen sensing mechanism in Arabidopsis [69]. Furthermore, the presence of molecular oxygen affected the stability of ERF-VII proteins, implying that the role of oxygen sensing was mediated via the N-end rule protein degradation pathway [70]. The N-end rule pathway, in particular, stagnates the stress response in plants, enabling optimal growth and development. Meanwhile, RAP2.12 was also an ethylene responsive TF, as ethylene is known to influence the growth and development of leaves and can be independent of or dependent on its interaction with other hormones [71,72]. We observed malformed curled leaves in CsWUS-OE lines, which showed similar phenotype behavior to those plants grown in the presence of polar auxin transport inhibitors. It is known that auxin regulates the cell division phase during leaf expansion [40]. We may speculate that the response of leaf growth to ethylene is likely to be auxin-dependent or auxin-independent [33]. In addition, overexpression of stabilized RAP2.12 alters the leaf phenotype by regulating non-hypoxic target genes in Arabidopsis [73]. Therefore, we suggest that CsRAP2.12 binds with promoter of CsWUS, and this interaction may be involved in hormonal signaling that regulates leaf development.
ZF-HD belongs to a subfamily of homeodomain TFs that has not been well characterized functionally in plants [68]. In Arabidopsis, 14 members of the ZF-HD family have been identified and predominantly expressed in flower tissues [74]. The CsHB22 ortholog in Arabidopsis (AT4G24660) was highly expressed in flower tissues [68,74]. However, its specific function needs further investigation in Arabidopsis. Recently, a member from the ZF-HD family (OsZHD2) was described that promotes root and meristem activity by biosynthesis of ethylene in rice [75]. Besides the regulatory role in the growth and development, some ZF-HD genes may play a vital role in response to abiotic stresses such as drought, heat, cold, and salt [76]. Recent studies indicate that ZF-HD is induced by cold stress, NaCl, and PEG in wheat [77]. Consistent with all previous reports, we also found the flower gynoecium in CsWUS-OE was smaller than that in WT. Thus, overexpression of CsWUS gene in sweet orange is expected to retard the ability to develop female gametophytes. In addition, tissue specific expression showed that CsRAP2.12, CsHB22, and CsCYCD3 were expressed in the flower, leaf, stem, and apical meristem. Taken together, we therefore hypothesized that the interaction of CsWUS-CsCYCD3, ProCsWUS-CsRAP2.12, and ProCsWUS-CsHB22 may be involved in development of sweet orange leaf, flower, stem, and apical meristem. However, further studies are required to confirm these initial findings.

Plant Materials
Tissues samples of leaf, flower, root, apical meristem, stem, and young fruit were obtained from sweet orange grown at the National Citrus Research Breeding Center, Huazhong Agriculture University, Wuhan, China. Flower and fruit samples were collected in flowering and fruiting seasons, respectively. The samples were collected for CsWOX expression analysis and promptly placed in liquid nitrogen and further preserved at -80 • C for RNA extraction and expression analysis.

Identification of WOX Genes in Sweet Orange
To recognize all putative WOX genes of sweet orange, a local BLAST search using the citrus database (http://citrus.hzau.edu.cn/orange/) was done [48]. Two BLASTP approaches were implemented to search the sweet orange WOX genes. First, 15 known Arabidopsis WOX genes were downloaded from the Arabidopsis Information Resource (TAIR) database and used to query the citrus database, and candidate genes were recognized by BLASTP search scores of ≥100 and e-value of ≤e −10 . Second, the same procedure was done using all known rice genes downloaded from the NCBI database, which were used as query sequences [78]. The 11 identified WOX genes were designated as CsWOX genes.

Phylogenetic Analysis and Gene Structure of CsWOX Family
To explore the phylogenetic relationship between sweet orange, Arabidopsis, and rice, a phylogenetic tree was constructed using the Clustal Omega program (Guide Tree) (www.ebi.ac.uk/Tools/msa/clustalo/) follow the default parameters. The corresponding phylogenetic tree data was downloaded, and the phylogenetic tree was drawn using Interactive Tree of Life (IToL) v. 4 (https://itol.embl.de/), scale bars correspond to 0.1 substitution [79,80]. In this study, Arabidopsis and rice CsWOX proteins were used as the out group. The protein sequences of Arabidopsis and rice were downloaded from the TAIR and NCBI databases, respectively. Gene Structure Display Server 2.0 was used to evaluate the exon/intron structures of CsWOX genes [81].

Analysis of Conserved Motif and Predicted Subcellular Localization
To study the structural modification of CsWOX genes, the conserved motifs in the encoded CsWOX proteins were analyzed by Multiple Expectation Maximization for Motif Elicitation (MEME) v. 5.3 (https://meme-suite.org/meme/tools/meme) with default parameters [50]. Conserved motifs were identified with the motif widths of 6-50 residues. The online program WOLF PSORT II (http://www.genescript.com/wolfpsort.html) was used to predicate subcellular positions of CsWOX genes with default parameters, the organism type selected plants [82].

Analysis of Protein Structures
The biochemical properties of CsWOX proteins such as amino acid composition, molecular weight (MW), theoretical pI, instability index, aliphatic index, and grand average of hydropathicity (GRAVY) were obtained by an online tool on the Bioinformatics Resource Portal ExPASy server (http://web.expasy.org/protparam/) with default parameters [83]. The protein structure of CsWOX proteins was predicted by the SWISS-MODEL (https: //swissmodel.expasy.org) online programs with default parameters [82].

Expression Profile of CsWOX Genes in Sweet Orange
Spatial expression of CsWOX genes in various organs was analyzed by qRT-PCR. Total RNA was isolated from young fruits, flowers, healthy leaves, roots, stems, and apical meristems of sweet orange. cDNA synthesis was conducted with a PrimeScript ® RT reagent kit (Takara, Dalian, China) following the manufacturer's instructions. The synthesized cDNA was diluted 1:10 as a template for qRT-PCR. The primers used are listed in Supplementary Table S1. Real-time PCR was performed with Hieff ® qPCR SYBR ® Green Master Mix (Yeasen Biotech Co., Ltd., Shanghai, China) on an ABI PRISM 7000 system (Applied Biosystems) using 1 µL of cDNA template, 8 µL of double distilled water, and 0.5 µL of forward and reverse primers (10 µM) in the following PCR condition: 95 • C for 5 min, 95 • C for 10 s, 55-60 • C for 20 s, 72 • C for 20 s, and 40 cycles. The relative expression level of target genes was calculated using the 2 −∆∆Ct method by normalizing CsActin as described previously [87].

Co-Expression Evaluation of CsWOX Genes during Induction of Floral Water Deficit
The specific role of CsWOX genes was evaluated through published transcriptome data on lemon bud under floral inductive water deficit conditions [51]. Co-expression network of CsWOX genes was created, and the relationships between two genes above 0.85 were retained and then visualized as a comparison in Cytoscape. Co-expression gene networks were evaluated through Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) using the clusterProfiler and pathview packages in R, as defined in [88].

Subcellular Localization, Motif Binding, and Transactivation Activity Assay of CsWUS
The coding sequence of CsWUS without stop codon was fused into PRI101 vector containing GFP gene under control of CaMV 35S promoter to form 35Ss:CsWUS-GFP. The 35S:CsWUS-GFP and 35S:GFP constructs were transformed into Agrobacterium tumefaciens GV3101. The 35S:CsWUS-GFP and 35S:GFP constructs were transiently expressed in Nicotiana benthamiana leaf epidermal cells. After 36 h of incubation, leaf cells were examined by confocal microscopy as described previously [89]. The full-length coding sequence of CsWUS was inserted into pGBKT7 as bait according to the manufacturer's instructions. The bait clone with the empty prey pGADT7 and co-transformed cells on SD/-Ade/-His/-Leu/-Trp medium was transformed into yeast AH109 strain. Transactivation activity was determined by a previously reported method [19]. CsWUS DNA binding activity was investigated by the Matchmaker TM Gold Yeast-one Hybrid Library screening system according to the user manual (Clontech, Mountain View, CA, USA). Three tandem repeats of the predicted TAATTCA motif were inserted into the pAbAi vector and transformed in the yeast-one hybrid system according to the manufacturer's instruction. The yeast clones expressing pAbAi and pGADT7-CsWUS were grown normally on SD/-Leu medium with AbA ng/mL [90]. The experiment was repeated three times. The primers used are listed in Supplementary Table S2.

Construct p35s-CsWUS Preparation and Overexpression in Tobacco Plants
The full-length coding sequence of CsWUS was used for overexpression analysis. For this analysis, CDS of CsWUS was cloned into pBI121 vector by replacing GUS gene (Supplementary Table S1). The sequence was then inserted into the pBI121 vector. The vector was transformed into Agrobacterium tumefaciens GV3101 by the heat shock method. A previously described method of Agrobacterium mediated transformation was used in tobacco [91]. The transgenic plants T 0 were confirmed by PCR amplification.

Vector Construction and VIGS of CsWUS in Lemon
The tobacco rattle virus (TRV) system (pTRV-RNA1 and pTRV-RNA2) was used for VIGS analysis. A 428 bp gene fragment of CsWUS was cloned into the pTRV-RNA2 vector to produce the pTRV2: CsWUS construct. pTRV1, TRV2 (negative control), and pTRV2-CsWUS were transformed into Agrobacterium tumefaciens strain GV3101 [89]. Agroinfiltration proceeded by dipping germinating lemon seeds with a shoot length of around 1 cm in a bacterial suspension in a vacuum chamber. The plants were dried with filter paper after vacuum infiltration and grown in the dark for 3 days, then sown in soil containers under a growth chamber at 25 • C, 16 h light/8 h dark [89]. After 2 weeks, DNA of each seedling was extracted with a DNeasy Plant Mini kit (Qiagen, Hilden, Germany) and subjected to genomic PCR using one pair of primers for detection of positive plants (Table S1), and qRT-PCR was done for each positive plant to measure the transcript level of CsWUS [92]. The primers used are listed in Supplementary Table S1.

Yeast Two-Hybrid Screening
The yeast two-hybrid library screening was completed with the Matchmaker Gold Yeast Two-Hybrid system (Takara Bio, Beijing, China), using yeast strain AH109. Conserved domain of CsWUS was inserted into pGBKT7 vector and used as bait. Single yeast clones were selected on solid high-stringency SD/-Ade/-His/-Leu/-Trp medium and grown on liquid SD/-His/-Leu medium. Recombinant pGADT7 plasmids with cDNA inserts were selected, re-transformed, and checked by X-α-Gal filter-lift assay before sequencing. The selected query was reconfirmed through NCBI, phytozome, and citrus genome databases for further identification of corresponding genes. After library screening, these recombinant plasmids were transformed again on yeast high-stringency SD/-Ade/-His/-Leu/-Trp X-α-Gal medium [13]. The primers used are listed in Supplementary Table S1.

Bimolecular Florescence Complementation Assay (BiFC)
To investigate the bimolecular florescence complementation of CsWUS and interacting proteins, their ORF sequences without stop codon were amplified and cloned into pUC-SPYNE (nYFP) and pUC-SPYCE (cYFP) vectors. Full-length CDS of CsCYCD3 was inserted into pFGC-nYFP vector to generate N-terminal in-frame fusions with N-YFP, while CsWUS coding sequences were cloned into pFGC-cYFP vector to form C-terminal in-frame fusions with C-YFP. All plasmids were transformed into Agraobacterium GV3101, and infiltration of tobacco leaves was performed following a previously described method [13]. After 36 h, tobacco leaves were observed under a Leica confocal laser scanning electron microscope.

Yeast One-Hybrid Screen and Assay
The 0.5 kb CsWUS promoter was amplified and cloned into pAbAi and used as bait, and pGADT7 library was used as a prey. Yeast one-hybrid library screening assay was performed using the Matchmaker TM Gold Yeast-one Hybrid Library system (Clontech, Mountain View, CA, USA) according to the manufacturer's instructions. Protein-DNA interaction was revealed by the growth ability of co-transformed yeast cells on high-stringency SD/-leu medium supplemented with AbA following the manufacturer's protocol [19].

Dual Luciferase Reporter Assay
The CDS of CsHB22 and CsRAP2.12 was cloned into pGreenII 62-SK vector using the ClonExpressTM II One Step Cloning Kit (Vazyme Biotech Co., Ltd., Nanjing, China), and this construct was utilized as effector plasmid. The CsWUS promoter with specific binding motifs was cloned and inserted into pGreen0800-LUC vector plasmid to attain the reporter plasmid. The empty pGreenII 62-SK vector was used as control effector and CsRAP2.12 and CsHB22 were used as treatment effectors. Agrobacterium pSoup-19 was used for transformation of the effectors and promoters. For transient gene expression analysis, the reporter and effector recombinant plasmid constructs were co-transformed into leaves of Nicotiana benthamiana. After 2-3 days of infiltration, a Dual-Luciferase ® Reporter Assay System (Promega Biotech Co., Ltd., Beijing, China) was utilized to qualify LUC and REN activity according to the manufacturer's instructions. At least 6 biological replicates were organized for each co-transformation.

Conclusions
This study provides a genomic framework for the citrus WOX gene family and its phylogenetic relation with rice and Arabidopsis. A total of 11 CsWOX genes were identified from the citrus genome. Bioinformatics analysis, including gene structure, conserved motif, cis-regulatory elements, protein physiochemical properties, predicted structure, and subcellular localization was performed, providing a framework for further study of this gene family in citrus. Comprehensive analysis of spatial expression patterns and co-expression analysis suggest that the WOX gene family is probably involved in distinct developmental mechanisms and the response to water deficit conditions in citrus. Functional analysis of CsWUS shows that CsRAP2.12 and CsHB22 regulate CsWUS expression. Further, CsWUS was found to interact with CsCYCD3. Ectopic overexpression of CsWUS is involved in aberrant cell development, malformed leaves, and defective gynoecium and ovary development. CsWUS gene silencing displayed radially symmetric thorn-like outgrowth in lemon. Thus, we assumed that CsRAP2.12 and CsHB22 regulates CsWUS expression, which further interacts with CsCYCD3, which may be involved in stem cell activity and other intrinsic development in citrus. The functional characterization of CsWUS genes the necessary foundation for follow-up more research to analyze the role of CsWOX genes in citrus development.

Conflicts of Interest:
The authors declare no conflict of interest.