The association population used in this study consisted of 333 diverse rapeseed accessions (20 are winter type, 308 are semi-winter type and 5 are unknown) collected as part of a recently published study (Sun F. M. et al., 2016). The rapeseed accessions were grown in two replicates over the course of 4 years (2012–2015) in Wuhan, China (E 114.32°, N 30.52°). Individuals from germplasm populations were genotyped using the 60 K Brassica Infinium SNP array. Accessions were grown in plots of 2 rows, with 10–15 plants in each row. Between eight and ten plants were selected from the center of each row to measure the three traits at maturity: PH, BIH, and BN.

The length of the plant from the base of the stem to the tip of the main inflorescence was noted as PH. The length from the base of the stem to the first primary branch base was recorded as BIH. We also measured the number of primary branches arising from the main shoot (BN). Correlation analyses of PH, BIH, and BN for the association panel were performed in R3.3.0 (Ihaka and Gentleman, 1996).

The probe sequences of 52,157 SNPs were used to perform a BLASTN search through the B. napus genome database (Chalhoub et al., 2014). Only top blast-hits with an e-value threshold of e⁻¹⁵ were used. Some SNPs were excluded due to a match in BLAST with multiple loci with identical scores. SNPs that showed a minor allele frequency (MAF) of <5% or that showed a maximal missing frequency of >2% were excluded (Brachi et al., 2011), remaining 32,297 SNPs identified, covering all 19 chromosomes of the B. napus genome and providing approximately one SNP per 25 kb. An estimate of the linkage disequilibrium (LD) was calculated using squared allele frequency correlations (r²) between all pairs of SNPs using TASSEL 5.2.28 (Bradbury et al., 2007).

All the SNPs were used to calculate the principal component analysis (PCA) matrix using the GCTA tool, and a subset of 3,571 SNPs (MAF ≥ 0.2) evenly distributed across the whole genome (every 150 k) were selected to perform population structure (Q) and relative kinship analysis (K). STRUCTURE v2.3.4 was used to calculate a Q matrix, the putative number of genetic groups (k-value) setting from 1 to 10 with five independent runs (Pritchard et al., 2000). The length of the burn-in period and the number of Markov Chain Monte Carlo (MCMC) replications after burn-in were set to 50,000 and 100,000 iterations, respectively. The relative kinship K matrix was performed by the SPAGeDi software package (Hardy and Vekemans, 2002; Yang et al., 2011). The Δk method was used to determine the most likely number of groups or subpopulations, as described by Evanno et al. (2005).

Association analysis was performed using TASSEL 5.2.28 (Bradbury et al., 2007), using a mixed linear model (MLM) to calculate the association in all analyses, incorporating Q matrix/PCA and kinship data (K) (Zhao et al., 2011). The MLM was applied using default settings (P3D for variance component analysis and compression set to optimum level). For MLM (Q + K), the significance threshold for significantly associated markers was set to p ≤ 4.06 × 10⁻⁴ [−log₁₀ (p-value) = 3.39].

To identify candidate genes, local LD decay was calculated using TASSEL5.2.28 to capture flanking regions of up to 450 kb on either side of significant SNPs, with a cut-off value of r² = 0.2. Gene sequences that correspond to putative orthologs in A. thaliana were based on GO annotations (http://www.arabidopsis.org/index.jsp). The genes with GO terms for auxin, GA, IAA, SL, CK, and flowering time were highlighted, and the closest one of these genes was considered as the most likely candidate.

Extensive phenotypic variations for PH, BIH, and BN were observed in the 333 accessions that were grown over the course of 4 years (2012–2015). PH varied from 86.2 to 206.0 cm, with 1.6 to 2.4-fold variations across the 4 years (Figure 1A, Table 1). BIH varied from 6.9 to 122.0 cm, with 5.9 to 17.3-fold variations across the 4 years (Figure 1B, Table 1). BN varied from 1.7 to 15.0, with 3.5 to 8.6-fold variations across the 4 years (Figure 1C, Table 1).

In addition, we assessed these three traits for any significant correlation. PH and BIH were significantly correlated across the 4 years with a Pearson's correlation coefficient of 0.3 to 0.8 (Table S1). This suggests that these two traits may show genetic linkage or that some genes have pleiotropic roles in controlling these phenotypes. However, no other significant correlations were found in our accession panel, between the PH and BN, or between BN and BIH (Table S1).

SNP genotyping was performed using the Brassica 60 K SNP array. We blasted the SNP probe sequences to the B. napus genome database (http://www.Genoscope.cns.fr/brassicanapus/) and a total of 34,292 SNP markers (65.7%) were selected to genotype the panel of 333 accessions (Table 2). The C04 linkage group had the most SNPs (3,182), with a marker density of one per 20.89 kb, whereas the C05 linkage group had the least SNPs (1,085), with a marker density of one per 39.8 kb (Table 2). Altogether, the overall mapping results demonstrate the high quality of the genotyping in this study.

We calculated LD in 333 accessions using the parameter r² with the 34,292 SNP markers. In this analysis, three chromosomes (A08, C02, and C07) showed strong LD, with distances ranging from 1,300 to 1,710 kb; while others exhibited modest LD, with distances ranging from 100 to 825 kb when r² = 0.2 (Table 2). The average LD decay of chromosomes A and C were of 190 and 700 kb respectively, while the average LD decay of the whole genome (A + C chromosomes) was of 450 kb when r² = 0.2 (Sun F. M. et al., 2016; Figure 2A).

We performed GWAS with mixed linear modeling (MLM) across 4 years. For the groups-populations, we created a Q matrix (k = 2) for all 333 accessions, in order to obtain the highest Δ k-value (Figure 2B). In addition, minor differences were observed among the PCA, Q, K, Q + K, and PCA + K models, and the Q + K model is an appropriate model (Figure 3). Thus, the Q + K model was used for GWAS. As in other GWAS studies, we set a threshold of −log₁₀ P = 3.39 as a significant association (Figure S1; Benjamini and Hochberg, 1995; Wang et al., 2016). Furthermore, QQ plots of PH in 2012 and 2013 were not linear. As a result, the SNPs showing the 50 highest −log₁₀ P-values in each year were compiled for the PH trait for these 2 years to identify potentially significant SNP clusters (Figure 4; Ueda et al., 2015). Overall, for all traits analyzed, a total of 158 SNPs associated with PH, BIH, and BN were identified (Table S2). Furthermore, LD blocks (450 kb) harboring the repeating SNPs across the 4 years (and detected in at least two of those years) were identified as regions containing putative candidate loci. This led to the selection of 38 SNPs (Table 3), which in turn could be merged into seven, four, and five loci for PH, BIH, and BN, respectively (Figures 4–6, Table 3, Figures S2–S4).

For PH, we detected nine GWAS peak SNPs in the seven loci located on chromosomes A02, A06, C03, and A07, and five of which were detected in at least 2 years (Figure 4, Table 3). These seven peak SNPs explained 35.56% of the total phenotypic variance, with the largest contribution (R² ~9.79%) from SNP Bn-A02-p9505646 (Table 3). For BIH, four peak SNPs were detected on chromosomes A02, A07, A08, and A09, and three of these SNPs were detected in at least 2 years (Figure 5, Table 3). These three SNPs explained 27.96% of the phenotypic variance (Table 3). SNP Bn-A02-p9505646 was detected in both PH and BIH traits and it accounted for 9.36% of the total phenotypic variance in BIH (Table 3). For BN, five peak SNPs were detected on chromosomes A01, A03, A07, C04, and C09 (Figure 6), two of which were detected in at least 3 years, with the largest contribution (R² ~11.54%) from SNP Bn-A07-p21413042 (Table 3).

The GWAS identified several significant SNP markers in the regulation of three traits (Table 3). To find useful SNP markers for marker-based breeding, we estimated the allelic effect of the peak SNPs across these three traits. The marker Bn-A02-p9505646 (G/A) of the GG allele showed the largest contribution (R² ~9.79%) to PH regulation, measuring 21.5 cm more than those of the AA allele (P < 0.05; Figure 7A). Apart from the Bn-A02-p9505646 marker, individuals with the CC allele marker Bn-A09-p9669847 showed the largest contribution to BIH regulation (R² ~9.87%) and measured 11.4 cm more than those of the AA allele (P < 0.05; Figure 7B). The average BN of individuals with the AA allele with marker Bn-A03-p7326202 showed the higher contribution (R²~8.61%) and was significantly greater than GG alleles on chromosome A3 (Figure 7C). These SNP markers may permit marker-based breeding for rapeseed plant architecture.

We searched for candidate genes from the genomic regions with significant SNPs according to GO annotations from the Arabidopsis database (http://www.arabidopsis.org/index.jsp). Notably, 31, 15, and 17 candidate genes for PH, BIH, and BN were located in each region, respectively (Table S3).

For PH, the main QTL was detected on chromosome A02, located at 6.0~6.7 M. From the peak SNP Bn-A02-p9607756 [−log₁₀ (p) = 5.47], two candidate genes associated with GA signaling and flowering time repression were located at 77.2 and 29.4 kb, respectively (Table 4). They encode proteins GAI-DELLA (BnaA02g12260D) and FT (BnaA02g12130D), which may be involved in PH formation (Wang and Li, 2008; Salas Fernandez et al., 2009). A BR biosynthesis gene CYP450 (BnaA02g02150D), a gibberellin-regulated protein GASA4 (BnaA02g02560D), and two cell wall biosynthesis proteins (BnaA02g02700D and BnaA02g03310D) were identified up- and down-stream of the peak SNP of Bn-A02-p3760177 [−log₁₀ (p) = 4.90] (Table S3). In addition, we also identified four candidate genes for PH regulation on chromosomes A06 and A07 (Table 4).

For BIH, the C2H2 zinc finger protein, TCP 12, and SOB five-like 2 (BnaA02g13870D, BnaA02g14010D, and BnaA02g14060D) were also identified in the candidate region around the peak SNP Bn-A02-p9505646 [−log₁₀ (p) = 4.70] on chromosome A02 (Table 4). On chromosome A07, a CLAVATA3 protein and a LOF2 transcription factor (BnaA07g24210D and BnaA07g24240D) were identified downstream from the peak SNP Bn-A07-p15758978) [−log₁₀ (p) = 4.88] (Table 4). Their homologs in Arabidopsis are also involved in meristem initiation and maintenance (Lee et al., 2009). In addition, on chromosome A09, we found a peak SNP Bn-A09-p9669847 with the highest probability of marker-trait association, with −log₁₀ (p) = 5.04, and with the largest contribution to BIH variation, with R² = 9.87%. Notably, a MATE efflux family protein (BnaA09g14730D) involved in determining the rate of organ initiation was identified close to the peak SNP (Table 4).

For BN, most candidate genes around peak SNPs were associated with Auxin, CK, and SL (Table S2). Five candidate genes, ARR11, TCP12, ARF8, LOF2, and CUC3, were identified up- and downstream of the peak SNP Bn-A07-p16996070 [−log₁₀ (p) = 6.08] (Table 4). Another peak SNP Bn-A01-p4070059 [−log₁₀ (p) = 3.93] contained two candidate genes (YUCC8 and IAA11). On chromosome C09, a WOX5 gene (BnaC09g49730D) was located 0.45 M from the peak SNP Bn-scaff_28053_1-p60542 [−log₁₀ (p) = 6.18] (Table 4).

GWAS was used to detect QTLs of PH, BIH, and BN, which are main factors regulating plant architecture in rapeseed. According to the marker-trait association SNPs, we identified seven QTLs for PH, four for BIH, and five for BN (Table 3). In the PH candidate loci, these QTLs were distributed on six chromosomes and some QTL regions were similar to those identified in previous studies. Among them, the QTLs at 1.2 and 6.4 M on chromosome A02 overlapped regions detected by Udall et al. (2006); Li X. N. et al. (2013), and Sun et al. (2016b). The QTL at 18.0 M on chromosome A07 was close to the QTLs detected by Udall et al. (2006) and Li et al. (2016). The QTL near the top of chromosome C03 was also detected by Ding et al. (2012) and Shi et al. (2009). Additionally, the QTL at 6.2 M on chromosome A03 was found by Mei et al. (2009) and Li X. N. et al. (2013), and also corresponds closely to the findings of Li et al. (2016). The QTL at 45.8 M on chromosome C02 corresponds closely to the QTL on the bottom of C02 detected by Udall et al. (2006). These results indicate that these QTLs were stable as they were detected by different map populations and in different environments using distinct methods of analysis, and they may also have been selected in rapeseed breeding. Notably, the QTL at 19.7 M of A06, which was detected in 2 years has not been found in previous studies, and may be a novel QTL.

For the BN candidate loci, only one QTL was located at 6.6 M on chromosome A03, as shown previously in other studies (Shi et al., 2011; Li et al., 2016). Interestingly, we found that all the other four QTLs were detected in a rapeseed accession that has a stable BN but increased branching. In our study, most of the accessions are varieties cultivated from different areas, we therefore deduced that the four QTLs might have converged or have been artificially selected in rapeseed breeding. However, several QTLs of PH and BN regulation reported in other studies were not detected in this study. This may be ascribed to insufficient numbers of rapeseed accessions in GWAS or to diversity of climate. In addition, all QTLs might not have been detected due to an insufficient density of SNPs in some genomic regions.

Previously, only one study identified QTLs associated with the height of the lowest primary effective branch (HPB) using two populations in rapeseed, and it identified 10 QTLs distributed on chromosomes A02, A07, A08, C04, C06, and C07 (Chen et al., 2007). In our study, BIH was measured from the ground level up to the base of the lowest primary effective branch, and four QTLs were detected in our natural population. One of these QTLs was located at 18.0 M on chromosome A08, which corresponds to the QTL detected on the bottom of the same chromosome by Chen et al. (2007). Two more QTLs were found on corresponding chromosomes A02 and A07. However, their specific regions did not overlap. To verify the reliability of these QTLs, further investigations of this trait must be undertaken.

In our study, the QTLs at 6.3 M on chromosome A02 and at 18.0 M on chromosome A07 influenced PH and BIH, and a significant correlation was found between these two traits during our correlation assay. In the A02 QTL regions, PH and BIH traits shared three SNPs and the two traits were also co-localized, which suggests pleiotropic regulation of a single gene or the existence of closely linked genes. In contrast, in the A07 QTLs, the SNP associated with the PH QTL was separated from the related SNP associated with the BIH QTL by about 700 kb, while the LD (r²) between the two SNPs was 0.27, indicating that they are separate QTLs. In addition, QTLs related to PH on chromosome A03 and A07 were found close to the QTLs related to BN on the same chromosomes, at a distance of 500 kb and 1 M, respectively. The QTL at 17.7 M on chromosome A07 related to BIH regulation was close to the QTL related to BN, at ~1.3 M. These distances are also beyond the average LD decay and the LD (r²) between the two peak SNPs for each group were relatively low, suggesting that they are separate QTLs.

In Arabidopsis, rice, and other plant species, BR, IAA, GA, and SL biosynthesis and signaling pathways are known to regulate PH, while mutations in most of these genes cause dwarf phenotypes (Wang and Li, 2008; Sun et al., 2010; Clouse, 2011). Additionally, genes involved in flowering time also play a role in PH regulation, such as DTH8 and Ghd7 in rice (Xue et al., 2008; Wei et al., 2010), and Flt-2L in wheat (Chen et al., 2009). In rapeseed, Mei et al. (2009) found an overlap between a flowering time locus and PH QTLs. Moreover, in some genes regulating cell wall formation, like cellulose synthase genes, mutations affecting cell elongation lead to dwarf phenotypes (Tanaka et al., 2003). Therefore, we defined PH regulation candidate genes in rapeseed based on similarity to orthologs of the QTL regions. Using the average LD decay, 5 genes (16.7%) are related to BR biosynthesis or signaling pathways, 6 genes (20%) are related to GA, 2 genes (6.7%) are associated with Auxin, 4 genes (13.3%) are involved in cell wall formation, 4 genes (13.3%) are related to flowering time regulation, and 2 gene (6.7%) are related to trichome branching development in Arabidopsis (Bischoff et al., 2010; Table S3). Notably, the RGA gene, encoding a DELLA protein, displays a signal close to Bn-A02-p9505646, which itself is close to the Bn-A02-p9610453 marker identified in Sun et al. (2016b). Furthermore, we found that some flowering time QTLs (BnFT) are also located in the same region of chromosome A02 (Xu et al., 2016). In Arabidopsis and rice, most studies showed that FLOWERING LOCUS T (FT) and the rice FT homolog HEADINGDATE3a (Hd3a) gene affected flowering and PH (Kardailsky et al., 1999; Kobayashi et al., 1999; Tamaki et al., 2007). In addition, overexpressed the FT in tobacco also exhibited the dwarf phenotype (Lewis and Kernodle, 2009). Thus, we deduce that BnRGA (BnaA02g12260D) and BnFT (BnaA02g12130D) are main candidate genes on chromosome A02 that are involved in PH regulation.

Several transcription factors (TF) and hormones (Auxin, CKs, and SLs) also affect shoot outgrowth in plant architecture (Yang and Jiao, 2016). We compared these orthologs to our main candidate genes for BN. Four genes (23.5%) are related to CK biosynthesis or signaling pathways, 6 genes (35.3%) are related to IAA, 3 genes (17.6%) are associated with flowering or flower identity formation, and one gene is involved in the SL signaling pathway. Many recent studies report that SL is functionally important in BN and PH regulation; however, only one candidate gene was identified in our GWAS. This may be due to a low density of SNP markers, or due to the candidate genes being farther away from the QTL locus than the average LD decay, or due to the confounding of the population structure. In Arabidopsis, LATERAL ORGAN FUSION2 (LOF2) encodes a MYB-domain transcription factor that functions in both lateral organ separation and axillary meristem formation, in part through interaction with CUC2, CUC3, and STM (Raman et al., 2008; Lee et al., 2009). In novel loci on chromosome A07, genes BnLOF2 (BnaA07g26170D) and BnCUC3 (BnaA07g28210D) were identified as the most likely candidate genes and may play roles in both organ separation and axillary meristem formation.

Although several candidate genes were identified following the GWAS, the function of these genes remains ambiguous. In the future, it is necessary to isolate target genes by cloning these putative genes or by developing mapping populations for the QTLs, and to illustrate their functional role in relative trait regulation by transformation experiments.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.