Regulatory Mechanisms of Branch Angle in Brassica rapa L. in Xizang Based on RNA-seq and Targeted-Metabolome Analyses

1 Introduction

Rapeseed is an important oil crop worldwide. The primary objectives of rapeseed breeding and cultivation in China are to increase seed yield and achieve mechanized production [1]. Ideal rapeseed plant characteristics include a compact plant type, reduced height, short branches, an increased number of branches, short silique wall length, oblique growth, tolerance to dense planting, strong resilience, resistance to lodging, suitability for mechanical harvesting, an appropriate yield structure, a high harvest index, and high-quality seeds [2,3,4,5]. Among these, the development of a compact plant type is the key focus for genetic improvement. By modifying the branch angle of rapeseed, growth and development can be optimized under dense planting conditions, which is the foundation for increasing yield under high-density planting conditions [6]. The compact branch type is crucial for the mechanized production of rapeseed. Therefore, understanding the mechanisms controlling rapeseed branch angles is vital for breeding ideal plants.

Plant branching significantly influences plant morphogenesis, as it involves the transformation of axillary buds in the leaf axils into branches. The axillary meristem (AM) in the leaf axils initiates and forms the axillary buds [7]. Rapid growth and development of axillary buds into branches are crucial for seed yield, as they are driven by the initiation and growth of axillary buds. The angle of the branches is crucial for shaping the compactness of the rapeseed canopy configuration [8,9]. A moderately compact plant form can enhance planting density, minimize shading between branches, and improve ventilation and light transmission in the middle and lower layers of the plant population. Furthermore, reducing the branch angle can minimize branch entanglement during mechanical harvesting, thereby facilitating machine operations in the field [6]. This reduction could also decrease the likelihood of silique cracking and dislodgement caused by entanglement. Interactions between various environmental factors and endogenous substances influence the complex spatiotemporal regulation of multiple genes during branch formation [10,11,12]. Environmental factors, including light intensity, light quality, and photoperiod, play significant roles in this process. Endogenous factors such as hormones, sugars, and nutrients also exert critical regulatory effects on branching development [13,14,15]. The branching angle of rapeseed is primarily determined by the gravitropism of the aboveground tissues. Abnormal gravitropism can result in changes in branch angle. Gravitropism in plants can be divided into four stages: gravity perception, gravity signal transduction, asymmetric auxin distribution, and bending of gravitropic organs [16,17]. The asymmetric distribution of auxin on both sides of a gravitropic organ is a key step in plant gravitropism. After gravitational stimulation, a higher auxin concentration on the side of the flower cluster and hypocotyl closer to the ground promotes upward bending, whereas the roots were inhibited and bent downward [18]. Spraying auxin transport inhibitors disrupts the asymmetric auxin distribution, leading to abnormal plant gravitropism [19,20,21]. Auxin synthesis, signal transduction, metabolism, and transport pathways influence the establishment of auxin gradients [22]. Rapeseed plants are tall with numerous branches. If the branch angles are wide, the branches of the entire plant loosen, causing individual plant branches to become entangled and occupy a large space that hinders air circulation and light capture [23]. The photosynthetic efficiency per unit area of the plant decreases and the possibility of pest and disease outbreaks increases. This is not conducive to mechanized harvesting, which severely affects crop yield [24]. The loose plant architecture of rapeseed also increases the load on stems and roots, making the crops more susceptible to lodging, thereby correspondingly reducing yield [1]. Reducing the branching angle in rapeseed significantly enhances planting density and improves photosynthetic efficiency and light energy utilization, which are crucial for achieving high and stable rapeseed yields. The SNP loci affecting the branch angle were examined using the 60 K Infinium SNP chip from rapeseed (Brassica napus). In previous genome-wide association studies on B. napus, genes associated with these loci include gravity response genes (SGR2/4/9, BnaA04g09380, BnaC04g31610, and BnaA10g26980) located at A3 and A10 [20,21]. The homologous gene BnaA10g19550 of LAZY1 is located at A10. The auxin synthesis gene TAR2, receptor gene AFB5, signal pathway gene ARF10, and transport-related gene PIN3/7 were also identified [23,24]. The Xizang Autonomous Region, located on the Qinghai–Tibet Plateau, is characterized by diverse climatic conditions, uneven water resource distribution, and marked variations in soil fertility. Rapeseed yield is constrained by natural environmental factors and the inherent traits of the varieties. Rapeseed (Brassica rapa L.) is the primary oil crop cultivated in Xizang and a key source of edible vegetable oil for the local population. Rapeseed exhibits strong resistance and tolerance to poor soil conditions during early maturity. However, it lacks lodging resistance and is unsuitable for mechanical harvesting because of its relatively thin stems, poor stem toughness, and irregular structure. Local rapeseed exhibits tolerance to poor soil conditions and strong resistance and is a rich germplasm resource. However, its yield is relatively low, the plant shape is poor, and the branches are loose, which is not conducive to actual production needs. Relevant research on branch angles in rapeseed is lacking. Therefore, it is essential to increase yield and address the research gap by improving plant type. In this study, the PBH of compact rapeseed (with a branch angle of 25°) was collected from previous projects, providing a material foundation for research on compact plant types and mechanized harvesting of rapeseed in Xizang. Based on compact rapeseed PBH, the endogenous hormone content of the main stem and different stem and branch sites was detected during the branch formation period. RNA-seq was used to identify changes in the expression of differentially expressed genes (DEGs) at different stem and branch sites in the main stem and branches. The results provide a theoretical basis for studying endogenous hormone contents at different stem and branch sites of the main stem and branches of various rapeseed plant types. It also provides guidance for the mechanized harvesting and breeding of high-altitude rapeseed.

2 Materials and Methods

2.1 Plant Growth and Materials Collection

Two distinct B. rapa L. cultivars (compact and dispersed plant types) were collected from a field (No. 3 experimental field, Agricultural Research Institute, Xizang Academy of Agriculture and Animal Husbandry Sciences, Xizang Autonomous Region, China). The S cultivar (PBH, compact plant type) and the F cultivar (dispersed plant type) were used in this study. Two rapeseed varieties were planted, and at the branching stage, different stem materials were collected for testing. Here, F1 and S1, F2 and S2, F3 and S3, F4 and S4, and F5 and S5 represent the top of the main stem, the midpoint between the top of the main stem and the junction of the first branch and main stem, the top of the first branch, the midpoint of the first branch, and the junction of the first branch and main stem in compact and dispersed plant types, respectively. Three biological replicates were collected, frozen in liquid nitrogen, and then stored at −80°C.

2.2 Detection of Endogenous Hormone Content

The endogenous hormone content was analyzed using Agilent 1260 HPLC system (Agilent Technologies Inc., Palo Alto, CA, USA) [25]. For HPLC analysis, the standards (GA (gibberellin), ABA (Abscisic Acid), IAA, CTK (cytokinin), SA (salicylic acid), and JA (jasmonic acid)) were purchased from Sigma–Aldrich (USA), and an Agilent C18 column (250 mm × 46 mm, 5 μm) was used with the Agilent 1260 HPLC system. The mobile phase consisted of methanol and water (40:60, containing 1% acetic acid) at a flow rate of 1.0 mL/min. The column temperature was set at 30°C, the injection volume was 10 μL, and the detection wavelength was set at 254 nm. The analytes were quantified using an external standard.

2.3 Metabolite Profiling by LC–MS

The samples were placed in a freeze-dryer (Scientz-100F, Scientz, Ningbo, China) for vacuum freeze-drying and were ground using a grinder (MM400, Retsch, Munich, Germany). After the samples were dissolved in methanol extract, the supernatant was centrifuged and used for SHIMADZU Nexera X2 UPLC-MS/MS (Shimadzu Corporation, Kyoto, Japan) analysis. The data acquisition instrument system primarily included UPLC and tandem mass spectrometry (MS/MS) (AB Sciex 4500 QTRAP, AB Sciex, Massachusetts, USA). Based on the database (MWDB), qualitative substance analysis was conducted using secondary mass spectrometry, and metabolites were quantified using the multi-reaction monitoring mode of triple quadrupole mass spectrometry. Spectral analysis data for metabolic substances across various samples were collected, followed by peak area integration and the necessary corrections. Quality control samples, prepared by combining extracts from rapeseed plant samples, were utilized to assess the repeatability of the samples under identical treatment conditions. During instrumental analysis, one quality control sample was included after every ten detection and analysis samples to monitor analytical repeatability. Multivariate statistical analysis was used to retain the original information to the greatest possible extent. A digital model was established after simplifying and reducing the data dimensions, standard curves, retention times, LOD/LOQ, recovery rates, repeatability, and internal standards were in File S1. The built-in statistical prcomp function of R software was used to compare analyses between groups and differences between sample groups. A heatmap was generated using the R software pheatmap, and hierarchical cluster analysis was conducted on metabolite accumulation in various samples. OPLS-DA was applied to extract the components of the independent variables X and Y, followed by differential variable screening. Based on the OPLS-DA results, the VIP of the multivariate analysis OPLS-DA model was obtained, and differential metabolites were further screened by combining the p-value and fold change. Permutation tests (200 permutations) were performed to avoid overfitting. The significantly regulated metabolites in each group were screened for differential metabolites with VIP ≥ 1, fold change ≥ 2, and fold change ≤ 0.5. Metabolite compounds were identified using the KEGG database, and metabolite annotations were mapped to the KEGG Pathway database [26]. Pathways with significant metabolite regulation were subjected to metabolite set enrichment analysis, and their significance was determined by the p-value of the hypergeometric test.

2.4 RNA Extraction, Illumina Sequencing, and Data Analysis

Approximately 0.3 g of each sample was subjected to total RNA extraction using TRIzol reagent (Invitrogen, Carlsbad, CA, USA). The quality of the isolated RNA was evaluated on an Agilent 2100 Bioanalyzer, employing oligo(dT) beads in accordance with the manufacturer’s protocol. For RNA library construction, the NEBNext Ultra RNA Library Prep Kit for Illumina (NEB #7530, New England Biolabs, Springfield, MA, USA) was utilized, yielding a total cDNA amount of no less than 0.1 μg. Subsequently, the prepared libraries were sequenced on the Illumina NovaSeq 6000 platform by Gene Denovo Biotechnology Co., which generated paired-end reads of 100 bp. Raw data were processed to acquire clean reads by filtering with Fastp (version 0.19.3). During this step, any reads containing ambiguous bases (N) or those deemed to be of low quality were discarded. A genomic index was constructed based on the B. rapa L. reference genome (NCBI accession: GCF_000309985.2) [27]. The clean paired-end reads were then aligned to this reference using the HISAT-3N aligner, which facilitated the acquisition of positional information and unique sequence features for subsequent assembly. Transcript assembly for each sample was performed with StringTie (v1.3.4) based on the reference genome. Gene expression levels were quantified using the featureCounts package (v1.6.2) within the R environment. To determine expression abundances, the FPKM (fragments per kilobase of transcript per million mapped reads) metric was calculated using RSEM, a software tool designed for precise transcript quantification from RNA-seq data. Identification of DEGs between experimental groups was carried out using the DESeq2 package. Statistical significance was determined using a threshold of FDR ≤ 0.01, |log2 fold change| ≥ 2, and FPKM ≥ 2. Functional annotation of the identified DEGs was performed via BLAST version 2.28.2+ searches against several public databases, including KEGG, GO, and KOG.

2.5 Combined Analysis of Transcriptome and Metabolome

To explore the interplay between transcriptional and metabolic profiles, DEMs and DEGs derived from the same tissue samples were integrated and mapped onto KEGG pathways. This approach aimed to uncover potential regulatory relationships between specific genes and metabolites. To visually compare the pathway enrichment patterns, bar charts were generated based on the enrichment analysis results of both DEMs and DEGs, highlighting differential enrichment levels between the two omics layers. For each defined subgroup, pairwise correlation analysis between the DEG and DEM datasets was conducted. Specifically, Pearson correlation coefficients were calculated for all gene–metabolite pairs using the cor function implemented in the R programming environment. To further filter tightly correlated pairs, nine-quadrant plots were employed to visualize the relationship between gene expression/log2FC and metabolite/log2FC. Only those gene–metabolite pairs exhibiting an absolute Pearson correlation coefficient greater than 0.8 within each comparative group were retained for further analysis. An integrative analysis was performed by constructing a two-way orthogonal projection to latent structures (O2PLS) model incorporating all identified DEMs and DEGs. The resulting loading plots were then examined to pinpoint variables exhibiting significant correlation and high weight across the transcriptomic and metabolomic datasets, thereby facilitating the selection of key features with substantial inter-omic influence.

2.6 qRT-PCR

To verify the reliability of genes encoding DEGs, qRT-PCR was performed in a 96-well plate on a CFX96 Touch Real-Time PCR system (Bio-Rad, Hercules, CA, USA). Actin2 was used as the internal standard for qRT-PCR. Thermal cycling was performed at 95°C for 10 s, followed by 40 cycles at 95°C for 10 s, 60°C for 30 s, and 75°C for 15 s. Gene-specific primers (https://www.ncbi.nlm.nih.gov/tools/primer-blast) were designed. The relative expression levels were calculated using the 2−∆∆CT method.

2.7 Cloning of LAZY1 Promoters and Prediction of Cis-Acting Elements Using PlantCARE

The sequences of the LAZY1 promoter were extracted from the reference genome using SPDEv1.2 [28]. The cis-acting elements of these promoters were predicted using the PlantCARE website. Based on the published nucleotide sequence of the LAZY1 (LOC103846512) promoter (Table S1) in the reference genome, primers were designed using Primer Premier 6.0 by focusing on the highest homology to the B. rapa L. sequence to amplify the promoter sequence. The cDNAs were used as templates for RT-PCR, with the following PCR amplification program: pre-denaturation at 98°C for 3 min; denaturation at 98°C for 15 s; annealing at 60°C for 20 s; extension at 72°C for 35 s for30 cycles; and a final extension at 72°C for 5 min for 30 cycles. RT-PCR products were stored at 4°C until further analysis. Subsequently, the products were purified and sent to AUG Biotech (AUG Biotech. Ltd., Shanxi, China) for sequencing.

2.8 Statistical Analysis

Statistical analyses and data visualization were conducted using Microsoft Office 2021 (Microsoft, Washington, DC, USA) and GraphPad Prism 10.1 (GraphPad Software Inc., San Diego, CA, USA).

3 Results

3.1 Changes in Branch Angle with Germplasm Resources (B. rapa L.) in Xizang

The branch angles of the 54 germplasm resources collected from Xizang (Table S2) were determined. Most rapeseed resources exhibited large branch angles (≥35°), with PBH (S) being a compact type with a branch angle of 25° (Fig. 1b, Fig. S1A). In contrast, the branch angle for 188028 (F) was 55°, classifying it as a loose line (Fig. 1a, Fig. S1B). Therefore, PBH provides important breeding material for developing compact-type rapeseed (B. rapa L.).

Figure 1: Changes in branch angle in germplasm resources (Brassica rapa L.). (a–g) The first and second branch angles in germplasm resources (188028, PBH, Deqingyou, Laziyou, Sa’eryou, Mouzhuyou, and 2021541172).

3.2 Differences in Endogenous Hormone Contents in Different Branch Sites of F and S Plant Types

In this study, the contents of six endogenous hormones (IAA, CK, ET, ABA, GA, and JA) (Fig. 2B) were determined in the S and F varieties at five locations on the stem (top of the main stem; midpoint between the top of the main stem and the junction of the first branch and main stem; top of the first branch; midpoint of the first branch; and junction of the first branch and main stem) (Fig. 2A). Except for ET, in S and F, the trends in the levels of the five endogenous hormones (IAA, CK, ABA, GA, and JA) were consistent across various stem and branch locations, peaking at the top of the main stem and reaching the lowest levels at the junction of the first branch and main stem. In addition, the levels in S were significantly higher than those in F. This suggests that endogenous hormones may influence the branching angle of rapeseed. IAA, GA, and ABA are crucial in determining the branch angle of rapeseed. The contents of IAA, GA, and ABA varied greatly between S and F, which led us to focus on the expression levels of DEGs and transcription factors related to the synthesis, transport, and degradation of IAA, GA, and ABA in transcriptome data analysis.

Figure 2: Detection of branch sites and endogenous hormone content. ①: Top of main stem; ②: Midpoint between the top of main stem and the junction of the first branch and main stem; ③: Top of the first branch; ④: Midpoint of the first branch; ⑤: Junction of first branch and main stem.

3.3 Differentially Expressed Metabolites (DEMs) and Enrichment at Different Branch Sites of F and S Plant Types

PCA can reveal the overall metabolic differences between samples across various groups, as well as the extent of variation within each group. OPLS-DA (Fig. 3A) and correlation analysis (Fig. 3B) were performed on the F and S plant types, revealing that the R2 values among biological replicates of F and S were greater than 0.60, which indicates significant differences in gene expression and shows high biological reproducibility for each sample (Fig. 3A). Based on the targeted metabolome of endogenous hormones, VIP and fold changes were used to screen DEMs. The number of DEMs in the different comparison groups (F1 vs. S1, F2 vs. S2, F3 vs. S3, F4 vs. S4, and F5 vs. S5) was analyzed (Fig. 3C, Table 1). The number of upregulated DEMs in F4 vs. S4 and F5 vs. S5 was much greater than that of the downregulated DEMs, indicating that, in the S plant type, branch angle was regulated by the upregulated DEMs. By assessing the up- and downregulation dynamics of DEMs, it was observed that the levels of the auxin hormone (IAA-glu-diME, IAA-Asp, and ICAld) and SA (SAG) contents (Fig. 3D) were significantly increased, whereas the levels of CK (DHZR) and JA (H2JA) contents (Fig. 3D) were significantly decreased in S. Combined with the endogenous hormone content (Fig. 2B) and the significance of DEMs (Fig. 3), it was evident that auxin might play a role in regulating branch angle. These results provide the basis for RNA-seq analysis.

Figure 3: Basic analysis of targeted-metabolome data. (A) The principal component analysis (PCA) score chart of RNA-seq data for all quality control (QC) samples. (B) Correlation heat map. Pearson’s correlation coefficient (R2) > 0.7 between the three biological replicate samples. (C) Venn diagram of differentially expressed metabolites (DEMs). Abscissa indicates different samples. Ordinate indicates the logarithmic values of the sample content log2Fold Change. (D) Heat map of DEMs (Key endogenous hormones). IAA-Glu-diMe: Indole-3-acetyl-L-glutamic acid dimethyl ester, IAA-Asp: Indole-3-acetyl-L-aspartic acid, ICAld: Indole-3-carboxaldehyde, DHZR: Dihydrozeatin ribonucleoside, H2JA: Dihydrojasmonic acid, SAG: Salicylic acid 2-O-β-glucoside.

3.4 RNA-seq Analysis of Different Branch Sites in F and S Plant Types

RNA-seq analysis was performed on samples from five stem and branch sites (top of main stem; midpoint between the top of the main stem and the junction of the first branch and main stem; top of the first branch; midpoint of the first branch; and junction of first branch and main stem) of S and F. After filtering the raw data, assessing the sequencing error rate, establishing the GC content distribution, 179.80 Gb of clean data were obtained. More than 5.4 Gb of clean reads with Q30 > 91% (Table S3) were obtained for each sample. The proportion of total mapped reads showed at least 93% similarity to the B. rapa L. genome, and uniquely mapped reads exceeded 75%. Thus, the RNA-seq data obtained were accurate and met the requirements for further analysis. PCA of the S and F stems (Fig. 4A) revealed consistency among the three biological replicates. The correlation between the biological replicates was >0.75. Normalized FPKM of the DEGs were extracted for hierarchical clustering analysis, and clustering heat maps of each group were generated. The results indicated significant differences in gene expression among various stem samples (Fig. 4B). These results can be used to identify the DEGs.

Following the analysis of DEGs (FDR ≤ 0.01, |log2Fold Change| ≥ 2, and FPKM ≥ 2) using DESeq2, the total number of DEGs, including those that were upregulated and downregulated (Table 2), was determined for the different rapeseed stem comparisons (F1 vs. S1, F2 vs. S2, F3 vs. S3, F4 vs. S4, F5 vs. S5). The UpSet plot reflects a relatively sufficient number of DEGs across different stem and branch locations of the S and F stems (Fig. 4C, Fig. S2). The number of DEGs in stems from identical parts of S and F plants was examined. The number of DEGs in F1 vs. S1 was 3999 (top of main stem, 1451 up- and 2548 downregulated); the number of DEGs in F2 vs. S2 was 4232 (midpoint between the top of the main stem and the junction of the first branch and main stem, 1612 up- and 2620 down-regulated); the number of DEGs in F3 vs. S3 was 4440 (top of the first branch, 1521 up- and 2919 downregulated); and the number of DEGs in F4 vs. S4 was 3896 (midpoint of the first branch, 1556 up- and 2340 downregulated). F5 vs. S5 showed 4070 DEGs (junction of the first branch and main stem, 1559 up- and 2511 downregulated). When combined with the targeted metabolome results of endogenous hormones, the transcriptome data were analyzed for hormone content variations at different stem and branch sites of the two plant materials. This analysis yielded three trend profiles of DEGs that aligned with the trends of the DEMs (3, 5, and 15) (Fig. 3D). These findings enhance the reliability of screening for DEGs related to branch angle formation.

Based on the trend analysis, each comparison group was annotated using GO (Fig. S3) and KEGG annotations (Fig. S4). GO terms were mainly enriched in CC (cellular component) (GO:0009986, F1 vs. S1:1658 DEGs, 487 up- and 1171 downregulated; F2 vs. S2:1910 DEGs, 674 up- and 1236 downregulated; F3 vs. S3:1994 DEGs, 649 up- and 1345 downregulated; F4 vs. S4:1771 DEGs, 648 up- and 1123 downregulated; F5 vs. S5:1833 DEGs, 623 up- and 1210 downregulated). The BP (biological process) (GO:0008150, F1 vs. S1:1789 DEGs, 612 up- and 1177 downregulated; F2 vs. S2:1916 DEGs, 686 up- and 1230 downregulated; F3 vs. S3:2013 DEGs, 670 up- and 1343 downregulated; F4 vs. S4:1806 DEGs, 669 up- and 1127 downregulated; F5 vs. S5:1831 DEGs, 641 up- and 1190 downregulated) and other 26 significantly enriched terms were also identified. The MF (molecular function) category was significantly enriched in 13 terms, such as binding (GO:0005488) and catalytic activity (GO:0003824), and cellular component was significantly enriched in 20 terms including cell part (GO:0044464) and cell (GO:0005623).

The biological function classification of DEGs related to branch angle was analyzed using KEGG. The top 20 KEGG enrichment analysis (selected pathway terms of DEGs >2, with each term ranged by −log10 p-value) showed that, at the same stem sites, the MAPK signaling pathway (ko04016), phosphatidylinositol signaling system (ko04070), and plant hormone signal transduction (ko04075) were significantly enriched KEGG terms (Table S4). The number of DEGs with plant hormone signal transduction and carbon metabolism was significantly higher than that in other pathways (Table 3), indicating that hormones play a role in regulating the branch angle in rapeseed.

Figure 4: Basic analysis and the trends of gene expression from RNA-seq. (A) PCA score chart of RNA-seq data for all QC samples; (B) Correlation heatmap; (C) UpSet plot of different groups; (D,E) Enrichment scatter diagrams of S and F. Black lines represent the trends of gene expression of S and F at different branch sites. Asterisks indicate the statistical significance based on an independent samples t-test: *p < 0.05, **p < 0.01, ***p < 0.001.

3.5 The Expression Levels of DEGs in Pathways Related to Auxin Metabolism

By examining the endogenous hormone content, KEGG enrichment, and DEMs, DEGs associated with hormone biosynthesis, transport, receptors, and degradation were scrutinized. The DEGs linked to auxin transport (Fig. 5A), including the PIN family genes (LOC103857147, LOC103857314, LOC103830450, LOC103872803, LOC103855472, LOC103873246, LOC103839184, LOC103827806, LOC103850554, LOC103830548, and LOC103860563), auxin responsive genes (Fig. 5D), SGR family genes (LOC103871580, and LOC103861824) were identified at the top of the main stem, the midpoint between the top of main stem and the junction of the first branch and main stem, the top of first branch, the midpoint of the first branch, and the junction of the first branch and main stem. The DEGs were all significantly higher in S plant types than in F. Conversely, the DEGs related to IAA transport and response (Fig. 5B), SAUR family genes (auxin-responsive, involved in cell expansion, Fig. 4B, LOC103852693, LOC103835461, LOC103831754, and LOC103875218), and LA family genes (Fig. 5C, LOC103850600, and LOC103836603) were significantly downregulated DEGs in S compared to F. It is evident that the expression of IAA synthesis and receptor genes in S was upregulated at the top of the main stem, the midpoint between the top of the main stem and the junction of the first branch and main stem, the top of the first branch, the midpoint of the first branch, and junction of the first branch and main stem, resulting in a significantly higher accumulation of IAA in these parts of S than in F. Similar findings have been documented in previous studies [14,29,30,31].

Figure 5: Fold changes of DEGs related to IAA metabolism. (A) Heat map of PIN family genes, (B) SAUR family genes, (C) LA family genes, and (D) SGR family genes.

3.6 DEGs in Pathways Involved in Regulation of Branch Angle

The expression levels of DEGs (LAZY1 (LOC103846512, LOC103828711, and LOC103873494)) associated with hormone biosynthesis, transport, receptors, and degradation were analyzed to understand their impact on the regulated branch angle (Fig. 6). This assessment revealed that in S, the LAZY1 (LOC103846512, LOC103828711, and LOC103873494) expression levels were significantly upregulated at top of the main stem, the midpoint between the top of the main stem and the junction of the first branch and main stem, the top of the first branch, the midpoint of the first branch, and the junction of the first branch and main stem compared to F. However, LOC103873494 did not show this increase at the junction of the first branch and main stem. The results indicate that LAZY1 plays a role in branch angle in rapeseed.

Figure 6: Fold changes of LAZY1 related to branch angle.

3.7 Combined Transcriptome and Metabolomics Analyses of the F and S Branch Sites

Based on the variations in IAA content, the expression of DEGs associated with IAA, and LAZY1 expression levels, we illustrated the pathway by which IAA regulates branch angle in rapeseed (Fig. 7). The upregulation of auxin transport DEGs (Fig. 4A; PIN family genes (LOC103857147, LOC103857314, LOC103830450, LOC103872803, LOC103855472, LOC103873246, LOC103839184, LOC103827806, LOC103850554, LOC103830548 and LOC103860563)) enabled S to accumulate more IAA. Conversely, the downregulated expression of DEGs related to IAA transport (SAUR family genes (Fig. 4B; LOC103852693, LOC103835461, LOC103831754 and LOC103875218)) and downregulation of LA family genes (Fig. 4C; LOC103850600 and LOC103836603) indicated that IAA accumulated at the top of stem (top of the branch) during synthesis. IAA promoted the upregulation of LAZY1 (Fig. 7, LOC103846512, LOC103828711 and LOC103873494) in controlling the branch angle. Modulation of the LAZY1 gene, influenced by genes involved in IAA and branch angle regulation, contributed to the reduction of the S branch angle in rapeseed, which is beneficial for mechanical harvesting and the breeding of ideal plant types.

Figure 7: Proposed model for the role of DEGs in the regulation of branch angle in S and F. In picture, heat map represent log2 Fold Change of DEGs, and the line graph shows the IAA content at different tissue locations. Different lowercase letters indicate significance at 0.05 levels.

3.8 Predicted Anlysis of Cis-Acting Elements in LAZY1 Promoters

Previous studies have shown that upregulation of LAZY1 decreases the branch angle of Arabidopsis thaliana and B. napus L., resulting in a more compact plant type. Based on the variations in IAA content, RNA-seq data, and targeted metabolome analysis, LAZY1 were significantly upregulated in S plants, with the branch angle being smaller than that of F. Compared to F, LAZY1 expression was significantly elevated in various parts of S, including the top of the main stem, the midpoint between the top of the main stem and the junction of the first branch and main stem, the top of the first branch, the midpoint of the first branch, and the junction of the first branch and main stem. Consequently, the CDS and promoter of LAZY1 were cloned. The results showed that in F, the LAZY1 CDS was shifted (deletion) forward by three bases, resulting in the absence of an amino acid sequence from the secondary structure (Fig. 8A,B). However, its tertiary structure (Fig. 8C,D) remained largely unaffected, and the missing base was not located in a conserved sequence or critical region, suggesting minimal impact on the expression of LAZY1. Therefore, the promoter sequences of LAZY1 in F and S were cloned and regulatory elements were predicted (Fig. 8E, File S2, and File S3). The promoter of LAZY1 in S was found to contain a cis-acting element AuxRR core (−112 to −106 bp), associated with IAA responsiveness, whereas base mutations in F altered the relevant element sequences (Fig. 8E, G to T), rendering LAZY1 unregulated by IAA. This insensitivity to IAA in F resulted in a larger branch angle, potentially explaining why the branch angle in S was smaller than that in F.

RNA-seq analysis and qRT-PCR were performed on randomly selected DEGs to determine the authenticity and reliability of the transcriptome data based on qRT-PCR. The qRT-PCR (Fig. S5) and RNA-seq results were consistent for nine (LOC103856327, LOC103873246, LOC103846512, LOC103828711, LOC103873494, LOC103850600, LOC103836603, LOC103871580, LOC103861824, LOC103836229, LOC117130206, LOC103829166, LOC103837423, LOC103851145 and LOC103842242) of the 15 validated genes. Hence, transcriptome sequencing was reliable.

Figure 8: (A,B) LAZY1 amino acid sequences in S and F; (C,D) tertiary structures in S and F; (E) prediction of cis-acting elements in the promoter.

4 Discussion

The mechanical harvesting of rapeseed represents a significant area of research. The optimal plant type of rapeseed is best suited for mechanical harvesting [32]. High branch position, narrow branch angle (compact plant type), and superior quality are important indicators of an ideal rapeseed plant type [33]. Research on branch angle has predominantly focused on mechanical harvesting of B. napus L., whereas research on the branch angle of B. rapa L. is relatively scarce. This is because the yield, quality, and other attributes of B. napus are far superior to those of B. rapa [1]. However, there are few varieties or resources with smaller branch angles in B. rapa, which poses challenges for related research and breeding. Despite this, rapeseed (B.rapa L.) offers notable advantages such as delayed sowing, a short growth period, tolerance to poor soil, and strong stress resistance, making it widely applicable [34]. Xizang is located on the Qinghai–Tibet Plateau at high altitudes, with an average elevation exceeding 4000 m, this region experiences diverse climatic conditions, uneven distribution of water resources, and considerable variation in soil fertility [35,36]. Rapeseed yield is not only limited by natural conditions but also by the characteristics of the variety itself. Local rapeseed is characterized by adaptation to poor soil, strong resistance, and diverse resources and has great application value [37]. However, the yield is relatively low, the plant type is suboptimal, and the branches are loose, which are not conducive to actual production. Relatively few studies have been conducted on the branching angles of rapeseed [5,10,23].

The pathways involved in the synthesis, signal transduction, metabolism, and transport of auxins in the main stem and branches affect the establishment of the auxin gradient, thereby affecting branch angle [38,39,40]. The asymmetric distribution of auxin is essential for branches to respond to gravitropism, which involves changes in branch angles. An increase in auxin concentration causes the main inflorescence on the side nearer to the ground in A. thaliana to develop smaller branch angles [41]. Tryptophan is a precursor of auxin synthesis. The TAA/TAR family catalyzes the conversion of tryptophan into indole-pyruvic acid, whereas the YUCCA family catalyzes the synthesis of IAA [42]. Tryptophan is an important precursor for auxin synthesis. Studies have shown that gain-of-function mutant yucca-1d in A. thaliana has a smaller branching angle, whereas the double mutant taa-tar has an increased branch angle. Auxin and its transport play key roles in regulating the gravitational responses of branch angles [43]. The PIN gene family encodes a group of auxin efflux carrier proteins. The asymmetric distribution of these proteins on the cell membrane results in the formation of auxin gradients, which subsequently lead to asymmetric cell growth and the development of growth angles [44]. The SAUR family encodes a class of small proteins that rapidly respond to auxins. A recent study showed that SAUR10 was specifically expressed at the distal end of the branch, whereas its expression was inhibited by FUL. The average branch angle of functional-defect mutants decreased by 6°, whereas the average branch angle increased by 10° when FUL was overexpressed [45]. The gene AtLA1 in A. thaliana, which is homologous to LA1 in rice, is involved in the gravity response signaling pathway, thereby regulating the branching angle [46]. In rapeseed, the gene homologous to LA1 is LAZY1. Another important branch angle regulatory gene in dicotyledonous plants, AtTAC1, was identified in a study of branch-angle regulation in peach trees [47]. In A. thaliana, the branch angle of the attac1 mutant was significantly reduced, indicating that, compared with LA1, AtTAC1 promotes an increase in the branch angle of A. thaliana. Low-level expression of PpeTAC1 is associated with upright growth of peach tree branches [48]. In A. thaliana, because of the absence of endoderm in the hypocotyl and stem of the mutants sgr1 and sgr7, the branch angle is relatively small.

Many auxin response genes (AUX1, IAA, GH3, and ARF) were enriched in the compact line [49]. BnaA0639380D is a homolog of AtYUCCA6. Sequence comparison of BnaA0639380D from oilseed rapeseed with small and large branch angles revealed six SNPs and four amino acid variations in the promoter and coding regions, respectively. Four known candidate genes (BnaA02g16500D, BnaA03g10430D, BnaC03g06250D, and BnaC06g20640D) were identified through both GWAS and RNA-seq, all of which play roles in regulating asymmetric auxin distribution [50]. IAA7 contributes to yield heterosis by improving plant architecture and could be beneficial for breeding superior rapeseed hybrid cultivars [51,52]. Such a mutation may increase the yield of other Brassica species. TAC1 and LA1 belong to different evolutionary branches of the insulin-like growth factor family [53]. Structurally, AtTAC1 lacks a conserved EAR class motif with a transcriptional inhibitory function in the C-terminal V domain, unlike LA1. It is hypothesized that TAC1 may have evolved from truncation of the C-terminus of LA1 [54]. TAC1 and LA1 might function antagonistically in auxin transport and signal transduction processes [53]; the absence of the C-terminal sequence of TAC1 or its loss of binding ability to microfilaments may regulate branch angle by inhibiting LA1-induced gravitropism [55]. However, there is currently no evidence linking TAC1 and auxin signaling.

This study measured the top of the main stem and the midpoint between the top of the main stem and the junction of the first branch and main stem in S plants. Both the top of the first branch and midpoint of the first branch were significantly higher than those in F. Additionally, targeted metabolomics showed that IAA (IAA-Glu-DIME, IAA-ASP, and ICAld) levels were significantly higher than those in F. Transcriptome analysis revealed that PIN family genes were related to auxin transport (LOC103857147, LOC103857314, LOC103830450, LOC103872803, LOC103855472, LOC103873246, LOC103839184, LOC103827806, LOC103850554, LOC103830548 and LOC103860563), as well as IAA12 (LOC103870470) and LAZY1, associated with branch angle regulation, were significantly upregulated. The increased expression of these genes resulted in a smaller branch angle in S compared to F. The promoter of LAZY1 in S had a cis-acting AuxRR core element associated with IAA. The presence of the AuxRR core renders S more sensitive to IAA, promotes the upregulation of LAZY1, and results in a smaller branch angle in S than in F.

5 Conclusions

Rapeseed is a primary oil crop cultivated in Xizang. Investigating ideal rapeseed plant types is beneficial for industrial development. In this study, rapeseed plants with small branch angles (25°, PBH, or S-type) were collected from previous projects. Measurements were taken at the top of the main stem, the midpoint between the top of the main stem and the junction of the first branch and main stem, and the top of the first branch. The levels of IAA, ABA, and GA were estimated in five parts, including the midpoint of the first branch and the junction of the first branch and main stem, where IAA (IAA-Glu-DIME, IAA-ASP) and ICAld were significantly higher than in the F plant type. DEGs associated with IAA, such as PIN family genes (LOC103857147, LOC103857314, LOC103830450, LOC103872803, LOC103855472, LOC103873246, LOC103839184, LOC103827806, LOC103850554, LOC103830548, and LOC103860563) and IAA12 (LOC103870470), were significantly upregulated in S and F; LAZY1, which influences branch angle regulation, was also significantly upregulated. The upregulated expression of these genes caused the branch angle of S to be smaller than that of F. The promoter of LAZY1 in S had a cis-acting element, the AuxRR core (a cis-acting regulatory element involved in auxin responsiveness), which is related to IAA. The presence of the AuxRR core makes S more sensitive to IAA, promotes upregulated expression of the LAZY1 gene, and results in a smaller branch angle in S than in the F plant type. The findings of this study elucidate the differences in hormone levels in stems and branches to develop compact plant types for superior yield and machine harvesting.