Genome-Wide Distribution, Expression and Function Analysis of the U-Box Gene Family in Brassica oleracea L.

The plant U-box (PUB) protein family plays an important role in plant growth and development. The U-box gene family has been well studied in Arabidopsis thaliana, Brassica rapa, rice, etc., but there have been no systematic studies in Brassica oleracea. In this study, we performed genome-wide identification and evolutionary analysis of the U-box protein family of B. oleracea. Firstly, based on the Brassica database (BRAD) and the Bolbase database, 99 Brassica oleracea PUB genes were identified and divided into seven groups (I–VII). The BoPUB genes are unevenly distributed on the nine chromosomes of B. oleracea, and there are tandem repeat genes, leading to family expansion from the A. thaliana genome to the B. oleracea genome. The protein interaction network, GO annotation, and KEGG pathway enrichment analysis indicated that the biological processes and specific functions of the BoPUB genes may mainly involve abiotic stress. RNA-seq transcriptome data of different pollination times revealed spatiotemporal expression specificity of the BoPUB genes. The differential expression profile was consistent with the results of RT-qPCR analysis. Additionally, a large number of pollen-specific cis-acting elements were found in promoters of differentially expressed genes (DEG), which verified that these significantly differentially expressed genes after self-pollination (SP) were likely to participate in the self-incompatibility (SI) process, including gene encoding ARC1, a well-known downstream protein of SI in B. oleracea. Our study provides valuable information indicating that the BoPUB genes participates not only in the abiotic stress response, but are also involved in pollination.


Introduction
Plant self-incompatibility exists mainly to prevent inbreeding and promote outcrossing [1], while the decline of self-inbreeding makes outcrossing more urgent. The signal transduction of sporphyric self-incompatibility responses in Brassica oleracea has been uncovered. When self-fertilized pollen falls on the stigma, the S-locus cysteine rich protein (SCR) in the pollen, combined with the extracellular domain of the S-locus receptor kinase (SRK) of the stigma mastoid cell, the SRK intracellular kinase domain, and the M-locus protein kinase (MLPK) are phosphorylated. At this time, arm repeat-containing protein 1 (ARC1) is located near the membrane of the nucleus [2]. The ARC1 arm repeat region interacts with the SRK kinase domain and MLPK to cause ARC1 activation, and its N-terminus then binds to Exo70A1. With the assistance of E1 and E2, Exo70A1 is ubiquitinated to a certain extent, guided by the U-box at the N-terminus of ARC1 to 26S proteasome, and degraded [3]. It

Retrieval and Identification of U-Box Family Genes in B. oleracea
In order to identify and classify the gene members of the U-box family of B. oleracea, 103 BoU-box gene family members were first searched for in the Brassica database (http://brassicadb.org/brad/) and the Bolbase database (http://ocri-genomics.org/bolbase/). Subsequently, after SMART (https: //smart.embl.de/) online analysis, it was found that four genes clearly did not have a U-box structure. For the accuracy of the results, we downloaded the genome-wide sequence from the Bolbase database and downloaded the conserved domain (PF04564) hidden Markov model (HMM) shared by U-box from the Pfam 32.0 database (http://pfam.xfam.org/). The protein sequences containing a conserved domain were obtained using the HMMER search tool (http://hmmer.org/) [20], and further SMART analysis was used to verify that four protein sequences without a U-box structure were deleted. Finally, 99 B. oleracea U-box genes were identified and the results were consistent with the first method (i.e., retrieval from the BRAD and Bolbase database). At the same time, 60 AtPUB genes homologous to the BoPUB genes were obtained by sequence alignment. Additionally, based on all the genome sequences of B. oleracea by transcriptome sequencing in our laboratory, we found all of these genes in the transcriptomes of pollinated B. oleracea A4 pistils by BLAST alignment (Table 1).

Phylogenetic, Conservative Motif, and Structure Analysis
The ClutsalW program was used for multiple protein sequence alignments, and MEGA X was used to construct the phylogenetic tree using the maximum likelihood (ML) method with 1000 bootstrap replications [21] (Figure 1). The constructed phylogenetic tree was imported into FigTree software (http://tree.bio.ed.ac.uk/) for presentation. MEME (Multiple EM for Motif Elicitation, http://meme-suite.org/tools/meme) was used to analyze the conserved motifs of proteins online; maximum motifs was set to 15, other parameters were set to defaults, and the analysis results were saved in a MEME.xml format [22]. The genomic-wide gff3 format file of B. oleracea was downloaded from EnsemblePlants, and the intron and exon structures of the U-box genome-wide were visualized with GSDS 2.0 (http://gsds.cbi.pku.edu.cn/) [23]. The other domains of the U-box protein were analyzed online using Batch CD-Search (https://www.ncbi.nlm.nih.gov/Structure/bwrpsb/bwrpsb.cgi) and the results were saved in the format of hitdate.txt [24]. The above results were visualized using TBtools [25].

Chromosome Localization, Genome Collinearity, and Ka/Ks Ratio Analysis
Under the Linux system, the chromosome location information of the U-box gene was extracted from the whole genome gff3 file of B. oleracea using the Perl language. The chromosome localization and collinearity visualization of the U-box gene were carried out using TBtools and Adobe Acrobat DC (Adobe, Beijing, China).The coding sequence (CDS) of U-box protein was downloaded from the Bolbase and TAIR (https://www.arabidopsis.org/index.jsp) databases, and compared using BLASTALL under the Linux system. The tandem repeat genes were screened according to the standard that the similarity between the two sequences was >80% and the alignment length of the two sequences was greater than 75% of the longer sequences [26]. Collinearity analysis was the completed using Circos software (http://circos.ca/) ( Figure 2). The ratio of non-synonymous substitution rate (Ka) to synonymous substitution rate (Ks) of each pair of tandem duplication genes was calculated using KaKs_Calculator2.0 [27].

Expression Analysis of U-Box Gene Family and RT-qPCR
Total RNA from the different pollination B. oleracea pistil was extracted using an RNA Isolation kit (Tiangen, Beijing, China), and the samples were sent to the Biomarker Technologies Corporation (Beijing, China) for transcriptome sequencing. Using the measured RNA-seq data, including the whole genome information of the U-box in the pistils of B. oleracea after 0 min pollination, self-pollination (SP), and cross-pollination (CP) for 15 min, 30 min, and 60 min, the expression levels of these genes were quantified by log 2 (FPKM) values, and the expression results were visualized by TBtools. Fold change ≥2 and FDR <0.01 were used as screening criteria [28]. Thus, 36 differentially expressed genes were screened and identified by DESeq [29], and the expression levels of these genes were then verified by RT-qPCR. Gene-specific primers were designed in the online tool qPrimerDB (https: //biodb.swu.edu.cn/qprimerdb/), with BoActin as the internal gene; related RT-qPCR primers are listed in Table S1. RT-qPCR was carried out with a total volume of 20 µL (comprising 10 µL of TB Green Premix Ex TaqTM II, 0.6 µL of each primer, 2 µL of template (×10 diluted complementary DNAs), ROX 0.4 µL, and 6.4 µL of sterile distilled water). Real-time fluorescence quantitative analysis was performed using the Bio-Rad cfx-96 system with the following steps: denaturation at 95 • C for 30 s, followed by 40 cycles of denaturation at 95 • C for 5 s, and annealing/extension at 60 • C for 30 s [30]. Each sample was run through three technical repeats, and the relative expression level of each gene was evaluated using the 2 −∆∆Ct algorithm [31].

Promoter Region Analysis of U-Box Differentially Expressed Genes
The 1500-bp sequences upstream of the translation initiation sites of 36 differentially expressed genes were extracted from the B. oleracea genome sequence. The cis-acting elements of the promoter region were analyzed and identified with the PLACE database [32]. The cis-acting elements of the promoter region were visualized by GSDS 2.0 (http://gsds.cbi.pku.edu.cn/).

Interaction Network of B. oleracea U-Box Protein
The amino acid sequences of 60 U-box genes in A. thaliana were extracted and the interaction network was constructed by the STRING (https://string-db.org/) program [33]. Based on the relationship between orthologous genes of B. oleracea and A. thaliana, the interaction of the U-box gene was presumed.

Gene Ontology Annotation and KEGG Pathway Enrichment of B. oleracea U-Box Gene
All B. oleracea U-box genes were mapped to Gene Ontology (GO, http://www.geneontology.org/) and the GOseq method was selected to calculate the number of genes for each GO term. The KOBAS online database was used to analyze the metabolic and signal transduction pathways involved in the U-box gene [34], and the Kyoto Encyclopedia of Genes and Genomes (KEGG, https://www.genome.jp/kegg/) was then used to verify its specific metabolic pathways.

Identification, Classification, and Phylogenetic Analysis of U-Box Gene Family
We used two methods (retrieval from the BRAD and Bolbase databases, and the HMMER search tool) to retrieve and identify members of the B. oleracea U-box gene family; the results were consistent. Four protein sequences (Bol019081, Bol018658, Bol011584, and Bol045810) that clearly did not contain the U-box structure were deleted, and finally 99 B. oleracea U-box genes were screened (Table 1). At the same time, 60 AtPUB genes homologous to BoPUB were obtained (Table 1), in which AT3G10370.1 and AT3G10370.2 indicate different transcripts of the A. thaliana AtPUB gene. Based on the results of transcriptome sequencing analysis in our laboratory, we matched the retrieved BoPUB genes with the transcriptome sequencing gene by sequence alignment and further confirmed that they were consistent ( Table 1). In addition, we retained the four transcriptome genes without a U-box structure excluded above for subsequent analyses.
Azevedo et al. clustered PUBs from A. thaliana based on conserved domains [35]. However, the domains of the 99 PUB genes in B. oleracea are somewhat complicated, so we considered combining phylogenetic trees and domains as a basis for grouping. According to the evolutionary relationship of the phylogenetic tree, conserved motif, gene structure, etc., the retrieved BoPUB genes were divided into seven groups (I-VII), and named BoPUB 1-BoPUB 99 ( Table 1). The data showed that most A. thaliana genes have three homologous genes in B. oleracea, such as AtPUB5, 9,22,24,29,36,46, and 56. AtPUB17 has four homologous genes in B. oleracea. This indicates that the U-box gene expansion and evolution between B. oleracea and A. thaliana are closely related.
In order to better understand and verify the phylogenetic relationships between members of the U-box gene family in B. oleracea and A. thaliana, 103 BoU-box and 60 AtU-box genes were constructed using MEGA X software with the maximum likelihood method and divided into seven groups ( Figure 1). The classification results were consistent with those in Table 1. Groups I, II, III, and IV of B. oleracea genes in Figure 3A are also clustered together in Figure 1 and included in each group. The clustering results of Figure 1 are consistent with those shown in Figure 3A; the slight difference is that the three groups (V, VI, and VII) come together to form a large branch ( Figure 1). This indicates that the classification and naming of members of the B. oleracea U-box gene family in Table 1 are reliable. Careful observation of Figure 1 shows that each A. thaliana U-box gene produces at least one copy of the B. oleracea U-box gene and up to four copies.
Compared with the study by Wang et al. [10]. The largest group consists of proteins containing the ARM domain (I), but in different numbers (41 proteins in B. rapa and 28 proteins in B. oleracea). We classify proteins containing a large number of S_TKC (N-terminal domain of the eukaryotic serine threonine kinase) domains into a second broad category (II), and almost all U-box gene domains in B. oleracea contain a RING-finger. It is worth noting that Wang et al. identified two WD40-containing proteins from the evolution of AT1G04510.1 and AT2G33340.1 in B. rapa. However, three WD40 repeats were identified in B. oleracea, and AT2G33340.1 evolved to produce two WD repeats (VI). Additionally, AT3G07370.1 evolved a protein containing the TPR domain in B. rapa, while two proteins containing the TPR domain were produced in B. oleracea (VI).

Collinear Analysis within the Genome of B. oleracea and between the Genomes of B. oleracea and A. thaliana
In order to reveal the gene duplication of the U-box family members in B. oleracea, the genome differences between B. oleracea and A. thaliana were analyzed using the Circos software and the Perl language under the Linux system to complete collinearity analysis. Figure 2A is an intra-genomic collinearity analysis map based on the location information of nine chromosomes in the B. oleracea genome and the selected 24 pairs of tandem repeat genes; the red line indicates the collinearity U-box gene in the genome and the different colors of the background are collinear blocks from collinear gene pairs of the entire B. oleracea genome. Figure 2B shows 39 pairs of tandem repeats between B. oleracea and A. thaliana, and the red line indicates the collinearity of the U-box gene between these two genomes. Schedule S-2 provides detailed information on tandem duplication genes (Table S2).
The results of collinear analysis indicated that both orthologous copies and paralogous copies were present between B. oleracea and A. thaliana. Bo1g005580.1, Bo3g064620.1, Bo3g154110.1, Bo4g036170.1, and Bo5g117170.1 each have two paralogous copies in the B. oleracea genome (Figure 2A). These paralogous copies may also be an important factor in the occurrence of gene clusters and diversity distribution in the U-box family of B. oleracea [36]. Figure 2B shows that multiple AtPUB genes (AT1G04510.1, AT1G20780.1, AT2G33340.1, AT3G11840.1, AT3G46510.1, AT5G1540.1, and AT5G65920.1) in A. thaliana have two orthologous copies, and the AT2G35930.1 gene has three orthologous copies in B. oleracea. Although these orthologous copies are located in different species (B. oleracea and A. thaliana), the BoPUB gene was derived from a common ancestor shared with A. thaliana, so their functions are generally similar. In addition to these polygene duplications, there are a large number of single gene duplications. Although these single gene duplications often result in gene loss [37], the final retained genes are of great significance in evolution [38]. In short, extensive whole genome duplication (WGD) is considered to be the main driver of species diversity and the primary mechanism for acquiring new genes [39,40].
The vast majority of gene duplicates may be silenced over time, but a small number of genes purifying selection remain. Although these may rarely evolve new functions, the stochastic silencing of such genes also plays an important role in the origin of new species [41]. To gain a deeper understanding of the selection of these genes after replication, we calculated the ratio of the non-synonymous substitution rate (Ka) to the synonymous substitution rate (Ks) ( Table S3-1 and Table S3-2). It was found that the Ka/Ks value of 24 pairs of tandem genes in the genome ranged from 0.09 to 0.65; the Ka/Ks value of 39 pairs of tandem genes between genomes ranged from 0.04 to 0.50; and all Ka/Ks values were less than 1, indicating that members of the U-box family in the B. oleracea genome and between the B. oleracea and A. thaliana genomes underwent purified selection during long-term evolution [42]. Gene duplication events are very important for studying the evolutionary mechanisms of plant genes, and tandem repeats have a significant impact on the amplification and evolution of gene families in plant genomes [43].  Table 2. Batch CD-Search was used to analyze the gene structure online ( Figure 3C). From comprehensive analysis of this information, it was finally determined that the 103 U-box genes could be roughly divided into the I-VII groups. Additionally, four genes (Bo8g003030.1, Bo9g018240.1, Bo7g118020.1, and Bo8g100600.1) were found that obviously lacked a U-box domain in the analysis of the gene structure. UniProt (https://www.uniprot.org/) online analysis confirmed this result.
In general, U-box proteins that are clustered in the same type have similar conserved motifs, indicating that these proteins have similar functions. In combination with Figure 3A,B, the U-box protein motif in the same group of phylogenetic trees is highly similar and conserved, and almost all of them have Motif 4, 1, 7, and 2 (except for Group VII); this indicates that these four motifs play an important role in the BoU-box gene family ( Figure 3B). Referring to EnsemblePlants and UniProt databases, the structure information of the U-box gene family in B. oleracea was determined ( Figure 3C), including CDS (exon), Ubox, RING-Ubox_PUB, Arm, WD40, RING-Ubox_WDSUB1_like, S_TKc, and RING-Ubox_CHIP. According to the existing research findings, the ligase (E3) recruits E2 and the target protein is the last step of ubiquitin (Ub) transfer to the substrate. These three types of E3 enzymes are: HECT (homologous to E6-AP C-terminus), RING (Really Interesting New Gene), and U-box [35]. Protein interaction modules (e.g., WD40 repeats, TPR, and ARM repeat domains) are involved when U-box and RING act as scaffold molecules to recruit and localize substrates [44,45]. Therefore, the domain of the U-box family, contains many RING-finger and U-box domains, which indicates that the U-box family gene of B. oleracea has abundant E3 ligase characteristics. This conclusion is demonstrated in the published work of Trujillo et al. [6].
According to the NCBI query, it is clear that RING-Ubox_CHIP, RING-Ubox_PUB, and RING-Ubox_WDSUB1_like belong to the RING_Ubox Superfamily. RING_Ubox contains the RING-finger and U-box domains. RING-finger is a special type of zinc finger with 40 to 60 residues. The U-box is a modified form of the RING-finger domain, lacking metal chelate cysteine and histidine. This consists of three β-sheets and an α-helix, which mainly relied on salt and hydrogen bonds to stabilize their structures, rather than chelated metal ions [46,47]. The results of the graph ( Figure 3C) indicate that almost all of the identified genes contain the U-box and RING-Ubox_PUB domains, and the U-box domain coincides with the RING-Ubox_CHIP/WDSUB1_like PUB, suggesting that they may have synergistic effects. RING-Ubox_PUB is a common PUB protein in plants, also known as U-box domain protein. Most of the family's AtPUB is called PUB protein containing the ARM domain, a type of improved RING-finger, which is not only related to the regulation of cell death and defense, but also plays an important role in other plant-specific pathways, such as controlling self-incompatibility and acting in abiotic stress. The two genes (Bo3g061190.1 and Bo5g141580.1) containing the RING-Ubox_CHIP domain belong to the Group VI, which is homologous to A. thaliana AtCHIP (AT3G07370.1) from the protein interaction network, indicating that these three genes may be functionally similar. The RING-Ubox_WDSUB1_like structure exists in Group II. WDSUB1 is a non-characterized protein containing seven WD40 repeats and one SAM domain in addition to the U-box. Group VI has three genes in the WD40 domain, and the WD-40 repetitive sequence, also known as the WD or β-transferrin repeat, is a motif with approximately 40 amino acids that normally terminates in the Trp-Asp dipeptide [48]. WD40-repeat proteins are a large family found in all eukaryotes and are involved in a variety of functions, including signal transduction, transcriptional regulation, and apoptosis [49]. Several WD40-containing proteins in A. thaliana are key regulators of plant-specific developmental events. Group II also has an N-terminal domain of the eukaryotic serine threonine kinase, S_TKc. Protein kinases are enzymes belonging to a very broad protein family that play a role in a variety of cellular processes, including division, proliferation, apoptosis, and differentiation. They have a common conserved catalysis, and homology of the N-terminal domain with the USP family indicates that the N-terminal domain is involved in ATP binding [50][51][52].
In addition, there are many proteins with ARM and RING-Ubox domains in groups I, V, and VI, which have the same domain as ARC1 and belong to the PUB-ARM Superfamily. ARM repeats are characteristic repetitive amino acid sequences with a length of approximately 40 residues found in proteins, and proteins containing these sequences typically contain several tandem repeat copies [53,54]. Studies have shown that the interaction between the C-terminal arm repeat region ARM of the ARC1 protein with functional E3 ubiquitin ligase and the S-locus receptor kinase domain is part of the signal transduction pathway, leading to self-incompatibility pollen rejection [2,55]. It has been suggested that the A. thaliana ARM repeat protein is a member of the U-box E3 ubiquitin ligase family, and AtPUB-ARM can interact with the target substrate through its ARM repeat or UND (the N-terminal region of U-box). Samuel et al., using AtPUB13, 14, 45 (similar to ARC1 structure, containing UND, U-box, and ARM domain) and S-locus receptor kinase (SRK) for yeast double-hybrid experiments, proved that they can interact well [18], but whether all genes with UND structure participate in SI signal transduction still needs further study. Mudgil et al., identified 17 PUB genes with an UND structure through protein phylogenetic analysis of A. thaliana [56]; whether the homologous genes of these PUB genes in B. oleracea also have an UND structure also needs to be further studied. In addition, not all AtPUB-ARMs contain UND. In vitro ubiquitination assays also indicate that UND only promotes binding to ARM repeat regions and does not affect E3 ligase activity, so the AtPUB-ARM family may regulate the substrate due to the diversity of the ARM repeat region, thus controlling a series of different cellular processes [56]. Thus, it can be seen that the specific function of the UND domain in self-incompatibility needs further research. Trujillo proposed that in order to have a deeper understanding of the function of PUB, structural analysis would be an important means of understanding the complex interactions between PUB domains, especially for the unique PUB composed of UND, U-box, and ARM repeats [6]. It is hoped that the structural analysis of this section will provide valuable information for the in-depth study of the PUB gene in B. oleracea.

Chromosome Localization of B. oleracea U-Box Gene
As shown in Figure 4, a total of 102 BoPUB genes are mapped to C1-C9 chromosomes. The Bo00615s220.1 gene is a novel gene identified by transcriptome sequencing, and is located on the chromosome scaffold and cannot be visualized. Different groups of genes are represented by different colors according to the aforementioned grouping. The number of genes on each chromosome is unevenly distributed, and up to 21 genes are distributed on the C3 chromosome (Figure 4), indicating the complex diversity of the U-box family genes in B. oleracea. In addition, there are more gene clusters in C3, C4, C5, C7, C8, and C9 chromosomes, indicating that these gene amplification products belong to a common ancestor. The homologous gene relationship between C3, C1, and C5 is the most obvious, since there are many tandem duplication events ( Figure 2A and Table S2).
Genes 2019, 10, x FOR PEER REVIEW 7 of 26 that in order to have a deeper understanding of the function of PUB, structural analysis would be an important means of understanding the complex interactions between PUB domains, especially for the unique PUB composed of UND, U-box, and ARM repeats [6]. It is hoped that the structural analysis of this section will provide valuable information for the in-depth study of the PUB gene in B. oleracea.

Chromosome Localization of B. oleracea U-box Gene
As shown in Figure 4, a total of 102 BoPUB genes are mapped to C1-C9 chromosomes. The Bo00615s220.1 gene is a novel gene identified by transcriptome sequencing, and is located on the chromosome scaffold and cannot be visualized. Different groups of genes are represented by different colors according to the aforementioned grouping. The number of genes on each chromosome is unevenly distributed, and up to 21 genes are distributed on the C3 chromosome (Figure 4), indicating the complex diversity of the U-box family genes in B. oleracea. In addition, there are more gene clusters in C3, C4, C5, C7, C8, and C9 chromosomes, indicating that these gene amplification products belong to a common ancestor. The homologous gene relationship between C3, C1, and C5 is the most obvious, since there are many tandem duplication events (Figure 2A and Table S2).

Interaction Network of U-box Protein
To analyze the protein interactions between the U-box genes in B. oleracea, a protein interaction network was constructed with 60 orthologous genes in A. thaliana using the STRING software ( Figure  5). Each A. thaliana gene corresponds to several orthologous B. oleracea genes, and a total of 23 A. thaliana genes and 45 B. oleracea genes are involved in the interaction network. These genes have catalytic activity and are related to the ubiquitination of proteins. The central node gene (AT5G15400) is Mutant snc1-enhancing 3 (MUSE3), an E4 ubiquitin ligase involved in polyubiquitination of NLRs, which seems to function downstream of CPR1 to facilitate the polyubiquitination and degradation of SNC1 and RPS2 [57]. In addition, the MUSE3 ortholog yeast E4 ubiquitin ligase UFD2 can promote the extension of ubiquitin chains [58]. This indicates that MUSE3 plays an important role in the ubiquitination pathway of proteins. Bo2g012540.1 and Bo9g164140.1 are the orthologous genes of MUSE3, and the evolutionary relationship is close and the gene structure is similar (Figure 3, Group VII), so it is likely that both have similar functions to MUSE3.

Interaction Network of U-Box Protein
To analyze the protein interactions between the U-box genes in B. oleracea, a protein interaction network was constructed with 60 orthologous genes in A. thaliana using the STRING software ( Figure 5). Each A. thaliana gene corresponds to several orthologous B. oleracea genes, and a total of 23 A. thaliana genes and 45 B. oleracea genes are involved in the interaction network. These genes have catalytic activity and are related to the ubiquitination of proteins. The central node gene (AT5G15400) is Mutant snc1-enhancing 3 (MUSE3), an E4 ubiquitin ligase involved in polyubiquitination of NLRs, which seems to function downstream of CPR1 to facilitate the polyubiquitination and degradation of SNC1 and RPS2 [57]. In addition, the MUSE3 ortholog yeast E4 ubiquitin ligase UFD2 can promote the extension of ubiquitin chains [58]. This indicates that MUSE3 plays an important role in the ubiquitination pathway of proteins. Bo2g012540.1 and Bo9g164140.1 are the orthologous genes of MUSE3, and the evolutionary relationship is close and the gene structure is similar (Figure 3, Group VII), so it is likely that both have similar functions to MUSE3. In this network, AtCHIP (C-terminal of Hsp70 interacting protein) contains a U-box domain and three tetratricopeptide repeats (TPR) domain that mediates its interaction with Hsp90 and Hsp70 chaperones. As an E3 ubiquitin ligase of the protein phosphatase 2A subunit, it can alter the response of plants to abscisic acid treatment and can also play a key role under temperature stress conditions [59][60][61]. Through the interaction network, it can reliably reveal the homology, co-occurrence, and coexpression of B. oleracea and A. thaliana U-box gene. Several genes (AtCHIP, AtPUB 17/18/22/23/43/14, and CMPG2) formed a relatively close interaction and showed homology and co-expression ( Figure  5). Among these, the interaction of AtPUB18/ 22/ 23 can participate in a variety of biological processes. AtPUB18 may mediate EXO70B1 ubiquitination and participate in the regulation of ABA-mediated stomatal movement, and AtPUB22 and AtPUB23 coordinate control of drought signaling pathways [62,63]. The above results indicate that a variety of AtPUB proteins can individually affect or interact with each other to regulate related proteins in different signaling pathways, which provides a strong rationale for the study of related pathways involved in the BoPUB genes.

Gene Ontology Annotation and KEGG Pathway Enrichment Analysis of U-box genes in B. oleracea
After studying the interaction of U-box proteins, we continued the GO annotation and KEGG pathway enrichment analysis to gain a deeper understanding and validation of the biological functions of these genes ( Figure 6, Table S4, and Table 3). The results showed that a total of 1689 GO terms were significantly enriched, including biological processes, cell components, and molecular functions ( Figure 6). Among these, biological processes are enriched with the most items, mainly cell processes, metabolic regulation, stimulation reactions, signal transmission, etc. The most common cellular components include cells, organelles, and cell membranes. The most abundant molecular functions are catalytic activity and molecular binding. The GO annotation is consistent with the results of the protein interaction network analysis (Table S4 and Figure 5). In this network, AtCHIP (C-terminal of Hsp70 interacting protein) contains a U-box domain and three tetratricopeptide repeats (TPR) domain that mediates its interaction with Hsp90 and Hsp70 chaperones. As an E3 ubiquitin ligase of the protein phosphatase 2A subunit, it can alter the response of plants to abscisic acid treatment and can also play a key role under temperature stress conditions [59][60][61]. Through the interaction network, it can reliably reveal the homology, co-occurrence, and co-expression of B. oleracea and A. thaliana U-box gene. Several genes (AtCHIP, AtPUB 17/18/22/23/43/14, and CMPG2) formed a relatively close interaction and showed homology and co-expression ( Figure 5). Among these, the interaction of AtPUB18/ 22/ 23 can participate in a variety of biological processes. AtPUB18 may mediate EXO70B1 ubiquitination and participate in the regulation of ABA-mediated stomatal movement, and AtPUB22 and AtPUB23 coordinate control of drought signaling pathways [62,63]. The above results indicate that a variety of AtPUB proteins can individually affect or interact with each other to regulate related proteins in different signaling pathways, which provides a strong rationale for the study of related pathways involved in the BoPUB genes.

Gene Ontology Annotation and KEGG Pathway Enrichment Analysis of U-Box genes in B. oleracea
After studying the interaction of U-box proteins, we continued the GO annotation and KEGG pathway enrichment analysis to gain a deeper understanding and validation of the biological functions of these genes ( Figure 6, Table S4, and Table 3). The results showed that a total of 1689 GO terms were significantly enriched, including biological processes, cell components, and molecular functions ( Figure 6). Among these, biological processes are enriched with the most items, mainly cell processes, metabolic regulation, stimulation reactions, signal transmission, etc. The most common cellular components include cells, organelles, and cell membranes. The most abundant molecular functions are catalytic activity and molecular binding. The GO annotation is consistent with the results of the protein interaction network analysis (Table S4 and Figure 5).

SNEV, UBOX4, hPSO4
Bo4g043800.  The metabolic or signal transduction pathways involved in U-box genes were analyzed online using the KOBAS database, and the specific metabolic pathways were verified by KEGG. It was found that some genes are mainly involved in three pathways: ubiquitin-mediated proteolysis, splicing, and protein processing in the endoplasmic reticulum (Table 3).

Expression Analysis and RT-qPCR
In order to understand the expression of the U-box whole gene family members of B. oleracea in different pollination conditions and times, the expression levels of these 103 genes were quantified by FPKM (fragments per kilobase of transcript per million fragments mapped) values of RNA-seq data [64]. The heat map was drawn according to the log2 (FPKM) values (Table S5) of pollinated 0 min, CP-15/30/60 min, and SP-15/30/60 min, and clustered into A-F based on the expression profile characteristics of these genes (Figure 7 and Table S5).  The metabolic or signal transduction pathways involved in U-box genes were analyzed online using the KOBAS database, and the specific metabolic pathways were verified by KEGG. It was found that some genes are mainly involved in three pathways: ubiquitin-mediated proteolysis, splicing, and protein processing in the endoplasmic reticulum (Table 3).

Expression Analysis and RT-qPCR
In order to understand the expression of the U-box whole gene family members of B. oleracea in different pollination conditions and times, the expression levels of these 103 genes were quantified by FPKM (fragments per kilobase of transcript per million fragments mapped) values of RNA-seq data [64]. The heat map was drawn according to the log 2 (FPKM) values (Table S5) of pollinated 0 min, CP-15/30/60 min, and SP-15/30/60 min, and clustered into A-F based on the expression profile characteristics of these genes (Figure 7 and Table S5).  (Table S5), and the color code represents the transformation value of log2(FPKM). Bright green represents low expression and dark red represents high expression. According to the trend of gene expression in different pollination conditions and times, the heat map is divided into six groups (A-F).
The five genes contained were significantly up-regulated after CP for 30 min ( Figure 7A). As show in Figure 7B, with the exception of Bo4g025080.1, which was expressed at a higher level after CP-60 min, most were expressed at 0 min pollination. Figure 7C shows that most genes were significantly expressed after SP-15 min, and decreased after 30 min and 60 min, and weakly expressed at each time point under the condition of CP. Most of the genes were highly expressed after SP-15 and SP-30 min or CP-15 and CP-30 min, but the expression level was extremely low at 0 min  (Table S5), and the color code represents the transformation value of log 2 (FPKM). Bright green represents low expression and dark red represents high expression. According to the trend of gene expression in different pollination conditions and times, the heat map is divided into six groups (A-F).
The five genes contained were significantly up-regulated after CP for 30 min ( Figure 7A). As show in Figure 7B, with the exception of Bo4g025080.1, which was expressed at a higher level after CP-60 min, most were expressed at 0 min pollination. Figure 7C shows that most genes were significantly expressed after SP-15 min, and decreased after 30 min and 60 min, and weakly expressed at each time point under the condition of CP. Most of the genes were highly expressed after SP-15 and SP-30 min or CP-15 and CP-30 min, but the expression level was extremely low at 0 min pollination ( Figure 7D). The results of Figure 7E,F are more intuitive; the genes in Figure 7E were significantly expressed at SP-60 min, while the expression in Figure 7F is extremely high after CP-60 min.
Gene expression is temporally and spatially specific. Genes with significant differences in expression levels under different conditions are called differentially expressed genes (DEG). A total of 36 differentially expressed genes were screened out using fold-change ≥2 and FDR <0.01 as screening criteria in the process of differentially expressed gene detection [29]. The expression levels of 36 differentially expressed genes were detected by RT-qPCR and are plotted as a line graph to verify the up-and -down-regulation of pollination at 0 min, CP-15/30/60 min, and SP-15/30/60 min (Figure 8 and Table S1). The abscissa in Figure 8 indicates the time after pollination, the ordinate indicates the relative expression level, and the error bars indicate the standard deviation (±SD) of the three biological replicates under different pollination conditions. RT-qPCR results showed that some genes were significantly up-regulated after SP-15 min or SP-30 min, such as Bo1g147370.1, Bo3g029820.1, Bo4g125390.1, Bo5g052370.1, and Bo9g130630.1. It is notable that the up-regulation of these genes under the condition of CP is not obvious, indicating that these genes are highly likely to be involved in the SI-related signaling transduction pathways in B. oleracea. We further analyzed the relationship between these genes and the results of bioinformatics analysis, and found that the conserved motifs and structures of these genes are similar (Figure 3). Additionally, Bo1g147370.1, Bo3g029820.1, and Bo4g125390.1 formed a relatively close interaction and showed homology and co-expression ( Figure 5). These proteins have ubiquitin-protein transferase activities [56]. Bo1g147370.1 may be involved in the abscisic acid-mediated signaling pathway [18], while Bo3g029820.1 and Bo4g125390.1 coordinate control of drought signaling pathways [62]. Interestingly, from the chromosomal localization of highly significant differentially expressed genes (Figure 4), such genes are mostly located at C1, C3, and C5. Combined with tandem replication events (Figure 2A and Table S2), there is significant gene duplication among the three chromosomes. Whether there is a connection between the genes on these chromosomes is worthy of consideration. expression ( Figure 5). These proteins have ubiquitin-protein transferase activities [56]. Bo1g147370.1 may be involved in the abscisic acid-mediated signaling pathway [18], while Bo3g029820.1 and Bo4g125390.1 coordinate control of drought signaling pathways [62]. Interestingly, from the chromosomal localization of highly significant differentially expressed genes (Figure 4), such genes are mostly located at C1, C3, and C5. Combined with tandem replication events (Figure 2A and Table  S2), there is significant gene duplication among the three chromosomes. Whether there is a connection between the genes on these chromosomes is worthy of consideration.

Cis-Regulatory Element Analysis in the Promoter Regions of U-Box Differentially Expressed Genes
To further elucidate the differential expression results of 36 genes and their putative role in the pollination process, PLACE and GSDS 2.0 were used to analyze and visualize the cis-regulatory elements in the promoter region [32]. The main regulatory elements include: hormone-related regulatory elements (ABRE, ERE), pollen-specific cis-regulators (POLLEN_lat52), dehydration-related action elements (MYC), and pressure and defense related regulatory elements (W-box) (Figure 9 and Table S6).
(SA)-induced WRKY DNA binding proteins. WRKY proteins are DNA-binding proteins that recognize the W-box elements found in the promoters of a large number of plant defense-related genes [71]. Arabidopsis transcription factor AtWRKY18 enhances developmental regulation of plant defense responses [72]. Analysis of the above cis-acting elements indicated that these differentially expressed B. oleracea U-box genes can not only respond to different biotic or abiotic stresses, but are also closely related to the pollen development of B. oleracea.  The ABRE regulatory element can transcriptionally regulate the abscisic acid (ABA) response gene [65], and a novel cis-acting element ABRE_ERD1 has been shown to be up-regulated in response to water stress and yellowing in A. thaliana, indicating that EDR1 elements play an important role in dehydration induction and stress [66]. MYC_CONSENSUSAT is the promoter region recognition site of the dehydration reaction gene rd22 in A. thaliana, and MYC is a cis-acting element in drought-induced rd22 gene expression [67]. Follow-up studies have shown that the AtMYC2 protein plays the role of transcriptional activator in ABA inducible gene expression under plant drought stress [68]. These differentially expressed U-box genes contain a large number of pollen-specific cis-acting elements of POLLEN1LELAT52 (AGAAA) (Figure 9). Studies have shown that this regulatory element is involved in the precise pollen development and tissue-specific expression of tomato lat52 and can affect the transcriptional activation of lat52 [69]. LeMAN5 endo-β-mannanase, which is related to pollen germination and pollen tube elongation in tomato, has a 5'-upstream region containing four copies of POLLEN1LELAT52, which can be used for pollen fertility control [70]. "W-box" found in promoter of A. thaliana NPR1 gene, they were recognized specifically by salicylic acid (SA)-induced WRKY DNA binding proteins. WRKY proteins are DNA-binding proteins that recognize the W-box elements found in the promoters of a large number of plant defense-related genes [71]. Arabidopsis transcription factor AtWRKY18 enhances developmental regulation of plant defense responses [72]. Analysis of the above cis-acting elements indicated that these differentially expressed B. oleracea U-box genes can not only respond to different biotic or abiotic stresses, but are also closely related to the pollen development of B. oleracea.

Conclusions
In this study, four PUB family genes without U-box characteristics were excluded from the B. oleracea genome database, and a total of 99 BoPUB genes were identified. Through the comprehensive analysis of the phylogeny, conservative motifs, and gene structure of the BoPUB genes, verified by their homologous evolutionary relationship and structure of the AtPUB gene, these genes were grouped and named. At the same time, the evolutionary relationship, gene expansion, and functional diversity of the PUB gene family in B. oleracea were also revealed. Chromosome localization and collinear analysis showed that members of the U-box family were purified and selected during long-term evolution, and these genes may be products of amplification from the same ancestral gene. Protein interaction network analysis indicated that the BoPUB protein is widely involved in the ubiquitination pathway and regulates downstream proteins, which was verified by GO annotation and KEGG pathway enrichment analysis. We analyzed the expression characteristics of the BoPUB genes under different pollination conditions, and the related genes in the U-box gene family of B. oleracea likely to be involved in SI reaction were identified by RT-qPCR. A large number of cis-acting elements related to pollen specificity, hormone regulation, and dehydration reactions were found in the promoters of differentially expressed genes, which further proved that these genes not only play an important role in the regulation of plant development and abiotic stress, but are also involved in pollination. These genes should be selected for follow-up study of SI reaction in B. oleracea. This study not only systematically analyzed the characteristics and evolution of the U-box family of B. oleracea, but also provided candidate genes for SI reaction and signal transduction in B. oleracea, thereby providing a reference for revealing how U-box genes participate in SI.

Supplementary Materials:
The following are available online at http://www.mdpi.com/2073-4425/10/12/1000/s1. Table S1: Primers for RT-qPCR analysis of selected U-box genes in B. oleracea, Table S2: Tandem replication gene pairs of the U-box family members in the B. oleracea genome and between the B. oleracea and A. thaliana genomes,  Table S4: GO Annotation Analysis of U-box genes in B. oleracea, Table S5: Expression levels of the U-box gene in B. oleracea after different pollination conditions and time, Table S6: Cis-Acting elements analysis in the 1.5 kb promoters of U-box genes in B. oleracea.