A panel of 472 rapeseed accessions collected worldwide and stored at the National Mid-term Gene Bank for Oil Crops of China was used for association analysis in this study. The informations of inbred lines about their origin and germplasm type have listed in a previous report (Li et al., 2014).

Field experiments were implemented in six environments across three growing seasons. During the 2013/2014 growing season, the association population was grown at Yangluo (114.50°E, 30.38°N) in Hubei province, which is referred to as E1; during the 2014/2015 growing season, the experiment was conducted at Wuhan (113.68°E, 30.58°N) and Yangluo, which are both in Hubei province, and are referred to as E2 and E3, respectively; and during the 2015/2016 growing season, the association panel was cultivated at Wuhan, Yangluo, and Changsha (113.00°E, 28.22°N, in Hunan province), and are referred to as E4, E5, and E6, respectively. A randomized complete block design with three replicates was adopted in each environment. Each plot contained two rows and 12–15 plants in each row, the distance between plants was 0.2 m within each row, and the space between rows was 0.3 m.

Forty randomly-selected accessions were considered as a sub-panel for observing the positional variation trend of branch angle and for determining the appropriate measurement region on individual plants. In the field, five typical plants in each plot were selected to identify the branch angle at 6 weeks after pollination. The branch angle was defined as the angle between the main stem and its branch and measured by a digital protractor. In each plot of the sub-panel, the branch angle value was obtained by measuring all of the branches of five individual plants. In association panel, the values of the middle three branches of a plant were recorded as the individual plant branch angle value, and the average value of five plants in a plot represents the phenotypic data of a line in this plot.

The broad-sense heritability was estimated according to the following equation: H² = δ_g²/(δ_g² + δ_ge²/n + δ_e²/nr), where δ_g², δ_ge², δ_e², n, and r represent the genetic variance, the interaction variance between genotypes and environments, the error variance, the number of years/locations, and the number of replicates within each environment, respectively. For the branch angle trait, the variance components and best linear unbiased predictors (BLUP) of the multi-environment for each line were estimated using the lme4 package in R software based on a linear model (Merk et al., 2012). The final trait values for association analysis included the BLUP-value and single environment phenotypic data of each accession. The frequency distribution, correlation analysis, and comparative analysis were performed using R software.

In previous reports, detailed descriptions about the process of SNP genotyping and mapping are provided, as are analyses of population structure and linkage disequilibrium (LD) (Li et al., 2014; Wang et al., 2016a).

In brief, the raw SNP data generated from the Brassica 60K Illumina® Infinium SNP array were clustered and automatically called using Illumina BeadStudio genotyping software. Subsequently, 26,841 high-quality SNPs with minor allele frequency (MAF) of more than 0.05 were retained for further analysis. In order to mapping the SNP to an exact position of the reference genome, a BLAST search (Altschul et al., 1990) was performed against B. napus genome sequences (Chalhoub et al., 2014) using the SNP sequences. Only the top and unique blast-hits were reserved.

Eventually, 19,945 SNPs were selected for principal component analysis (PCA), and a relative kinship and population structure analysis. The GCTA tool was used to construct a P matrix of PCA (Yang et al., 2011), SPAGedi software was served to build a K matrix of relative kinship (Hardy and Vekemans, 2002), STRUCTURE v2.3.4 was employed to infer a Q matrix of population structure (Pritchard et al., 2000) and TASSEL 5.0 was used to calculate LD (Bradbury et al., 2007).

The haplotype block structure across 472 rapeseed accessions with the ultimately selected 19,945 SNPs was evaluated using Haploview v4.2 software (Barrett et al., 2005). The analysis referred to the definition of “strong LD” by Gabriel et al. (2002), i.e., the upper minimum of the confidence interval was 0.98 and the lower was 0.7. Furthermore, the fraction of strong LD in informative comparisons must be at least 0.95. Since Haploview software would ignore pairwise comparisons if the distance between markers was beyond 500 kb following the default parameter, to estimate all marker pairs, especially for a strong LD between markers above 500 kb, this setting was adjusted to zero.

The GWAS was implemented using two methods: a MLM (Yu et al., 2006) and a MRMLM (Wang et al., 2016b). The Q+K model, one of the MLMs, including both a fixed effect as the population structure matrix (Q) and a random effect as the kinship matrix (K) was adopted as the optimal model and was performed using TASSEL 5.0 software (Bradbury et al., 2007). An MLM can be described by the following matrix notation: y = Xβ + Zu + e, in which y is the phenotype; X is the genotype; β is a vector containing the fixed effects, including the genetic marker and the population structure (Q); Z is the relative kinship matrix; u is a vector of random additive genetic effects; and e is the unobserved vector of the random residual. The threshold of significant association between a trait and the SNPs in the MLM was p < 1.0 × 10⁻³ [i.e., −log₁₀(p) = 3.0], which has been broadly adopted in the literature (Cai et al., 2014; Hatzig et al., 2015; Raman et al., 2015). The GWAS results were visualized with Manhattan and quantile-quantile (Q-Q) plots that were yielded from the qqman package in R software (Turner, 2014).

An MRMLM, which would improve the power and accuracy of the GWAS, was employed using the R package mrMLM, and the critical log of odds (LOD) score was set as 2.5 (Wang et al., 2016b).

The total phenotypic variation that was explained by the significant SNPs in the best fitting multiple regression model was estimated using the “stepAIC” function from the MASS package in R (Ihaka and Gentleman, 1996).

Two methods were performed to ascertain the region where the potential candidate gene was situated. The first method was based on a definition of the QTL: a locus showing true marker-trait association should harbor at least two SNPs with p-values above the threshold in a 1.5 Mb region (Wang et al., 2016a). Then, for the SNPs that could not be assigned to any QTL, the LD blocks where the associated SNPs were located, in which flanking markers had strong LD (r² > 0.4), were regarded as the candidate gene regions. If the SNPs were also not located in the LD blocks, the 100 kb region centered on an unassigned SNP was considered as the potential candidate gene interval. An LD block analysis was performed using Haploview v4.2 with the default settings (Barrett et al., 2005).

To predict the function of candidate genes, a functional annotation was implemented. First, the protein sequences coded by candidate genes within the definitive region were extracted by referring to the annotation information for the B. napus “Darmor-Bzh” genome (http://www.genoscope.cns.fr/brassicanapus, Chalhoub et al., 2014). Later, the BlastP program was run against Arabidopsis protein sequences with the E-value ≤ 1E-10; then, the candidate genes were functionally annotated using the top hit of Arabidopsis homologous genes. Based on the research progress of branch angles, the orthologous genes involved in gravitropism and auxin transport were focused on.

To validate the expression level of candidate genes between extremely large and small branch angle accessions, four identified candidate genes were randomly selected to perform the quantitative real-time PCR (qRT-PCR) analysis. Gene-specific primers were designed using Primer-BLAST (https://www.ncbi.nlm.nih.gov/tools/primer-blast). For each accession, three middle branches at three weeks after pollination were harvested to extract total RNA. The procedure of total RNA extraction, cDNA synthesis, qRT-PCR amplification, and candidate genes expression analysis were as previously described (Yan et al., 2016). Each sample was examined in three independent biological replicates with three technical replicates.

The obvious tendency for branch angle was observed in different branch positions of individual plants. Generally, the branch angle would increase with elevating branch physical position regardless of plant architecture (Table 1, Figure 1). For example, in the case of accession 1218, representing loose plant architecture, the angle of the first branch (bottom branch) was 35.52° and the angles of the third, fifth, seventh, ninth and eleventh branch rose gradually to 38.96°, 46.03°, 52.00°, 53.40°, and 56.80°, respectively (Table 1, Figure 1). Another line, 3304, representing compact plant architecture, exhibited a similar branch bending pattern (Table 1, Figure 1). Therefore, the linear growth of branch angle may be affected by the dose-response of auxin concentration. To determine the appropriate measurement region on individual rapeseed plants, a t-test of significant differences in different branch regions was conducted. The results demonstrated that there were distinct differences between the upper plant branch angle and the whole plant branch angle in the majority of the rapeseed accessions. The lower plant branch angle showed a similar result, but the middle plant branch angle showed no marked difference compared to the whole plant branch angle (Table S1). For example, in accession 1218, the mean values of the upper, middle, lower and whole plant branch angle were 51.71°, 44.84°, 37.54°, and 44.79°, respectively, and it is evident that the middle plant branch angle was much closer to the whole plant branch angle. Thus, the most representative phenotypic data could be obtained by just measuring the middle branches.

A distinct phenotype variation in branch angle, ranging from 17.5° to 53.6°, was found across the 472 rapeseed accessions in the six environments (Figure 2, Table 2). The maxima in the observed phenotype data were 1.8–2.5 times the minima, varying from 17.5° to 43.9° in E1, 28.0° to 49.3° in E2, 26.4° to 48.3° in E3, 20.6° to 46.9° in E4, 22.1° to 50.4° in E5 and 24.1° to 53.6° in E6. Moreover, two adjacent locations, Wuhan and Yangluo, exhibited analogous phenotypic variation both in 2015 (E2 and E3) and in 2016 (E4 and E5), indicating that the branch angle is a relative stably inherited trait (Table 2). The correlation coefficient of branch angle, ranging from 0.49 to 0.64, indicated that the branch angles in six environments have a significantly positive correlation (p < 0.05, Table 3).The broad-sense heritability (H²) of branch angle among the rapeseed panel was 76.06% (Table S2), suggesting that environmental factors had limited influence on the branch angle, which exhibited a fairly stable manner.

A total of 2,423 conserved haplotype blocks were detected after estimating 19,945 high-quality SNPs distributing on the whole-genome in 472 rapeseed accessions, which spanned 181.53 Mb and covered 28.28% of the assembled B. napus genome (Figure 3, Table S3). The haplotype block position, length, and SNP number within each block is also provided in Table S3. The average number of haplotype blocks in the A-subgenome chromosomes was 143.6 (ranging from 62 to 246) with an average block size 34.10 kb (ranging from 19.30 to 78.97 kb), and a haplotype block coverage percentage ranging from 7.23 to 27.52%, with a mean percentage of 21.00% (Table 4, Figures 4A–C). In the C-subgenome chromosomes, the haplotype block number varied from 57 to 209 (average number of 109.7) with a fairly larger haplotype block size ranging from 80.86 to 274.85 kb (average size 134.30 Kb), and with a block coverage proportion varing from 14.87 to 51.49% (average proportion of 32.80%, Table 4, Figures 4A–C). In the A-subgenome, haplotype blocks less than 30 kb in size were roughly four-fifths of all blocks (80.3%), and, this proportion was more than one-half (53.2%) in the C-subgenome (Figure 4D).

There were 23 haplotype blocks whose size was more than 1 Mb, accumulatively accounting for more than one-third of the total block size (Table 5). These blocks were distributed on 12 rapeseed chromosomes and most of them (18/23) were located on the C-subgenome (Table 5). Approximately half of the blocks (11/23) were across or in the vicinity of their corresponding centromere (Table 5 and Figure 3, Mason et al., 2016), meaning that stronger LD existed in these blocks; furthermore, if the SNPs located on the blocks are excluded, the LD decay will depress sharply (Qian et al., 2014; Sun et al., 2016b).

A total of 144 and 69 significantly-associated SNPs of branch angle were detected, which could explain up to 62.2 and 66.2% of the cumulative phenotypic variation, using the BLUP value and individual environment in MLM and MRMLM, respectively (Figure 5, Table S4). The significantly-associated SNPs corresponded to 46 and 38 loci in MLM and MRMLM, respectively (Figure 3, Table S5), and 21 loci among them partly shared the region between the two models, which accounted for 45.65% (21/46) and 55.26% (21/38) of the total identified loci (Table 6). For instance, two vicinity loci on A6 were both repeatedly detected in the two models. In MLM, the first locus on A6 was crossed from 19,416,065 to 20,815,553 with a peak SNP (highest significant) of Bn-A06-p18028879, which contributed to 4.14% of the phenotypic variance. In comparison, the homologous locus on A6 in MRMLM spanned from 19,630,281 to 20,815,553, with a peak SNP Bn-A06-p18246821, which explained 5.01% of the phenotypic variance. The second locus on A6 in MLM was crossed from 23,211,156 to 23,499,018, with a peak SNP Bn-A06-p24551529 which contributed to 3.53% of the phenotypic variance. In MRMLM, the corresponding locus on A6 spanned from 23,362,162 to 23,495,968, with a peak SNP Bn-A06-p24544753, which explained 3.75% of the phenotypic variance. In a word, the two loci on A6 in MLM shared 84.7 and 46.5% of their regions with the homologous loci in MRMLM (Table 6).

Furthermore, 45.65% (21/46) and 52.63% (20/38) of the identified loci in MLM and MRMLM were verified in at least two environments, illustrating that our association results were credible and reproducible (Table S5). Furthermore, these loci were distributed on all the chromosomes in both models except for A2 in MLM (Figure 3, Table S5). Approximately three-fifths (27/46) of the loci were located on the A sub-genome in MLM, and the average number of loci in each chromosome was 2.5 (ranging from 1 to 4; Figure 3, Table S5). In MRMLM, a similar proportion of loci (23/38) were located on the A sub-genome; the average number of loci in each chromosome was 2.0 (ranging from 1 to 4; Figure 3, Table S5).

Using the A. thaliana orthologous genes and published literatures about branch angle as a reference, altogether 73 and 65 candidate genes, corresponding to 32 and 28 loci, were identified in MLM and MRMLM, and 43 candidate genes were commonly detected in the two models (Figure 3, Table S5).

LAZY1 and TAC1 are well-known genes modulating branch angle, whose mutants display converse branch angle morphology in plants (Yu et al., 2007; Yoshihara et al., 2013). We identified two orthologs of LAZY1, BnaA10g19550D and BnaC03g06250D, at 13.9 Mb on A10 and 3 Mb on C3, which are 2,213 kb and 2,202 kb from the peak SNPs of Bn-A10-p10252741 and Bn-scaff_16614_1-p1291979, respectively (Table 7, Table S5A). Furthermore, the TAC1 ortholog, BnaC04g00780D, was detected at 0.7 Mb on C4, which is 183 kb from the peak SNP Bn-scaff_16935_1-p98166 (Table 7, Table S5A).

The mutants of sgr1–sgr7 showed a defective in gravitropic response resulting in the alteration of normal branch angle in Arabidopsis (Fukaki et al., 1996; Yamauchi et al., 1997). Two orthologs of SGR3 in rapeseed, BnaA06g35880D and BnaC09g19750D, were identified at 23.5 Mb on A6 and 16.8 Mb on C9, which were 5 kb from the peak SNP Bn-A06-p24551529 and 551 kb away from the peak SNP Bn-scaff_15650_1-p624000, respectively (Table 7, Figure 6A, Table S5A).The haplotype of the peak SNP Bn-A06-p24551529 for SGR3 was analyzed, and 472 rapeseed genotypes were classified into four haplotype groups (Figure 7A). Haplotype 3 (Hap3, n = 250) was the largest group, Hap1 (n = 137) and Hap3 (n = 74) were the second and third largest group, and Hap4 (n = 3) was a minor group comprising a few rapeseed lines. Statistically, accessions with Hap1 and Hap2 had a significantly lower branch angle than those with Hap3 (P = 2.86 × 10⁻⁶ and 4.08 × 10⁻⁴, respectively, Figure 7A).

The ortholog of SGR1, BnaC08g25070D, was located at 26.9 Mb on C8, which is 43 kb from the peak SNP Bn-scaff_16770_1-p4296727 (Table 7, Figure 6B, Table S5A). The results of the haplotype effect analysis for Bn-scaff_16770_1-p4296727 illustrates that the average branch angle of individuals with the Hap1 allele was prominently lower than that with Hap4 (P = 1.17 × 10⁻³, Figure 7B). The SGR5 ortholog BnaA06g34390D was detected at 22.7 Mb on A6, which is 779 kb upper from the peak SNP Bn-A06-p24551529, which was shared with SGR3 (Table S5A). In addition, the ortholog of SGR7, BnaA08g15740D, was characterized at 13.1 Mb on A8, which is 205 kb from the SNP Bn-A08-p15792942 (Table 7, Table S5A).

The branch curvature growth results from auxin asymmetry accumulation between the upper and bottom portion of this organ; the genes involved in auxin homeostasis are required in this process (Roychoudhry and Kepinski, 2015). This category gene was also characterized in our study, for example, the ortholog of PIN3 BnaA07g23670D, the member of the auxin efflux carrier family, was detected at 17.8 Mb on A7, which is 1,085 kb down from the peak SNP Bn-A07-p14798978 (Table 7, Table S5A). More information about the other genes identified in the present study for branch angle are available in the Table S5A; the above-mentioned orthologs in rapeseed were derived from the MLM, and the homologous genes from MRMLM are listed in Table S5B.

Four identified candidate genes, i.e., TAC1, SGR1, SGR3, and SGR5, were selected to validate the gene expression level between extremely large branch angle lines (1218 and 3078, with average branch angle of 44.12° and 47.53°, respectively) and extremely small branch angle lines (2874 and 3304, with average branch angle of 26.84° and 28.33°, respectively). Gene-specific primers were listed in Table S6. As shown in Figure 8, the expression patterns of the four candidate genes detected by qRT-PCR showed significant difference between extremely large and small branch angle lines, confirming the reliability of the association mapping results. For instance, the expression levels of TAC1 in line 1218 and line 3078 were significantly higher than that in line 2874 and line 3304 (P < 0.05, Figure 8A). And the expression levels of SGR1, SGR3, and SGR5 in line 1218 and line 3078 were significantly lower than that in line 2874 and line 3304 (P < 0.05, Figures 8B–D).

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.