ORIGINAL RESEARCH article
Genome-Wide Association Mapping Reveals the Genetic Control Underlying Branch Angle in Rapeseed (Brassica napus L.)
- 1Oil Crops Research Institute of the Chinese Academy of Agricultural Sciences, Key Laboratory of Biology and Genetic Improvement of Oil Crops, Ministry of Agriculture, Wuhan, China
- 2National Key Lab of Crop Genetic Improvement, National Center of Crop molecular Breeding, National Center of Oil Crop Improvement, College of Plant Science and Technology, Huazhong Agricultural University, Wuhan, China
Plant architecture is vital not only for crop yield, but also for field management, such as mechanical harvesting. The branch angle is one of the key factors determining plant architecture. With the aim of revealing the genetic control underlying branch angle in rapeseed (Brassica napus L.), the positional variation of branch angles on individual plants was evaluated, and the branch angle increased with the elevation of branch position. Furthermore, three middle branches of individual plants were selected to measure the branch angle because they exhibited the most representative phenotypic values. An association panel with 472 diverse accessions was estimated for branch angle trait in six environments and genotyped with a 60K Brassica Infinium® SNP array. As a result of association mapping, 46 and 38 significantly-associated loci were detected using a mixed linear model (MLM) and a multi-locus random-SNP-effect mixed linear model (MRMLM), which explained up to 62.2 and 66.2% of the cumulative phenotypic variation, respectively. Numerous highly-promising candidate genes were identified by annotating against Arabidopsis thaliana homologous, including some first found in rapeseed, such as TAC1, SGR1, SGR3, and SGR5. These findings reveal the genetic control underlying branch angle and provide insight into genetic improvements that are possible in the plant architecture of rapeseed.
In nature, a particular plant specializes its architecture and corresponding function. For crops, such as rapeseed (Brassica napus L.), the desirable architecture is able to produce high grain yields (Wang and Li, 2008). Shoot branching, such as branch angle (BA), is a principal factor in plant architecture (Ariyaratne et al., 2009). Plant density is a vital environmental factor influencing the plant architecture (Diepenbrock, 2000), and results in the capacity to bend the branches to suitable angles for increasing the photosynthetic efficiency. Mechanical harvesting, which is affected by many aspects especially shoot branching, is an inevitable option for the rapeseed industry in the future because of the resulting decreases in required labor resources. A higher plant density with decreased branching angle would produce the highest mechanical seed yields in rapeseed (Kuai et al., 2015).
Essentially, branch angle is one of the ways to adapt to diverse environmental conditions through gravitropism in plants, whereas gravitropism is achieved through asymmetric distribution of the auxin concentration (Roychoudhry and Kepinski, 2015). Considerable genes modulating branch angle have been identified in plants. In rice, lazy1 was insensitive to gravity stimuli and exhibited a prostrate morphology due to impaired polar auxin transport (Li et al., 2007), and then, the orthologs of LAZY1 in Arabidopsis and maize were characterized and cloned (Dong et al., 2013; Yoshihara et al., 2013). In contrast to the lazy1 mutant, the tac1 mutant had an almost vertical tiller angle in rice (Yu et al., 2007), and TAC1 played an antagonistic role to LAZY1, although they belong to the same IGT gene family (Dardick et al., 2013). Furthermore, the TAC1 orthologs in maize and Miscanthus were reported to regulate leaf angles (Ku et al., 2011; Zhao et al., 2014). A series of sgr mutants have been identified and manifested as a defective in gravitropic response that lead to a discrepant branch growth angle in Arabidopsis (Fukaki et al., 1996; Yamauchi et al., 1997). For instance, the normal endodermis, where gravity-sensing cells are located, was absent in the hypocotyls and inflorescence stems of sgr1 and sgr7 mutants (Fukaki et al., 1998). Furthermore, aberrant vacuoles affect amyloplast accumulation in the tissues of sgr2–5 mutants (Hashiguchi et al., 2013). The other auxin homeostasis genes were also reported as participating in branch angle regulation, for example, gravity-induced PIN3 polarization diverts the auxin flow to mediate the asymmetric distribution of auxin for shoot bending (Rakusova et al., 2011).
A genome-wide association study (GWAS) has emerged as a powerful approach for dissecting important causal loci that correlated with complex traits. An excellent opportunity for insight into the genetic basis of agronomic traits at the DNA level for rapeseed is provided by the development of a 60K Brassica single nucleotide polymorphism (SNP) Infinium array (Edwards et al., 2013) and the completion of B. napus genome sequencing (Chalhoub et al., 2014). Thus, association mapping has been widely implemented for numerous traits in rapeseed in recent years, such as the seed weight, plant height, oil content, and clubroot resistance (Cai et al., 2014; Li et al., 2014, 2016a,b). Nevertheless, an association analysis for branch angle in rapeseed has not been well elucidated. Liu et al. (2016) detected 25 significantly associated quantitative trait loci (QTLs) and identified three candidate genes, including LAZY1, for branch angle across 143 rapeseed accessions. As reported by Sun et al. (2016a), 56 loci significantly associated with branch angle among 520 rapeseed accessions were confirmed, and many candidate othologs were detected, such as LAZY1, SGR2, and PIN3. However, some vital potential genes for branch angle, such as TAC1, SGR1, SGR3, and SGR5, have not been detected and remain to be further mined.
In this study, a massive phenotypic identification of branch angle was conducted in the association mapping of a population of 472 diverse rapeseed accessions in six different environments; the association mapping population was genotyped with a high-through 60K SNP array. Genome-wide association analysis was performed using two models, mixed linear model (MLM) and multi-locus random-SNP-effect mixed linear model (MRMLM), and 46 and 38 loci significantly associated with branch angle were mined, respectively. Subsequently, considerable highly-promising candidate genes were identified by annotating against Arabidopsis thaliana homologous. These findings will sharpen our understanding of genetic mechanisms underlying branch angle and will provide insight into genetic improvements that are possible for the plant architecture of rapeseed.
Materials and Methods
Plant Materials and Field Experiments
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.
Trait Measurements and Statistical Analysis
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: H2 = δg2/(δg2 + δge2/n + δe2/nr), where δg2, δge2, δe2, 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.
Genotype Data Acquisition
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).
Haplotype Block Structure Analysis
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.
Genome-Wide Association Study
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−3 [i.e., −log10(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).
Candidate Genes Identification
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 (r2 > 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.
Quantitative Real-Time PCR Analysis
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.
Positional Variation of the Branch Angle on Individual Plants
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.
Figure 1. Models of branch angle for loose and compact plant architecture in rapeseed. (A) A model for branch angle of loose plant architecture in accession 1218 with a large branch angle and the angle tendency of different branch positions. (B) A model for branch angle of compact plant architecture in accession 3304 with a small branch angle and the angle tendency of different branch positions.
Phenotypic Variation of Branch Angle in an Association Mapping Population
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 (H2) 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.
Figure 2. The distribution of branch angle across 472 rapeseed accessions in six environments. (A) The environments of E1, E2, and E3. (B) The environments of E4, E5, and E6.
Table 2. Phenotypic variation in branch angle for a rapeseed association population in six environments.
Haplotype Block Structure Study
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).
Figure 3. Mapping of haplotype blocks and association loci for branch angle on the Brassica napus genome. The black bands in the outer circle indicate the haplotype blocks; the location, and size of centromeres are marked by the rectangles attached to the outer circle. The bands in the middle and inner circle indicate the association loci for branch angle that are identified in the MLM and MRMLM model, respectively. Three types of loci (QTL, LD block, and no block) are colored in red, blue and violet, respectively. The identified genes are listed on the outside of the middle and inner circles, and the genes of interest are colored in red.
Figure 4. Comparative analysis of haplotype block structure in the A-subgenome and C-subgenome of rapeseed. (A) A comparison of the number of haplotype blocks on 19 rapeseed chromosomes. (B) A comparison of the average size of haplotype blocks on 19 rapeseed chromosomes. (C) A comparison of the haplotype block coverage percentage on 19 rapeseed chromosomes. (D) A comparison of the haplotype block size range distributions on the A- and C-subgenome.
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).
Genome-Wide Association Analysis
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).
Figure 5. Manhattan plots of association analysis for branch angle using the Q+K model in six environments. (A) E1. (B) E2. (C) E3. (D) E4. (E) E5. (F) E6. The horizontal red line indicates the significance threshold [–log10 (p) = 3.0].
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).
Candidate Genes Identification
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−6 and 4.08 × 10−4, respectively, Figure 7A).
Figure 6. Significantly associated SNPs and corresponding candidate genes identified for branch angle in MLM. (A) Manhattan plot of the A6 chromosomal region around the candidate BnaA06g35880D (Bna.SGR3) in E2 and LD plot with a peak SNP Bn-A06-p24551529 and both flanking SNPs. (B) Manhattan plot of C8 chromosomal region around candidate BnaC08g25070D (Bna.SGR1) in E1 and LD plot with a peak SNP Bn-scaff_16770_1-p4296727 and both flanking SNPs. The SNP written in bold and marked with an asterisk indicates the peak SNP.
Figure 7. Haplotypes of associated SNPs among rapeseed natural variations. (A) Haplotype analysis of the peak SNP Bn-A06-p24551529 for Bna.SGR3 in MLM. (B) Haplotype analysis of the peak SNP Bn-scaff_16770_1-p4296727 for Bna.SGR1 in MLM. n denotes the number of genotypes belonging to each haplotype group, and the genotypes less than three are not shown. Statistical significance was determined with a t-test, different letters represent significant a difference at 5% level. The branch angle distribution of each haplotype group is displayed with a box plot.
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−3, 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 S4B.
Candidate Genes Validation
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).
Figure 8. Expression levels of four candidate genes between extremely large and small branch angle accessions. (A) Bna.TAC1. (B) Bna.SGR1. (C) Bna.SGR3. (D) Bna.SGR5. Error bars, s.d.; statistical significance was determined with a t-test, different letters above the bar represent significant a difference at 5% level.
The MLM that accounts for population structure (Q) and kinship (K), namely, the Q+K model, is a popular and powerful method used for GWASs, and it could reasonably resolve the spurious association between traits and markers caused by population structure (Yu et al., 2006; Bradbury et al., 2007). In this study, the same conclusion, i.e., the Q+K model being selected as the first-rank model, was drawn by comparing the different models (Figure S2). The Bonferroni correction is one of the typical multiple test corrections used for the threshold value of a significance test. However, it is often too conservative, such that many important loci may not pass the stringent criterion of significance test. A similar situation existed in the present study: when a GWAS was performed using the BLUP values in an MLM based on a modified Bonferroni threshold of p < 5.0 × 10−5 [−log10(p) = 4.3, 1/19,945], only one significant SNP on the A5 chromosome was discovered (Figure S1). Thus, to detect as many association signals as possible for use in further research, the significance threshold of association analysis in the MLM was dropped to a less stringent value (i.e., p < 1.0 × 10−3, −log10(p) = 3.0, Figure 5, Figure S1), which has been widely used in association mapping in rapeseed (Cai et al., 2014; Hatzig et al., 2015; Raman et al., 2015).
To prove the association results produced by the MLM and to take more advantage of the phenotypic and genotypic information obtained from an enormous amount of accessions and SNPs in this study, another model, a so-called MRMLM, was employed for a GWAS (Wang et al., 2016b). As a result, an additional 38 significance loci were identified using MRMLM, in which more than 55% of the loci overlapped part or most of the region with those obtained using MLM (Table S5), demonstrating the reliability of association analysis consequences and the practicality of combining MLM and MRMLM to improve the power and robustness of association analysis. Nevertheless, there are two prominent features in MRMLM compared to MLM. First, the MRMLM method treats marker effects as random. One advantage of this approach is that the model will shrink the effects of markers that are independent of target traits toward zero, leading to a maximum correlation between the observed and predicted phenotypic values (Goddard et al., 2009). Second, multiple test correction is not required due to the multi-locus and shrinkage nature. The MLM method is a single-locus analysis approach, in which only one marker is tested at a time. Thus, a Bonferroni correction for multiple tests is required to control the experimental error. In particular, when the number of markers is extremely large, the Bonferroni correction will be so stringent that many false-negative loci are introduced, which are significantly associated with traits in fact. Therefore, the MRMLM provides an alternative to GWASs in virtue of the power in QTL detection and the precision in locus effect estimation.
Two or more tightly related SNPs in strong LD were assigned to haplotype blocks, which were separated by recombination regions and defined the genetic variation across the genome. The block structure analysis will provide insight into the vital functional genomic regions in the course of selection and evolution (Qian et al., 2014). Therefore, genome-wide sweeping across the association panel using a high-throughput SNP chip was implemented for haplotype block structure analysis. One of the important conclusions was given based on our analysis: the large haplotype blocks were mostly distributed on the C-subgenome and were enriched around the centromere regions (Figure 3, Table 5), which is consistent with previous articles (Qian et al., 2014; Sun et al., 2016b). Here, we intend to give a plausible explanation of this phenomenon. First, the superficial reason is that the considerably stronger retention of LD leads to more long-range haplotype blocks on the C-subgenome (Qian et al., 2014). However, the ultimate contributor is the lack of genetic diversity in the C-subgenome. During Chinese B. napus breeding, the interspecific hybridization with B. rapa improves the genetic recombination and genetic diversity of the A-subgenome (Qian et al., 2006; Chen et al., 2007). However, the efforts to diversify the C-subgenome genetic component through B. napus × B. oleracea crosses were constrained due to cross-incompatibility (Bennett et al., 2008). Second, the transposon-rich regions often represent the recombination-poor (Gorelick, 2003). Recently, (Mason et al., 2016) observed the peak in transposable element density and the troughs in gene density in the centromere regions (Mason et al., 2016), which implies that lower frequency recombination events have occurred in centromere regions. Furthermore, considerably greater expansion of transposable elements was found in the C-subgenome of rapeseed (Chalhoub et al., 2014).
As depicted in previous reports, GWASs have been employed for branch angle research in rapeseed (Liu et al., 2016; Sun et al., 2016a). Hence, we compared the association consequences in this study with previous works. Unfortunately, the alignment results indicate that no identical SNP was found among them. But, encouragingly, there were nine SNPs detected by Sun et al. (2016a) that were within or proximate to loci detected in our study (Table S7). For example, two SNPs identified in the published literature on A7, Bn-A07-p15007983, and Bn-A07-p15505090, were within the locus on A7 with a peak SNP Bn-A07-p14798978 in the present paper. However, some loci detected by previous studies were still not discovered in our study, which may be affected by environmental factors, such as the location and year. In this study, the broad-sense heritability (H2) of branch angle was 76.06% (Table S2), hinting that environmental factors have a certain extent influence on the branch angle variation. Furthermore, the population size also has an important impact on the detection power of loci in GWASs, especially for rare alleles (Huang et al., 2012; Huang and Han, 2014; Li et al., 2016a). For branch angle, extensive variations are mainly caused by the cumulative effects of numerous polygenes with small effect (Sun et al., 2016a); the alleles with large effects may become rare, even extinct, in the gene pools of modern cultivars because of intensive artificial selection during domestication and modern breeding (Huang and Han, 2014). In addition, models based on a discrepant algorithm will depress the consistency of the results in GWASs. For example, approximately thirty percent of genes identified in one model could not be detected in another model in the present study.
Branch angle is regulated mainly by shoot gravitropism, which is a complex multistep process including the perception of gravity, transduction of the gravity signal into a biochemical signal, transport of the biochemical signal to a response site, and organ curvature (Sang et al., 2014). In the recent decade, many genes controlling the branch angle have been identified. LAZY1 plays a negative role in polar auxin transport and regulates the shoot gravitropism by which the rice tiller angle is controlled (Li et al., 2007). TAC1, a major gene involved in branch (tiller) angle and leaf angle control in plants, has been extensively studied (Yu et al., 2007; Ku et al., 2011; Dardick et al., 2013; Zhao et al., 2014). TAC1 and LAZY1 are both part of the same IGT gene family, but the gene structures of TAC1 and LAZY1 differ due to the presence of an additional EAR repression motif, which has the function of transcriptional repression, in the LAZY1 gene (Dardick et al., 2013). The difference in gene structure between TAC1 and LAZY1 may result in a discrepancy in molecular function; for example, there is no evidence that TAC1 plays a role in polar auxin transport thus far, leading to tac1 mutants with more vertical branch (tiller) angle in rice and Arabidopsis (Yu et al., 2007; Dardick et al., 2013). In the present study, the expression levels of TAC1 in large branch angle lines were significantly higher (Figure 8A), suggesting that the gene functions universally to promote the horizontal growth of branches. A series of Arabidopsis sgr mutants have been shown to exhibit disturbed shoot gravitropism. For example, loss-of-function of SGR5 in Arabidopsis and its ortholog in rice LPA1 displays less vertical branch (tiller) angle, in which the distribution of auxin was affected through regulation of auxin biosynthesis and transport (Cui et al., 2013; Wu et al., 2013). In the present study, the expression levels of SGR5 in small branch angle lines were significantly higher (Figure 8C), meaning that the gene is contrary to TAC1 in the function of branch angle regulation. The polarization of PIN-mediated auxin transport leads to changes in branch angle in the Arabidopsis and rice (Rakusova et al., 2011; Chen et al., 2012), demonstrating the central role of auxin and auxin transport in branch growth angle control.
The genes involved in branch angle control play an important role in modulating plant architecture mostly through auxin-dependent gravitropism. In this paper, many genes for branch angle, including TAC1, SGR1, SGR3, and SGR5, were first identified in rapeseed. Although extensive studies in branch angle genes have been done, the modulation basis underlying branch angle formation and maintenance is still elusive. Digby and Firn (1995) put forward the concept of gravitropic set-point angle (GSA), defined as the growth angle with respect to gravity. Recently, Roychoudhry et al. (2013) proposed a model for GSA maintenance based on the antagonistic interaction of auxin-dependent gravitropism and the anti-gravitropic offset component (AGO), the magnitude of which is regulated by gravity sensing cells in the shoot via Aux/IAA-TIR1-ARF-dependent auxin signaling. The model provided a conceptual framework for understanding GSA variation. However, the gravitropic-AGO model may not be a case of another class of growth angles, where the organ in question is not being actively maintained relative to gravity, such as the higher order secondary branches in peach trees (Dardick et al., 2013). In particular, in rice, there is no clear evidence indicating that the already cloned genes control the tiller angle through the gravity response, except for LAZY1 and LPA1 (Wu et al., 2016). It is suggested that there may be some other patterns regulating branch growth angle in plants. Therefore, more thorough research is required to elucidate the molecular mechanism underlying branch angle.
HL, ZL, and XW conceived and designed the study. BC, KX, and GG organized the implementation of field trials. LZ, FZ, HL, and TZ performed the phenotyping measurements. HL wrote the paper, JH, ZL, and XW modified the manuscript. All the authors have read and approved the publication of the manuscript.
This work was supported by the National Key Program for Research and Development (2016YFD0100202) and the Germplasm Resources Protection Project in China (NB2011-2130135-33).
Conflict of Interest Statement
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.
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/article/10.3389/fpls.2017.01054/full#supplementary-material
Figure S1. Manhattan plot of association analysis for branch angle using BLUP value in Q+K model. The horizontal blue line indicates the suggestive threshold [Bonferroni-corrected threshold −log10(p) = 4.3].The horizontal red line represents the significance threshold [−log10 (p) = 3.0].
Figure S2. Quantile–quantile plot of estimated −log10(p) from association analysis using six methods for branch angle. The black line represents the expected p-values with no association existed.
Table S1. Significant difference test in different branch region.
Table S2. Variance components and broad-sense heritability of branch angle.
Table S3. Haplotype block structure in rapeseed genome across 472 association population.
Table S4A. Significant associated SNPs for branch angle in MLM.
Table S4B. Significant associated SNPs for branch angle in MRMLM.
Table S5A. Significant associated loci and candidate genes for branch angle in MLM.
Table S5B. Significant associated loci and candidate genes for branch angle in MRMLM.
Table S6. Primers used for qRT-PCR.
Table S7. A comparison of the loci in the present paper and the SNPs in reference.
Bennett, R. A., Thiagarajah, M. R., King, J. R., and Rahman, M. H. (2008). Interspecific cross of Brassica oleracea var. alboglabra and B. napus: effects of growth condition and silique age on the efficiency of hybrid production, and inheritance of erucic acid in the self-pollinated backcross generation. Euphytica 164, 593–601. doi: 10.1007/s10681-008-9788-0
Bradbury, P. J., Zhang, Z., Kroon, D. E., Casstevens, T. M., Ramdoss, Y., and Buckler, E. S. (2007). TASSEL: software for association mapping of complex traits in diverse samples. Bioinformatics 23, 2633–2635. doi: 10.1093/bioinformatics/btm308
Cai, D., Xiao, Y., Yang, W., Ye, W., Wang, B., Younas, M., et al. (2014). Association mapping of six yield-related traits in rapeseed (Brassica napus L.). Theor. Appl. Genet. 127, 85–96. doi: 10.1007/s00122-013-2203-9
Chalhoub, B., Denoeud, F., Liu, S., Parkin, I. A., Tang, H., Wang, X., et al. (2014). Early allopolyploid evolution in the post-neolithic Brassica napus oilseed genome. Science 345, 950–953. doi: 10.1126/science.1253435
Chen, S., Nelson, M., Ghamkhar, K., Fu, T., and Cowling, W. (2007). Divergent patterns of allelic diversity from similar origins: the case of oilseed rape (Brassica napus L.) in China and Australia. Genome 51, 1–10. doi: 10.1139/g07-095
Chen, Y., Fan, X., Song, W., Zhang, Y., and Xu, G. (2012). Over-expression of OsPIN2 leads to increased tiller numbers, angle and shorter plant height through suppression of OsLAZY1. Plant Biotechnol. J. 10, 139–149. doi: 10.1111/j.1467-7652.2011.00637.x
Cui, D., Zhao, J., Jing, Y., Fan, M., Liu, J., Wang, Z., et al. (2013). The arabidopsis IDD14, IDD15, and IDD16 cooperatively regulate lateral organ morphogenesis and gravitropism by promoting auxin biosynthesis and transport. PLoS Genet. 9:e1003759. doi: 10.1371/journal.pgen.1003759
Dardick, C., Callahan, A., Horn, R., Ruiz, K. B., Zhebentyayeva, T., Hollender, C., et al. (2013). PpeTAC1 promotes the horizontal growth of branches in peach trees and is a member of a functionally conserved gene family found in diverse plants species. Plant J. 75, 618–630. doi: 10.1111/tpj.12234
Digby, J., and Firn, R. D. (1995). The gravitropic set-point angle (GSA): the identification of an important developmentally controlled variable governing plant architecture. Plant Cell Environ. 18, 1434–1440. doi: 10.1111/j.1365-3040.1995.tb00205.x
Dong, Z., Jiang, C., Chen, X., Zhang, T., Ding, L., Song, W., et al. (2013). Maize LAZY1 mediates shoot gravitropism and inflorescence development through regulating auxin transport, auxin signaling, and light response. Plant Physiol. 163, 1306–1322. doi: 10.1104/pp.113.227314
Fukaki, H., Wysocka-Diller, J., Kato, T., Fujisawa, H., Benfey, P. N., and Tasaka, M. (1998). Genetic evidence that the endodermis is essential for shoot gravitropism in Arabidopsis thaliana. Plant J. 14, 425–430. doi: 10.1046/j.1365-313X.1998.00137.x
Gabriel, S. B., Schaffner, S. F., Nguyen, H., Moore, J. M., Roy, J., Blumenstiel, B., et al. (2002). The structure of haplotype blocks in the human genome. Science 296, 2225–2229. doi: 10.1126/science.1069424
Gorelick, R. (2003). Transposable elements suppress recombination in all meiotic eukaryotes, including automictic ancient asexuals: a reply to Schon and Martens. J. Nat. Hist. 37, 903–909. doi: 10.1080/0022293021000007705
Hardy, O. J., and Vekemans, X. (2002). SPAGEDi: a versatile computer program to analyse spatial genetic structure at the individual or population levels. Mol. Ecol. Notes 2, 618–620. doi: 10.1046/j.1471-8286.2002.00305.x
Hatzig, S. V., Frisch, M., Breuer, F., Nesi, N., Ducournau, S., Wagner, M. H., et al. (2015). Genome-wide association mapping unravels the genetic control of seed germination and vigor in Brassica napus. Front. Plant Sci. 6:221. doi: 10.3389/fpls.2015.00221
Huang, X., Zhao, Y., Wei, X., Li, C., Wang, A., Zhao, Q., et al. (2012). Genome-wide association study of flowering time and grain yield traits in a worldwide collection of rice germplasm. Nat. Genet. 44, U32–U53. doi: 10.1038/ng.1018
Ku, L., Wei, X., Zhang, S., Zhang, J., Guo, S., and Chen, Y. (2011). Cloning and characterization of a putative TAC1 ortholog associated with leaf angle in maize (Zea mays L.). PLoS ONE 6:e20621. doi: 10.1371/journal.pone.0020621
Kuai, J., Sun, Y., Zuo, Q., Huang, H., Liao, Q., Wu, C., et al. (2015). The yield of mechanically harvested rapeseed (Brassica napus L.) can be increased by optimum plant density and row spacing. Sci Rep. 5:18835. doi: 10.1038/srep18835
Li, F., Chen, B., Xu, K., Gao, G., Yan, G., Qiao, J., et al. (2016a). A genome-wide association study of plant height and primary branch number in rapeseed (Brassica napus). Plant Sci. 242, 169–177. doi: 10.1016/j.plantsci.2015.05.012
Li, F., Chen, B., Xu, K., Wu, J., Song, W., Bancroft, I., et al. (2014). Genome-wide association study dissects the genetic architecture of seed weight and seed quality in rapeseed (Brassica napus L.). DNA Res. 21, 355–367. doi: 10.1093/dnares/dsu002
Li, L., Luo, Y., Chen, B., Xu, K., Zhang, F., Li, H., et al. (2016b). A genome-wide association study reveals new loci for resistance to clubroot disease in Brassica napus. Front. Plant Sci. 7:1483. doi: 10.3389/fpls.2016.01483
Liu, J., Wang, W., Mei, D., Wang, H., Fu, L., Liu, D., et al. (2016). Characterizing variation of branch angle and genome-wide association mapping in rapeseed (Brassica napus L.). Front. Plant Sci. 7:21. doi: 10.3389/fpls.2016.00021
Mason, A. S., Rousseau-Gueutin, M., Morice, J., Bayer, P. E., Besharat, N., Cousin, A., et al. (2016). Centromere locations in Brassica A and C genomes revealed through half-tetrad analysis. Genetics 202, 513–523. doi: 10.1534/genetics.115.183210
Merk, H. L., Yarnes, S. C., Van Deynze, A., Tong, N., Menda, N., Mueller, L. A., et al. (2012). Trait diversity and potential for selection indices based on variation among regionally adapted processing tomato germplasm. J. Am. Soc. Hortic. Sci. 137, 427–437.
Qian, L., Qian, W., and Snowdon, R. J. (2014). Sub-genomic selection patterns as a signature of breeding in the allopolyploid Brassica napus genome. BMC Genomics 15:1170. doi: 10.1186/1471-2164-15-1170
Qian, W., Meng, J., Li, M., Frauen, M., Sass, O., Noack, J., et al. (2006). Introgression of genomic components from Chinese Brassica rapa contributes to widening the genetic diversity in rapeseed (B. napus L.), with emphasis on the evolution of Chinese rapeseed. Theor. Appl. Genet. 113, 49–54. doi: 10.1007/s00122-006-0269-3
Rakusova, H., Gallego-Bartolome, J., Vanstraelen, M., Robert, H. S., Alabadi, D., Blazquez, M. A., et al. (2011). Polarization of PIN3-dependent auxin transport for hypocotyl gravitropic response in Arabidopsis thaliana. Plant J. 67, 817–826. doi: 10.1111/j.1365-313X.2011.04636.x
Raman, H., Raman, R., Coombes, N., Song, J., Prangnell, R., Bandaranayake, C., et al. (2015). Genome-wide association analyses reveal complex genetic architecture underlying natural variation for flowering time in canola. Plant Cell Environ. 39, 1228–1239. doi: 10.1111/pce.12644
Roychoudhry, S., Del Bianco, M., Kieffer, M., and Kepinski, S. (2013). Auxin controls gravitropic setpoint angle in higher plant lateral branches. Curr. Biol. 23, 1497–1504. doi: 10.1016/j.cub.2013.06.034
Sang, D., Chen, D., Liu, G., Liang, Y., Huang, L., Meng, X., et al. (2014). Strigolactones regulate rice tiller angle by attenuating shoot gravitropism through inhibiting auxin biosynthesis. Proc. Natl. Acad. Sci. U.S.A. 111, 11199–11204. doi: 10.1073/pnas.1411859111
Sun, C., Wang, B., Wang, X., Hu, K., Li, K., Li, Z., et al. (2016a). Genome-wide association study dissecting the genetic architecture underlying the branch angle trait in rapeseed (Brassica napus L.). Sci. Rep. 6:33673. doi: 10.1038/srep33673
Sun, C., Wang, B., Yan, L., Hu, K., Liu, S., Zhou, Y., et al. (2016b). Genome-wide association study provides insight into the genetic control of plant height in rapeseed (Brassica napus L.). Front Plant Sci. 7:1102. doi: 10.3389/fpls.2016.01102
Wang, N., Chen, B., Xu, K., Gao, G., Li, F., Qiao, J., et al. (2016a). Association mapping of flowering time QTLs and insight into their contributions to rapeseed growth habits. Front Plant Sci. 7:338. doi: 10.3359/fpls.2016.00338
Wang, S., Feng, J., Ren, W., Huang, B., Zhou, L., Wen, Y., et al. (2016b). Improving power and accuracy of genome-wide association studies via a multi-locus mixed linear model methodology. Sci. Rep. 6:19444. doi: 10.1038/srep19444
Wu, X., Tang, D., Li, M., Wang, K., and Cheng, Z. (2013). Loose Plant Architecture1, an INDETERMINATE DOMAIN protein involved in shoot gravitropism, regulates plant architecture in rice. Plant Physiol. 161, 317–329. doi: 10.1104/pp.112.208496
Yamauchi, Y., Fukaki, H., Fujisawa, H., and Tasaka, M. (1997). Mutations in the SGR4, SGR5 and SGR6 loci of Arabidopsis thaliana alter the shoot gravitropism. Plant Cell Physiol. 38, 530–535. doi: 10.1093/oxfordjournals.pcp.a029201
Yan, G., Lv, X., Gao, G., Li, F., Li, J., and Qiao, J., et al (2016). Identification and characterization of a glyoxalase I gene in a rapeseed cultivar with seed thermotolerance. Front. Plant Sci. 7:150. doi: 10.3389/fpls.2016.00150
Yu, J. M., Pressoir, G., Briggs, W. H., Bi, I. V., Yamasaki, M., Doebley, J. F., et al. (2006). A unified mixed-model method for association mapping that accounts for multiple levels of relatedness. Nat. Genet. 38, 203–208. doi: 10.1038/ng1702
Keywords: Brassica napus L., branch angle, plant architecture, association mapping, candidate-genes
Citation: Li H, Zhang L, Hu J, Zhang F, Chen B, Xu K, Gao G, Li H, Zhang T, Li Z and Wu X (2017) Genome-Wide Association Mapping Reveals the Genetic Control Underlying Branch Angle in Rapeseed (Brassica napus L.). Front. Plant Sci. 8:1054. doi: 10.3389/fpls.2017.01054
Received: 20 March 2017; Accepted: 31 May 2017;
Published: 19 June 2017.
Edited by:Maoteng Li, Huazhong University of Science and Technology, China
Reviewed by:Shengwu Hu, Northwest A&F University, China
Katarzyna Gacek, Plant Breeding and Acclimatization Institute, Poland
Copyright © 2017 Li, Zhang, Hu, Zhang, Chen, Xu, Gao, Li, Zhang, Li and Wu. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.