Extensive phenotypic variations of plant height were observed in 230 inbred accessions over two consecutive years ( Table 1 ). The plant height ranged from 149.23–230.59 cm in 2017–2018 and from 115.12–189.28 cm in 2018–2019, suggesting that the environment factors had a great impact on plant height ( Supplementary Figure 1A and Table 1 ). The plant heights in 2017–2018, 2018–2019, and the BLUP of the 230 rapeseed accessions displayed normal distributions ( Supplementary Figure 1A ). The coefficient of variation in 2018–2019 was 9.10%, which was higher than that in 2017–2018 (7.38%) and BLUP (7.69%) ( Table 1 ). Nevertheless, no significant difference was observed between the phenotype of 2017-2018 and 2018-2019, as shown by the correlation analysis (R² > 0.70) ( Supplementary Figure 1B ). These analyses revealed that the phenotype of 230 rapeseed accessions were reliable and feasible for association analysis.

A total of 230 rapeseed accessions, consisting of 25 spring-, 33 winter-, and 172 semi-winter ecotypes, were employed for WGS ( Supplementary Table 1 ). Approximately 1,097.37 Gb data were generated, with an average size of 4.77 Gb and an average depth of 6.46 × depth per accession ( Supplementary Table 1 ). The average coverage of B. napus reference genome was 82.19% ( Supplementary Table 1 ). A total of 1,707,194 informative SNPs were acquired with an average of 94,844 SNPs on each chromosome ( Table 2 ). The density of SNPs on different chromosomes ranged from 1.01 to 4.86 SNP/kb, with chromosome C09 having the lowest density and A10 the highest ( Table 2 ). These results suggested the reliability of SNP information and could be used for further analyses.

According to the Q+K model, the associated population could be divided into nine subgroups ( Supplementary Figure 2A, B ). More than 90% of the relative kinship coefficients among these accessions were found to be lower than 0.1, suggesting that most accessions in this population lacked or had weak genetic relatedness ( Supplementary Figure 2C ). The average linkage disequilibrium (LD) decay of the A and C sub-genomes were 4.1 and 120.3 kb, respectively. It was 33.4 kb for the whole genome (A + C) when r² decayed to its half ( Figure 1D ).

To dissect the genetic control of plant height in B. napus, we performed GWAS in two consecutive years. A significant locus on chromosome A01 was simultaneously identified using GLM, MLM, and BLINK models ( Figures 1A–C , and Supplementary Table 3 ). Quantile-quantile plots showed obvious deviations between the observed and expected values, indicating the selected models were correct and suitable for GWAS ( Supplementary Figure 3 ). Within the significance interval, 806 SNPs were repeatedly identified in different environments and models (GLM and MLM) ( Supplementary Table 3 ). According to the MLM model in BLUP and LD decay of A sub-genome ( Figure 1 ), three genes (BnaA01g09530D, BnaA01g09540D, and BnaA01g09550D) near the significant SNPs were strongly associated ( Figure 2A ). Based on the annotation of the B. napus reference genome, BnaA01g09530D, encoding BRASSINOSTEROID-INSENSITIVE 2 (BIN2) and involving in the brassinosteroid signaling pathway, may be a candidate gene controlling plant height in B. napus. In addition, the position of co-identified SNP by BLINK model in different environments was 4,772,232 ( Supplementary Table 3 ), which was far away from the co-identified significant SNPs in GLM and MLM, due to the different algorithm principle of BLINK. In the application of BLINK, a bin, containing all the linked SNPs in a region, is taken as a unit, rather than a single SNP as a unit like GLM and MLM (Huang et al., 2019), suggesting that BnaA01g09530D was also identified in the BLINK models. A total of ten SNPs variations were observed in the sequence of BnaA01.BIN2. Haplotype analysis of these ten SNPs revealed favorable allelic variation (Hap_II), conferring a significant reduction in plant height ( Figure 2B ).

The expression pattern of BnaA01.BIN2 showed that it was highly expressed in leaves, buds, and roots, followed by SAM, suggesting that it plays an important role in plant development. Subcellular localization, as indicated by green fluorescent protein (GFP), showed that BnaA01.BIN2 was localized in nucleus and cytoplasm ( Figure 2D ).

Candidate gene BnaA01.BIN2 belongs to the glycogen synthase kinase 3 (GSK3) gene family. Using protein sequences of AtGSK3s as the query of BLAST, a total of 38 BnGSK3s were identified in “Darmor-bzh” rapeseed genome, and 31 BnGSK3s with Pkinase domain were finally extracted ( Table 3 ). Of these BnGSK3s, 16% (5) resulted from dispersed duplications and 84% (26) originated from whole-genome duplication (WGD) or segmental duplication ( Table 3 ). Most AtGSK3s have several syntenic genes in B. napus, among which BnBIN2 contains six homologous genes and is the largest member of BnGSK3s ( Supplementary Table 2 ). However, there were no homologous genes for AtBIL2 and AtSK42 in B. napus ( Supplementary Table 2 ). We identified 15 BrGSK3s and 16 BoGSK3s according to the Brassica Database (BRAD) (http://brassicadb.cn/) and no homologous genes of AtBIL2 and AtSK42 were identified in the reference genomes of B. rapa (Brara_Chiifu_V3.5) and B. oleracea (Braol_JZS_V2.0) ( Supplementary Table 2 ). These results suggested that the GSK3s family is highly conserved in Brassicaceae, whereas the loss of BIL2s and SK42s may occur prior to Brassicaceae speciation.

The 31 BnGSK3s were unevenly distributed in 13 chromosomes and four random chromosomes, 16 and 15 BnGSK3s were located on the A and C sub-genomes, respectively ( Figure 3 and Table 3 ).

To explore the phylogenetic relationship of GSK3s family, we constructed a phylogenetic tree using GSK3s protein sequences from Arabidopsis and B. napus. The 10 AtGSK3s and 31 BnGSK3s were divided into four groups: Group I (SK11, SK12, and SK13), Group II (BIN2, BIL1, and BIL2), Group III (SK31 and SK32), and Group IV (SK41 and SK42) ( Figure 4A and Supplementary Table 2 ). Group I had 22 GSK3s, including 11 BnGSK3s, 5 BrGSK3s, and 6 BoGSK3s, accounting for the largest group. Group IV was the smallest, with only eight GSK3s ( Figure 4A and Supplementary Table 2 ). This suggests that Group I of GSK3s was more expanded compared to that of Group IV. Within each group, BnGSK3s belonging to the A and C sub-genomes in B. napus, along with the AtGSK3s in Arabidopsis, clustered into a small clade ( Figure 4A ), suggesting that the phylogenetic relationship of GSK3 was consistent with the evolution of rapeseed. The syntenic analysis between AtGSK3s and BnGSK3s showed that most of AtGSK3s have over two syntenic genes in B. napus ( Figure 4B ), which is consistent with phylogenetic relationship of GSK3s.

To explore the conservation of BnGSK3s, gene structure and protein motifs were analyzed ( Figure 5 ). In general, gene structures of the 31 BnGSK3s differed obviously between different groups. Among them, the syntenic genes showed relatively similar gene structures ( Figures 5A, C ). The gene structures of approximately 74% of BnGSK3s (23) exhibited 5’ and 3’ untranslated regions (UTR) ( Figure 5C ). Six BnGSK3s (BnA05g31460D, BnA03g05700D, BnA09g38810D, BnA09g51790D, BnCnng52760D, and BnAnng35300D) only had 5-’ or 3-UTR ( Figure 5C ). In addition, all BnGSK3s contained exons and introns ( Figure 5C ). As for motif analysis, except BnaA05g11700D possessed seven conserved motifs, the remaining BnGSK3s had ten conserved motifs ( Figure 5B ). These results suggested that the core sequences of the BnGSK3s were conserved.

Based on published transcriptome data (Li et al., 2019b; Dong et al., 2021), the expression patterns of the 31 BnGSK3s in 13 tissues of ZS11 were analyzed, which showed that the BnGSK3s were expressed in different tissues ( Figure 6A ). However, a set of homologous genes, including BnSK31, BnBIL1, and BnBIN2, showed similar expression patterns, suggesting a potential redundancy of function ( Figure 6A ). Different expression patterns were observed within the same group, suggesting functional divergence in BnGSK3s. For example, BnBIL1 was highly expressed in SAM, whereas BnBIN2 was highly expressed in roots ( Figure 6A ). In addition, BnSK13s were prone to express in pistils and buds, suggesting that these genes may be involved in flower development ( Figure 6A ). Thirteen BnGSK3s were highly expressed in SAM, indicating that BnGSK3s have a certain effect on the development of plant architecture. In addition, we selected eight BnGSK3s from different groups to perform qRT-PCR in six tissues, which suggested that the expression pattern was consistent with the RNA-seq data ( Figure 6B ).

Since the “Green Revolution” in the 1960s, researchers have carried out extensive research to come up with ideal plant architecture models for many crops (Teichmann and Muhr, 2015; Liu et al., 2020; Pearce, 2021). Several novel genes controlling aboveground plant architecture have been identified and their regulatory mechanisms have been expounded, laying the foundation for breeding new high-yielding varieties of rice (Monna et al., 2002; Sasaki et al., 2002; Wang and Li, 2008; Wang et al., 2020b), wheat (Chai et al., 2022; Xiong et al., 2022), maize (Phillips et al., 2011; Li et al., 2020a), and soybean (Guo et al., 2020; Chen et al., 2021b). Many researchers have proposed models of ideal rapeseed plant architecture (Liu et al., 2022; Zheng et al., 2022). However, these are only concepts and cannot solve actual problems in production. There are several difficulties in studying the ideal plant architecture of rapeseed: 1) the lack of materials with good plant architecture materials; 2) the uncertainty of proper index used for the research of rapeseed plant architecture; 3) severe environmental impact on plant architecture-related traits; 4) the lack of clear genetic basis. For many crops, such as rice and wheat, plant height has been used as a breakthrough point to study plant architecture (Peng et al., 1999; Hedden, 2003). Therefore, plant height is essential in shaping the ideal plant architecture of crops. Although research progress has been made in the study of plant height traits of rapeseed (Liu et al., 2010; Wang et al., 2016a; Wang et al., 2016b), the genetic basis of rapeseed plant height remains unclear.

Currently, a single genetic resource cannot effectively improve the present plant architecture of rapeseed. GWAS is often used as an effective method to unravel the genetic architecture of complex agronomic traits in crops. Combined with association analysis and linkage analysis, 61 SNPs significantly associated with low zinc tolerance and 15 QTLs were identified in maize. Expression and haplotype analyses were used to mine the favorable allele conferring low zinc tolerance (Xu et al., 2022). Similar study could be found in Guo et al. (2021), in which 63 loci related to stem strength and yield were identified and favorable alleles for both high stem strength and high yield were discovered using 524 rice germplasm resources and 193 recombinant inbred lines (Guo et al., 2021). Based on GWAS and a transcriptome-wide association study, 15 stable QTLs and 1,854 candidate genes were detected in B. napus, which were significantly associated with seed glucosinolate content. Haplotype analysis showed that seed low glucosinolate was mainly resulted by the co-action of multiple favorable alleles (Tan et al., 2022). In this study, GWAS was performed on plant height of 230 B. napus accessions using three models (GLM, MLM, and BLINK). An unreported gene, BnaA01.BIN2, was simultaneously identified by all three models ( Figures 1 , 2 ), which increased the confidence of the results. However, no other loci or reported genes were co-identified, probably due to population structure constraints. These results provide insights for subsequent adjustment of population structure to more effectively detect available loci, genes, or favorable alleles.

GSK3 is a group of highly conserved cytoplasmic serine/threonine protein kinases that are widely present in animal and plant cells. These proteins perform their functions mainly by phosphorylating key substrate proteins of different signaling pathways. GSK3 is regulated by a variety of post-translational modification mechanisms. BnaA01.BIN2, identified by GWAS in this study ( Figures 1 , 2 ), encodes BRASSINOSTEROID-INSENSITIVE 2 (BIN2) and, belongs to the BnGSK3 family. AtBIN2 plays a role in the crosstalk between auxin and brassinosteroid signaling pathways (https://www.arabidopsis.org/index.jsp). In B. napus, BnaC04.BIL1, which has been isolated from the dwarf mutant Bndwarf2 (Yang et al., 2021), encodes BIN2-LIKE 1, and is also a member of GSK3s.

As a core member of GSK3s, BIN2 is a constitutively active kinase in plants, whose activity is affected by various regulatory mechanisms, including nucleocytoplasmic distribution, protein-protein interaction strength, phosphorylation and dephosphorylation, acetylation, and ubiquitination (Mao and Li, 2020). The direct function of BIN2 is to participate in the signal transduction pathway of brassinolide, which plays an important role in plant development (Anne et al., 2015). BIN2 directly controls the transcriptional regulatory complex composed of WEREWOLF (WER), transcription factor GLABRA3 (GL3), and WD40 repeat protein TRANSPARENT TESTA GLABRA1 (TTG1), It can phosphorylates GL3 and TTG1 in the WER-GL3-TTG1 complex to inhibit their transcriptional activity, thereby regulating root hair development (Cheng et al., 2014). BIN2 participates in photomorphogenesis by interacting with HY5, an important transcription factor for photomorphogenesis (Li et al., 2020b). In addition, BIN2 is involved in osmotic stress and adverse effects, and can promote lateral root development by phosphorylating auxin-responsive factor ARF7 (Cho et al., 2014). BIN2 is also involved in abscisic acid signal transduction to regulate the osmotic stress response (Wang et al., 2018) and enhances plant drought tolerance by phosphorylating RSPONSIVE OT DESICCATION 26, NAC family transcription factor (Jiang et al., 2019).

GSK3 is involved in the regulation of plant growth and development. However, only one gene has been reported to be related with plant height in rapeseed (Yang et al., 2021). In this study, to determine the relationship between BnGSK3s and plant height in allotetraploid rapeseed, we investigated the BnGSK3s family, which consists of 16 homologs in A sub-genome and 15 in C sub-genome ( Figure 3 and Table 3 ). Based on the transcriptome data of 13 tissues in ZS11, the expression pattern of BnGSK3s were found to show obvious expression preference difference in rapeseed ( Figure 6 ), in which 16 genes were highly expressed in SAM and three were highly expressed in the pistil ( Figure 6A ). Moreover, we identified favorable allelic variations in BnaA01.BIN2 among 230 B. napus accessions, whereas we failed to detect any SNP variation in the corresponding syntenic gene BnaC01.BIN2 (BnaC01g11150D). This could be caused by the limited numbers of accessions used for GWAS in this study. As such, more rapeseed genetic resources should be collected to dissect more favorable allelic variations in BnGSK3s for plant height and plant architecture.

A total of 230 rapeseed cultivars or inbred lines ( Supplementary Table 1 ) were collected worldwide, representing the genetic diversity of B. napus for GWAS of plant height. Field trials were conducted by a randomized design with three replications. For each accession, 45 individuals were grown in a 2.0 × 1.0 m² plot with three rows in each environment (2017–2018, 2018–2019, winter-spring growing season) in the Yangluo experimental field of Oil Crops Research Institute of the Chinese Academy of Agricultural Sciences, Wuhan, China. The R script lme4 (CRAN-Package lme4 (r-project.org)) and lsmeans were used to calculate the best linear unbiased prediction (BLUP) of each inbred line in the natural population (Zhao et al., 2022).

At the mature stage, 10 plants with good growth and development were randomly selected from each plot for phenotype investigation. The length from the cotyledon node to the apical position of the whole plant was measured and recorded as plant height. The statistics of the phenotype variation and frequency distribution were calculated using SPSS 22 (IBM SPSS, Armonk, NY, United States) (Zhao et al., 2022). The Pearson’s product-moment correlation analysis of plant height between 2017-2018 and 2018-2019 was carried out by R Package.

Total genomic DNA from fresh young leaf tissue of each inbred line (230 accessions) was extracted using a Hi-DNAsecure Plant Kit (TIANGEN, Beijing). DNA libraries were constructed with high-quality genomic DNA and whole-genome resequencing (WGS) was performed using the Illumina NovaSeq 6000 system. Clean data (clean reads) were obtained by filtering the raw data. All clean reads were mapped to the B. napus reference genome (Darmor-bzh V5, https://www.genoscope.cns.fr/brassicanapus/data/) using the Burrows-Wheeler Aligner software (Li and Durbin, 2009; Chalhoub et al., 2014). SAMTools (parameter: -q 30; http://samtools.sourceforge.net/) and Sentieon Genomics (parameter: –algo Dedup –rmdup) software were used to filter alignment duplications (Li et al., 2009; Freed et al., 2017). GATK (version 4.1.4.0) and vcftools (version 4.2) were used for SNP identification and filtration (parameters: MQ < 50.0 || QD < 2.0, -min-alleles 2 -max-alleles 2 -maf 0.05 -max-missing 0.9, and -cluster-size 3 -cluster-window size 10) (McKenna et al., 2010; Danecek et al., 2011). At last, a total of 1,707,194 informative SNPs were acquired, and the original SNPs were obtained from published data of our lab (Tang, 2019; Ding et al., 2020).

To analyze the natural population structure and linkage disequilibrium (LD) decay, ADMIXTURE (Version 1.3.0) (Alexander et al., 2009), Q+K model, and PopLDdecay (Zhang et al., 2018) were performed according to detailed descriptions from previous studies (Zhao et al., 2022). Three software and models were used, including the general linear model (GLM) in trait analysis by association, evolution, and linkage (TASSEL, Version 5.0) (http://www.maizegenetics.net/tassel); mixed linear model (MLM) in Efficient Mixed-Model Association eXpedited (EMMAX); and Bayesian information and Linkage-disequilibrium Iteratively Nested Keyway (BLINK) (Huang et al., 2019). TASSEL was used to calculate the kinship of 230 B. napus accessions (Yu et al., 2006). The LD block was displayed using LDBlockShow software (Dong et al., 2020).

Complete coding sequence of BnaA01.GSK3 was amplified from the B. napus cv. Zhongshuang11 (ZS11). The purified DNA fragment was fused with green fluorescent protein (GFP) in the backbone vector pBWA(V)HS-gfp, resulting in the plasmid 35S:BnaA01.BIN2-GFP via the ClonExpressMultiS One Step Cloning Kit C113-01 (Vazyme). The 35S:GFP plasmid was used as the mock control. These plasmids were transiently transformed into Arabidopsis protoplast cells using the Agrobacterium-mediated method. The subcellular localization of BnaA01.BIN2 was determined by observing GFP using a Nikon C2-ER confocal microscope (Nikon, Japan) (Zhao et al., 2021). The primers used for amplification of BnaA01.BIN2 are listed in Supplementary Table 4 .

The amino acid sequences of AtGSK3s family were obtained from the database “The Arabidopsis Information Resource (TAIR; https://www.arabidopsis.org/),” which were used to build a Hidden Markov Model, and HMMER3.0 was used to search the annotation and genome information of B. napus “Darmor-bzh” in the Brassicaceae Database (BRAD) (Mistry et al., 2013; Chalhoub et al., 2014; Chen et al., 2021a). The isoelectric point (pI) and molecular weight (MW) of BnGSK3s proteins were predicted using ProtParam online software (https://web.expasy.org/protparam/).

The National Center for Biotechnology Information (NCBI) Conserved Domain Database (https://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi) and the SMART database (http://smart.embl.de/) were performed to verify the candidate BnGSK3s (Lu et al., 2020; Letunic et al., 2021). Chromosomal locations of the candidate BnGSK3s were visualized via MapGene2Chromosome V2 (MG2C, http://mg2c.iask.in/mg2c_v2.0/).

The sequences of BnaGSK3s were downloaded from BRAD and the gene structures were displayed by Tbtools. The conserved motifs were analyzed by Multiple Expectation Maximization for Motif Elicitation (MEME, http://meme-suite.org) (Bailey et al., 2015; Chen et al., 2020).

The alignment of the amino acid sequences of AtGSK3s and BnGSK3s was performed by ClustalW (Larkin et al., 2007). Phylogenetic tree was constructed and visualized using the neighbor-joining (NJ) method in MEGA11 software with 1,000 bootstrap replications (Tamura et al., 2021), and visualized by Evolview7 software (He et al., 2016). The syntenic analysis of GSK3s between AtGSK3s and BnGSK3s were obtained from the BRAD database.

The RNA-seq data generated from 13 tissues of ZS11, including roots, SAM, stems, leaves, buds, siliques, stamens, pistils, blossomy petals, wilting petals, sepals, ovules, and pericarps, was previously published in our lab (Sequence Read Archive accession: PRJNA474576 in NCBI and CNP0001630 in China National GeneBank DataBase) were used for the expression pattern analysis of the BnGSK3s family (Li et al., 2019b; Dong et al., 2021). FastPure Plant Total RNA Isolation Kit RC401 (Vazyme) was used to extracted total RNA from three biological replicates of different tissues of ZS11, including the roots, SAM, stems, leaves, buds, and siliques using the. First-strand cDNA was generated using a HiScript III 1^st Strand cDNA Synthesis Kit (+gDNA wiper) R312 (Vazyme). Quantitative real-time PCR (qRT-PCR) was performed according to a previously described protocol, and the BnActin gene was used as an internal control to quantify the relative expression levels of target genes (Zhao et al., 2021). Gene-specific primers for BnGSK3s used for qRT-PCR were obtained from the qPrimerDB qPCR Primer Database (Lu et al., 2018) and the corresponding sequences were listed in Supplementary Table 4 . The heatmap is illustrated using OmicShare Tools (https://www.omicshare.com/tools/).