Genome-Wide Association Analyses Reveal Candidate Genes Controlling Harvest Index and Related Agronomic Traits in Brassica napus L.

Abstract

Harvest index (HI) is a complex and vital agronomic trait that is closely related to the economic benefits of rapeseed. In this study, we measured the HI and 13 HI-related agronomic traits of 104 core breeding lines of rapeseed during 3 years and sequenced the populations using the Bnapus50K array. The phenotypic analyses showed the complex connections among HI and other traits. A total of 212 significant SNPs related to the traits and 22 stable SNPs were identified. Four SNPs, A01_1783685 (PH and SYP), C06_26638717 (PH and NSS), C03_4731660 (MIL and MINS), and C09_36899682 (PH and BYP), were identified as potential pleiotropic loci. Compared to previous reports, 49 consensus loci were obtained that were related to PH, TSW, NSP, BAI, NSS, SL, BN, MINS, SYP, and BYP. Twelve stable SNPs were detected as promising novel loci related to BN (A05_19368584 and A05_19764389), SL (A06_23598999, A06_23608274, and C07_38735522), PH (C04_47349279, C04_47585236, and C09_36899680), MINS (C05_6251826), NSS (C06_22559430 and C06_22570315), and HI (C05_6554451). In addition, 39 putative genes were identified in the candidate intervals. This study provides novel insights into the genetic mechanisms of HI and HI-related traits, and lays a foundation for molecular marker development and casual gene cloning to improve the harvest index of rapeseed.

1. Introduction

Harvest index (HI), i.e., co-efficient of economics, is the ratio of economic yield to the total biomass production above the ground. For grain-producing crops, the “sink” capacity of the grains is the basis for the transfer and storage of assimilation products. The efficiency of the transformation and distribution of the assimilation product “source” to the seed represents the unobstructed “flow” [1,2]. Harvest index is always used for evaluating the “flow” in the “source-sink-translocation” theory.

Rapeseed (Brassica napus L., AACC, 2n = 38) is one of the most important oilseed crops worldwide. With the improvement of varieties and cultivation techniques, the crop harvest index continues to increase. Despite this, the HI of rapeseed is generally between 10% and 30%, which is much lower than that of tobacco (over 60%), rice (50–60%), peanut (50–60%), maize (40–50%), and wheat (40–50%) [3,4,5]. The increase in the output of rice, wheat, and barley mainly benefits from the increased harvest index after the Green Revolution [6], so the yield of rapeseed can be enhanced by increasing the HI as well. Although the complexity of HI can be affected by plant structures, various environmental factors, and nutrient distribution [3,7], we can increase HI by adjusting related factors.

In recent years, an increasing number of reports of B. napus have been published on the major agronomic traits, such as plant height (PH), siliqua length (SL), thousand-seed weight (STW), branch angle (BA), and so on [8,9,10]. Nevertheless, the studies about the harvest index or the relationship with the agronomic traits were denumerable. Nine SNPs were identified to be significantly associated with HI, explaining 3.42% of the phenotypic variance [11]. Lu et al. identified 294 single nucleotide polymorphisms (SNPs) significantly associated with HI and four HI-related traits by GWAS using 520 B. napus accessions, and 33 functional candidate genes were predicted in the intervals [12]. A segregation DH population containing 348 lines from the cross between “KenC-8” and “N53-2” was used to detect quantitative trait loci (QTL) for HI, seed yield, biomass yield, and PH, and 160 QTL, 163 epistatic loci pairs for the studied traits [2]. Luo et al. used a segregating population of B. napus, identified 1904 consensus QTL accounting for 22 traits, and formed a network containing developmental traits, seed quality traits, seed yield, and seed-yield components [13]. These results not only show us the feasibility of GWAS to locate quantitative trait loci (QTL) for complex agronomic traits, but improve our knowledge of interactions among complex traits and their related traits.

In this study, we focused on 104 core breeding inbred lines of B. napus, sequenced by the Bnapus50K array, and cultivated them across three cropping seasons. HI and 13 HI-related traits were measured for GWAS to associate significant SNPs. A trait-SNP network was constructed, and promising candidate genes were forecast. Our study will provide an important basement for the genetic mechanism improving the harvest index of rapeseed.

2. Materials and Methods

2.1. Material Planting and Field Management

The 104 rapeseed germplasms are the core breeding inbred lines from the rapeseed germplasm resource lab of Northwest A&F University, and each line has at least one excellent agronomic trait. All of these materials were planted and measured in the Caoxinzuang experimental base (34°30′ N, 108°09′ E) of Northwest A&F University in Yangling during the 2017–2018 cropping season, 2018–2019 cropping season, and 2020–2021 cropping season. A completely randomized design was adopted with row width of 2 m, row spacing of 35 cm, and plant spacing of 10 cm, respectively. Each material was planted in two rows. Fresh and tender leaves were sampled at the budding stage and kept at −80 °C until further use.

2.2. Phenotype Collection

Five representative plants were selected randomly from each inbred line and labeled in order separately. The plant traits of each individual were investigated. We directly measured the plant height (PH, cm), main inflorescence length (MIL, cm), the numbers of primary valid branch (BN), silique length (SL, cm), silique width (SW, cm), number of effective siliques on main inflorescence (MINS), and number of seeds per silique (NSS). The branch angles (BA) were photographed and the values were acquired by AutoCAD software, including top branch angle (TBA, °), middle branch angle (MBA, °), and basal branch angle (BBA, °). The individual above-ground portions were harvested, placed into net bags, and dried under the sunlight. The shoot dry weight (biomass yield per plant, BYP, g) and the weight of the seeds (seed yield per plant, SYP, g) were weighed. Thousand-seed weight (TSW, g) was measured by a Wansheng Automatic seed test analyzer. Then, the branch angle index (BAI) was calculated by TBA, MBA, and BBA as follows:

$$\mathrm{BAI}=\sum \left(xi-BA{i}_{\min }\right)/(BA{i}_{\max }-BA{i}_{\min })/3$$

(1)

where x is the measured value, and i represents the TBA, MBA, and BBA, respectively. Silique density of the main inflorescence (MISD) was calculated by MIL/MINS, the number of siliques per plant (NSP) was calculated by SYP/(NSS × TSW), and harvest index (HI, %) was calculated by SYP/BYP.

2.3. DNA Isolation, Sequencing, and Quality Controlling

The genomic DNA was isolated from the leaves using the Tiangen DNAsecure Plant kit (Cat. #DP320-03) following the manufacturer’s protocol. The quality of DNA was checked on 1% agarose gel and the DNA concentrations were diluted to the range of 50–100 ng/uL. DNA sequencing was performed according to the Infinium assay standard protocol (Infinium HD Assay Ultra Protocol Guide, http://www.illumina.com/, accessed on 7 April 2020) using the Bnapus50K array [14] at the Greenfafa Institute (Wuhan, China). The SNPs were eliminated according to the rules with an absence rate of more than 65% and the smallest allele frequency less than 1%. The missing genotypes were imputed in Beagle v.21Apr21.304 with default parameter settings [15].

2.4. Population Structure Analysis

Approximately 26 K SNPs from these 104 B. napus lines were used for phylogenetic tree construction. The phylogenetic tree was constructed with Tassel 5.0 using a neighbor-joining algorithm [16] and displayed by the iTOL web tool (https://itol.embl.de/, accessed on 28 August 2021). Linkage disequilibrium between any two SNPs on one chromosome was estimated with r² using PLINK [17]. The parameters with a sliding window of 50 kb (in steps of 10 SNPs) were used, and SNPs related to other SNPs in the window with r² > 0.2 were removed. The population structure was determined by ADMIXTURE v.1.3.0 [18] and STRUCTURE 2.3.4 [19]. For ADMIXTURE five runs were performed with default parameters for each number of populations (K) set from 1 to 10, and the most likely K value was determined by the log probability of the data. For STRUCTURE, five independent runs were performed under the admixture model with K-value from 1 to 10, burn-in period of 10,000 iterations, and 100,000 Markov chain Monte Carlo. The principal component score of the samples was also computed by the principal component analysis (PCA).

2.5. Genome-Wide Associated Loci

All the analyses were performed using TASSEL5 [16] and the mrMLM v4.0.2 [20]. A mix linear model (MLM) was used by TASSEL5 with the covariate matrix from PCA and the kinship matrix. The p-value = 10 × 10⁻⁴ was used to establish a significance threshold for the result of TASSEL. For the mrMLM, five algorithms were applied to associate significant SNPs, including FASTmrMLM, ISIS EM-BLASSO, mrMLM, FASTmrEMMA, and pLARmEB [20]. The analysis process was referred to the protocol provided by the mrMLM and the PCA matrix was used as the covariance for GWAS by the default parameters. The previously obtained genotypic data of these 104 plants, along with their BLUP values, were used for genome-wide association. The significant SNPs detected by any two of the algorithms or environments above were regarded as reliable and stable loci for the trait.

2.6. Network Construction

The traits, environments, and associated SNPs were retained to build the phenotype-SNP network. The same traits among different environments were linked by Pearson’s correlation coefficients, and the different traits were connected by partial correlation coefficients. The width of edges connecting the SNPs to the traits represented the SNP phenotypic variation. We also linked the physical proximity SNPs. Network construction and visualization were performed using Cytoscape v.3.5.1 [21].

2.7. Candidate Gene Identification

Linkage disequilibrium (LD) analysis and visualization were performed by the software LDBlockShow with the default parameters [22]. All genes within the same LD block (r² > 0.6) of the significant SNPs were assessed as potential candidates. In addition, candidate genes outside the LD blocks but within 100 KB flanking were also predicted as possible candidate genes [23,24,25]. Gene functional annotations were obtained by protein BLAST in the Kyoto Encyclopedia of Genes and Genomes (KEGG, https://www.genome.jp/kegg/, accessed on 4 December 2021), TAIR (https://www.arabidopsis.org/, accessed on 4 December 2021), and Brassicaceae Database (BRAD, http://39.100.233.196/, accessed on 4 December 2021). Then, gene expressions during five developmental stages, including bolting, initial flowering, full-bloom, podding, and maturation, were extracted from the BrassicaEDB v.1.0 (https://brassica.biodb.org/, accessed on 7 December 2021) [26].

2.8. Data Analyses and Visualization

The phenotypic data were analyzed by Excel, and visualized by ggpubr in R, TBtools [27] and the online webtool “Hiplot” (https://hiplot.com.cn, accessed on 9 January 2022). The partial correlation coefficient is a coefficient that describes the relationship between X and Y when removing the effects of the other control variables under multi-factor conditions [28]. The correlation and partial correlation analyses were performed in the R environment (R v.4.0.2, http://www.r-project.org/, accessed on 9 January 2022) by the R packages PerformanceAnalytics v.2.0.4 (https://github.com/braverock/PerformanceAnalytics, accessed on 2 September 2021) and pltPcorrelation v.0.1.0 (https://gitee.com/wqssf102/pltpcorrelation, accessed on 2 September 2021), respectively. The heritability (h²) was computed as follows [29]:

$$h^2={\sigma }_G^2/({\sigma }_{GE}^2/n+{\sigma }_G^2+{\sigma }_e^2/nr)$$

(2)

where ${\sigma }_G^2$ is the genotypic variance, ${\sigma }_{GE}^2$ is the genotype × environment variance, ${\sigma }_e^2$ is the error variance, n is the number of locations, and r is the number of replications. The estimates of ${\sigma }_{GE}^2$, ${\sigma }_G^2$, and ${\sigma }_e^2$ were analyzed by the analysis of variance using the lmer function of the R/lme4. The best linear unbiased prediction (BLUP) was obtained by fitting the mixed linear model in the R/lem4 for the inbred lines using Y~1 + (1|line) + (1|year), where Y is the trait data, the parenthesis represents the random effects, “1|” represents groups, “line” represents all experimental lines used, and “year” represents the 3 years.

The locations of SNP/QTL on chromosomes were displayed via MapChart, and the Manhattan plots and quantile-quantile (QQ) plots were drawn by CMplot in R. The significance between the two groups of data was determined by the Wilcoxon test.

3. Results

3.1. Phenotypic Variation for Different Agronomy Traits

The phenotypic data of the agronomy traits were recorded for 3 years. The statistical parameters (average, maximum, minimum, standard deviation, coefficient of variation (CV), and h²) related to the traits are summarized in Table S1. The violin plots display the distribution patterns of phenotypic values (Figure 1). Variations of all the traits were widely and continuously distributed, and the correlations for the trait among different years were middle to high (Figure 1). The PH (CV = 9.32%) and SW (CV = 8.4%) were the two least CVs, while the highest average variation coefficients were SYP (CV = 31.71%), NSP (CV = 29.22%), and BYP (CV = 26.61%). The HI (CV = 16.36%) ranged from 0.18 to 0.49 in 2019 and 0.08 to 0.39 in 2021. Moreover, these traits (SYP, BYP, NSP, and HI) showed weak or no correlation among different years (Figure 1) and the heritabilities were also the lowest with the h² from 0.43 to 0.57 (Table S1). It suggested that these four traits were easily affected by environmental conditions and difficult to measure accurately.

We measured the three-part branch angles (BA) for each plant including TBA, MBA, and BBA. However, there were large differences between the angles of different parts, and the correlations among them were not strong (Figure S1). Thus, we calculated BAI to integrate these angles for further studies. For the rest of the traits, the statistical parameters were displayed in Table S1 as well. The correlations for all the traits among different years were middle to high (Figure 1). The heritabilities (h²) of them were highly repeatable, ranging from 0.43 to 0.90. The most stable inherited traits were SL and NSS with a heritability of 0.90 and 0.88, respectively. Several significant partial correlations were observed between all traits (Figure S2). Notably, HI only showed significant correlations with SYP (0.46 ***) and BYP (−0.58 ***), suggesting it was determined by SYP and BYP. For SYP, it was extremely significantly related to PH (0.22 **), BN (0.29 ***), NSS (0.46 ***), NSP (0.35 ***), and TSW (0.29 ***). For BYP, extremely significantly correlations were determined with MIL (−0.27 ***), MINS (0.24 ***), and MISD (0.23 ***). In addition, PH had a positive contribution to SYP (0.22 *) and BYP (0.15 *), and MISD showed the opposite contributions to SYP (−0.18 *) and BYP (0.23 **). Thus, although HI was a complex trait that was difficult to measure accurately and repeatedly, we can improve it indirectly by regulating its key factors.

3.2. Population Structure and Genome-Wide Association Analysis

All of the pairwise genetic distances among the 104 rapeseed lines were determined from the SNP genotypes, and 3476 independent SNPs were used for calculating population structure. Population structure was assessed for K values ranging from 1 to 10 with the lowest cross-validation error using the ADMIXTURE software. As a result, the LD decay changed continuously, and when K = 5, the cross-validation error was the lowest (Figure S3). For STRUCTURE, variation of Delta K suggested that the population could be assigned to seven groups (Figure S3C). Furthermore, a neighbor-joining tree determined from the SNP genotypes showed that the 104 rapeseed lines could be classified into five divergent groups (Figure S3D). Therefore, PCA = 5 was applied as the covariance for GWAS.

To ensure the reliability of GWAS results, genome-wide association mapping was executed with six algorithms by Tassel5 [16] and mrMLM [20]. Most of the ideal models were presented in accordance with the observed −log10(p-value) vs. the expected −log10(p-value) in the QQ plots (Figures S4 and S5). For different algorithms (i.e., SNPs detected by the same trait in different environments), there were differences in the detections of QTL of the six methods (Figure 2A), FASTmrEMMA (70), FASTmrMLM (318), ISIS EM-BLASSO (236), mrMLM (188), pLARmEB (300), and MLM (163). Among them, 183 SNPs detected by at least two algorithms were considered as credible loci. For different environments (i.e., SNPs detected by the same trait with different algorithms), a total of 279, 260, 305, and 431 loci were detected in 2018, 2019, 2021, and BLUP, respectively (Figure S2B). One hundred and seventeen significant SNPs (the same SNPs for different traits were treated as the repeated SNPs) associated with a certain phenotype were identified by any two environments of 2018, 2019, 2021, and BLUP, which were considered as credible loci as well. The SNPs C06_22559430 and A03_22103527 were simultaneously detected by the four environments of NSS and TSW, respectively (Table S2). Nine significant SNPs were identified across three different environments related to PH (4/9), SL (3/9), NSS (1/9), and MINS (1/9). These SNPs repeatedly detected in multiple environments can be regarded as reliable loci for controlling traits.

We combined the above credible SNPs detected by different algorithms and environments and excluded the minor QTLs with PVE of less than 1%. Three hundred and twenty-four QTLs involved in 212 SNPs above were used for further analysis. The traits mostly detected at QTLs were BAI (39/324), MISD (37/324), NSS (35/324), SL (34/324), TSW (34/324), PH (33/324), and MIL (26/324) (Table S3). These QTL were widely distributed on 19 chromosomes. The most distributed chromosomes were A01, A03, A05, C05, C06, and C09. Chromosome A08 contained only one SNP associated with MISD at A08_16931989. Five SNPs for the harvest index were detected on A05, C05, C08, and C09, three for BYP on A06, C02, and C05, and fourteen for SYP, mainly on A01, C04, A03_random, and A06_random.

3.3. Comparison of the Stable SNPs

Appropriately reducing plant height can not only improve lodging tolerance but also increase harvest index. For PH, four SNPs, namely A02_random_355708, A03_25488780, C09_36899682, and C04_47349279, were co-detected by three environments with PVE 5.2–35.6%. These were stable SNPs for controlling plant PH. Since multi-environmental phenotypic data were used in GWAS, we examined the association of the SNPs with BLUP values. As expected, all of them showed extremely significant associations with PH (Figure 3A). Interestingly, the majority of them showed a significant correlation with SYP, such as A01_1783685 and A03_25488780, while the A allele of SNP C09_36899682 was also associated with low BYP. Moreover, C09_36899682 explained for over 27.4% of the phenotypic variation of PH, and others only explained for 5.5–18.1% (Table S2). C09_36899682 could be an effective locus to control PH to decrease BYP without reducing SYP.

Silique and seed characteristics are important factors of seed yield. For NSS (Figure 3B), C06_22559430 and A07_random_1623623 were co-detected by at least three environments explaining the phenotypic variation from 4.5 to 21.2%. The T allele of C06_22559430 was associated with high NSS, but the A allele of A07_random_1623623 did not show significant association with NSS distribution. For TSW, A03_22103527 was detected by TSW across four environments, and the G allele showed higher TSW than the A allele (Figure 3C).

For pleiotropic SNPs, A01_1783685 was identified by PH_2018 and SYP_2018, and C06_26638717 was detected by PH_2019, PH_2021, NSS_2021, and NSS_BLUP (Figure 4 and Table S2). Notably, A01_1783685 and C06_26638717 were located in the interval related to PH as shown in previous studies [13,30,31] (Figure 4). However, they have not been detected by SYP or NSS. Meanwhile, C03_4731660, co-detected by MIL_2019 and MINS_2019 (Figure 5 and Table S2), has not been reported to be associated with MIL or MINS. These results implied that these loci might play multiple roles in plant growth.

For harvest index, C05_6554451 was significantly associated with HI_2019 and HI_BLUP, explaining for 13.4% and 4.6% phenotypic variation, respectively. C05_6554451 was detected to be significantly associated with HI_BLUP (p = 0.0044), and the A allele was high harvest index (Figure 3D). In addition, we also analyzed the significance between other agronomic traits and the allele of stable detection loci of traits (Figure 3). These results provided favorable alleles to regulate HI directly and indirectly.

3.4. Network of Significant SNPs Associated with Phenotypes

Based on the partial correlated relationships among the traits and the significant SNPs for the traits (Figure 5), we found that the 14 traits in this study tended to be correlated relationships, suggesting that they might be genetically coregulated. To dissect the correlations across different traits clearly, a phenotype-SNP network was constructed. At the trait level, all the traits were linked directly or indirectly by partial correlated relationships. For trait to SNP, they were associated with the PVE of significant SNP. At the SNP level, SNPs related to more than one trait or environment or adjacent to other significant SNPs were marked. Based on the networks, it was easily to find that some SNPs were associated with the certain trait across multiple years, and some with multiple phenotypes, such as A03_22153527, C06_26638717, and A06_23598999. Meanwhile, the adjacent SNPs were also displayed, and these genomic intervals may control multiple phenotypes. Overall, the network provided a visualized map for hub stable or pleiotropic SNPs, and suggested that complex relationships existed in HI and the other traits.

3.5. Identification of Candidate Genes for Stable Loci

According to the phenotype-SNP network of the GWAS results (Figure 5), the stable significant SNPs and pleiotropic SNPs were determined as stable loci. In addition, significant SNPs for HI and significant SNPs adjacent to the stable loci were also listed as stable loci. Twenty-two SNPs associated with the agronomic traits were filtered as promising QTLs. Potential candidate genes were identified in the intervals of significant marker-trait associations by investigating all genes within shared LD blocks or with physical proximity within 100 KB upstream or downstream of the SNPs (Figure S6). The genomic regions from 100 KB of both flanks of the significant SNPs were inspected for putative candidate genes for each agronomic trait. The candidate genes of different QTLs were selected by spatial and temporal expression and functional annotations. A total of 39 putative candidate genes were identified in these intervals (

Table 1. Summary of candidate genes predicted for agronomic traits.

Gene Id	SNP	Trait	Trait & SNP	Gene Name	Protein Name
BnaA01g03380D	A01_1783685	PH; SYP	PH-SYP.A01_1783685	KUP5	K+ uptake permease 5
BnaA01g03490D	A01_1783685	PH; SYP	PH-SYP.A01_1783685	MYB69	MYB domain protein 69
BnaA01g03590D	A01_1783685	PH; SYP	PH-SYP.A01_1783685	CDC20	Cell division cycle 20
BnaA01g03600D	A01_1783685	PH; SYP	PH-SYP.A01_1783685	TIF3K1	Eukaryotic translation initiation factor
BnaA01g03630D	A01_1783685	PH; SYP	PH-SYP.A01_1783685	SLOMO	Slow motion
BnaA02g35610D	A02_random_355708	PH	PH.A02_random_355708	BCP	Blue copper protein
BnaA02g35620D	A02_random_355708	PH	PH.A02_random_355708	ATL54	RING-H2 finger protein ATL54
BnaA03g43740D	A03_22103527	TSW	TSW.A03_22103527	TRANS11	Translocase 11
BnaA03g43820D	A03_22103527	TSW	TSW.A03_22103527	AG1	Floral homeotic protein AGAMOUS
BnaA03g43840D	A03_22103527	TSW	TSW.A03_22103527	MYA1	MYOSIN 1
BnaA03g49330D	A03_25488780	PH	PH.A03_25488780	CFP	Cotton fiber protein
BnaA05g26410D	A05_19368584; A05_19764389	BN	BN.A05_19368584; BN.A05_19764389	DCN1L	Defective in cullin neddylation protein
BnaA05g26840D	A05_19368584; A05_19764389	BN	BN.A05_19368584; BN.A05_19764389	FTM4.1	Leucine-rich repeat protein
BnaA05g26860D	A05_19368584; A05_19764389	BN	BN.A05_19368584; BN.A05_19764389	FTM4.2	Leucine-rich repeat protein
BnaA06g35950D	A06_23598999; A06_23608274	SL	SL.A06_23598999; SL.A06_23608274	MSCT	Man1-Src1p-carboxy-terminal domain protein
BnaA06g35970D	A06_23598999; A06_23608274	SL	SL.A06_23598999; SL.A06_23608274	RPL32A	60S ribosomal protein L32-1
BnaA07g38220D	A07_random_1623623	NSS	NSS.A07_random_1623623	ABH	alpha/beta-Hydrolases superfamily protein
BnaA07g38350D	A07_random_1623623	NSS	NSS.A07_random_1623623	PPR596	Pentatricopeptide repeat 596
BnaA07g38370D	A07_random_1623623	NSS	NSS.A07_random_1623623	KTNA1	Katanin p60 ATPase-containing subunit A1
BnaA09g39690D	A09_28188769	SL	SL.A09_28188769	PGF10	Polygalaturonase clade F 10
BnaA09g39760D	A09_28188769	SL	SL.A09_28188769	UBQ5	Ubiquitin 5
BnaA10g18690D	A10_13564961; A10_13678989	BAI	BAI.A10_13564961; BAI.A10_13678989	BnaA10g18690D	Hypothetical protein
BnaA10g18940D	A10_13564961; A10_13678989	BAI	BAI.A10_13564961; BAI.A10_13678989	BNQ2	Banquo 2
BnaC03g09950D	C03_4731660	MIL; MINS	MIL-MINS.C03_4731660	CBR2L	NADH-cytochrome b5 reductase-like protein
BnaC03g09960D	C03_4731660	MIL; MINS	MIL-MINS.C03_4731660	MPC1	Mitochondrial pyruvate carrier 1
BnaC04g49260D	C04_47349279; C04_47585236	PH	PH.C04_47349279; PH.C04_47585236	JAL24.1	Jacalin-related lectin 24
BnaC04g49270D	C04_47349279; C04_47585236	PH	PH.C04_47349279; PH.C04_47585236	JAL24.2	Jacalin-related lectin 24
BnaC05g10890D	C05_6251826	MINS	MINS.C05_6251826	HRGP	Hydroxyproline-rich glycoprotein family protein
BnaC05g11300D	C05_6554451	HI	HI.C05_6554451	SCP50.1	Serine carboxypeptidase-like 50
BnaC05g11310D	C05_6554451	HI	HI.C05_6554451	SCP50.2	Serine carboxypeptidase-like 50
BnaC05g11350D	C05_6554451	HI	HI.C05_6554451	CRP	Cysteine-rich peptide family protein
BnaC06g20430D	C06_22559430; C06_22570315	NSS	NSS.C06_22559430; NSS.C06_22570315	GSTU20	Glutathione S-transferase U20
BnaC06g20510D	C06_22559430; C06_22570315	NSS	NSS.C06_22559430; NSS.C06_22570315	QUA2	Pectin methyltransferase QUA2
BnaC06g25100D	C06_26638717	PH; NSS	PH-NSS.C06_26638717	PK	Protein kinase superfamily protein
BnaC06g25110D	C06_26638717	PH; NSS	PH-NSS.C06_26638717	ACL	Actin cross-linking protein
BnaC07g36830D	C07_38735522	SL	SL.C07_38735522	PPR334	Pentatricopeptide repeat 334
BnaC07g36960D	C07_38735522	SL	SL.C07_38735522	BnaC07g36960D	Hypothetical protein
BnaC07g37120D	C07_38735522	SL	SL.C07_38735522	EBS	Early Bolting in Short Days
BnaC09g33520D	C09_36899682	PH	PH.C09_36899682	RTN21	Reticulon protein 21

), and the gene expressions during bolting, initial flowering, full-bloom, podding, and maturation stage in the BrassicaEDB were displayed (Figure 6 and Table S4). All of the candidates were highly expressed in certain tissues of the bud, stem, carpel, inflorescence tip, silique pericarp, embryo, and so on.

For HI, BnaC05g11300D (SCP50.1, Serine carboxypeptidase-like 50), BnaC05g11310D (SCP50.2), and BnaC05g11350D (CRP, Cysteine-rich peptide family) were located in the region of C05_6554451. They showed higher expression in inflorescence tip, embryo, and root.

For TSW, three genes including BnaA03g43740D (TRANS11, Translocase 11), BnaA03g43820D (AG1, Floral homeotic protein AGAMOUS), and BnaA03g43840D (MYA1, MYOSIN 1) were located in the candidate region of A03_22103527. For SL, BnaA09g39690D (PGF10, Polygalaturonase clade F 10) and BnaA09g39760D (UBQ5, Ubiquitin 5) were identified in the A09_28188769 interval, and BnaC07g36830D (PPR334, Pentatricopeptide repeat 334) and BnaC07g36960D (Hypothetical protein) were located in the interval of C07_38735522 and C07_38598079. For the two SNPs only significantly associated with NSS, five genes were predicted as candidate genes. BnaC06g20430D (GSTU20, Glutathione S-transferase U20) and BnaC06g20510D (QUA2, Pectin methyltransferase QUA2) were identified in the intervals of C06_22559430 and C06_22570315. The other three, viz. ABH (alpha/beta-Hydrolasese superfamily), PPR596 (Pentatricopeptide repeat 596), and KTNA1 (Katanin p60 ATPase-containing subunit A1), were located in the regions of A07_random_1623623.

For BAI, two putative candidate genes were identified on A10 (A10_13564961 and A10_13678989), including BnaA10g18690D (Hypothetical protein) and BnaA10g18940D (BNQ2, Banquo 2). For the significant SNPs of BN (A05_19368584 and A05_19764389), three promising genes were identified including BnaA05g26410D (DCN1L, Defective in cullin neddylation protein), BnaA05g26840D (FTM4.1, Floral transition at the meristem 4), and BnaA05g26860D (FTM4.2). For PH, 13 candidate genes were identified on three chromosomes, viz. A01, A03, C04, C06, C09, and A02_random for the significant SNPs. Among them, six genes, BnaA02g35610D (BCP, Blue copper protein), BnaA03g49330D (CFP, Cotton fiber protein), BnaC04g49260D (JAL24.1, Jacalin-related lectin 24), BnaC04g49270D (JAL24.2), BnaA02g35620D (ATL54, RING-H2 finger protein ATL54), and BnaC09g33520D (RTN21, Reticulon protein 21), were predicted for controlling plant height. Two genes, BnaC06g25100D (PK, Protein kinase superfamily) and BnaC06g25110D (ACL, Actin cross-linking protein), were identified in the region of C06_26638717, which were associated with PH and NSS. Five genes were located in the interval of A01_1783685 co-detected by PH and SYP, including KUP5 (K+ uptake permease 5), MYB69 (MYB domain protein 69), CDC20 (Cell division cycle 20), TIF3K1 (Eukaryotic translation initiation factor), and SLOMO (Slow motion). These genes were highly expressed in the stem during bolting and flowering stages.

For MINS, BnaC05g10890D (HRGP, Hydroxyproline-rich glycoprotein family protein), BnaC03g09950D (CBR2L, NADH-cytochrome b5 reductase-like protein), and BnaC03g09960D (MPC1, Mitochondrial pyruvate carrier 1) were the promising candidate genes located in C05_6251826 or C03_4731660. For C03_4731660 (MINS and MIL), the genes BnaC03g09950D (CBR2L, NADH-cytochrome b5 reductase-like protein) and BnaC03g09960D (MPC1, Mitochondrial pyruvate carrier 1) were located in the region of C03_4731660.

4. Discussion

HI was found to be one of the crucial factors for enhancing biomass and seed yield [32,33]. However, the harvest index shows a big difference among different crops and is sensitive to environmental factors. In the present study, the correlation coefficient and heritability of HI in different years are relatively low, consistent with previous reports [2,11].

Population genetic structure refers to a non-random distribution of genetic variation in a species or population. In previous studies, B. napus accessions were generally divided into three subgroups, which were mainly contributed to their ecotypes [34,35,36]. There were also reports that GWAS groups were divided into other numbers of subgroups, such as five subgroups [31]. In this study, the ADMIXTURE and STRUCTURE software were used for determining the population structure, results of which were more than three subgroups, probably because 104 inbred lines were mainly semi-winter ecotypes. Combined with the population evolution tree and PCA, the top five principal components were used for the correction of the population structure to control the false-positive results.

For a complex trait, it should take associated traits into consideration, which could provide a better understanding of complex traits for crop [13]. In this study, a total of 14 major agronomic traits were measured that involved plant architecture and economic traits. The partial correlations among HI and the other traits indicated that HI was negatively correlated with BYP and positively correlated with SYP, while some specific traits could connect to HI indirectly by interacting with other traits (Figure S2). In addition, we constructed the association networks across different traits and SNPs (Figure 4). It revealed the complex genetic connections among the agronomic traits. In breeding practice, to improve target traits, favorable alleles of other traits in recipient varieties should be maintained as possible. The network could help to establish a strategy for variety development. For instance, NSS and PH showed the highest connectivity with the most edges linked to other traits, indicating they should be vital traits. If HI is improved by NSS, this should increase SYP and avoid the BYP benefited to break the heritable covariations. Thus, C06_2259430 and A03_10522683 could be the best ones chosen. The same goes for the other traits.

Three hundred and twenty-four SNPs related to 212 SNPs were identified to be involved in the 14 traits. Since the increase of the HI of rice is due to the integration of semi-dwarf genes [6,37], breeders have been working to reduce the plant height of rapeseed to increase the HI of rapeseed. However, the PH usually shows positive correlations with the seed yield and biomass [13,33], indicating PH may be genetically linked to seed yield and biomass. We identified one stable SNPs C09_36899682 detected by PH, which was also related to BYP but had no correlation with SYP. This result provides a novel locus for dwarfing the plant to improve HI.

Plant architecture improvement plays an important role in the increase of crop yield. Several related QTLs of agronomic traits on rapeseed have been reported. Compared to previous reports [10,12,13,30,31,36,38,39,40,41,42,43,44], 49 of 212 loci shared an overlapped region for the same traits based on the reference genome of Darmor-bzh [23] (Figure 5), including PH (11 loci), TSW (9 loci), BAI (7 loci), NSP (6 loci), NSS (5 loci), SL (4 loci), BN (3 loci), MINS (2 loci), SYP (1 locus), and BYP (1 locus). It implied that the significant SNPs of our study were dependable compared with previous results. From the location of the significant SNP, some SNPs for different traits were linked in a certain segment, such as 12.4–13.9 Mb on A10 and 26.6–28.2 Mb on C06. Due to the linkage of the trait locus and the locus pleiotropism, it is challenging to decipher the genetic network of agronomic traits. There are consensus or physical close loci among previous studies and this study, such as the interval of SNPs of A01_1783865 for PH, A03_6291348 for BN, A07_10401876, and A07_11506350 for TSW. These results demonstrate the reliability of this report and provide potential loci for controlling agronomic traits. Some pleiotropic SNPs were also detected, including A01_1783685 (PH and SYP), C06_26638717 (PH and NSS), and C03_4731660 (MIL and MINS), which might be the loci playing multiple roles during plant development. Besides, 12 promising novel significant SNPs were detected that are related to BN (A05_19368584 and A05_19764389), SL (A06_23598999, A06_23608274, and C07_38735522), PH (C04_47349279, C04_47585236, and C09_36899680), MINS (C05_6251826), NSS (C06_22559430 and C06_22570315), and HI (C05_6554451).

HI is mainly related to plant height, inflorescences, flower, pod, branches, and flowers, which is also influenced by environmental conditions such as light, temperature, and nutrient status [45]. Gene functions and expressions help predict candidate genes controlling the traits. In this study, 39 candidate genes were predicted in 22 significant SNP intervals, all of which were highly expressed in relative tissues and involved in the development of cell wall, flower, xylem, and whole plant. The important candidate genes were discussed as follows.

Within the interval of C05_6554451 for HI, BnaC05g11350D codes the CRP protein, while CRPs function as regulators of cell–cell communication in plants, and participate in plant growth, pollen tube growth, and fertilization process [46,47,48].

Plant development is the process of cell division, growth and differentiation. Among these genes, CDC20, identified by A01_1783685 and A01_1565817 for PH and SYP, encodes a cell division cycle 20 protein that is indispensable for normal plant development and fertility. In this region, SLOMO coding an F-box protein is considered as a candidate gene. In Arabidopsis thaliana, the rate of inflorescence meristems organ formation of the slomo mutant is significantly reduced, indicating SLOMO regulates organ initiation at the shoot meristem [49]. KTNA1 was detected on A07_random for NSS, whose orthologs in Arabidopsis thaliana were shown to control cell plate/daughter wall formation [50]. HRGP located on C05 was identified for MINS, encoding a hydroxyproline-rich glycoprotein family protein that functions in early leaf and root vascular differentiation [51].

Floral organogenesis plays a vital role in plant height, branch number, branch angel, main inflorescence length, silique length and width, and number of seeds per silique. In this study, five genes, AG1, BNQ2, TRANS11, EBS, FTM4.1, and FTM4.2, were predicted as the requirements for floral organ growth in Arabidopsis thaliana [52,53,54,55,56]. For example, AG1 is a floral homeotic gene encoding a MADS domain transcription factor, suppression of which favors a reversion of floral meristems from determinate to indeterminate development [52]. FTM4, encoding a leucine-rich repeat protein, is induced in the early inflorescence meristem, and mutations of FTM4 delay flowering [53]. The phenotypes of bnq2 and bnq3 mutants show pale-green sepals and carpels, and purple inflorescence stems and siliques [54], and the double or high order mutants show dwarf vegetative growth and reduced fertility with unequal genetic redundancy [55]. Moreover, we identified some genes functioning in pollen germination and tube growth, including PGF10 and RPL32A. All of these genes showed high expressions in the inflorescence tips, silique pericarp, and stems (Figure 6).

Three candidate genes were related to fiber and pectin involved in SL, PH, and NSS. UBQ5 and QUA2, encoding pectin methyl esterase and pectin methyltransferase, respectively, participate in cell wall formation [57,58,59]. qua2 showed a 50% reduction in homogalacturonan content compared with the wild type [59]. CFP was predicted to encode a cotton fiber protein, but this has not been verified till now.

In addition, KUP5 was located on A01 for PH and SYP, while the uptake and transport of K⁺ were the key to plant growth and responses to the environment. Several genes have not been functionally annotated or their functions are still indeterminate. However, their expressions showed highly specific in certain tissues. Thus, these genes were also considered as putative candidates for HI and HI-related traits.

These B. napus lines are core breeding materials with at least one outstanding trait, from which some excellent rapeseed varieties have been cultivated. We phenotyped the traits over three-year periods and approached GWAS using six algorithms to ensure the accuracy of our GWAS results. However, the number of the GWAS population is not high enough, which might make it prone to false positive results. Hence, significant SNPs are needed to validate the repeatability, and the candidate genes need to be verified as well.

5. Conclusions

In this study, HI and 13 HI-related traits were investigated in 104 core breeding lines of rapeseed. HI showed a complex network with other traits. We performed genome-wide association analyses using the Bnapus50K array. A total of 212 significant SNPs involved in 324 QTL related to the studied traits, and 22 stable SNPs were identified. Three pleiotropic SNPs were obtained, namely A01_1783685 (PH and SYP), C06_26638717 (PH and NSS), and C03_4731660 (MIL and MINS). Twelve promising novel SNPs were detected that are related to BN (A05_19368584 and A05_19764389), SL (A06_23598999, A06_23608274, and C07_38735522), PH (C04_47349279, C04_47585236, and C09_36899680), MINS (C05_6251826), NSS (C06_22559430 and C06_22570315), and HI (C05_6554451). Meanwhile, 39 candidate genes were predicted in stable SNP intervals by their gene expressions and function annotations. The obtained significant SNPs and candidate genes facilitate the development of molecular markers and casual gene cloning. These results provide further support for the improvement of the harvest index in rapeseed.