Marker Metabolite‐Based Multi‐Omics Analysis Identifies New Loci Controlling Thousand Seed Weight in Brassica Napus
Long Li, Zhitao Tian, Xiaowei Wu, Zengdong Tan, Xu Han, Yuyan Xiang, Jie Chen, Hu Zhao, Wei Chen, Liang Guo, Xuan Yao
TL;DR
This study uses multi-omics data to uncover new genetic factors influencing seed weight in Brassica napus, identifying a key gene that negatively affects seed weight.
Contribution
The first marker metabolite-based multi-omics analysis to unravel the genetic basis of thousand seed weight in Brassica napus.
Findings
137 marker metabolites significantly correlate with thousand seed weight (TSW) in Brassica napus.
BnaA03.TGA6 is identified as a novel gene negatively influencing TSW through activation of E3 ubiquitin-protein ligase DA2.
Abstract
The thousand seed weight (TSW) of Brassica napus (B. napus) is a crucial agronomic trait, and metabolites serve as a bridge between genotype and phenotype. Leveraging previously released metabolome data from a 388 B. napus population, 137 marker metabolites significantly correlated with TSW are identified, and 21 markers significantly correlate with both TSW and seed oil content. To delve deeper into the relationships between selected marker metabolites and phenotypes, a comprehensive multi‐omics analysis is conducted, which unveils 734 metabolite‐QTLs, 5,225 expression‐QTLs, and 9,077 transcriptome‐wide significantly associated genes. Employing weighted correlation network analysis, a triple‐link network is constructed involving 133 metabolites, 731 QTLs, and 3,321 genes. The multi‐omics analysis highlights the impact of the transcriptional factor TGACG motif‐binding factor 6…
Genes, proteins, chemicals, diseases, species, mutations and cell lines named across the full text — each resolved to its canonical identifier and authoritative record.
Click any figure to enlarge with its caption.
Figure 1
Figure 2
Figure 3
Figure 4
Figure 5
Figure 6
Figure 7| Metabolite | r_SOC | r_TSW | Name | Class |
|---|---|---|---|---|
| mr1058 | −0.36 | −0.15 | Tricin | Flavonoids |
| mr1076 | −0.32 | −0.15 | Chrysoeriol 5‐O‐hexoside | Flavonoids |
| mr1077 | −0.41 | −0.14 | Chrysoeriol C‐hexoside | Flavonoids |
| mr1082 | −0.37 | −0.14 | Selgin 5‐O‐hexoside | Flavonoids |
| mr1083 | −0.40 | −0.14 | Tricin 5‐O‐hexoside | Flavonoids |
| mr1089 | −0.26 | −0.16 | Chrysoeriol O‐malonylhexoside | Flavonoids |
| mr1093 | −0.38 | −0.14 | Tricin O‐malonylhexoside | Flavonoids |
| mr1204 | −0.35 | −0.16 | Chrysoeriol 7‐O‐hexoside | Flavonoids |
| mr1275 | −0.38 | −0.13 | Tricin 7‐O‐hexoside | Flavonoids |
| mr2090 | −0.28 | −0.18 | Phloretin | Flavonoids |
| S21_L‐0942 | −0.36 | −0.12 | Isorhamnetin 3‐galactoside | Flavonoids |
| n00282 | −0.38 | −0.13 | #N/A | Unknown |
| S21_L‐3824 | −0.20 | −0.12 | #N/A | Unknown |
| S21_L‐4082 | −0.21 | −0.22 | #N/A | Unknown |
| S21_L‐4831 | −0.23 | −0.10 | #N/A | Unknown |
| S21_S‐0678 | −0.22 | 0.24 | #N/A | Unknown |
| S21_S‐1014 | −0.24 | −0.11 | #N/A | Unknown |
| S21_S‐1266 | 0.22 | 0.19 | #N/A | Unknown |
| S21_S‐4234 | −0.21 | 0.17 | #N/A | Unknown |
| wm0031 | −0.23 | 0.11 | #N/A | Unknown |
| mr1170 | 0.27 | 0.11 | Caffeoyl shikimic acid | Phenylpropanoids |
- —National Key Research and Development Plan of China
- —National Science Fund for Distinguished Young Scholars10.13039/501100014219
- —Basic Research Project of Yazhouwan National Laboratory
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Taxonomy
TopicsLipid metabolism and biosynthesis · Genetic Mapping and Diversity in Plants and Animals · Nitrogen and Sulfur Effects on Brassica
Introduction
1
More than 200000 types of metabolites are estimated to exist in plants.^[^ 1, 2 ^]^ Typically, variations in plant growth or responses to stresses are accompanied by the decrease or increase in specific metabolites. These compounds serve as valuable markers or indicators for distinct biological processes or physical traits.^[^ 3, 4, 5, 6, 7 ^]^ Plant metabolites play a pivotal role in environmental defense and adaptation, essential for plant growth and development.^[^ 8, 9 ^]^ Numerous studies have linked certain metabolites strongly to physical traits,^[^ 6, 7 ^]^ such as their role in regulating maize lignin content to increase plant height and biomass,^[^ 10 ^]^ the discovery of amino acid metabolites as markers for improving soybean seed oil content (SOC) or seed protein content,^[^ 11 ^]^ and the identification of flavonoid metabolites as markers for SOC in Brassica napus (B. napus).^[^ 12 ^]^ Given that secondary metabolites are more diverse compared to primary metabolites,^[^ 13, 14 ^]^ their utilization, along with genome and transcriptome data, has led to an increasing number of studies showcasing how metabolomics advances the exploration of crop metabolic sciences.^[^ 15, 16, 17, 18 ^]^
QTL mapping is widely used in the agronomic loci and thousand‐seed weight (TSW) related gene discovery.^[^ 19, 20 ^]^ Such as Auxin‐response factor 18 (BnaA09.ARF18) is mapped at a major QTL on A09,^[^ 21 ^]^ P450 monooxygenase gene BnaA9.CYP78A9 also on A09 was mapped by a transposable element region controlling seed weight and silique length in B. napus.^[^ 22 ^]^ And shaggy‐like protein kinase 41 (SK41) is found at QTL qTGW3, controlling seed size and weight in rice.^[^ 23 ^]^ With the development of de novo sequencing, genome‐wide association studies (GWAS) have been identifying more novel loci.^[^ 24, 25 ^]^ For example, GmWRI14‐like genes were identified to control soybean seed weight,^[^ 26 ^]^ and flavin‐binding monooxygenase family proteins in peanut are associated with seed weight.^[^ 27 ^]^ Metabolome‐based genome‐wide association studies (mGWAS) have successfully unearthed numerous loci linked to pivotal agronomic traits in major crops, such as rice, peach, sesame, and lettuce.^[^ 16, 17, 28, 29, 30, 31, 32, 33, 34 ^]^ In tandem with these identified loci, metabolome‐based transcriptome‐wide association studies (mTWAS) aim to elucidate gene‐metabolite relationships through comprehensive association analysis.^[^ 12, 35, 36, 37, 38, 39 ^]^ The integration of transcriptome and genome data in expression genome‐wide association studies (eGWAS) is designed to pinpoint genomic loci influencing gene transcript levels.^[^ 12, 40 ^]^ Consequently, a multi‐omics analysis incorporating genome, transcriptome, and metabolome data holds significant promise for advancing the discovery and utilization of loci and genes associated with crop traits.
TSW is one of the three key yield determinants for B. napus,^[^ 41 ^]^ which is a global secondary oil crop providing ≈13% of edible oil (USDA ERS, 2022). Recent investigations have unveiled several pathways governing seed weight, including the involvement of the ubiquitin receptor DA2 in the ubiquitin‐proteasome pathway,^[^ 42 ^]^ Arabidopsis G protein gamma subunit 3 (AGG3) in G‐protein signaling,^[^ 43, 44 ^]^ and mitogen‐activated protein kinase phosphatase1 (MKP1), mitogen‐activated protein kinase 6 (MKP6), and mitogen‐activated protein kinase kinase 5 (MKK5) in MAPK signaling.^[^ 45, 46, 47, 48, 49 ^]^ DNA METHYLTRANSFERASE1 (MET1) and DECREASE IN DNA METHYLATION1 (DDM1) exhibit parent‐of‐origin effects.^[^ 50 ^]^ Notably, downregulating the ubiquitin receptor BnaDA1 has been demonstrated to increase B. napus seed size.^[^ 51 ^]^
Numerous transcription factors (TFs) influencing seed weight have been reported.^[^ 52, 53 ^]^ In Arabidopsis, the chromatin status of TRANSPARENT TESTA GLABRA2 (TTG2), regulated by TOPOISOMERASE Iα (TOP1‐α), is instrumental in determining seed size.^[^ 54, 55 ^]^ Over‐expression of Leafy Cotyledon 1 (LEC1) and WRINKLED 1 (WRI1) from Pistacia chinensis in Arabidopsis has been shown to enhance seed weight.^[^ 56 ^]^ HECT E3 ubiquitin ligase (BnaC03.UPL3) has been demonstrated to directly impact the protein level of Leafy Cotyledon 2 (LEC2), exerting a negative influence on both TSW and SOC.^[^ 57 ^]^ In soybeans, the down‐regulation of BIG SEEDS1 (BS1), a plant‐specific TF, significantly increases seed size.^[^ 58 ^]^ DNA binding with one finger (Dof) and GmWRI14‐like genes also play roles in influencing seed weight in soybean.^[^ 26, 59 ^]^ On the other hand, TGA6, a transcription factor possessing a basic region domain leucine zipper domain (b‐ZIP) and belonging to the second branch of the Arabidopsis TGA family,^[^ 60 ^]^ has been reported to directly bind to the promoter sequence of PATHOGENESIS‐RELATED1.^[^ 61 ^]^ It participates in activating the salicylic acid response gene alongside TGA2 and TGA5.^[^ 62 ^]^ TGA6 inhibits immune response through physical interaction in the nucleus.^[^ 63 ^]^ However, its impact on TSW has not been previously reported in plants.
A comprehensive analysis quantified 2172 metabolites in mature seeds of B. napus in our previous work.^[^ 12 ^]^ This is the first metabolome‐coupled multi‐omics cohort study that reveals 137 marker metabolites significantly associated with TSW by computing correlation coefficients between TSW and the metabolome. Utilizing multi‐omics analysis and weighted correlation network analysis (WGCNA), we construct a ternary relationship network linking metabolites, genes, and QTLs to TSW. Our findings identify 30 seed size‐related candidate genes, including Auxin response factor 2 and AP2‐like ethylene‐responsive transcription factor ANT. Experimental validation confirms the feasibility and effectiveness of our approach. Notably, BnaA03.TGA6 emerges as a novel transcription factor influencing TSW in B. napus.
Results
2
Metabolic Profiling Identifies 137 Marker Metabolites of TSW in B. Napus
2.1
In our prior investigation, a comprehensive metabolomics approach quantified a total of 2172 metabolites in mature seeds from 388 B. napus inbred lines.^[^ 12 ^]^ To explore marker metabolites associated with TSW in B. napus, we have conducted a corresponding calculation of the correlation between TSW and metabolites. The correlation indices are computed separately for two consecutive years, revealing 46 consistently negatively correlated TSW marker metabolites (Figure 1a) and 91 consistently positively correlated ones across two years (Figure 1b). Consequently, 137 metabolites across two years (p < 0.05) are identified as TSW marker metabolites (Table S1, Supporting Information). These marker metabolites are classified based on NPClassifier,^[^ 64 ^]^ including 14 flavonoid metabolites, 10 amino acid derivatives, 4 phenylpropanoid metabolites, and 109 unknown metabolites. The correlation indices among TSW marker metabolites over the two years demonstrate good repeatability (Figure 1c). Furthermore, correlation analysis for the 137 TSW marker metabolites reveals high inter‐class correlation among metabolites within the same class (Figure 1d). The average heritability value of TSW marker metabolites is 0.67, surpassing the average value of all detected metabolites (Figure S1, Supporting Information).
B. napus TSW and metabolome data in two years. a) The Venn diagram of metabolites negatively correlated with TSW in two years. b) The Venn diagram of metabolites positively correlated with TSW in two years. c) The correlation index of 137 TSW‐correlated metabolites in two years. d) The heatmap of different classes of 137 TSW‐correlated metabolites in the natural population. The color of each cell represents the correlation levels with each TSW‐correlated metabolite.
Multi‐Omics Analysis Reveals QTLs Potentially Regulating TSW
2.2
mGWAS analysis has been conducted on the 137 TSW marker metabolites,^[^ 12 ^]^ resulting in the detection of 1154 significantly associated variations for 133 TSW marker metabolites (Figure 2a; Table S2, Supporting Information). Notably, 111 significant signals are identified for 13 flavonoids, 27 for 4 phenylpropanoid metabolites, 183 for 9 amino acid derivatives, and 1046 for 107 unknown metabolites (Figure 2a; Table S2, Supporting Information). Utilizing a linkage disequilibrium (LD) threshold of 100 kb, these significantly associated loci are categorized into 734 metabolite‐QTLs (mQTLs). These mQTLs colocalization with more than 30 signals are distributed on chromosomes A03, C03, C05, C06, C07, and C09 (Figure 2b; Table S2, Supporting Information). These results indicate that there are obvious mQTL hotspots of TSW marker metabolites.
mGWAS and module‐trait association studies of TSW and TSW‐correlated metabolites. a) Chromosomal distribution of all lead variations associated with TSW‐correlated metabolites (the results of 2017 and 2018 are combined). b) hromosomal distribution of all mQTLs associated with TSW‐correlated metabolites (the results of 2017 and 2018 are combined). The interval size is 100 kb. c) Heatmap of correlation among TSW‐correlated metabolites, TSW, and eigenvalues of modules obtained by WGCNA. Each eigenvalue represents the profile of the genes in the correlated module. The vertical axis shows the modules. d) The genome‐wide distribution of eQTLs. The X‐axis represents the physical position of the eQTLs, and the Y‐axis represents the number of genes associated with the eQTL. The orange point indicates that the eQTL colocalizes with the mQTL detected by mGWAS of TSW‐correlated metabolites.
In our prior research, transcriptome data from seeds developed 40 days after flowering (DAF) were obtained.^[^ 37 ^]^ To enhance our understanding of the relationship between transcripts and TSW marker metabolites, mTWAS has been performed on the 40 DAF transcriptome datasets, revealing 9077 genes significantly associated with TSW marker metabolites (p < 1.41 × 10^−05^; Table S3, Supporting Information). WGCNA has been executed using the transcripts from the 40 DAF transcriptome data, clustering them into 139 modules. Correlation coefficients of these modules with both TSW and TSW marker metabolites have been calculated (Figure 2c; Table S4, Supporting Information). Among these modules, 19 exhibit significant correlations with TSW (p < 0.01), with module sienna3 displaying the strongest correlation (Table S5, Supporting Information). GO enrichment analysis reveals that module sienna3 is enriched in cell fate determination and seed coat development (Figure S2a, Supporting Information), while other modules, such as blue, cyan, and turquoise are enriched in pathways such as lipid metabolism and seed development (Figure S2b–d, Supporting Information). These findings underscore the utility of gene co‐expression networks in unveiling candidate gene sets correlated to TSW.
To elucidate the connection between genomic variation and gene expression, we have conducted eGWAS on the 9077 genes significantly correlated with TSW marker metabolites. The eGWAS reveals 59121 variations significantly associated with these genes. Employing a linkage disequilibrium (LD) threshold of 100 kb, we have categorized these variations into 5225 expression‐QTLs (eQTLs) (Table S6, Supporting Information). Notably, eQTL hotspot analysis highlights specific genomic regions, particularly on chromosomes A02, A03, A05, A09, C01, C02, C03, C06, C07, and C09, where eQTLs are associated with over 100 genes (Figure 2d). This observation implies the presence of transcription regulatory factors in these eQTL hotspots, potentially manipulating a cluster of genes that directly or indirectly influence TSW marker metabolite content.
A Seed Development‐Related Gene BnaC09.DDM1 is Significantly Associated with TSW
2.3
DECREASED DNA METHYLATION 1 (DDM1) plays a pivotal role as an ATP‐dependent chromatin remodeler crucial for maintaining DNA methylation in plants, a process essential for seed development.^[^ 50 ^]^ In our study, BnaC09.DDM1 (BnaC09g07810D) is predicted to be significantly associated with TSW‐correlated metabolites. Six mQTLs are identified to be associated with S21–2486 (Figure 3a), and haplotype analysis within mQTL4855 (Table s7, Supporting Information) reveals that, at the population level, inbred lines with haplotype C exhibit higher S21–2486 content in seeds (p = 8.1 × 10^−27^; Figure 3b). Since it's a TSW correlated metabolite QTL hotspots, all variations, position and genomic region information has been summarized in Table S7 (Supporting Information). mTWAS analysis demonstrates a significant association between BnaC09.DDM1 and S21–2486 content (Figure 3c). Furthermore, the expression levels of BnaC09.DDM1 show a significant correlation with S21–2486 (r = 0.22, p = 8.2 × 10^−04^; Figure 3d). eGWAS results indicate that the expression level of BnaC09.DDM1 is significantly associated with 1 cis‐eQTL and 7 trans‐eQTLs (Figure 3e). Haploid analysis reveals that haplotype A represents a high expression level of BnaC09.DDM1 (p = 8.9 × 10^−04^, Figure 3f; Table S6, Supporting Information). Notably, BnaC09.DDM1 is associated with 14 metabolites whose 18 mQTLs are co‐located on the C09 chromosome, suggesting the identification of a mQTL hotspot with potential influence on TSW (Table S2, Supporting Information). Additionally, mTWAS results highlight the significant association of BnaC09.DDM1 with three metabolites: wm0034 (p = 7.23 × 10^−06^), S21–2486 (p = 1.13 × 10^−06^), and S21_L‐3616 (p = 3.08 × 10^−06^) (Table S3, Supporting Information). We have annotated wm0034 as 4‐indolecarbaldehyde and provided high‐resolution spectra for the other two unknown metabolites (Figure S3, Supporting Information). These findings collectively position BnaC09.DDM1 as a promising candidate gene influencing TSW.
Multi‐omics identifies BnaC09.DDM1 as TSW candidate gene. a) Manhattan plot of mGWAS for S21‐2486. b) Haplotype analysis of the lead variation in mQTL4855 for S21–2486 content. c) Manhattan plot of mTWAS for S21–2486. d) Correlation analysis between S21 and 2486 content and BnaC09.DDM1 expression level. e) Manhattan plot of eGWAS for BnaC09.DDM1. f) Haplotype analysis of the lead variation in eQTL4855 for BnaC09.DDM1 expression level.
Ternary Relationship Networks Boost the Understanding of TSW and SOC
2.4
To elucidate our network‐building approach, we have outlined a comprehensive flowchart (Figure 4a). The outcomes of mGWAS and transcriptome analysis have been amalgamated into a ternary relationship network (gene‐QTL‐metabolite) (Figure 4b; Table S7, Supporting Information), encompassing 731 QTLs, representing both mQTLs and eQTLs, 3321 genes, and 133 TSW marker metabolites. Based on the amalgamated results of TSW, marker metabolites, mQTLs, eQTLs, and mTWAS, we have compiled a list of 30 seed size‐related candidate genes influencing B. napus TSW (Table S8, Supporting Information).
TSW‐correlated metabolites, genes, and QTLs network. a) Flow chart of network building. b) Network built based on the correlation among TSW‐correlated metabolites, genes, and QTLs. Candidate genes, BnaA10.ARF4, BnaA04.MET2, BnaC01.CLM, BnaA01.RPT2‐α, BnaC09.BG1, BnaA06.SRE1, and BnaC09.ARF2 are shown as circles with a distinct color per co‐expression module. All modules significantly correlated with TSW are shown in this network. The QTLs are detected simultaneously by the significant associated genes (eQTL) and metabolites (mQTL) in this triple relationship network.
Additionally, through a comparative analysis of marker metabolites for SOC and TSW, we have identified 21 marker metabolites shared between SOC and TSW. These include 11 flavonoid metabolites, 1 phenylpropanoid metabolite, and 9 unknown metabolites (Table 1; Figure S4, Supporting Information). Among the flavonoid metabolites, 11 exhibiting negative correlations with both SOC and TSW, are associated with 31 mQTLs (Table S2, Supporting Information). Three unknown metabolites S21_S‐0678, S21_S‐4234, and wm0031, displaying negative correlations with SOC and positive correlations with TSW, are associated with 34 mQTLs (Table S2, Supporting Information). The mQTLs mapped based on these metabolites may unveil candidate genes potentially influencing both TSW and SOC. Integrating SOC and TSW results, we observe that a flavonoid metabolite, tricin 5‐O‐hexoside (mr1083, Table 1), exhibits a correlation of r = −0.40 with SOC and r = −0.14 with TSW, while a phenylpropanoid metabolite, caffeoyl shikimic acid (mr1170, Table 1), shows a correlation of r = 0.27 with SOC and r = 0.11 with TSW.
Referring to our ternary relationship network for SOC and TSW (Figure S4, Supporting Information), it's noteworthy that 26 mQTLs corresponding to 43 metabolites, including mr1083 and mr1170, are co‐localized on the C09 chromosome (Table S8, Supporting Information). Notably, eGWAS results for Auxin response factor 2 (BnaA09.ARF2, BnaA09g05840D) coincide with mQTL1635 for mr1083, and BnaC09.ARF2 aligns with mQTL4857 for mr1170 (Table S8, Supporting Information). The mTWAS results reveal significant associations of BnaA09.ARF2 with mr1421 (p = 6.22 × 10^−06^) and S21_S‐0305 (p = 5.97 × 10^−06^); BnaC09.ARF2 with wm0034 (p = 3.73 × 10^−06^) and S21–2486 (p = 7.68 × 10^−07^) (Table S3, Supporting Information). Furthermore, eGWAS of BnaC09.ARF2 indicates several cis‐eQTLs on the C09 chromosome and trans‐eQTLs on A02, A03, A07, A09, A10, C02, C07, and C09. Similar observations are made for BnaA09.ARF2 (Table S6, Supporting Information). ARF2, known for specifically binding to auxin response promoter elements 5′‐TGTCTC‐3′, plays a role in auxin signaling, cell division, and seed size regulation.^[^ 65, 66, 67 ^]^ BnaARF2s may impact both SOC and TSW by influencing the content of phenylpropanoid metabolites (mr1170 and mr1083), making them potential candidates for future improvements in seed traits of B. napus. In summary, we have identified 67 candidate genes associated with 21 marker metabolites for both SOC and TSW, presenting promising avenues for improving both B. napus TSW and SOC breeding (Table S9, Supporting Information).
BnaA03. TGA6 is a Key Transcription Factor Gene Negatively Regulating TSW in B. Napus
2.5
A notable mQTL hotspot, designated as mQTL555 on A03, is associated with a total of 11 metabolites. All variations, position and genomic region information about the TSW correlated mQTL555 hotspots have been summarized in Table S7 (Supporting Information). Among these, 10 metabolites S21–2077, S21–1527, S21_L‐2853, mr1169 (roseoside), S21_L‐4193, S21_S‐5484, S21–4031, S21–4360, S21–0694, and S21_S‐0841, are significantly positively correlated with TSW, while S21_S‐1322 is significantly negatively correlated with TSW (Figure 5a–c; Figures S5 and S6, Supporting Information). High‐resolution spectra for these 10 unknown metabolites and mr1169 are provided (Figure S7, Supporting Information). Within these 11 unknown metabolites, three intriguing marker metabolites S21–0694, S21–1527, and S21–4360 positively correlate with TSW (Figure 5d–f). Haplotype analysis reveals that haplotype T of lead variations for S21–0694 represents higher content (Figure 5g, p = 6.6 × 10^−22^), haplotype TAGA represents higher content of S21–1527 (Figure 5h, p = 6.3 × 10^−08^), and haplotype A represents higher content of S21–4360 (Figure 5i, p = 3.2 × 10^−18^). In the mTWAS results, BnaA03.TGA6 is associated with three metabolites (Figure 5j; S21–4360, p = 2.84 × 10^−08^; S21–4031, p = 4.93 × 10^−08^; S21–0694, p = 2.93 × 10^−11^; Table S3, Supporting Information). eGWAS analysis of BnaA03.TGA6 reveals its association with one cis‐eQTL and two trans‐eQTLs (Figure 5k). A local Manhattan map and haplotype analysis indicate significant variation in this gene, with a p‐value of 7.8 × 10^−08^ (Figure 5l–n). Additionally, WGCNA analysis assigns BnaA03.TGA6 to the turquoise module (Table S4, Supporting Information). Based on multi‐omics analysis, we posit that BnaA03.TGA6 is a novel transcription factor affecting TSW in B. napus.
Multi‐omics analysis predicts BnaA03.TGA6 as a regulator of TSW. a) Manhattan plot of mGWAS results (2017) of S21–0694, b) S21–1527, c) S21–4360. d) Correlation analysis of TSW and S21–0694, e) S21–1527, f) S21–4360. g) Haplotype analysis of the lead variation in mQTL555 for S21–0694, h) S21–1527, i) S21–4360. j) Manhattan plot of mTWAS results of S21–4360. k) Manhattan plot of eGWAS results of BnaA03.TGA6. l,m) Local Manhattan plot of eGWAS results of eQTL555. n) Haplotype analysis of the lead variation in eQTL555 for BnaA03.TGA6.
Examining the expression levels of BnaTGA6 genes across various developmental stages of ZS11 seeds (https://yanglab.hzau.edu.cn/BnIR) reveals that BnaTGA6 encompasses 6 homologous genes in ZS11. Notably, the copies on A01, A03, A05, C01, C03, and C05 exhibit similar expression patterns (Figure S8, Supporting Information). Delving into the function of BnaA03.TGA6, we have generated single‐gene mutants, which display no significant alterations in TSW (Figure S9, Supporting Information). Subsequently, we have created hexamutants by crossing different mutant lines of the 6 BnaTGA6 genes in B. napus Westar (Table S10, Supporting Information). The TSW of two independent BnaTGA6 hexamutant lines, CR1 and CR2, is 27.6% and 29.9% higher than that of the wild type (WT), respectively (Figure 6a,b). It indicates that there is functional redundancy between BnaA03.TGA6 and other TGA6 genes in the genome of B. napus. Metabolomic analysis has unveiled significantly elevated levels of metabolites, including S21–0694 (Figure 6c), S21–1527 (Figure 6d), and S21–4360 (Figure 6e) in the mutants compared to WT. Metabolome analysis of TGA6 mutants and WT indicates that 108 metabolites are down‐regulated, while 58 metabolites are up‐regulated (Figure 6f). Specifically, the most down‐regulated metabolites are flavonoids and lipids, while the most up‐regulated metabolites are amino acid derivatives (Figure 6g,h).
Functional identification of BnA03.TGA6 as a negative regulator of TSW in B. napus. a) Seeds of CRISPR/Cas9 TGA6 mutants (CR1 and CR2). Bar = 1 cm. b) TSW in TGA6 mutants seeds and WT seeds. c) Peak area of S21–0694 d) S21–1527, e) S21–4360 in WT and TGA6 mutants seeds. f) The volcano plot for differentially accumulated metabolites in TGA6 mutants vs WT (FC > 1 and adjusted p < 0.05). The red and blue points indicate up‐ and down‐regulated metabolites, respectively. g,h) Number of up‐ and downregulated metabolites and their classification. i) The volcano plot for differentially expressed genes in TGA6 mutants vs WT of developing seed at 40 DAF (FC > 2 and adjusted p < 0.05). j,k) Reported TSW‐related gene expression correlation analysis with BnaA03.TGA6 and S21–0694. The horizontal coordinates indicate the correlation between gene expression and S21–0694 at 20 and 40 DAF, and the vertical coordinates indicate the correlation between gene expression and BnaA03.TGA6. Blue dots indicate significantly correlated genes, and orange dots indicate reported TSW‐related genes. Values are means ± s.e.m., n = 3. Statistical analysis is using Student's t‐test (, p < 0.05; *, p < 0.01).
Furthermore, we have conducted RNA‐seq on developing seeds (40 DAF) of both WT and CR1. The transcriptome analysis has revealed 3715 differentially expressed genes (DEGs, log_2_FC > 2, P‐adjust<0.05) in CR1 compared to WT, with 1595 genes up‐regulated and 2120 down‐regulated and 6 genes related to seed size regulation exhibit significantly altered expression levels and 6 TSW related DEGs (Figure 6i). Interestingly there are two Polycomb group protein FERTILIZATION‐INDEPENDENT ENDOSPERM (BnaA05.FIE and BnaAnn.FIE) involved in endosperm development, are significantly decreased (Figure 6i).^[^ 68 ^]^ And two gene cyclin‐D7‐1 (BnaC03.CYCD7) and cyclin‐D3‐1 (BnaA01.CYCD3) involved in early endosperm and embryo development (Figure 6i).^[^ 69, 70 ^]^ We also found another DDM1 (BnaC07.DDM1) is significantly increased in the tga6 mutant. Interestingly, DA2 ^[^ 42, 71 ^]^ participate in the ubiquitin pathway among the 6 TSW related genes (Figure 6i). Go enrichment of all 3715 DEGs we found that they are enriched in processes related to photosynthesis, photosystem I, and photosystem II (Figure S10, Supporting Information), indicating that the tga6 mutant influences photosynthesis in developing seeds.
We have proved that S21–0694, S21–1527, and S21–4360 significantly increased in the tga6 mutants (Figure 6c–e). To explore the BnaTAG6 metabolome level mechanism, we have also provided predicted structures for these three TSW marker metabolites (Figure S11, Supporting Information). Importantly, S21_0694, the highest correlated TSW marker metabolite, has been further validated its representativeness for TSW through correlation analysis with both the population expression of BnaA03.TGA6 and the population metabolite level of S21_0694 using transcriptome data at 20 DAF and 40 DAF. Notably, 10 reported TSW genes stand out in the 20 DAF transcriptome (Figure 6j), and 12 genes in the 40 DAF (Figure 6k). This suggests that S21_0694 is validated as a major marker for TSW, as its concentration increased in parallel with the TSW gain observed in the tga6 mutants and associated with many TSW‐correlated genes at a population level.
BnaA03. TGA6 Regulates TSW by Inhibiting BnaA07.DA2
2.6
To further explore BnaTGA6 mediated regulation, we have performed a genome‐wide prediction of BnaA03.TGA6 binding motifs among the significantly altered genes and identify 1313 potential BnaA03.TGA6 interacting targets (Table S11, Supporting Information). We have found that among 6 TSW related DEGs (Figure 6i), only BnaA07.DA2 harbored predicted BnaA03.TGA6 binding motif in its promoter region (Table S11, Supporting Information), and it also shows the significant expression fold change among the 6 DEGs (Figure 7a). To detect BnaA03.TGA6 enrichment at the endogenous BnaA07.DA2 promoter under near‐physiological chromatin conditions in vivo, we have performed Cleavage Under Targets and Tagmentation (CUT&Tag) using protoplasts isolated from rapeseed leaves. We have provided annotated sequence data in the Table S12 (Supporting Information), and there are 16390 hits in total and 8720 hits are sequenced at genes’ promoter region within 2kb. There are 437 genes are overlapped between 3715 tga6 mutant transcriptome and 8720 BnaA03.TGA6 CUT&Tag hit genes (Figure 7b). Go enrichment analysis indicates that BnaA03.TGA6 significantly effect photosystem (Figure S12, Supporting Information). Sequencing binding profile indicates significant peaks in the BnaA07.DA2 promoter region P1 (Figure 7c). BnaA03.TGA6 is predicted to have two types of putative binding motifs.^[^ 72 ^]^ Among these, the TGACGT motif, P1, is found within the BnaA07.DA2 promoter (Figure 7d). In addition, we have adopted a targeted and sensitive CUT&Tag‐qPCR approach to quantify BnaA03.TGA6 binding at the P1 region, suggesting that BnaA03.TGA6 is significantly enriched at Bna07.DA2 P1 promoter region (Figure 7e). To validate whether BnaA03.TGA6 regulates BnaA07.DA2 expression, we conduct a transient dual‐luciferase assay in Arabidopsis protoplasts. The results demonstrate that Luciferase (LUC) expression is driven by the BnaA07.DA2 promoter is significantly upregulated when co‐transformed with the BnaA03.TGA6 effector (Figure 7f,g), suggesting that BnaA03.TGA6 may activate BnaA07.DA2 transcription. Complementarily, electrophoretic mobility shift assays (EMSAs) using purified His‐BnaA03.TGA6 protein and synthetic promoter fragments confirmed direct binding to the P1 regions of the BnaA07.DA2 promoter in vitro (Figure 7h). Collectively, the in vivo and in vitro results demonstrate that BnaA03.TGA6 directly binds to the promoter region of BnaA07.DA2 in the ubiquitin pathway and activates its expression, thereby negatively regulating TSW in rapeseed (Figure 7i).
Regulatory mechanism of BnaA03.TGA6. a) Expression level of BnaA07.DA2 in WT and tga6 mutant developing seed (40DAF). b) Venn gram of tga6 mutant transcriptome (log2FC > 2, P.adj < 0.05) and BnaA03.TGA6 CUT&Tag hits (promoter region within 2kb) c) BnaA03.TGA6 CUT&Tag binding profile and corresponding expression levels near the BnaA07.DA2 region. Peaks are normalized to counts per million (CPM). Blue peaks represent BnaA03.TGA6 binding sites on the BnaA07.DA2 promoter. Grey and red peaks represent expression data for BnaA07.DA2; red is not visible due to low expression. Orange rectangles indicate the promoter region targeted in CUT&Tag‐qPCR, shown in f. d) Enrichment of the “TGACGT” motif within the BnaA07.DA2 promoter region targeted by BnaA03.TGA6. e) CUT&Tag‐qPCR results showing in vivo binding of BnaA03.TGA6 to the BnaA07.DA2 promoter region. f) Schematic representation of the constructs used in the transient dual‐luciferase assay. The effector plasmid contains the BnaA03.TGA6 coding sequence driven by the CaMV 35S promoter. The reporter plasmid includes the REN gene under the control of the 35S promoter and the LUC gene driven by the BnaA07.DA2 promoter (ProBnaA07.DA2). g) Bar graph showing the relative LUC/REN ratio from the dual‐luciferase assay, indicating the transcriptional activation of BnaA07.DA2 by BnaA03.TGA6. h) EMSA results demonstrating that the His‐BnaA03.TGA6 fusion protein directly binds to specific motif within the BnaA07.DA2 promoter. i) Proposed model of transcriptional regulation: BnaTGA6 activates BnaA07.DA2, contributing to TSW regulation via the ubiquitin pathway. Values are means ± s.e.m. (n = 3). Statistical significance was determined by Student's t‐test (, p < 0.05; *, p < 0.01).
Discussion
3
In recent years, as omics data become increasingly available, association analysis, including GWAS and TWAS, has emerged as a crucial tool for exploring genetic variations in crops.^[^ 37, 38, 39, 73, 74, 75 ^]^ Linking marker metabolites to multiple agronomic traits has identified potential breeding loci for various crops.^[^ 6, 7, 76 ^]^ In this study, we have identified 137 marker metabolites correlated with TSW using two years of B. napus mature seed metabolome and TSW data. Through association analysis (Figure 2a,b,d), we have uncovered numerous potential loci influencing B. napus TSW. In contrast to previous reports on genes and loci affecting B. napus TSW,^[^ 51, 77, 78, 79, 80 ^]^ our findings reveal positive correlations between amino acid derivatives and phenylpropanoid metabolites with TSW, while flavonoid metabolites exhibit a negative correlation. These results suggest the potential for improving B. napus TSW through metabolic pathway engineering.
Compared to limited loci from single agronomic trait linkage or association analysis, high‐throughput metabolomic data provide a rich resource for potential breeding loci.^[^ 78, 81, 82, 83 ^]^ Previous studies have identified B. napus other TSW loci on the A03 chromosome.^[^ 78, 84 ^]^ Such as fine mapping of cqSW.A03‐2 narrowed the locus to a 61.6‐kb region, in which the histidine kinase gene BnaA03G37960D was identified as the candidate gene,^[^ 78 ^]^ and GWAS discovered a SNP, rs6515 on A03 (22079–22,2 kb) is associated with B. napus TSW.^[^ 84 ^]^ The discovery of mQTL555 in this study underscores the power of combining metabolomic analysis to unveil new B. napus TSW loci (Figure 2a,b). But we also noticed that metabolites S21–0709 (putative flavonoid), S21–4140 (putative amino acid derivative), and S21–1134 (other) are co‐localized at mQTL1534 (Position: A08 12758480), which is co‐localized with a previous SNP (rs17921, A08, 12464–12777 kb).^[^ 84 ^]^ mQTL3177 (S21_L‐4233, S21‐4731, S21_L‐1836, S21‐1134; position C03 52264 kb) is co‐localized with rs17613 (Position C03 51969–52357 kb).^[^ 84 ^]^ These co‐localized QTLs suggest potential metabolite associations for those loci. We provided a total of 734 mQTLs in this study, most of which are not reported by earlier TSW GWAS, which provides guide information for the B. napus TSW breeding and needs further validation. Furthermore, TWAS offers unique advantages over GWAS by associating single genes with traits rather than linkage intervals.^[^ 35, 36, 37, 38, 39 ^]^ Our comprehensive mTWAS of B. napus TSW correlated metabolites furnishes extensive genetic evidence and candidate genes potentially influencing B. napus TSW (Table S3, Supporting Information). Notably, based on multi‐omics data, we have constructed the first TSW triple relationship network (metabolite‐gene‐QTL) for B. napus (Figure 4), laying a metabonomic foundation for future B. napus TSW breeding.
The identification of BnaC09.DDM1, an ATP‐dependent chromatin remodeler essential for maintaining DNA methylation, suggests a potential mechanistic link between epigenetic regulation and metabolic control during seed development. Our findings imply that it may also influence seed metabolite composition, thereby indirectly contributing to agronomic traits such as TSW. On the other hand, through multi‐omics analysis, we have identified a novel gene, BnaA03.TGA6. Considering gene function redundancy, we have generated hexamutants of BnaTGA6 (Figure S9, Supporting Information). The seeds of the TGA6 mutant materials are significantly larger than the WT (Figure 5a,b). Interestingly, three out of eleven metabolites positively correlated with TSW are associated with the BnaA03.TGA6 locus, and their metabolite levels changed accordingly (Figure 6c–e). This suggests that alterations in B. napus TSW are accompanied by changes in metabolite content. Whole‐genome scanning reveals more metabolites associated with mQTL555, with the majority being annotated as flavonoids and phenylpropanoid metabolites (Table S12, Supporting Information). Given the tight association between mQTL and metabolites, those metabolites having the same mQTL may share structural similarities. It's possible that BnaA03.TGA6 may also participate in flavonoid and phenylpropanoid biosynthesis regulation, or there are more genes clustered at mQTL555 involved in the biosynthesis pathways. We have also analyzed the correlation between 20 DAF and 40 DAF with the expression profile of BnaA03.TGA6 (Figure 6j,k) and emphasize the potential of S21_0694 been a valuable TSW marker metabolite. For BnaA03.TGA6 regulation validation, we have performed dual‐luciferase assay, CUT&Tag, and EMSA with BnaA07.DA2. Combining TGA6 mutant transcriptional data, we find that BnaA07.DA2 is the downstream target of BnaA03.TGA6, activated by it (Figure 7). In addition, we have fully investigated the tga6 40DAF transcriptome data, find that BnaA09.MET1 ^[^ 50 ^]^ is significantly increased in the tga6 mutant (Figure S13a, Supporting Information) and predicted with 1 binding motif of BnaA03.TGA6 (Figure S13b, Supporting Information). However, dual‐luciferase assay shows that BnaA09.MET1 is activated by BnaA03.TGA6 (Figure S13c, Supporting Information) and not identified by CUT&Tag. It suggests that BnaA09.MET1 is not the downstream target of BnaA03.TGA6. We also thoroughly investigate the BnaA03.TGA6 CUT&Tag sequencing data, find that BnaA03.TGA6 is enriched at the promoters’ region of BnaA01.TTG2 ^[^ 54 ^]^ and BnaC01.TOP1‐α ^[^ 54 ^]^ within 1kb (Table S12, Supporting Information). And BnaA01.TTG2 is predicted with 1 binding motif of BnaA03.TGA6 and BnaC01.TOP1‐α with 4 (Figure S13b, Supporting Information). Dual‐luciferase assay shows that BnaC01.TOP1‐α is repressed by BnaA03.TGA6 and BnaA01.TTG2 is activated by BnaA03.TGA6 (Figure S13d,e, Supporting Information). Interestingly, BnaA03.TGA6 is enriched at the promoter region of 4 BnaDA1 ^[^ 42 ^]^ homologue genes (Table S12, Supporting Information) but not significantly changed at tga6 40DAF transcriptome. BnaA08.DA1 and BnaA08.DA1 are both predicted with 5 binding motifs of BnaA03.TGA6 and BnaA06.DA1 and BnaC05.DA1 don't have BnaA03.TGA6 binding motifs (Figure S13b, Supporting Information). Then, we have chosen BnaA08.DA1 to perform validation experiment by CUT&Tag‐qPCR and EMSA. Only two BnaA03.TGA6 binding motifs at the promoter region of BnaA08.DA1 are showing interaction (Figure S13f,g, Supporting Information). We believe that BnaA03.TGA6 is interacting with multiple downstream targets within different pathways, especially BnaA08.DA1 and BnaA07.DA2 of ubiquitin pathway at different seed developing stages. And BnaA03.TGA6 is binding BnaA07.DA2 promoter region P1 to control TSW at 40DAF.
We consider this an exemplary instance of discovering novel genes for crop agronomic traits, building on our earlier work on SOC.^[^ 12 ^]^ Through multi‐omics analysis, we have identified and confirmed the impact of BnaA03.TGA6 on TSW, demonstrating its ability to activate BnaA07.DA2 expression. Our findings contribute to the understanding of the regulatory mechanisms governing seed size in plants. BnaA03.TGA6, a transcription factor of the b‐ZIP superfamily and part of the second branch of TGA in Arabidopsis, is known for its involvement in salicylic acid reactions and plant immunity.^[^ 60, 61, 62, 63 ^]^ Prior research on TGA6 has primarily focused on plant immune responses; this study discovered for the first time its new function in rapeseed seed development, suggesting the differentiation of gene function in evolution.
Due to the extraction method, our study has concentrated on water‐soluble metabolites, known for their high confidence.^[^ 16, 17, 18 ^]^ While our analysis spanned two years, the significant variance in metabolites correlated to TSW between 2017 and 2018 led us to intersect the data for subsequent analysis (Figure 1a,b). The 137 TSW‐correlated metabolites demonstrated good repeatability (Figure 1c), mitigating systematic errors to a certain extent. Most of these markers are unknown metabolites, necessitating further investigation. Techniques involving genetic analysis and mass spectrum ion information can aid in inferring the structure of unknown metabolites.^[^ 16 ^]^ Additionally, algorithm analysis, including machine learning and artificial intelligence, has been proposed to assist in metabolite annotation and structure prediction, potentially expediting metabolomics research.^[^ 85, 86, 87 ^]^ Although the functional annotation of unknown metabolites still needs to be further explored, the metabolite‐QTL network provided by this study lays a foundation for the development of molecular markers.
Experimental Section
4
Plant Materials, Collection of Phenotype Data
Based on the previous study, 388 accessions were selected from a population containing 505 accessions that were grown in Wuhan (2017, 2018).^[^ 37 ^]^ The WSeen SC‐G TSW analysis system (Hangzhou WSeen Detection Technology Co., Ltd., China) was used to acquire TSW phenotype data, 6 replicates for each accession. Transgenic B. napus materials were arranged in the transgenic crop field at Huazhong Agricultural University. Randomized complete block designs were applied for the field trial with three replications. Seed metabolome data used in this study were obtained from Li et al. (2023),^[^ 12 ^]^ in which 2172 metabolites were obtained. All metabolite extraction and data analysis procedures were performed as described in the original publication.^[^ 12 ^]^
mGWAS, mTWAS, eGWAS, and Co‐Expression Analysis
The variants used for mGWAS and eGWAS were called by a meta‐analysis of reads.^[^ 37, 88 ^]^ B. napus genome was used as the reference genome (B. napus v4.1, http://www.genoscope.cns.fr/brassicanapus/). Variations filtered with minor allele frequency > 5% were collected with PLINK.^[^ 89 ^]^ TSW‐correlated metabolites content in 2017 and 2018 (8274830), and gene expressions (8258336) were performed using mGWAS and eGWAS by FaST‐LMM software.^[^ 90 ^]^ The threshold was set at Bonferroni correction of p = 1.2 × 10^−07^, and Manhattan plots were generated using the R package. Similarly, 274 accessions transcriptome data at 40 DAF (2017, Wuhan) were used as described before.^[^ 37 ^]^ Using linear regression, mTWAS was performed within 2017 TSW‐correlated metabolites and gene expression, and the threshold was set at 1/70781.^[^ 91 ^]^ Then, these genes were then used to conduct co‐expression analysis with an R package (WGCNA).^[^ 92 ^]^ Finally, 139 modules were clustered, and the scale‐free topology of the determined networks was approximated based on a soft‐thresholding power (β = 6).
Network Building
For each QTL, the edges between TSW‐correlated metabolites were generated based on the significant signals of mQTLs; likewise, the edges between genes were based on the significant signals of eQTLs. In detail, all the significant signals over the threshold for 137 TSW‐correlated metabolites were clustered into mQTL LD blocks (100 kb). On the other hand, 9077 genes related to the 137 TSW‐correlated metabolites (p < 1.41 × 10^−5^) were selected for conducting eGWAS.^[^ 91 ^]^ All the significant signals over the threshold were clustered into eQTL LD blocks (100 kb). Finally, both mQTLs and eQTLs distance was calculated; if the distance was within 100 kb, it was considerd that there was a connection among TSW‐correlated metabolite, gene, and QTL (representing both mQTL and eQTL).
Vector Construction and Plant Transformation
The CRISPR‐Cas9 genome editing system^[^ 93 ^]^ was used to create TGA6 mutant materials. Golden Gate Assembly was used to conduct the transformation pKSE401 vector with the PCR fragment. Then, the recombinant vector was used to create TGA6 mutant materials. For detailed methods, please refer to Dai et al.^[^ 94 ^]^ BnaTGA6 hexamutants were created by crossing different mutant lines, and sequencing primers for TGA6 mutant materials were listed, which were used to confirm the mutant types (Table S14, Supporting Information).
Metabolite Extraction
A total of 0.1 g of mature B. napus seeds for WT and tga6 mutant was collected for each sample and homogenized in 1 mL of 70% methanol containing 0.1 mg L^−1^ acyclovir (internal standard) using a TissueLyser II (Qiagen, Germany) at 29 Hz for 60 seconds. The homogenate was vortexed and incubated at 4 °C for 10 h to extract metabolites. Following centrifugation at 9000 × g for 10 min, the supernatants were collected, pooled, and filtered through a 0.22 µm organic membrane filter (SCAA‐104; ANPEL, Shanghai, China) prior to LC‐MS analysis.
Dual‐Luciferase Assay
The Dual‐Luciferase Reporter Assay System (Promega, Madison, WI, USA) was utilized to perform dual‐luciferase assays. To construct the reporter plasmid, the promoters for BnaA09.MET1, BnaC01.TOP1‐α, BnaA01.TTG2 and BnaA07.DA2 was amplified and inserted into the pGreenII0800‐LUC vector. Similarly, the effecter plasmids were created by amplifying and inserting the open reading frame (ORF) of BnaA03.TGA6 into the pM999‐YFP vector. The experimental procedure was conducted in Arabidopsis protoplasts as described by Yoo et al.^[^ 95 ^]^ After 16 h, the transformed Arabidopsis protoplasts were lysed in 100 µL of passive lysate. Following a 30 s incubation, luminescence measurements were carried out using the SPARKRMULTIMODE MICROPLATE (TECAN, Swiss) with the subsequent steps: 50 µL of the firefly luciferase reagent (LARII) was added to the test sample, followed by a 10 s equilibration period and luminescence measurement (10 s integration time). This was followed by the addition of 50 µL of the REN reagent and firefly quenching (Stop and Glow TM buffer), another 10 s equilibration period, and luminescence measurement (10 s integration time). The data were then presented as the ratio of firefly to Renilla luciferase activity (Fluc/Rluc). Each data point was obtained from a minimum of three biological replicates, and a total of 15 repeats were performed for each assay to ensure statistical reliability and accuracy. Primers were designed using NCBI Primer‐BLAST and were listed in Table S14 (Supporting Information).
RNA Extraction and RNA‐seq
Developing seeds from wild‐type (WT) and mutant plants grown in Wuhan were subjected to RNA‐seq analysis. Seeds at 40 days after flowering (DAF) were collected, with three biological replicates per genotype. Detailed sampling procedures were described previously.^[^ 37 ^]^ Total RNA was extracted using the RNAprep Pure Plant Kit (TIANGEN, China). RNA‐seq libraries were prepared with the NEBNext Ultra RNA Library Prep Kit for Illumina (NEB, USA; Cat. #E7530L) and sequenced on an Illumina HiSeq 4000 platform.
Quality control of sequencing data was performed using MultiQC,^[^ 96 ^]^ and transcript quantification was conducted with STAR^[^ 97 ^]^ against the Darmor v5 reference genome annotation. Gene expression levels were estimated using the RSEM.^[^ 98 ^]^ Differential gene expression analysis was performed with the R package DESeq2.^[^ 99 ^]^
Protoplast Isolation, CUT&Tag Assay, and CUT&Tag‐qPCR
Transient transformation of rapeseed leaf protoplasts was performed for CUT&Tag assays using the Hyperactive Universal CUT&Tag Assay Kit (Vazyme, TD904). Rapeseed protoplast isolation was adapted from the established protocol for Arabidopsis leaf protoplast extraction.^[^ 95 ^]^ The full‐length coding sequence of BnaA03.TGA6 was cloned into the pM999‐GFP expression vector under the control of the CaMV 35S promoter via the XbaI restriction site. The recombinant plasmid was introduced into rapeseed protoplasts by polyethylene glycol (PEG)‐mediated transformation. After incubation for 16 h, protoplasts exhibiting GFP fluorescence were collected.
The CUT&Tag sequencing data analysis method is adapted from a previous paper.^[^ 100 ^]^ In short, the successfully constructed libraries were subjected to high‐throughput sequencing on the Illumina NovaSeq platform by Novogene (Beijing, China). Clean reads were aligned to the ZS11 v0 reference genome using BWA (v0.7.17) with the MEM algorithm.^[^ 101 ^]^ Peak calling was performed with high confidence using MACS2 (v2.2.7.1) under the parameters “–gsize 1e9 –tsize 150 ‐B –nomodel –shift 100 –extsize 200 –qvalue 0.01 ‐f BAMPE”.^[^ 102 ^]^ Peak regions were visualized with pyGenomeTracks.^[^ 103 ^]^ Motif analysis was conducted using FIMO and meme‐chip from the MEME Suite.^[^ 104 ^]^ Immunoprecipitated DNA was subjected to CUT&Tag‐qPCR to validate the enrichment of target regions, using SYBR Premix TransStart Green qPCR SuperMix (TransGen, China) on a Bio‐Rad CFX Connect Real‐Time PCR System. A spike‐in DNA was included as a negative control. Primers for CUT&Tag‐qPCR validation were designed based on identified motifs and were listed in Table S14 (Supporting Information).
RT‐qPCR
RT‐qPCR was performed using SYBR Green chemistry on a CFX Connect Real‐Time PCR System (Bio‐Rad, USA). Gene expression levels were normalized to ACTIN7 (BnaA06G0089200WE). Reactions were carried out in 96‐well Hard‐Shell PCR plates (Bio‐Rad, USA) using 0.4 µm gene‐specific primers and SYBR Premix TransStart Green qPCR SuperMix (TransGen, China) in a final volume of 15 µL. Each sample included three biological replicates and three technical replicates. The thermal cycling conditions were: initial denaturation at 95 °C for 1 min, followed by 45 cycles of 95 °C for 15 s, 60 °C for 15 s, and 72 °C for 20 s. Primers were designed using NCBI Primer‐BLAST and were listed in Table S14 (Supporting Information).
Electrophoretic Mobility Shift Assay (EMSA)
The full‐length coding sequence of BnaA03.TGA6 was amplified from ZS11 cDNA and cloned into the pET15D vector. The resulting His‐BnaA03.TGA6 fusion protein was expressed in Escherichia coli Rosetta cells induced with 0.4 mM IPTG at 16 °C for 16 h. The fusion protein was purified by NTA‐Ni.
Promoter fragments of BnaA08.DA1 and BnaA07.DA2, containing predicted binding motifs, was synthesized and labeled with FAM at the 5′ end. Forward and reverse oligonucleotides were annealed by heating at 98 °C for 10 min, followed by gradual cooling to room temperature to form double‐stranded probes. Both FAM‐labeled and unlabeled probes were incubated with purified His‐BnaA03.TGA6 protein in EMSA/Gel‐shift Binding Buffer (Beyotime, GS005) at 23 °C for 30 min.
Binding reactions were resolved on a 6% (w/v) native polyacrylamide gel at 80 V for 1 h in 0.5× TBE buffer (45 mm Tris‐base, 45 mm boric acid, 0.5 mm EDTA, pH 8.3) at 4 °C in the dark. Fluorescent signals from FAM‐labeled probes were detected using an Amersham Typhoon imaging system (Cytiva). All oligonucleotide sequences used for EMSA were listed in Table S14 (Supporting Information).
Conflict of Interest
The authors declare no conflict of interest.
Author Contributions
L.L., Z.T., and X.W. contributed equally to this work. L.G., W.C., and X.Y. designed and supervised this study. L.L., X.W., and J.C. performed the experiments. Z.T., Z.T., H.Z., X.H., Y.X., and L.L. analyzed the data. L.L. and Z.T. wrote the manuscript. L.G., W.C., and X.Y. revised the manuscript. All authors read and approved the manuscript.
Supporting information
Supporting Information
Supplemental Table 1
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1A. R. Fernie , R. N. Trethewey , A. J. Krotzky , L. Willmitzer , Nat. Rev. Mol. Cell Biol. 2004, 5, 763.15340383 10.1038/nrm 1451 · doi ↗ · pubmed ↗
- 2J. Venegas‐Molina , F. J. Molina‐Hidalgo , E. Clicque , A. Goossens , Trends Plant Sci. 2021, 26, 472.33478816 10.1016/j.tplants.2020.12.008 · doi ↗ · pubmed ↗
- 3R. C. Meyer , M. Steinfath , J. Lisec , M. Becher , H. Witucka‐Wall , O. Törjék , O. Fiehn , Ä. Eckardt , L. Willmitzer , J. Selbig , T. Altmann , Proc. Natl. Acad. Sci. USA 2007, 104, 4759.17360597 10.1073/pnas.0609709104 PMC 1810331 · doi ↗ · pubmed ↗
- 4R. Sulpice , E.‐T. Pyl , H. Ishihara , S. Trenkamp , M. Steinfath , H. Witucka‐Wall , Y. Gibon , B. Usadel , F. Poree , M. C. O. Piques , M. Von Korff , M. C. Steinhauser , J. J. B. Keurentjes , M. Guenther , M. Hoehne , J. Selbig , A. R. Fernie , T. Altmann , M. Stitt , Proc. Natl. Acad. Sci. USA 2009, 106, 10348.19506259 10.1073/pnas.0903478106 PMC 2693182 · doi ↗ · pubmed ↗
- 5S. Bijlsma , I. Bobeldijk , E. R. Verheij , R. Ramaker , S. Kochhar , I. A. Macdonald , B. van Ommen , A. K. Smilde , Anal. Chem. 2006, 78, 567.16408941 10.1021/ac 051495 j · doi ↗ · pubmed ↗
- 6N. Carreno‐Quintero , A. Acharjee , C. Maliepaard , C. W. B. Bachem , R. Mumm , H. Bouwmeester , R. G. F. Visser , J. J. B. Keurentjes , Plant Physiol. 2012, 158, 1306.22223596 10.1104/pp.111.188441 PMC 3291263 · doi ↗ · pubmed ↗
- 7P. T. Do , M. Prudent , R. Sulpice , M. Causse , A. R. Fernie , Plant Physiol. 2010, 154, 1128.20841452 10.1104/pp.110.163030 PMC 2971594 · doi ↗ · pubmed ↗
- 8T. Obata , A. R. Fernie , Cell. Mol. Life Sci. 2012, 69, 3225.22885821 10.1007/s 00018-012-1091-5PMC 3437017 · doi ↗ · pubmed ↗