The B. napus population used in this study comprises 2,311 B. napus accessions, including 1,259 from Asia, 929 from Europe, 60 from North America, two from South America, 38 from Oceania and four from Africa (Table 1) (Lu et al., 2019; Song et al., 2020; Tang et al., 2021; Wu et al., 2019). Three ecotypes, including spring (354 accessions), winter (756 accessions), and semi-winter (1,122 accessions) were included in the population.

The re-sequencing reads with low-quality were firstly filtered out by Trimmomatic (v0.36) (Bolger et al., 2014). Next, the remaining reads were mapped to ZS11 and Darmor-bzh v4 (Chalhoub et al., 2014) reference genome using BWA-MEM (v0.7.17) (Li and Durbin 2009) with default parameters, respectively. Reads with MAPQ lower than 10 were removed by SAMtools (v1.9) (Li et al., 2009), and reads of PCR duplication were removed by Picard (v2.17.6, http://broadinstitute.github.io/picard/). Then, SNP and InDel calling was performed with Sentieon (v201704.01). To obtain high-quality variations, raw SNPs and InDels were filtered using the VariantFiltration module of GATK (v 3.6-0) (McKenna et al., 2010) with the parameter "QUAL < 30.0 || MQ < 50.0 || QD < 2". The SNPs and InDels with a minor allele frequency (MAF) lower than 0.01 and a missing ratio higher than 30% were filtered out by VCFtools (v0.1.16) (Danecek et al., 2011), and the SNPs and InDels with heterozygosity lower than 0.5 in the population were also removed. After that, genotype imputation for the SNP set was performed using beagle (v5.1) (Browning and Browning 2007). SV identification was also conducted in the population. Finally, 8,709,516 SNPs, 1,287,483 InDels, and 93,562 SVs were obtained based on ZS11 reference. The SVs included 61,002 insertions, 31,552 deletions, 647 inversions and 361 duplications. A total of 5,781,762 SNPs and 906,714 InDels were obtained based on Darmor-bzh reference.

Gene haplotype analysis was conducted by mainly referring to the method of Yano et al (Yano et al., 2016). Firstly, gene haplotype identification based on polymorphisms localized on the gene regions was performed by PHP codes written in house. Then, frequencies of haplotypes were estimated using PLINK (v1.90b4.4) (Purcell et al., 2007) with the parameter "--hap-freq". Haplotypes with a frequency lower than 0.01 were filtered.

SNPs with strong linkage disequilibrium (LD) were discarded with PLINK (v1.90b4.4) (Purcell et al., 2007) with parameter "--indep-pairwise 100 50 0.8" for phylogenetic tree construction using TreeBeST (v1.9.2) (Vilella et al., 2009) with the "treebest nj -b 1000" parameter. The population structure analysis of all accessions was performed using fastStructure (v1.0) (Raj et al., 2014) with the number of ancestry kinships (K) from 2 to 10, and K = 3 was used. Principal component analysis (PCA) was performed using GCTA (v1.92.4 beta2) (Yang et al., 2011). Linkage disequilibrium (LD) between pairs of SNPs on each chromosome was estimated as the correlation coefficient (r²) in PLINK (v1.90b4.4) (Purcell et al., 2007) and visualized using LDheatmap R package (Shin et al., 2006).

To explore the candidate regions potentially effected by selection, nucleotide diversity (π) and population fixation statistics (). π was calculated with expected heterozygosity per site derived from the average number of sequence differences in three subpopulations (winter, spring and semi-winter) with VCFtools (v0.1.13) (Danecek et al., 2011) in a 500-kb sliding window with a step size of 50 kb. , which indicates genomic differentiation between populations, was calculated for each pair of subpopulations using VCFtools (v0.1.13) (Danecek et al., 2011) in a same windows size as that of π. The outliers of were determined as windows with the top 5% of values, indicating divergent selection of the regions. To explore the selection pressure of genes in B. napus genome during domestication, K_a, K_s and K_a/K_s ratio were calculated for orthologous gene pairs between B. napus and the ancestors, B. rapa and B. oleracea, with K_a/K_s_Calculator (v 2.0) (Wang et al., 2010). Selection signals were calculated based on ZS11 and Darmor-bzh reference genome, respectively.

RNA-seq data of seeds from 309 B. napus accessions at 20 and 40 days after flowering (DAF) were obtained from a previous study (Tang et al., 2021). For RNA-seq analysis, adaptors and low-quality read pairs with read quality < 80 were removed using Trimmomatic (v0.38) (Bolger et al., 2014). Then, the filtered reads were mapped to the ZS11 and Darmor-bzh v4 (Chalhoub et al., 2014) reference genome using Hisat2 (v2-2.1.0) (Kim et al., 2015) with default parameters, respectively. Finally, FPKM of genes were calculated with StringTie (v1.3.3b) (Pertea et al., 2015). Moreover, RNA-seq data of multiple tissues, including root, flower, bud, leaf, silique, silique wall, were retrieved from BnTIR (Liu et al., 2021). For each tissue, samples were merged and then about 25 Gb of reads were randomly extracted for visualization in JBrowse.

Eighteen traits of 1,703 B. napus accessions from multiple planting areas and years were collected from published studies (Liu et al., 2016; Song et al., 2020; Tang et al., 2021; Wu et al., 2019; Zhao et al., 2019) (Table 2). The best linear unbiased prediction (BLUP) values of all planting areas and years were used as final phenotypic values for further analysis.

The BLUP results of Eighteen traits and the SNPs and InDels with a minor allele frequency (MAF) higher than 0.05 were used for genome-wide association study (GWAS). GWAS was performed for the Eighteen traits with GEMMA (Zhou and Stephens 2012) (v0.98.1). The population structure was controlled by the first two principal components as covariates, as well as an IBS kinship matrix derived from SNPs and InDels calculated by GEMMA (Zhou and Stephens 2012) (v0.98.1). The cutoff for determining significant associations was set as -lg(1/n), where n represents the total number of variations.

The BnVIR was built with the software of MySQL (https://www.mysql.com), PHP 7.4 (http://www.php.net), ThinkPHP 5.0 (http://www.thinkphp.cn/) and Apache 2.4 (http://www.apache.org), and all procedures were run on a Linux Ubuntu 7.5 (https://ubuntu.com/) operation system. To construct a user-friendly web interface, JQuery (http://jquery.com) and Bootstrap (https://getbootstrap.com) were used in all database pages. A variety of graphs, such as geomap, polygenetic tree, gene structure, heatmap and violin plot, were drawn by multiple plugins (Table 3). We used the Jbrowser (http://www.jbrowse.org) (Buels et al., 2016) as a genome browser for the visualization of gene model, expression and variation in B. napus genome.

Gene index of B. napus genomes was obtained from BnTIR (http://yanglab.hzau.edu.cn/BnTIR) (Liu et al., 2021) for convenience in variation search. Reference genome, gene annotation, TE annotation of ZS11 (Song et al., 2020) and Darmor-bzh (Chalhoub et al., 2014) were used for visualization of genetic variations genome in JBrowse.

Differences in phenotypes and gene expression of accessions with different genotypes were determined by Wilcoxon test and Student's t test in Variation module.