Genome sequence and annotation data of B. napus cv. ‘Darmor-bzh’ (v_5.0) were obtained from the Genoscope database (http://www.genoscope.cns.fr/brassicanapus/) (Chalhoub et al., 2014). B. rapa ‘Chiifu’ (v3.0) datasets were obtained from the database (https://bigd.big.ac.cn/gwh/Assembly/134/show) (Zhang et al., 2018) and B. oleracea ‘HDEM’ (broccoli) datasets were acquired from Bolbase database (http://www.ocri-genomics.org/bolbase/index.html) (Belser et al., 2018). The ubiquitin-conjugating enzyme (UQ_con) domain (PF00179) annotation file obtained from the Pfam database (https://www.ebi.ac.uk/interpro/) was applied as queries to perform HMM search against protein sequences of annotated genes using HMMER version 3.1 (http://hmmer.org/) (Eddy, 2011) with E-value cutoff of 1e-5 (Finn et al., 2011). The amino acid sequences of UBCs (predicted above) were utilized as queries to carry out BLASTP searches against the full-length protein sequences of 48 AtUBCs from The Arabidopsis Information Resource (TAIR) database (http://www.arabidopsis.org/) and 39 OsUBCs from Rice Genome Annotation Project (http://rice.plantbiology.msu.edu/downloads_gad.shtml) with E-value < 1e-5. The putative BnUBCs with hits of both AtUBCs and OsUBCs were deployed for further confirming the existence of the UQ_con domain using the Pfam (https://www.ebi.ac.uk/interpro/), SMART (http://smart.embl-heidelberg.de/) and CDD (https://www.ncbi.nlm.nih.gov/Structure/cdd) databases. The molecular weight (MW), isoelectric point (PI), instability index, aliphatic index and grand average of hydropathicity (GRAVY) of UBC proteins, were calculated using the ProtParam tool ExPASy (https://web.expasy.org/protparam/). The subcellular location of BnUBCs was predicted using Plant-mPLoc (http://www.csbio.sjtu.edu.cn/bioinf/plant-multi/#) (Chou and Shen, 2010).

The chromosomal location, coding sequence (CDS) and amino acid sequences were determined based on the genome annotation of B. napus in the Genoscope database (http://www.genoscope.cns.fr/brassicanapus/). The chromosomal distribution of BnUBCs was graphically represented using the RIdeogram package of the R software (https://github.com/TickingClock1992/RIdeogram) (Hao et al., 2020). Multiple alignments of the BnUBC protein sequences were performed using CLUSTAL v2.1, with default parameters. A schematic diagram of BnUBCs gene structure was conducted by Gene Structure Display Server 2.0 (http://gsds.cbi.pku.edu.cn/) (Hu et al., 2015). Conserved motifs in the BnUBC proteins were investigated using the online MEME server (http://meme-suite.org/tools/meme) with the following parameters settings: the maximum number of motifs, 10; minimum width of motifs, 6; maximum width of motifs, 100 aa; and E-value < 1e-10 (Bailey et al., 2009). Moreover, the 2-kb upstream sequence of the BnUBCs gene sequences was extracted and submitted to PlantCARE (Lescot et al., 2002) for the detection of the cis-elements.

To better understand the evolutionary relationships among the UBC of B. rapa, B. oleracea, B. napus and A. thaliana, a phylogenetic analysis was performed. A phylogenetic tree was constructed applying MEGA 5.2 software (Tamura et al., 2011), based on the neighbor-joining (NJ) method and 1,000 bootstrap replications. The online software Interactive Tree Of Life (iTOL, http://itol.embl.de/) (Letunic and Bork, 2016) was used to decorate this phylogenetic tree. Genes orthologous between B. napus and its ancestors (B. rapa, B. oleracea, and A. thaliana) were identified by BLASTn searches of their CDSs based on two criteria: coverage of sequence length > 80%, and identity of aligned regions > 80% (Kong et al., 2013). In addition, DupGen_finder (https://github.com/qiao-xin/DupGen_finder) (Qiao et al., 2019) was employed to identify the modes of gene duplication of UBC paralogous genes in A. thaliana, B. napus, B. rapa and B. oleracea. The syntenic relationship among these paralogs was presented using the fmsb package of the R software.

The KaKs calculator (Wang et al., 2010) was applied to estimate the divergence between pairwise nonsynonymous substitution rates (Ka) and synonymous substitution rates (Ks). The evolutionary constraint (Ka/Ks) of UBC orthologous genes between B. napus and the other three species (A. thaliana, B. rapa and B. oleracea) was calculated according to their CDSs. In addition, to minimize errors, only the gene pairs with Ks < 1 were remained for further analysis (Guo et al., 2017; Wang T, et al., 2020).

The psRNATarget Server (Dai et al., 2018) was used to predict the miRNAs potentially targeting BnUBCs. The gene sequences of BnUBCs were submitted as candidates to search against the available sequences of reference B. napus miRNA with default arguments. A network of interactions between the predicted miRNAs and their target UBC genes in B. napus was visualized using the Cytoscape software (Shannon et al., 2003).

The publicly available RNA-seq datasets of five different B. napus tissues (callus, root, leaf, bud and silique) were collected from the NCBI SRA database (accession no. SRP136038) (Yao et al., 2020). The expression levels of BnUBCs were quantified based on their fragments per kilobase of exon per million reads mapped (FPKM) values using Cufflinks with default parameters (Trapnell et al., 2012). Furthermore, the expression levels of BnUBC genes were used to construct clustered heatmaps with TBtools (Chen et al., 2020). Then, GO enrichment analysis of BnUBC genes was performed using the clusterProfiler package of the R software.

The agronomic traits of 324 worldwide accessions comprising a natural population of B. napus were applied to detect the natural sequence variations of BnUBCs (Tang, 2019; Xie et al., 2022). Single nucleotide polymorphisms (SNPs) in BnUBCs were annotated using the SnpEff software (Cingolani et al., 2012). Moreover, agricultural traits involving the primary flowering time (PFT), full flowering time (FFT1), final flowering time (FFT2), plant height (PH), thousand seed weight (TSW), main inflorescence silique density (MISD), main inflorescence silique number (MISN), oil content (OC), protein content (PC) and main inflorescence length (MIL) were used to conduct associated mapping (Tang, 2019). The EMMAX software was applied to perform family-based association mapping analysis with a mixed linear model (Kang et al., 2010). Then, the visualization of haplotype blocks and linkage disequilibrium was conducted using the LDBlockShow software (Dong et al., 2021b). The interaction networks of B. napus proteins were obtained from the STRING database (Damian et al., 2020).

To identify the members in the UBC family, the annotation file of UBC domain obtained from the Pfam database (http://pfam.xfam.org/) was applied as a query for searching against the protein dataset of B. napus. The peptides of assumed BnUBCs showing hits of both AtUBCs and OsUBCs were further employed to confirm the presence of the UBC domains by searching in the Pfam, SMART and CDD databases. In total, 200 BnUBCs were identified in B. napus ( Supplementary Table 1 ). The full-length transcript sequences of BnUBCs ranged from 168 bp (BnaA08g09170D) to 5,656 bp (BnaC07g05960D) with the deduced amino acid sequences varying from 55 aa (BnaA08g09170D) to 1,649 aa (BnaA07g03330D) ( Supplementary Table 1 ). The predicted MW of the 200 deduced BnUBC proteins ranged from 6.10 kDa (BnaA08g09170D) to 184.41 kDa (BnaA07g03330D), and their GRAVY and pIs ranged from -0.904 (BnaC04g10080D) to 0.112 (BnaA08g13930D) and from 4.19 (BnaA04g11090D) to 9.57 (BnaA09g09660D), respectively ( Supplementary Table 1 ). As previously reported in other plant species (Die et al., 2018; Zhou et al., 2018), a wide range of pIs suggested that BnUBCs function in different subcellular environments. Furthermore, the instability index of 38 of the 200 BnUBC proteins was less than 39, and these proteins were classified as stable. In addition, the prediction of their subcellular localization indicated that 172 BnUBCs were located in the nucleus, 27 in the cytoplasm and the other one in the cell membrane. In addition, to explore the evolutionary relationships among members of this family across Brassicaceae crops, 93 BraUBCs and 95 BolUBCs were identified in the reference genomes of B. rapa and B. oleracea respectively, following the same identification pipeline ( Supplementary Table 2 ). Details of these BraUBCs and BolUBCs were listed in Supplementary Table 2 .

The chromosomal positions of detected BnUBCs were drafted to corresponding chromosomes with the RIdeogram package of R software. The 200 BnUBCs were unevenly distributed across all 19 chromosomes with 97 and 103 genes located in the A_n and C_n subgenomes, respectively ( Figure 1 ). The A_n subgenome carried, on average 9.7 BnUBCs on 10 chromosomes; A10 had the lowest number of BnUBCs (6), while A03 and A05 harbored the highest of BnUBCs (11 each). The average number of BnUBCs in the C_n subgenome was 11.4, with the lowest number (5) on C09 and the highest number on C03. Hence, no biased tendency was observed between the two subgenomes. Furthermore, on each chromosome, BnUBCs were unevenly distributed, especially, since the genes on A06 and A08 were located at the ends.

To explore the evolutionary relationships among the UBC genes in A. thaliana, B. rapa, B. oleracea, and B. napus, a tree of phylogenetic was constructed based on the multiple sequence alignments of 48 AtUBCs, 93 BraUBCs, 95 BolUBCs, and 200 BnUBCs. All 200 BnUBCs could be categorized into 18 subfamilies and the number of groups was consistent with a previous report (Jia et al., 2019) ( Figure 2 ). Among them, Group XVIII possessed the highest number (37) of the BnUBCss. Furthermore, group III (two BnUBCs) and group VIII (ten BnUBCs) members were involved in RUB and SUMO conjugation pathways, respectively, and other three subgroups (XVI, XV, VI) contained the members lacking the Cys active-site (UEVs).

Concerning the evolution of multigene families, the variety of gene structures provided an important resource (Ischoff et al., 2002; Yunpeng et al., 2016; Pellicer et al., 2018). For exploring the multiplicity of the gene structure of BnUBCs among the different groups, the intron-exon structure of BnUBCs was compared based on the phylogenetic relationships of these genes ( Figure 3A , Supplementary Table 3 ). The results indicated that the gene architecture of BnUBCs was relatively complex ( Figure 3B ). Only five genes (BnaA10g10530D, BnaA05g17560D, BnaC04g10080D, BnaC07g46220D and BnaA08g09170D) contained a single intron, while the others harbored multiple introns. BnaA07g03330D contained the most introns and exons. Besides, 33 of 200 BnUBCs lacked untranslated regions (UTR) at both ends or one end because of inaccurate annotation. The number of introns among the BnUBCs ranged from 1 to 26, and the length of introns ranged from 14 to 6,632 bp. Moreover, exon numbers inBnUBCs also varied widely across distinct subfamilies, ranging from 2 to 26. However, the majority of the members within the same subfamily tended to show similar intron-exon distribution patterns. For example, members of Group X generally possessed five exons, whereas those of Group XV, harbored eight exons. Therefore, the organization of introns and exons presented valuable evidence for the phylogenetic relationship among members of this gene family.

Furthermore, a total of 10 conserved motifs were recognized in 200 BnUBC genes ( Figure 3C ). Conserved motifs ranged in length from 11 to 74 aa, and among 10 motifs motif 3 as well as motif 1 were the most abundant. Moreover, the motif 1-4 were significantly matched to the UQ_con domain (Pfam: PF00179). Furthermore, a motif sequence search in the CDD database revealed that motifs 5 and 6 belonged to UBCc family. In addition, motifs 7-8 and motif 10 were matched to UBQ-conjugating_enzyme/RWD based on the search in the InterPro database. Motif composition varied across distinct subfamilies, whereas the conserved character of motif distributions in the same subcategory was similar, which emphasized their phylogenetic relationship. For example, all Group IX members possessed motif 1, motif 4-6 and motif 10, whereas in Group X, members generally contained motifs 1-6 and motif 10. In addition, three structural properties were detected and refined in BnUBC proteins ( Supplementary Figure 1 ).

The cis-acting elements in the promoter region of a gene can modulate the initiation and efficiency of gene expression through the binding of specific transcription factors, associated with plant development, growth and stress response (Ibraheem et al., 2010). To better understand the potential functions and transcriptional regulation of BnUBCs, we isolated the sequences 2-kb upstream of the BnUBC CDSs from the B. napus genome sequence to search for cis-elements ( Supplementary Table 4 ) (Chalhoub et al., 2014). A total of 115 functional cis-elements were detected in the promoter regions of BnUBCs ( Supplementary Figure 2 , Supplementary Table 5 ). Among these cis-elements, many were light-responsive and hormone-related related elements. Besides, majority of BnUBCs (197/200, 98.50%) contained both CAAT-box and TATA-box elements, which usually existed in eukaryon. Moreover, the MYB element, which participates in the regulation of phenylpropanoid secondary metabolism in plants was identified in 192 out of 200 BnUBCs. The MYC element was identified in 183 of the 200, (91.50%) BnUBCs, which were involved in the regulation of plant growth and development as well as resistance to environmental stress and secondary metabolite synthesis. Furthermore, many BnUBCs contained light-responsive elements such as Box 4 element (part of a conserved DNA module involved in the light response), G-box element (cis-acting regulatory element involved in the light response), GT1-motif (light responsive element), and TCT-motif (part of a light responsive element) ( Supplementary Figure 2 , Supplementary Table 5 ). Additionally, the stress-related elements, including the LTR element and TC-rich repeats were also detected in BnUBCs, suggesting that the BnUBCs played an important role in plant development and growth as well as in response to biotic and abiotic stresses.

Gene family expansion occurred mainly via the duplication of genes through whole-genome duplication (WGD), tandem duplication as well as segmental duplication events (Bodt et al., 2005; Mun et al., 2009). Moreover, the duplication of ubiquitin members through WGD events plays important role in plants (Hua et al., 2018). To further explore the evolutionary history of BnUBC gene family expansion we performed a syntenic comparison of genome sequences between A. thaliana, B. napus, B. rapa and B. oleracea. A. thaliana is the ancestor of Brassica species, whose structural genes have been identified and functionally annotated, therefore it is served as a prominent model system for the investigation of the evolutionary history of Brassica. homologous genes from inter-species comparison and paralogous genes from intra-species comparison were identified to trace the duplicated gene pairs of BnUBCs. We identified 268, 104, 86 and 29 paralogous UBC gene pairs of within B. napus, B. rapa, B. oleracea and A. thaliana, respectively ( Figure 4B ). Out of 286 paralogous UBC gene pairs in B. napus, 50 and 44 pairs were detected in the A_n subgenome, and C_n subgenome, respectively, while the remaining 174 pairs were identified across the A_n and C_n subgenomes. Moreover, we classified these paralogous gene pairs into five types: dispersed, proximal, tandem, WGD and transposed ( Table 1 , Supplementary Table 6 ). In B. napus, 195 genes, accounting for 97.5% of all BnUBCs, were derived from gene duplication, of which WGD (64.5%) and transposed (17.4%) gene duplication events were primarily responsible for gene family expansion ( Table 1 , Supplementary Table 6 ).

Furthermore, synteny comparison between B. napus and the other three species revealed that 180 BnUBCs were collinear to the orthologous genomic regions in one or more of the other three species, 164 orthologous gene pairs were identified between B. napus and B. rapa, 163 between B. napus and B. oleracea as well as 115 between B. napus and A. thaliana ( Figures 4A, B , Supplementary Table 7 ). It was obvious that majority of the BnUBCs were inherited from their progenitors. In addition, Brassica species experienced an extra whole-genome triplication (WGT) event after divergence from A. thaliana (Feng et al., 2014). Therefore, a single A. thaliana gene was inherited as three copies in B. oleracea as well as B. rapa and as six copies in B. napus, if no gene loss occurred after the WGT event. We found that only 25.9% of the AtUBCs corresponded to six copies in B. napus, suggesting that a substantial gene loss occurred during the process of polyploidization. To estimate the selection pressure on BnUBCs, the non-synonymous to synonymous substitution ratios (Ka/Ks) of BnUBC orthologous genes to BraUBCs, BolUBCs, and AtUBCs were calculated. The Ka/Ks less than one represents purifying selection, whereas Ka/Ks equal one indicates neutral evolution and Ka/Ks more than one represents positive selection (Nekrutenko, 2002). Thus, the Ka/Ks ratio can predict the pressure of selection on each duplicated pair throughout evolution as well as their divergence time. The Ka/Ks ratios of all the BnUBCs paralogous gene pairs, except two duplicated gene pairs (BnaA07g03430D/BnaC07g05830D, BnaC02g02910D/BnaC03g25550D) were less than one, suggesting that these genes were subjected to strong purifying selection ( Supplementary Table 8 ). Moreover, the duplicated gene pair BnaA07g03430D/BnaC07g05830D had Ka/Ks ratio considerably larger than 1, suggesting these genes were subjected to strong evolutionary pressure under positive selection. Furthermore, the Ka and Ks values of the identified orthologous gene pairs between A. thaliana and B. napus were calculated to estimate the selection pressure on orthologous gene pairs and the divergence time of the two species. The Ka/Ks ratios of all the orthologous gene pairs ranged from 0.001 to 0.25 with an average of 0.074, which indicated these genes were under purifying selection ( Supplementary Table 9 ). The orthologous gene pairs among A. thaliana and B. napus diverged around 16 million years ago (MYA) based on the mutational rate estimate, R = 1.5 × 10^-8 synonymous substitutions per site per year (Koch et al., 1999; Koch et al., 2001). This result was consistent with the previously estimated divergence time (14–24 MYA) between the A. thaliana and B. napus lineages (Koch et al., 2000; Cheung et al., 2009). In addition, the Ks values of the orthologous genes between B. napus and its progenitors were also calculated ( Supplementary Table 9 ). The Ka/Ks ratio of the orthologous gene pairs between B. napus and B. oleracea was significantly higher than that between B. napus and B. rapa, showing that genes in the C_n subgenome underwent relatively weaker selection pressure in B. napus during the process of evolution ( Figure 5 ). In addition, the Ka/Ks values of the orthologous gene pairs between B. napus and Arabidopsis were significantly lower than those of the orthologous gene pairs identified between B. napus and B. oleracea and between B. napus and B. rapa, indicating that orthologous gene pairs between B. napus and A. thaliana underwent intense purifying selection. Moreover, the sequences of their proteins might maintain more consistent characteristics during evolution ( Figure 5 ).

Many miRNAs were involved in regulating the function of their targeted genes in plants (Jones-Rhoades and Bartel, 2004; Yu et al., 2017; Song et al., 2019). To explore the miRNAs regulating BnUBCs, 106 BnUBC targeting genes of 71 putative miRNAs were identified, and their relationship network was constructed by Cytoscape software ( Figure 6 , Supplementary Table 10 ). It was found that BnaC06g03670D was targeted by the B. napus miRNA169 family based on an interaction network. A recent study showedthat miR169 regulated the function of the anaphase-promoting complex/cyclosome (APC/C), an essential ubiquitin-protein ligase, by targeting DUO POLLEN1 (DUO1). DUO1 upregulated the expression level of APC/C, which stimulates the production of miR159 (Zheng et al., 2011). Moreover, members of the miRNA156 family were found to regulate most of the BnUBCs. In a recent study, overexpression of miR156 inhibited gibberellic acid (GA)-induced and ubiquitination-mediated degradation of DELLA (Jerome Jeyakumar et al., 2020). In addition, the remaining BnUBCs were presumably targeted by members of the other 26 miRNA families ( Figure 6 ).

Several studies reported that members of the UBC family mediated crosstalk among diverse signaling pathways for various abiotic stress responses (Gao et al., 2017; Jue et al., 2018; Jia et al., 2019). Besides, UBCs were predicted to play an important role in protein and ion binding (Liu et al., 2019). To predict the putative functions of BnUBCs, GO enrichment analysis was performed. The GO terms were divided into three categories: biological process (BP), molecular function (MF), and cellular component (CC) ( Supplementary Table 11 ). The majority of enriched GO terms, such as postreplication repair, belonged to the BP category ( Supplementary Figure 3 ). The remaining GO terms were broadly associated with vegetative growth (negative regulation of flower development, root epidermal cell differentiation), response to stress (response to iron ion, cellular response to water deprivation and response to gibberellin), cell growth and protein process. The CC category included UBC13−MMS2 complex and the perinuclear region of the cytoplasm. Moreover, various MF were detected in this analysis, which included acid-amino acid ligase activity, SUMO transferase activity, ubiquitin protein ligase binding and endopeptidase activity. The results of GO enrichment analysis suggested that BnUBCs play crucial roles in plant development, growth, cell protein quality regulation and response to stress, consistent with the findings of previous studies (Dreher and Callis, 2007; Sadanandom et al., 2012; Jue et al., 2018).

Moreover, the expression profiles of BnUBC genes were examined across five major tissues (callus, leaf, root, bud and silique) based on the previously published RNA-seq datasets of B. napus (Yao et al., 2020). The expression levels of BnUBCs were estimated with FPKM and displayed as a heatmap ( Figure 7A ). A total of 38 BnUBCs showed relatively weak expression (FPKM < 1) and 24 BnUBCs were silience in any of the five tissues. The remaining BnUBCs showed high expression levels (FPKM ≥ 1), the majority of which expressed at a specific organ or differentially expressed in various tissues ( Figure 7B ). BnaC06g16140D (ortholog of AT3G57870) presented the highest expression in all tissues, suggesting its critical roles in plant growth. Genes showing low expression levels (FPKM <1) across all five tissues could be regarded as pseudogenes ( Figure 7 , Supplementary Table 12 ). Notably, in our study, most of BnUBCs exhibited relatively lower expression levels in the leaf than in other tissues ( Supplementary Table 12 ).

Furthermore, the expression divergence between the duplicated genes in B. napus investigated. For example, BnaA06g10330D was silence across all tissues, whereas its paralog BnaA03g31080D presented a high expression level in various tissues. Similarly, BnaA04g02190D expressed in all tissues, however, its paralogs BnaA01g13200D and BnaC01g24100D didn’t express in any tissue. All these results could provide a clue for the functional divergence between duplicated gene pairs in BnUBCs.

SNPs in a natural population with 324 B. napus accessions collected from worldwide countries were identified (Tang, 2019) ( Supplementary Table 13 ) for investigating the genetic variants of BnUBCs in germplasm. Each BnUBC gene contained 34 SNPs on average, which was close to the 36 SNPs in each gene across the entire genome. Taking genome size into account, the number of SNP in each kilobase (kb) among the BnUBCs was calculated (17 SNPs/kb), whereas 13 SNPs/kb for the genes across entire genome. BnUBCs in the A_n subgenome contained slightly higher SNP density (20 SNPs/kb) than that in the C_n subgenome (14 SNPs/kb). Moreover, variations in SNP number between some paralogous gene pairs of BnUBCs were observed. For example, BnaA03g11080D possessed 65 SNPs, however, its paralog BnaC02g12440D contained no SNP. Moreover, the paralogs BnaA07g03330D/BnaC03g46750D, contained 235 and 33 SNPs, respectively. Based on the paired t-test, no significant difference was detected in SNP density between the BnUBC paralogs. Additionally, SNP annotation presented that 1,671 SNPs existed in the exon regions of BnUBCs and 600 SNPs lead to missense mutations.

Ubiquitination is a well-characterized post-translational that regulates plant growth and developmental processes contributing to diverse phenotypes, and affect agronomic traits (Dreher and Callis, 2007; Vierstra, 2009; Sadanandom et al., 2012; Linden and Callis, 2020; Yu and Hua, 2022). Association mapping analysis was performed to explore the impact of BnUBCs on agronomic phenotype. In total, 96 SNPs in 33 BnUBC genes ( Supplementary Table 13 ) were significantly associated with the investigated agronomic traits (p < 0.001). Moreover, BnaC02g25260D was significantly associated with plant height ( Figure 8A, B ). Furthermore, the two genotypes of accessions were divided according to the most significantly associated SNP of BnaC02g25260D, and the plant height between the two genotypes was statistically significant according to the t-test (p < 2.7 e–4) ( Figure 8C ). The protein interaction network of BnaC02g25260D was obtained from STRING (Damian et al., 2020) ( Figure 8D ) and the interacted genes were enriched in the regulation of embryonic development (GO:0045995), positive regulation of cell size (GO:0045793), cellulose biosynthetic process (GO:0030244), circadian rhythm (GO:0007623), positive regulation of flower development (GO:0009911), detection of visible light (GO:0009584), leaf development (GO:0048366), regulation of hormone levels (GO:0010817) and chlorophyll biosynthetic process (GO:0015995) ( Figure 8E ). Overall, the above categories were related to plant development and growth, the regulation of hormones and oil content, which affected the plant height and the cell size in seed formation (Wang X, et al., 2020; Dong et al., 2021a). Furthermore, SNPs distributed in BnaAnng34240D were significantly associated with flowering time and silique density ( Supplementary Figure 4A-D, G-H ). The protein interaction networks of BnaAnng34240D were also obtained from STRING ( Supplementary Figure 5A ), and its interacted proteins enriched in the protein modification process, regulation of transcription and so on ( Supplementary Figure 5B ). Besides, BnaA07g28330D was significantly associated with TSW, which affect the crop yield ( Supplementary Figure 4E, F ). Furthermore, we found that BnaA07g28330D exhibited a higher expression level in the silique of high oil content accession than in that in the silique of low oil content accession (https://bnaomics.ocri-genomics.net/tools/exp-view/index.php) (Tang, 2019).