A diversity panel of 490 B. napus accessions collected from Xu et al. (2016) was used in this study. This panel was cultivated in eight natural environments, i.e., Wuhan 2013 and 2014 (WH2013 and WH2014), Changsha 2013 and 2014 (CS2013 and CS2014), Nanjing 2013 and 2014 (NJ2013 and NJ2014), Ezhou 2013 (EZ2013), and Chongqing 2013 (CQ2013). Additionally, the Gangan and ZS11 cultivars for RNA-seq were planted in Wuhan 2018 at the experimental stations of Huazhong Agricultural University. The design of field trial of the above materials and the acquisition of phenotypic data were the same as those used in the previous study (Xu et al., 2016). The re-sequencing genome data were obtained from Tang et al. (2021). The B. napus genome (B. napus ZS11 v0) from BnPIR (Song et al., 2020; Song et al., 2021) (http://cbi.hzau.edu.cn/bnapus/index.php) was used as the reference genome.

By using the “lme4” R package (Bates et al., 2015), the best linear unbiased prediction (BLUP) model was fitted to each B. napus accession:

Taking into account the variations between eight environments as phenotypic variance derived from environmental factors, broad-sense heritability ( h B 2 ) was estimated using the following equation by treating populations as a random effect and the environments as an environment effect, where σ g 2 and σ e 2 is the variance derived from genetic and environmental effects, respectively (Knapp et al., 1985).

Climatic data for daily average temperature (TAVG, °F) and precipitation (PRCP, in) were retrieved from the National Oceanic and Atmospheric Administration (https://www.noaa.gov/weather). Due to the lack of climatic data for Ezhou, there were climatic datasets of seven environments in total, i.e., WH2013 and WH2014 (114.05°E, 30.60°N; Station ID: GHCND: CHM00057494), CS2013 and CS2014 (112.87°E, 28.23°N; GHCND: CHM00057687), NJ2013 and NJ2014 (118.90°E, 31.93°N; GHCND: CHM00058238), and CQ2013 (106.48°E, 29.58°N; GHCND: CHM00057516). Climatic data were obtained from the day after being planted to the 200 days after planting (DAP). For each window from a starting day (3 DAP) to an end day (41 DAP) during B. napus growth, the average value of the climatic index and their correlation with the environmental mean vector for flowering time was calculated by CERIS analytical package (Li et al., 2021; https://github.com/jmyu/CERIS_JGRA). The most relevant climatic index for flowering time was chosen according to the highest correlation between environmental means and climatic index with corresponding window. Reaction norms were calculated as described in Guo et al. (2020) and Liu et al. (2020), using environmental mean and environmental climatic index as x-axis and phenotype as y-axis. Each line represented an individual and was shown by fitted linear regression. The intercept and slope were used to perform GWAS further.

The intersection of the accessions in phenotypic and genotypic datasets, i.e., 490 accessions with 11,700,689 SNPs, were used for GWAS using 3VmrMLM (Li et al., 2022a) via software IIIVmrMLM (Li et al., 2022b). Flowering time QTNs were obtained from separate analyses of phenotypic data from eight environments and joint environmental analyses of these datasets. The reaction norms between flowering time and climatic index were also used to conduct GWAS by 3VmrMLM. QEIs for flowering time were obtained by joint environment analyses of the above phenotypic datasets in eight environments. Population structure and kinship matrix were considered in 3VmrMLM analysis, and the “svpal” parameter was set as 0.01. According to Tang et al. (2021), the population structure calculated as K=3 was used in the analysis. The threshold was set at 0.05/m for significant QTNs and QEIs and LOD score ≥ 3.0 for suggested QTNs and QEIs, where m is the number of markers (Li et al., 2022a; Li et al., 2022b). According to the LD interval estimated by Tang et al. (2021), stable QTNs were defined as QTNs identified in at least three environments within the 100-kb upstream and downstream regions.

To identify candidate genes for flowering-time-related QTNs and QEIs, genes within the 100 kb upstream and downstream regions of each QTN or QEI were extracted according to the LD interval estimated by Tang et al. (2021). Then, two strategies were employed. First, the B. napus homologs of Arabidopsis flowering time genes downloaded from FLOR-ID (http://www.flor-id.org) were selected and considered as known genes. Second, new candidate genes were identified using differentially expressed genes (DEGs) in two SWORs (Gangan and ZS11) before and after vernalization and in two SORs (Westar and No. 2127). The t-test was adopted in the hypothesis testing for haplotype analysis; p< 0.05, p< 0.01, and p< 0.001 indicated the significances at 0.05, 0.01, and 0.001 probability levels, respectively.

The leaves of Westar, No. 2127, Gangan, ZS11 at 24 and 147 DAP were collected for RNA-seq with two biological replicates. Total RNA was extracted using the TIANGEN RNAprep Pure Plant Kit. Sequencing libraries were generated using the NEBNext^® UltraTM RNA Library Prep Kit for Illumina^® (NEB, USA) and were sequenced on an Illumina Hiseq 4000 platform. The detailed processes were described in Tan et al. (2022). We used MultiQC (Ewels et al., 2016) to perform quality control and Salmon (Patro et al., 2017) to quantify the RNA-seq reads of annotated genes in the reference ZS11. DESeq2 was used for differential expression analysis (Love et al., 2014). The threshold for DEG is set as the absolute value of log₂FoldChange >1 and adjusted p< 0.05 (two-tailed Student’s t-test; Tan et al., 2022).

To detect the regions under selective sweeps between SWOR and SOR, XP-CLR (v1.1.1), a genome scan using the composite likelihood approach was performed in sub-populations (Chen et al., 2010). Each chromosome was analyzed using the XP-CLR command with the parameters “–ld 0.99 –phased –maxsnps 200 –size 100000 –step 10000.” Non-overlapping 20-kp windows within the top 20% XP-CLR scores were merged into one single region, and then, these regions in the top 1% of XP-CLR scores were considered as candidate selective regions (An et al., 2019).

According to the above RNA-seq datasets, Pearson correlation analysis was calculated between candidate genes and DEGs in SWORs and SORs, respectively. Significant genes were considered to be co-expressed when Pearson correlation coefficient was >0.80 and p-value was<0.05. Network visualization was implemented with the Cytoscape package (Shannon et al., 2003).

Complex flowering time variation was observed in diversity group of 490 B. napus oilseed rapes, including 49 SORs, 20 winter oilseed rapes, 326 SWORs, and 95 mixed type oilseed rapes, grown in eight natural environments ( Figure 1A ; Supplementary Table S1 ). The means plus standard deviations of the eight environments WH2013, WH2014, CS2013, CS2014, NJ2013, NJ2014, CQ2013, EZ2013, and BLUP values were 155.49 ± 3.80, 153.56 ± 9.61, 160.27 ± 4.29, 166.55 ± 5.44, 160.50 ± 5.49, 167.55 ± 6.38, 151.31 ± 7.82, 162.57 ± 5.30, and 159.68 ± 4.83 (DAP), respectively ( Figure 1B ). The correlation of each pair of environments ranged from 0.37 to 0.72 (0.50 ± 0.09). The coefficients of variation, skewness, and kurtosis of the trait in eight environments illustrated that flowering time is a typical quantitative trait ( Supplementary Table S2 ). The broad-sense heritability for flowering time is 0.86. More importantly, joint regression analysis modeled with environmental mean showed the presence of a significant phenotypic plasticity ( Figure 1C ).

Climate change is altering the environment in which all plants grow. To understand the effect of climatic index on flowering time plasticity, the correlation between environmental means and climatic index (TAVG and PRCP) for different growth windows was predicted by CERIS ( Supplementary Table S3 ). The results of the correlation pattern between TAVG and flowering time showed a positive correlation at early seedling stage and a negative trend after bolting stage, while the pattern of PRCP was exactly opposite ( Figure 1D ; Supplementary Figure S1A ; Supplementary Table S4 ). The windows with the highest negative (TAVG_135–144 and PRCP_3–41) and positive correlations (TAVG_10–19 and PRCP_133–169) were chosen as the most related climatic index for further analysis ( Figures 1E , F ; Supplementary Figures S2A –F ; Supplementary Table S4 ). TAVG_135–144 (r = −0.986) showed higher correlation with flowering time than PRCP_3–41 (r = −0.809). TAVG_10–19 (r = 0.922) showed higher correlation with flowering time than PRCP_133–169 (r = 0.901). It is noted that these windows are surrounded by other windows with slightly decreasing correlation values ( Figure 1D ; Supplementary Figure S1 ).

To detect QTNs for flowering time, the phenotypes in each of the eight environments were used to associate with 11,700,689 SNPs using 3VmrMLM under population structure and polygenic background control. As a result, 55, 57, 42, 49, 54, 50, 44, and 43 significant QTNs at the critical p-value of 4.27e−09 (=0.05/m, where m is the number of markers) and 10, 5, 14, 10, 8, 13, 13, and 13 suggested QTNs (with the LOD score ≥ 3.0 but the p > 0.05/m) were identified for WH2013, WH2014, CS2013, CS2014, NJ2013, NJ2014, CQ2013, and EZ2013, respectively ( Supplementary Table S5 ; Supplementary Figure S3 ). In addition, flowering phenotypes from eight environments were used to perform joint analysis by 3VmrMLM. Sixty-eight significant and 11 suggested QTNs were identified. Based on the above QTNs in single and multiple environments analyses, 19 stable QTNs were identified in at least three environments ( Figure 2A ; Table 1 ).

All the datasets in eight environments were used to conduct joint analysis for identifying QEIs using 3VmrMLM. As a result, 32 significant QEIs and 4 suggested QEIs were identified, including 10 significant QEIs overlapped with the above stable QTNs ( Supplementary Table S6 ). Among these significant QEIs, 20 were found to have the highest absolute value of additive-by-environment interaction effects in WH2014 than those in other environments ( Figures 2B, D ), e.g., BnvaC0967693730 has an additive-by-environment interaction effect of −1.85 in WH2014 than those in other environments ( Supplementary Table S6 ; LOD = 67.17; R ² = 1.07%). The two loci BnvaC0967693730 and BnvaA0406097547 have the highest R ² (LOD = 67.17; R² = 1.07% and LOD = 66.42; R ² = 1.06%, respectively).

To obtain reaction norms of flowering-time-related climatic index, joint regression analyses were performed on phenotypes and the above flowering-time-related climatic indexes (TAVG_135–144, PRCP_3–41, TAVG_10–19, and PRCP_133–169; Supplementary Table S4 ). The intercept and slope of reaction-norm parameters were used to detect QTNs for flowering-time-related climatic indexes using 3VmrMLM. As a result, 10 QTNs for reaction norm parameters of B. napus flowering time were commonly identified with the above stable QTNs or QEIs, including 5, 2, 1, and 2 for TAVG_135–144, PRCP_3–41, TAVG_10–19, and PRCP_133–169, respectively ( Figure 2C ; Supplementary Table S7 ).

To mine candidate genes among the above QTNs and QEIs, DEGs analysis was conducted before and after vernalization. A total of 5,511 DEGs were identified in two SORs before and after vernalization ( Figure 3A ; Supplementary Table S8 ), and 5,850 DEGs were identified in two SWORs before and after vernalization ( Figure 3A ; Supplementary Table S9 ). Then, according to Arabidopsis gene annotation, 12 candidate genes were found to be associated with flowering time in approximately above 19 stable QTNs, including 7 known flowering-time-related homologs in Arabidopsis and 5 newly discovered genes ( Table 1 ). Using the same methods, 14 candidate genes were identified to be located in the above 32 QEIs, including 9 homologs of known genes, in which their homologs are related to flowering time and environments in Arabidopsis and 5 newly identified genes ( Table 2 ). In addition, five candidate genes were found to be associated with flowering-time-related climatic index, including two genes (BnaC02.DDB2 and BnaA05.COL9) commonly identified in QEIs and three genes (BnaA02.VIN3, BnaC02.RUG1, and BnaA06.UBC2) commonly found to be associated with the flowering time QTNs ( Supplementary Table S7 ).

Among these candidate genes, BnaFTs, BnaA05.COL9, BnaA08.SRC2, and BnaA03.DREB1B were significantly upregulated before vernalization in both SWORs and spring SORs, while BnaFLCs, BnaA01.FSD1, BnaC02.RUG1, BnaC05.PLGG, and BnaC03.PSAK were significantly upregulated after vernalization ( Figures 3A , B ). Interestingly, BnaA02.VIN3 and BnaC02.FPGS were only significantly upregulated before vernalization in SOR, which may indicate different functions between eco-types.

To validate the above flowering time candidate genes, we conducted selective sweep, haplotype, and co-expression analysis. First, by performing XP-CLR between SWORs and SORs, 954 selective sweeps were detected ( Supplementary Table S10 ). Eleven candidate genes for flowering time were found in the selective sweep, e.g., BnaFLCs, BnaFTs, BnaC02.FPGS1, BnaA08.SRC2, BnaA01.FSD1, BnaA02.VIN3, and BnaC09.PRR7. Second, haplotype analyses were further conducted in these genes. For BnaA02.FT, BnaA10.FLC, BnaA02.VIN3, and BnaC09.PRR7, significant difference exists between each haplotype in different environments ( Figures 4A–D ; Supplementary Figures S4A–D ). Interestingly, the haplotype for early flowering tends to exist in SORs, while the haplotype for late flowering prefers to exist in SWORs. Moreover, the co-expression networks of BnaA02.VIN3 and BnaC09.PRR7 have been constructed using DEGs in SWORs and SORs, respectively ( Figures 5A, B ). The co-expressed genes of BnaA02.VIN3 mainly participated in the circadian clock, photoperiodism, light perception, and signaling. Eight genes are specific co-expressed in SORs, including BnaA07.ZEP and BnaC09.ABCG22 in response to water deprivation. Five genes are specific co-expressed in SWORs. On the other hand, the co-expressed genes of BnaC09.PRR7 mainly participated in the circadian clock and autonomous pathway. Five and one genes are specific co-expressed in SORs and SWORs, respectively. BnaCKA2s and BnaPKDM7s participated in epigenetic regulation.

It is well-known that the flowering time regulation of B. napus is in response to day length or vernalization (Reeves and Coupland, 2000). This study calculated the correlations between two climatic factors, TAVG and PRCP, and flowering time in seven environments. TAVG correlated positively with flowering time in vernalization and negatively with flowering time after the seedling stage ( Figure 1D ). In a previous study, a reduction in autumn or winter chilling delays floral transition in B. napus (O’Neill et al., 2019). An elevated growth temperature is equally efficient in inducing the flowering of Arabidopsis (Balasubramanian et al., 2006). However, the transition or critical point of these two stages is unclear. For PRCP, this study reported the relationships between PRCP and flowering time in B. napus for the first time. Although the correlation coefficients are lower than TAVG, PRCP was found to be correlated negatively with flowering time in early development and positively later ( Supplementary Figure S1 ). This result is consistent with a previous study in Arabidopsis that flowering time correlated negatively with fall and winter precipitations and positively with summer precipitation (Vidigal et al., 2016).

In this study, multi-environment joint GWAS improved the power on identifying more QTNs than single environment GWAS. We dissected the genetic basis for flowering time in the following three aspects. First, 12 flowering time candidate genes were mined in approximately 19 stable QTNs for flowering time. Seven genes are previously reported, e.g., BnaFLCs (Tadege et al., 2001), BnaFTs (Wang et al., 2009), BnaAGL24 (Yu et al., 2002), and BnaCRY2 (Sharma et al., 2022), whereas five genes are newly identified, which are differentially expressed before and after the vernalization of different ecotypes ( Figure 3 ; Supplementary Tables S8 , S9 ). Second, it is worth noting that this study focused on the mining of flowering time genes related to the environments. Fourteen candidate genes were identified around 32 QEIs, including 9 known flowering time genes related to environments. For example, BnaCOL9 and BnaUBP12 are regulated by the circadian clock in the photoperiod pathway (Cheng and Wang, 2005; Cui et al., 2013). BnaFLD is subjected to the direct regulation by brassinosteroids (Zhang et al., 2013b). It has been reported that the overexpression of BnaDREB1B not only delayed flowering but also responded to cold (Seo et al., 2009). BnaEBS functions in the chromatin-mediated repression of floral initiation by H3K4me3 (López-González et al., 2014). Finally, five genes were found to be associated with flowering-time-related climatic index. BnaC02.DDB2 and BnaA05.COL9 were commonly identified in QEIs, and BnaA02.VIN3, BnaC02.RUG1, and BnaA06.UBC2 were commonly found to be associated with the main effect flowering time QTNs.

In this study, the missing heritability exists, in which the total phenotypic variance explained of QEIs and QTNs is much less than the estimated broad-sense heritability. This can be explained in several ways. First, the population is not enough to detect rare variants. Second, allelic heterogeneity may be the reason for this phenomenon. Lastly, epigenetic variation is likely to be a source of missing heritability (Brachi et al., 2011). Moreover, some candidate genes for stable QTNs, e.g., BnaA02.FT and BnaA10.FLC, were found to be related to environments but were not identified in QEIs ( Figure 4 ). This result is explained by multiple facets, e.g., the difference in phenotypic data among environments, the diversity of population accessions, and the power of QEI detection. In the previous study, COL9 and FLD have been reported to regulate FT and FLC, respectively (Cheng and Wang, 2005; Jiang et al., 2009). BnA05.COL9 and BnaA03.FLD were found to be candidate genes for QEIs in this study. We hypothesized that QEI may be associated with direct environmental response upstream regulators due to the complexity of transcription and epigenetic regulations of flowering (Bouché et al., 2016).

In Arabidopsis, VIN3 acts together with PRC2 to repress histone marks at FLC in response to vernalization (Kim and Sung, 2013). PRR7 was reported to coordinate with PRR9 and PRR5 and regulate flowering time through the canonical CO-dependent photoperiodic pathway (Nakamichi et al., 2007). In this study, BnaA02.VIN3 and BnaC09.PRR7 have been shown to be crucial G by E genes for flowering time. There are three pieces of evidence. First, BnaA02.VIN3 is significantly associated with ChrA02-6152101 (LOD = 13.14) for flowering-time-related climatic factors and with ChrA02-6374324 (LOD = 12.17) for flowering time in WH2013. BnaC09.PRR7 is significantly associated with the QEI, ChrC09-67693730 (LOD = 67.17), by multi-environment GWAS. Second, BnaA02.VIN3 and BnaC09.PRR7 are DEGs before and after vernalization and in the selective sweep between SORs and SWORs ( Figure 2 ). Then, in these genes with significant haplotype differences, their haplotypes for early flowering tend to exist more in SORs ( Figures 4C , D ). Lastly, co-expression networks were constructed for BnaA02.VIN3 and BnaC09.PRR7. Some relationships have been proven, e.g., PRR7 with LHY (Liu et al., 2013), PRR7 with PRR5 (Nakamichi et al., 2007), and VIN3 with CCA1 and LHY (Kyung et al., 2022).

In summary, we dissected the G by E for flowering time for B. napus from different eco-types in eight environments. Four windows of flowering-time-related climatic index were identified. Stable QTNs and QEIs for flowering time and their candidate genes were identified. These findings provide valuable information that can be used to breed B. napus varieties with optimized flowering time by pyramiding favorable alleles. The candidate genes will also greatly promote the dissection of flowering time mechanisms in different eco-types.