Skip to main content
Download PDF

Genome-wide association study and transcriptome analysis reveal key genes affecting root growth dynamics in rapeseed

Biotechnology for Biofuels volume 14, Article number: 178 (2021) Cite this article

Abstract

Background

In terms of global demand, rapeseed is the third-largest oilseed crop after soybeans and palm, which produces vegetable oil for human consumption and biofuel for industrial production. Roots are vital organs for plant to absorb water and attain mineral nutrients, thus they are of great importance to plant productivity. However, the genetic mechanisms regulating root development in rapeseed remain unclear. In the present study, seven root-related traits and shoot biomass traits in 280 Brassica napus accessions at five continuous vegetative stages were measured to establish the genetic basis of root growth in rapeseed.

Results

The persistent and stage-specific genetic mechanisms were revealed by root dynamic analysis. Sixteen persistent and 32 stage-specific quantitative trait loci (QTL) clusters were identified through genome-wide association study (GWAS). Root samples with contrasting (slow and fast) growth rates throughout the investigated stages and those with obvious stage-specific changes in growth rates were subjected to transcriptome analysis. A total of 367 differentially expressed genes (DEGs) with persistent differential expressions throughout root development were identified, and these DEGs were significantly enriched in GO terms, such as energy metabolism and response to biotic or abiotic stress. Totally, 485 stage-specific DEGs with different expressions at specific stage were identified, and these DEGs were enriched in GO terms, such as nitrogen metabolism. Four candidate genes were identified as key persistent genetic factors and eight as stage-specific ones by integrating GWAS, weighted gene co-expression network analysis (WGCNA), and differential expression analysis. These candidate genes were speculated to regulate root system development, and they were less than 100 kb away from peak SNPs of QTL clusters. The homologs of three genes (BnaA03g52990D, BnaA06g37280D, and BnaA09g07580D) out of 12 candidate genes have been reported to regulate root development in previous studies.

Conclusions

Sixteen QTL clusters and four candidate genes controlling persistently root development, and 32 QTL clusters and eight candidate genes stage-specifically regulating root growth in rapeseed were detected in this study. Our results provide new insights into the temporal genetic mechanisms of root growth by identifying key candidate QTL/genes in rapeseed.

Introduction

Rapeseed (Brassica napus L.; Brassicaceae), a globally cultivated crop, is not only one of the essential vegetable oil sources, but also an important emerging biodiesel and biofuel sources for industrial production [1]. Currently, biodiesel is mainly made from the monounsaturated fatty acids from vegetable oils [2]. Rapeseed oil has the highest percentage of monounsaturated fatty acids among the plant oils. Biodiesel has been manufactured primarily from rapeseed oil in Europe [3]. Rapestraw can be used to produce liquid biofuel, particularly ethanol, since it contains abundant lignocellulosic material [4]. Therefore, it is necessary to boost rapeseed biomass and yield so as to satisfy the increasing demand for edible oil and fuel worldwide.

The root system architecture (RSA) usually denotes the spatial configuration of complex assembly of the root system, and root shape plays key role in healthy plant growth, since root system penetrates the soil in search for water and nutrients [5]. Therefore, plants rely on the modulation of RSA in response to a changing soil environment to increase yield potential and yield stability. The genetic improvement of root architecture, such as increasing lateral root (LR) number, facilitates resource bioavailability in plants and increases crop yield and stress tolerance [6, 7]. To breed the crops with better RSA, a large number of studies have focused on variations in root architecture in many crops, such as rice, wheat, maize, soybean, and rapeseed [8,9,10,11,12]. Several studies have identified hundreds of root QTL in controlled environments or in the field [11, 13,14,15]. Besides, several of these QTL have also been reported to influence such traits as yield, water/nutrient uptake, and abiotic stress tolerance [13, 16,17,18,19,20,21].

Genome-wide association study (GWAS) has been successfully used for the identification of the polymorphism sites and/or genes related to complex traits including root traits in crops, such as rice, wheat, maize, and rapeseed [10, 14, 15, 22, 23]. The Brassica 60 K Illumina single-nucleotide polymorphism (SNP) array has facilitated the genetic improvement of different traits including flowering time, seed oil content, and phosphate-efficiency to obtain desirable alleles in B. napus [14, 24, 25]. Moreover, brassinosteroid signaling kinase 3 (BSK3) was confirmed to regulate root elongation at the low-nitrogen condition in Arabidopsis by GWAS [26]. Transcriptome analysis has become an effective technique for detecting candidate genes. Many crucial differentially expressed genes (DEGs) related to root development have been identified by RNA sequencing in rice, maize, and B. napus [27,28,29]. Weighted gene co-expression network analysis (WGCNA) has been usually used to analyze the relationship and network between different genes. Functional candidate genes related to root development were identified at different developmental stages in crops by WGCNA, including DcMYB113, which was reported to regulate anthocyanin transport in carrot root [30], and three hub genes (GRMZM2G477658, GRMZM2G15536, and GRMZM2G072121) played a possible role in maize root formation and growth through the division and/or elongation of cells [31]. Recently, the combination of GWAS, transcriptome sequencing, and/or WGCNA has been turned out to be a rapid and efficient approach to identifying crucial candidate genes regulating root development [9, 15, 32, 33]. For example, OsNal1 and OsJAZ1 located in the peak SNPs have been confirmed to facilitate the root development in rice [33].

Root growth is a continuous and complex process with temporal dynamics and spatial patterning. A previous study has defined seven root growth types in a B. napus recombinant inbred line (RIL) population derived from two rapeseed cultivars (Zhongshuang11 and NO. 73290) with contrasting root systems, and identified two types of QTL (persistent and stage-specific) by the analysis of root traits in rapeseed [11]. To further identify the genetic factors controlling the dynamic root growth, we examined five continuous stages during root development in 280 natural accessions of B. napus. Sixteen persistent and 32 stage-specific QTL clusters further confirmed the existence of the two types of QTL controlling root development. In addition, we performed a transcriptome analysis of samples of four root growth types with extremely contrasting growth rates during the investigated timepoints. A total of 12 crucial candidate genes involved in root growth were identified via combining GWAS, WGCNA, and differential expression analysis, some of which have been reported to be related to root development in previous studies.

Results

Phenotypic analysis of 280 B. napus accessions reveals genetic stability of root development

To examine dynamic growth patterns of roots during the vegetative stage, the hydroponic system was used for evaluating root-related traits and shoot biomass traits of 280 B. napus accessions which were sampled at 13 days after sowing (DAS) from the germination device and at 10 days after transplanting (10 DAT, equal to 16 DAS), three expanding leaves (3 EL), 5 EL, and 7 EL from the growth device with three biological replications for each sample, respectively (Additional file 1: Figure S1a–e). The statistics of the seven root-related traits (root fresh weight (RFW), root dry weight (RDW), primary root length (PRL), total root length (TRL), total root surface (TSA), total root volume (TRV), and total number of roots (TNR)), and two shoot biomass traits (shoot fresh weight (SFW) and shoot dry weight (SDW)) from each replication at the five sampling timepoints were listed in Additional file 2: Table S1, and the mean values of three replications are presented in Table 1. All the investigated traits showed a normal distribution or approximate normal distribution (Additional file 1: Figure S2). The coefficient of variation (CV) ranged from 15.21% to 25.66%, indicating considerable phenotypic variations for all the traits in the population (Table 1). All the traits showed a high broad-sense heritability (H2) at each developmental stage, ranging from 0.68 to 0.94. Furthermore, the H2 of all the traits was also high, ranging from 0.59 to 0.92 during the five developmental stages. For root traits, RFW, TRL, TSA, and TRV had heritability slightly higher than TNR and RDW (Table 1). We discovered significant correlations of each trait among all five stages with r2 ranging from 0.33 to 0.87 (P < 0.0001). In general, the highest root correlations were observed between two adjacent stages in spite of the gradually decreased correlation with the increased sampling interval, indicating that the effects of environment on these traits increased with plant development (Fig. 1). The PCA results of the traits suggested that component 1 (X axis, 47.9%) and component 2 (Y axis, 13.0%) explained the majority of genetic variation in this population (Fig. 2). With exception of PRL (Group 1), all the other traits examined at early stages (13 DAS and 10 DAT) were clustered into Group 3, whereas the traits recorded at late stages (3 EL, 5 EL, and 7 EL) were clustered into another group (Group 2). The separation of PRL and the other traits on the X axis indicated the substantial differences between PRL and the other traits. Traits captured at early stages (Group 3) and late stages (Group 2) were separated by the Y axis, but mapped to the same position on the X axis, suggesting the high correlations; however, a degree of specificity between the traits at early and late stages (Fig. 2). As shown in Additional file 2: Table S2, all the traits were significantly correlated with r2 ranging from 0.24 to 0.74 (P < 0.001) among the three biological replications. The results suggested that early root development traits were positively correlated with late root traits, thus suggesting that root development was a continuous process influenced by early genetic factors.

Table 1 Nine trait statistics of 280 accessions collected at five continuous stages
Full size table
Fig. 1
figure 1

Correlations of each captured trait at five stages. Forward slash represents the correlation, for example, 13DAS/10DAT represent the correlation in the traits between 13 DAS and 10 DAT, P < 0.0001

Full size image
Fig. 2
figure 2

Principal component analysis (PCA) of captured traits at five stages in association population. Three obvious groups are indicated with blue circles. PRL trait was clustered into group 1. Other traits captured at 13 DAS and 10 DAT were clustered into group 3. The traits captured at 3 EL, 5 EL, and 7 EL were clustered into group 2

Full size image

Persistent and stage-specific QTL clusters related to the root system are identified by GWAS

After filtering, a total of 23,542 SNPs with known physical position in the B. napus Darmor-bzh reference genome were selected for GWAS [34]. The distribution of the 23,542 SNP markers and LD decay on each chromosome were presented in Additional file 2: Table S3. Approximately 58.1% of the kinship coefficients between individual accessions were equal to zero, and 97.6% were less than 0.2, suggesting a weak kinship for most accessions in the natural population (Additional file 1: Figure S3).

The 1,107 significant trait–SNP associations were detected (− log10P > 4.37, − log101/23,542) using the mixed linear model (MLM) for three repetitions (Additional file 2: Table S4). The manhattan plots were drawn using the best linear unbiased prediction (BLUP) values of three repetitions for all the traits to give visual GWAS results at various stages (Fig. 3). We termed the SNPs with close proximity (within 1 Mb) and an LD r2 > 0.2 as one cluster, since these SNPs were identified as the same QTL [24]. As a result, a total of 683 identified significant trait–SNP associations with 134 significant SNPs markers, and 747 suggestive trait–SNP associations (3.5 <  − log10P ≤ 4.37) were integrated into 48 valid QTL clusters (Fig. 4, Additional file 2: Table S5). Of these 48 clusters, 21 QTL clusters contained multiple SNPs and 27 QTL clusters harboured single SNPs. The maximum genetic variation explained by these clusters ranged from 7.55% to 16.15%. We detected 28, 19, 4, 23, and 26 QTL clusters at the 13 DAS, 10 DAT, 3 EL, 5 EL, and 7 EL stages, respectively. Except 8 QTL clusters S1, S2, #1, S6, S8, S21, S26, and S18, all other clusters (40 out of 48) were detected at two or more stages. Noteworthy, two significant SNPs displaying the vast majority of trait–SNP associations for all the investigated traits except PRL on chromosome C8 were detected at multiple stages (Fig. 3; Additional file 2: Table S4). This suggested the existence of genetic factors controlling multiple root-related traits at various stages.

Fig. 3
figure 3

Manhattan plot of genome-wide association analysis for root and shoot-related traits at five stages. a–h Manhattan plot of genome-wide association analysis for SFW, SDW, RFW, RDW, TRL, TRV, TNR, and TSA, respectively. The different colours represent the trait-related SNPs at 13 DAS, 10 DAT, 3 EL, 5 EL, and 7 EL, respectively. The horizontal black lines indicate the significance threshold of GWAS (− log101/23,542 = 4.37). The x-axis shows the 19 chromosomes (A1–A10 and C1–C9) in B. napus. Each chromosome is scaled by the physical chromosome length

Full size image
Fig. 4
figure 4

Genetic landscapes of QTL clusters obtained in this study. QTL clusters refer to significant SNPs with close proximity (within 1 Mb) and an LD of r2 > 0.2. Position, physical range, number, largest genetic variation explanation of significant SNPs of each cluster. # represents clusters with multiple SNPs. S indicates the clusters harbouring single SNPs. ES-specific clusters: early stage specific clusters, including 13 DAS and 10 DAT; LS-specific clusters: late stage specific clusters, including 3 EL, 5 EL, and 7 EL; ES–LS common clusters: early and late stage clusters, including 13 DAS, 10 DAT, 3 EL, 5 EL, and 7 EL

Full size image

To reveal the genetic basis of root traits at the multiple vegetative stages, these QTL clusters were divided into three categories based on their identification stages: ES-specific clusters (early stage, 13 DAS and 10 DAT), LS-specific clusters (later stage, 3 EL, 5 EL, and 7 EL), and ES–LS common clusters. Sixteen out of the 48 QTL clusters constituted ES–LS common clusters, indicating the existence of the persistent QTL controlling root development. In addition, 16 ES-specific clusters and 16 LS-specific clusters revealed genetic mechanism in the root system at specific stages (Fig. 2; Additional file 2: Table S5). The major QTL identified in this study could be applied for improving root system architecture in rapeseed.

Transcriptome analysis reveals dynamic root development

Clustering analysis of the 280 accessions was performed to examine the similarity and diversity of their root growth patterns. At the same the growth stage, the traits (except PRL) exhibited significant correlations (P < 0.0001) with each other (from 0.44 to 0.97, P < 0.0001; Additional file 2: Table S6), suggesting developmental relevance among these root-related traits. SFW was considered as the trait reflecting the plant growth status. RFW showed higher correlations with SFW (0.70–0.79) than with other root traits (Additional file 2: Table S6). The traits investigated at 13 DAS and 10 DAT were from different growth devices (germination device and growth device), so growth rate (GR) from 13 DAS to 10 DAT was not shown in this study. The normalized GRs were calculated by RFW to present the root dynamic growth patterns, the heatmap showed that the 280 accessions fell into seven growth types (Types 1–7) (Fig. 5a). The 38 accessions (accounting for 13.57%) belonged to growth type 1 with their GRs below the average GR from 10 DAT to 7 EL, and at least one GR less than 80% of the average. The 48 accessions (17.14%) belonged to growth type 2 with their GRs greater than the average GR from 10 DAT to 7 EL, and at least one GR greater than 120% of the average. The majority of the accessions (64, 22.86%) belonged to type 3 with a relatively stable GR ranging from 80 to 120% of the average GR from 10 DAT to 7 EL. Type 4 possessed 48 accessions (17.14%) whose GRs were below average GR from 10 DAT to 3 EL or from 10 DAT to 5 EL, but were above average from 3 to 7 EL or from 5 to 7 EL. Contrastive, type 5 consisted of 50 accessions (17.86%) whose GRs were above average GR from 10 DAT to 3 EL or from 10 DAT to 5 EL, then below average from 3 to 7 EL or from 5 to 7 EL (Fig. 5a, b). Type 6 contained 17 accessions (6.07%) with its GR changing from fast to slow, and then to fast again during the investigated stages. On the contrary, the GRs of type 7 consisting of 15 accessions (5.36%) were subjected to the change pattern of first slow, and then fast, followed by slow. Obviously, the majority of accessions fell into type 1 (with consistent slow GR) and type 2 (with consistent fast GR), and type 3 (with stable GR), indicating that genes expressed at an early stage might control root growth at the late stage with prolonged effects. In addition, type 4 and type 5 displayed obvious stage-specific changes in GRs, suggesting the existence of genes functioning at a specific stage.

Fig. 5
figure 5

Clustering analysis of 280 accessions. a Hierarchical heatmap of 7 dynamic growth types of 280 accessions based on normalized GR values calculated by RFW. Red represents the GRs above the average, and green indicates the GRs below the average. b Diagram of 7 growth types with six accessions in each growth type. c Diagram of 4 transcriptome-sequenced growth types (type 1, 2, 4, and 5) with 5 accessions in each type

Full size image

The GRs of type 3 were close to average GR. The GRs of type 6 and type 7 exhibit two reverse changes during root development. Considering this, we excluded type 3, 6, and 7 in subsequent transcriptome analysis. We selected the four growth types (type 1, type 2, type 4, and type 5) for subsequent transcriptome analysis, because type 1 and type 2 had contrasting (slow and fast) GRs throughout root development stage, type 4 and type 5 exhibited opposite changing GRs at the specific stage. Five accessions from each growth type were sampled at the 10 DAT, 3 EL, 5 EL, and 7 EL stages, respectively (Fig. 5c), and were subjected to transcriptome analysis to reveal the temporal molecular mechanisms of root development. A total of more than 41 million clean reads were obtained from each library after adaptor trimming, of which, 73.47–91.59% clean reads were uniquely matched to B. napus reference genome (Additional file 2: Table S7). The qRT-PCR of 20 genes was performed in all the samples (Additional file 2: Table S14). The results of qRT-PCR were highly consistent with those of RNA-Seq data, suggesting the reliability of the RNA-Seq data (Additional file 1: Figure S4). The PCA of the RNA-Seq data indicated that all the four root growth types displayed an obvious separation between group 1 (10 DAT and 3 EL stage) and group 2 (5 EL and 7 EL stage) on component 1 (Additional file 1: Figure S5), indicating a change in gene expression from 3 to 5 EL during root development.

Persistent and stage-specific mechanisms underlying root development are revealed by transcriptome analysis

To explore the persistent genetic factors during root development, the VENN analysis of the DEGs from growth type 1 vs type 2 at the four root development stages was performed. A total of 367 DEGs were found to be overlapped within the four stages (Additional file 2: Table S8). A K-means clustering analysis of these persistent DEGs showed that the expressions of these genes were stabilized among the four stages, but they exhibited significant difference between type 1 and type 2 (Fig. 6a). In addition, 35 persistent DEGs encoded transcription factors, belonging to the families of bHLH NAC, MYB, MYB_related, MADS-box, and E2F/DP. The genes BnaAnng09810D, BnaA01g20660D, BnaC03g61210D homologous to the MADS-box family member ARABIDOPSIS NITRATE REGULATED 1, ANR1, the NAC transcription factor family PEROXIDASE 34, PER34, and bHLH25, respectively, which have been reported to participate in root development [35,36,37], were expressed higher at all the stages of type 1 than type 2. The GO enrichment analysis showed that the 367 DEGs were enriched in GO terms related to energy metabolism (including acetyl-CoA biosynthetic process from pyruvate and glycolytic process) and biotic or abiotic stress (such as response to oxidative stress, hydrogen peroxide catabolic process and cold acclimation) (Fig. 6b; Additional file 2: Table S9).

Fig. 6
figure 6

Information of persistent and stage-specific DEGs during root development in rapeseed. a K-mean clustering analysis of persistent DEGs in type 1 and type 2. The x axis indicates the samples at the four stages of type 1 and type 2. The y axis denotes the log2(FPKM). error bars represent the min and max data. b GO terms enriched with the persistent DEGs in type 1 and type 2. c K-mean clustering analysis of the stage-specific expression genes in type 4 and type 5. The x axis indicates the samples at the four stages of type 4 and type 5, and the y axis denotes the log2(FPKM). Error bars represent the min and max data. d GO terms enriched with the stage-specific expression genes in type 4 and type 5

Full size image

Meanwhile, 485 stage-specific DEGs were identified, which exhibited lower or higher expressions at type 4 early stages than at late stages, and displayed opposite expression patterns at type 5 corresponding stages (Fig. 6c; Additional file 2: Table S10). The GO enrichment analysis revealed that these genes were significantly enriched in GO terms, such as nitrate metabolism (including nitrate transport, response to nitrate, nitrate assimilation, and cellular response to nitrogen starvation), plant-type cell wall organization, and glucosinolate catabolic process. (Fig. 6d; Additional file 2: Table S11). In this study, a total of 16 stage-specific DEGs were highly expressed at 10 DAT of type 5, and they encoded multiple transcription factors, including trihelix family genes (PETAL LOSS and PTL) and MYB_related genes (CAPRICE and CPC). These transcription factors have been reported to affect root development in Arabidopsis [38, 39]. These results suggested that the biological processes, energy metabolism, and response to biotic or abiotic stress might influence the persistent root development, whereas the nitrate metabolite process might function at a specific stage during root growth. The root development was also regulated by several important transcription factors.

Crucial candidate genes are identified by integrating GWAS, WGCNA, and differential expression analysis

To investigate the gene regulatory network during root development, 26,039 DEGs from the four root growth types were used to identify co-expression gene modules by WGCNA. A total of 30 modules were identified in the dendrogram according to the correlations of genes (Fig. 7a), and the relationships between modules and samples were presented in Fig. 7b. The purple module was associated with all the stages of type 1, whereas the green module was associated with all the stages of type 2. The darkorange, darkturquoise, white, and darkred modules were significantly associated with 10 DAT, 3 EL, 5 EL and 7 EL of growth type 4, respectively, and the red, lightyellow, saddlebrown and darkgrey modules exhibited high correlations with 10 DAT, 3 EL, 5 EL and 7 EL of type 5, respectively. The heatmaps showed that the genes within one module were highly expressed in the samples highly correlated with the module (Additional file 2: Figure S6). GWAS results indicated that 2,461 genes were located in the haplotype blocks of the 48 QTL clusters (Additional file 2: Table S12). Considering a high correlation of WGCNA genes with each module (r2 > 0.85), 9 persistent and 13 stage-specific candidate genes each including 3 DEGs were screened from GWAS and WGCNA overlapped genes (Tables 2, 3). In addition, 6 GWAS and DEGs overlapped genes with correlations to the modules < 0.85 were also screened as stage-specific candidate genes (Table 3).

Fig. 7
figure 7

WGCNA of gene expression matrix in rapeseed. a Hierarchical clustering tree (dendrogram) of genes based on co-expression network analysis. b Module–sample association. Each row corresponds to a module labeled with a color identical to that in a, and each column corresponds to a sample

Full size image
Table 2 Key persistent candidate genes identified by GWAS, WGCNA, and differential expression analysis
Full size table
Table 3 Crucial stage-specific candidate genes identified by GWAS, WGCNA, and differential expression analysis
Full size table

Among the nine persistent candidate genes, two and seven genes with high correlations to the purple and green modules were highly expressed at all the stages of type 1 and type 2, respectively (Table 2). Four genes in the green module were located less than 100 kb away from the peak SNPs, including BnaA05g03210D, BnaC02g10480D, BnaA03g52990D, and BnaC08g39040D which were 35.9 Kb, 50.6 Kb, 69.3 Kb, and 6.3 Kb apart from the peak SNPs of #4, #9, S20, and S26 (Table 2). BnaA03g52990D encodes the GATA transcription factor, whose homolog influences root development by affecting auxin level and cell division in Arabidopsis [40]. Two genes BnaC02g10710D and BnaA05g03210D exhibited high correlation with BnaA03g52990D (Fig. 8a). ROOT INITIATION DEFECTIVE 1, RID1 (the homolog of BnaC02g10710D) has been reported to function in root apical meristem and root morphogenesis in Arabidopsis [41]. Furthermore, three persistent DEGs, BnaC02g10480D, BnaC02g10710D, and BnaA05g03210D, displayed high correlations to each other in the green module (Fig. 8a).

Fig. 8
figure 8

Networks of genes in green and red modules. a, b Correlation networks in green and red modules, respectively. Yellow colour in the network indicates the candidate genes overlapped by GWAS, WGCNA and differential expression analysis, and the purple colour in the network indicates the candidate genes overlapped by GWAS and WGCNA

Full size image

Of the 19 stage-specific candidate genes, BnaA03g47900D in the white module and seven genes BnaA03g42930D, BnaA09g07840D, BnaC01g22700D, BnaC02g14450D, BnaC03g26110D, BnaA06g37280D, and BnaC02g14330D in the red module were located less than 100 kb away from the peak SNPs (Table 3). Especially, BnaC01g22700D encoding ferredoxin-3 protein was 0.06 kb apart from the peak SNP of S19 QTL cluster. Our data indicated that PROTEIN BASIC PENTACYSTEINE5, BPC5 (the homologs of BnaA06g37280D) was located 28.6 kb away from the peak SNP of #11, and in previous study, BPC5 has been found to promote lateral root growth in Arabidopsis [42]. Our WGCNA, stage-specific DEG analysis, and GWAS results indicated that BnaA03g42930D, BnaA09g07840D, and BnaC01g22700D were detected and were highly expressed at 10 DAT of type 5 (Table 3), and that these three genes were highly correlated with another three closely linked genes BnaA08g24190D, BnaC02g14450D, and BnaC03g26110D (Fig. 8b).

In the present study, four candidate genes were screened as crucial persistent genetic factors and eight as stage-specific genetic factors with less than 100 kb physical distances from the peak SNPs in B. napus. Furthermore, homologs of three candidate genes (BnaA03g52990D, BnaA06g37280D, and BnaA09g07580D) have been reported to regulate root development in previous studies. The results showed that the method of screening candidate genes by combining GWAS, WGCNA, and differential expression analysis was effective.

Discussion

Two types (persistent and stage-specific) of temporal genetic factors controlling root development in B. napus

Recent advances in high-resolution imaging of root growth have indicated that the root system was determined by continuous spatial and temporal growth [11, 43,44,45]. Consistent with the previous report [11], our phenotypic correlation analyses and root growth dynamics study revealed two types (persistent and stage-specific) of temporal genetic factors controlling root development in B. napus. Furthermore, the persistent and stage-specific genetic factors were verified by our identified QTL clusters and DEGs (Fig. 2; Additional file 2: Table S5). Our dynamic QTL analysis results were in line with the previous reports on several dynamic traits at different developmental stages in Arabidopsis, barley, wheat, upland cotton, maize, and B. napus [11, 14, 46,47,48,49,50,51,52,53,54,55]. For example, 35 dynamic conditional QTL which can enhance the number of roots were detected at different root development stages in upland cotton, suggesting the dynamic development of roots [46].

Furthermore, the peak SNPs of 18 QTL clusters in this study were co-localized in the identical haplotype blocks of the 27 previously reported significant SNPs related to root traits at low or sufficient phosphorus conditions (Additional file 2: Table S13) [14]. Clusters S10, S20, #15, and S25 were also co-localized with previously identified QTL (qcA09-1, qcC02-2, qcC02-2, and uqPRLC06) related to root surface area (RSA) trait, respectively [11, 14]. Our results provided useful QTL and the major QTL can be used for marker-assisted selection of root traits in rapeseed.

Possible regulatory pathways of persistent and stage-specific genetic factors related to root development

Root growth, as a complex process, is determined by the interaction of many genes. Some genes play a persistent role during root development, whereas others function at a specific stage. In this study, we identified 367 persistent DEGs from growth type 1 vs type 2 controlling root development in rapeseed. Three persistent DEGs enriched in acetyl-CoA biosynthetic process were homologs of PYRUVATE DEHYDROGENASE E1, ALPHA in A. thaliana affecting polar auxin transport during root development [56]. Oxidative stress response is a general response of living organisms to biotic or abiotic stress [57]. Ten out of the 14 persistent DEGs (Additional file 2: Table S8) enriched in oxidative stress encoded 10 proteins PEROXIDASE 34 (PRX34), CYP709B3, CYP87A3, CYP78A6, PEROXIDASE 3 (PER3), PER34, PER44, PER71, CATALASE-2 (CAT2), and ALPHA-DIOXYGENASE 1 (DOX1), and these genes have been reported to act as regulators in root development in Arabidopsis [58,59,60,61,62]. For example, PRX34 mediated H2O2 generation and increased Ca2+ flux from the cytosol of Atmpk6 root cells to inhibit root elongation [63]. CYP709B subfamily was involved in cytokinin metabolism and signaling in roots [58]. These results suggested that biological processes, such as energy metabolism and biotic or abiotic stress response, especially oxidative stress response might act as the major molecule mechanisms influencing persistent root development.

NO3− and nitrate metabolites can serve as regulatory signals to control root system architecture [64]. Three stage-specific DEGs (Additional file 2: Table S10) homologous to NRT1/NPF6.3 not only regulated auxin biosynthesis to promote LR primordia emergence, but also repressed LR development by promoting auxin transport at low nitrate in Arabidopsis [65]. High-affinity nitrate transporter 2.1 (NRT2.1) homologous to our four stage-specific DEGs has been reported to play an essential role in root nitrate uptake (Additional file 2: Table S10) [66]. The 485 stage-specific DEGs were found to be enriched in GO terms, such as nitrate transport, response to nitrate, nitrate assimilation, and glucosinolate catabolic process (Fig. 6; Additional file 2: Table S11). Furthermore, some previous studies have reported that the glucosinolate accumulation can restrain root growth and development [67,68,69,70,71]. Defense metabolite Allyl-glucosinolates (allyl-GSL) have been reported to affect Arabidopsis root development through three different catabolic products [72]. AtTGG4 and AtTGG5 homologous to two genes, BnaA08g01990D and BnaC06g08840D, enriching in the glucosinolate catabolic process have been reported to regulate root growth and play a part in flood tolerance in Arabidopsis [73]. The above results suggested that nitrate metabolism process and glucosinolate catabolic process might mainly regulate the stage-specific root development.

Efficient discovery of candidate genes by combining GWAS, WGCNA, and differential expression analysis

Combination of GWAS, WGCNA, and differential expression analysis has been reported as an efficient way to acquire crucial genes in maize, rice, soybean, carrot, and other crops [9, 15, 32, 33]. We identified four persistent and eight stage-specific crucial candidate genes related to root development by integrating GWAS, WGCNA and differential expression analysis in rapeseed.

Four crucial persistent candidate genes BnaA03g52990D, BnaA05g03210D, BnaC02g10480D, BnaC08g39040D in the green module displayed high correlations to each other (Fig. 8), and two genes, BnaA03g52990D and BnaC02g10710D in green module were homologous to ATGATA3 and ATRID1 which have been reported to function in root development in Arabidopsis [40, 41]. Our identified ATRID1 had similar function with SRD2 which affected LR morphogenesis by reducing the level of auxin efflux facilitator (PIN) in A. thaliana [41, 74]. Furthermore, the homologs of other genes in the green module, PROTEIN PHOSPHATASE 2C (AIP1), REPLICATION PROTEIN A SUBUNIT B (RPA1B), HISTIDINE KINASE 3 (AHK3), POLYADENYLATE-BINDING PROTEIN 2 (PAB2), have also been reported to affect root development by regulating phytohormone or promoting cell elongation [75,76,77,78]. All these results indicated that the crucial persistent candidate genes in the green module might have similar functions during root development.

Seven out of eight stage-specific crucial candidate genes were in the red module, of which BnaA06g37280D and BnaA09g07580D were homologous to BPC5 and RALFL34 reported to promote LR development by inhibiting the abscisic acid insensitive 4 expression and activating PIN1 level in Arabidopsis [42, 79]. Moreover, the red module included several function-known genes involved in root development, such as BnaA08g06170D and BnaC08g06550D which were homologous to AtSMAP1 reported to modulate root development by interacting with 2,4-Dichlorophenoxyacetic acid [80]. The results above further demonstrated that these seven candidate genes played significant roles in root growth.

The candidate genes and dynamic QTL identified in this study can serve as exploitable resources to broaden our research on molecular mechanism of root development. More studies are needed to further analyze these candidate genes and validate their functions.

Conclusions

Rapeseed provides not only edible vegetable oil for human consumption, but also an important source for biofuel production. To construct excellent root system by genetic improvement is conducive to improve rapeseed productivity. The seven dynamic patterns of root growth rates and 16 persistent and 32 stage-specific quantitative trait loci (QTL) clusters which were obtained by GWAS supported the existence of two types of QTL (persistent and stage-specific) controlling root growth at specific or multiple developmental stages, respectively. Total of 367 identified persistent DEGs were enriched in energy metabolism and biotic or abiotic stress. Whereas 485 stage-specific DEGs were enriched in nitrogen metabolism. By integrating GWAS, WGCNA, and differential expression analysis, we identified four candidate genes as crucial persistent genetic factors and eight as stage-specific genes. Among these, three candidate genes (BnaA03g52990D, BnaA06g37280D, and BnaA09g07580D) had been reported to regulate root development in previous studies, supporting the validity of this method to obtain candidate genes. Our results provide new insights and useful candidate QTLs/genes into the temporal genetic mechanisms of root growth in rapeseed.

Methods

Plant materials and growth conditions

The natural population used in this study consisted of 280 B. napus lines, including 156 semi-winter accessions, 86 spring accessions and 38 winter accessions. A total of the 280 rapeseed germplasm accessions were collected, including 222 accessions from the Yangtze River of China, 23 from northwestern China, 16 from Europe, 14 from Australia, and 5 from other places or unknown origins. All the accessions were strictly self-crossed.

The previously reported hydroponic system was used for the root-related trait evaluation of the 280 B. napus accessions [11]. Briefly, uniform and stout rapeseed seeds were placed on the medical gauze of the germination device for 2 days in the dark, then they grew in the light (180 μmol photons m−2 s−1) for 4 d in a greenhouse (60–80% relative humidity) under 16/8 h day/night cycles at 24 °C (Additional file 1: Figure S1f, g). A quarter of modified Hoagland’s nutrient solution was filled into germination device to retain moisture and provide nutrients for seed germination [81]. Six days after sowing, uniform seedlings were transferred to the growth device containing 1/4 Hoagland’s solution. The 1/4 solution was replaced with a 1/2 solution, and then with a 100% solution once a week until harvesting.

Phenotypic evaluation of association panel

The accessions from the natural population were completely randomly grown and evaluated with three replications. In each replication, three uniform plants per accession were collected from the germination device at 13 days after sowing (DAS). At 6 DAS, 24 plants of an accession were transplanted to one growth device (Additional file 1: Figure S1h). Then three plants per accession were sampled from the growth device at four timepoints, namely, 10 days after transplanting (10 DAT, equal to 16 DAS), three expanding leaves (3 EL), 5 EL, and 7 EL, respectively. In total, 12,600 plants (280 accessions × 3 replicates × 3 plants × 5 timepoints) were sampled. Once the plants were sampled, shoot fresh weight (SFW), root fresh weight (RFW), and primary root length (PRL) were measured manually. The intact roots in a transparent box full of water were scanned with the root scanner (EPSON, 11000XL). The obtained high-resolution root images were analyzed using WinRHIZO-Pro software (Regent Instruments, QC, Canada) to determine total root length (TRL), total root surface (TSA), total root volume (TRV), and total number of roots (TNR). Subsequently, shoot and root samples were dehydrated at 65 °C for a week to determine shoot dry weight (SDW) and root dry weight (RDW).

Data analysis

The variance and correlation analyses of the investigated traits were performed using the software SAS 9.2. The broad-sense heritability was calculated using the formula reported by Liu et al. [24]. The principal component analysis (PCA) of all the investigated traits were conducted by the software SAS 9.2. The first step of PCA was to obtain the correlation matrix between different traits, then the dimensionality reduction was performed to obtain eight principal components, and PC1 and PC2 were plotted by R. According to previous reported method [11], the growth rate (GR) of accession sample was calculated as the root fresh weight (RFW) value at the late stage minus that at the early stage, and then divided by the growing days. GR of an accession was normalized according to the following formula. Normalized GR = (GRg − GRp)/GRp In the formula, GRg represented the GR of a genotype, and GRp was the average GR of the population of 280 accessions which were clustered in terms of the normalized GR using MeV_4_9_0 software (http://mev.ro/en/).

Population structure, relative kinship, and association analysis

The Brassica 60 K Illumina® Infinium consortium SNP array [82] (http://www.illumina.com/technology/beadarray-technology/infinium-hd-assay.html) was used for accessions genotype. SNP data were analyzed using Illumina BeadStudio genotyping software (http://www.illumina.com/) with parameters set as a missing rate ≤ 0.2, heterozygous rate ≤ 0.2, and minor allele frequency (MAF) > 0.05. BLAST was performed to search the probe sequences of these SNPs against the B.napus Darmor-bzh reference genome [34] with an threshold of e−10. SNPs with merely one matched position in reference genome were used for further analysis. The population structure and relative kinship of the 280 B. napus accessions were analyzed using STRUCTURE v. 2.3.4 and SPAGeDi software, respectively [83]. The linkage disequilibrium (LD) decay between all SNPs was assessed by TASSEL 4.0 [84]. The trait–SNP association was analyzed using mixed linear model (MLM) for both the single repetition and the BLUP [85]. Marker haplotypes at each associated locus were identified using the four-gamete rule with Haploview software [86].

Transcriptome sequencing and analysis

Five accessions from each of the four growth types (type 1, type 2, type 4, and type 5, detailed information presented in “Result” section) with contrasting GRs were selected and replanted for further transcriptome analysis. Total roots of three plants for each accession were sampled at four timepoints (10 DAT, 3 EL, 5 EL, and 7 EL) with two biological replications. Samples of the five accessions within one growth type at each sampling timepoint with the same weight were mixed as a single sample. A total of 32 obtained samples were fully mixed for total RNA extraction with the IRIzol reagent (Invitrogen, USA). Sequencing library construction and Illumina sequencing were performed by the Oebiotech Company in Shanghai, China using an Illumina HiSeq™ 2500 platform. The raw reads with 150 paired-end base pair (bp) were filtered and aligned as previously reported [25]. The raw data were submitted into database of the National Center for Biotechnology Information Sequence Read Archive (SRA; http://www.ncbi.nlm.nih.gov/sra) (Accession No. PRJNA714285).

The clean reads were mapped to the B.napus Darmor-bzh reference genome [34] (http://www.genoscope.cns.fr/brassicanapus/data/) using Hisat2. The gene expression levels were expressed as FPKM (fragments per kilobase per million reads) value. The PCA of the gene expressions was performed using the PCAtools package in R. The WGCNA was conducted using the WGCNA package in R [87]. P ≤ 0.05 for the false discovery rate (FDR) and |log2ratio|≥ 1 were used as criteria to identify DEGs with the DESeq package in R. The k-mean clustering was performed by MeV_4_9_0 software. Gene ontology (GO) enrichment analysis was performed using the ClusterProfiler package in R.

Real-time reverse transcription PCR

Quantitative real-time PCR (qRT-PCR) of 20 genes randomly selected from the DEGs was performed to verify the accuracy of RNA-seq data. The primer sequences were presented in Additional file 2: Table S14. The SYBR qPCR Master Mix (Vazyme) was used for qRT-PCR analysis with the CFX96 (BIO-RAD). Three technical replications were performed for each sample. The B. napus ACTIN2 was used as an internal control to compute the relative expression of target genes by the 2−ΔΔCT method.

Availability of data and materials

The raw sequence data have been deposited in the National Center for Biotechnology Information Sequence Read Archive (http://www.ncbi.nlm.nih.gov/sra/) under Accession number PRJNA714285. All other relevant data during this study are included in the manuscript and additional files.

Abbreviations

13 DAS:

13 Days after sowing

10 DAT:

10 Days after transferring

3 EL:

Three expending leaves

5 EL:

Five expending leaves

7 EL:

Seven expending leaves

PRL:

Primary root length

RDW:

Root dry weight

RFW:

Root fresh weight

SDW:

Shoot dry weight

SFW:

Shoot fresh weight

TNR:

Total root number

TRL:

Total root length

TRV:

Total root volume

TSA:

Total root surface area

LR:

Lateral root

GR:

Growth rate

GWAS:

Genome-wide association study

SNP:

Single nucleotide polymorphism

MAF:

Minor allele frequency

DEGs:

Differentially expressed genes

QTL:

Quantitative trait loci

qRT-PCR:

Quantitative real-time PCR

RNA-seq:

RNA sequencing

FPKM:

Fragments per kilobase of exon per million reads mapped

CV:

Coefficient of variation

H2 :

Broad-sense heritability

BLUP:

The best linear unbiased prediction

References

  1. Meyer M. Rapeseed oil fuel—the crisis-proof home-made eco-fuel. Agrarforschung. 2009;16(7):262–7.

    Google Scholar 

  2. Emiliani D, Pistocchi M. Biodiesel production from vegetable oils. Ind Sacc Ital. 2006;99:91–8.

    CAS  Google Scholar 

  3. Pullen J, Saeed K. Investigation of the factors affecting the progress of base-catalyzed transesterification of rapeseed oil to biodiesel FAME. Fuel Process Technol. 2015;130:127–35.

    Article  CAS  Google Scholar 

  4. Wang X, Zheng M, Liu H, Zhang L, Hua W. Fine-mapping and transcriptome analysis of a candidate gene controlling plant height in Brassica napus L. Biotechnol Biofuels. 2020;13(42).

  5. Giehl RF, Gruber BD, von Wiren N. It’s time to make changes: modulation of root system architecture by nutrient signals. J Exp Bot. 2014;65(3):769–78.

    Article  CAS  PubMed  Google Scholar 

  6. Uga Y, Sugimoto K, Ogawa S, Rane J, Ishitani M, Hara N, et al. Control of root system architecture by DEEPER ROOTING 1 increases rice yield under drought conditions. Nat Genet. 2013;45(9):1097–102.

    Article  CAS  PubMed  Google Scholar 

  7. Chen YS, Lo SF, Sun PK, Lu CA, Ho TH, Yu SM. A late embryogenesis abundant protein HVA1 regulated by an inducible promoter enhances root growth and abiotic stress tolerance in rice without yield penalty. Plant Biotechnol J. 2015;13(1):105–16.

    Article  CAS  PubMed  Google Scholar 

  8. Burton AL, Johnson JM, Foerster JM, Hirsch CN, Buell CR, Hanlon MT, et al. QTL mapping and phenotypic variation for root architectural traits in maize (Zea mays L.). Theor Appl Genet. 2014;127(11):2293–311.

    Article  PubMed  Google Scholar 

  9. Guo J, Li C, Zhang X, Li Y, Zhang D, Shi Y, et al. Transcriptome and GWAS analyses reveal candidate gene for seminal root length of maize seedlings under drought stress. Plant Sci. 2020;292:110380.

    Article  CAS  PubMed  Google Scholar 

  10. Beyer S, Daba S, Tyagi P, Bockelman H, Brown-Guedira G, Mohammadi M. Loci and candidate genes controlling root traits in wheat seedlings—a wheat root GWAS. Funct Integr Genomics. 2019;19(1):91–107.

    Article  CAS  PubMed  Google Scholar 

  11. Wang J, Kuang L, Wang X, Liu G, Dun X, Wang H. Temporal genetic patterns of root growth in Brassica napus L. revealed by a low-cost, high-efficiency hydroponic system. Theor Appl Genet. 2019;132(8):2309–23.

    Article  CAS  PubMed  Google Scholar 

  12. Song L, Prince S, Valliyodan B, Joshi T, Maldonado dos Santos JV, Wang J, et al. Genome-wide transcriptome analysis of soybean primary root under varying water-deficit conditions. BMC Genomics. 2016;17:57.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  13. Song W, Wang B, Hauck AL, Dong X, Li J, Lai J. Genetic dissection of maize seedling root system architecture traits using an ultra-high density bin-map and a recombinant inbred line population. J Integr Plant Biol. 2016;58(3):266–79.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  14. Wang X, Chen Y, Thomas CL, Ding G, Xu P, Shi D, et al. Genetic variants associated with the root system architecture of oilseed rape (Brassica napus L.) under contrasting phosphate supply. DNA Res. 2017;24(4):407–17.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  15. Wang H, Wei J, Li P, Wang Y, Ge Z, Qian J, et al. Integrating GWAS and gene expression analysis identifies candidate genes for root morphology traits in maize at the seedling stage. Genes. 2019;10(10):773.

    Article  PubMed Central  CAS  Google Scholar 

  16. Hua Y, Zhang D, Zhou T, He M, Ding G, Shi L, et al. Transcriptomics-assisted quantitative trait locus fine mapping for the rapid identification of a nodulin 26-like intrinsic protein gene regulating boron efficiency in allotetraploid rapeseed. Plant Cell Environ. 2016;39(7):1601–18.

    Article  CAS  PubMed  Google Scholar 

  17. Bernardino KC, Pastina MM, Menezes CB, de Sousa SM, Maciel LS, Carvalho G, et al. The genetic architecture of phosphorus efficiency in sorghum involves pleiotropic QTL for root morphology and grain yield under low phosphorus availability in the soil. BMC Plant Biol. 2019;19(1):87.

    Article  PubMed  PubMed Central  Google Scholar 

  18. Fan X, Zhang W, Zhang N, Chen M, Zheng S, Zhao C, et al. Identification of QTL regions for seedling root traits and their effect on nitrogen use efficiency in wheat (Triticum aestivum L.). Theor Appl Genet. 2018;131(12):2677–98.

    Article  CAS  PubMed  Google Scholar 

  19. Gong X, McDonald G. QTL mapping of root traits in phosphorus-deficient soils reveals important genomic regions for improving NDVI and grain yield in barley. Theor Appl Genet. 2017;130(9):1885–902.

    Article  CAS  PubMed  Google Scholar 

  20. Islam A, Zhang Y, Anis G, Rani MH, Anley W, Yang Q, et al. Fine mapping and candidate gene analysis of qRN5a, a novel QTL promoting root number in rice under low potassium. Theor Appl Genet. 2021;134(1):213–27.

    Article  CAS  PubMed  Google Scholar 

  21. Schierholt A, Tietz T, Bienert GP, Gertz A, Miersch S, Becker HC. Root system size response of bzh semi-dwarf oilseed rape hybrids to different nitrogen levels in the field. Ann Bot. 2019;124(6):891–901.

    Article  CAS  PubMed  Google Scholar 

  22. Sanchez DL, Liu S, Ibrahim R, Blanco M, Lubberstedt T. Genome-wide association studies of doubled haploid exotic introgression lines for root system architecture traits in maize (Zea mays L.). Plant Sci. 2018;268:30–8.

    Article  CAS  PubMed  Google Scholar 

  23. Zhang H, San ML, Jang SG, Lee JH, Kim NE, Lee AR, et al. Genome-wide association study of root system development at seedling stage in rice. Genes. 2020;11(12):1395.

    Article  CAS  PubMed Central  Google Scholar 

  24. Liu S, Fan C, Li J, Cai G, Zhou Y. A genome-wide association study reveals novel elite allelic variations in seed oil content of Brassica napus. Theor Appl Genet. 2016;1(6):1203–15.

    Article  CAS  Google Scholar 

  25. Liu S, Huang H, Yi X, Zhang Y, Yang Q, Zhang C, et al. Dissection of genetic architecture for glucosinolate accumulations in leaves and seeds of Brassica napus by genome-wide association study. Plant Biotechnol J. 2020;18(6):1472–84.

    Article  CAS  PubMed  Google Scholar 

  26. Jia Z, Giehl RFH. Natural variation of BSK3 tunes brassinosteroid signaling to regulate root foraging under low nitrogen. Nat Commun. 2019;10(1):2378.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  27. Dun X, Tao Z, Wang J, Wang X, Liu G, Wang H. Comparative transcriptome analysis of primary roots of Brassica napus seedlings with extremely different primary root lengths using RNA sequencing. Front Plant Sci. 2016;7:1238.

    Article  PubMed  PubMed Central  Google Scholar 

  28. Xu Y, Zou J. RNA-Seq transcriptome analysis of rice primary roots reveals the role of flavonoids in regulating the rice primary root growth. Genes. 2019;10(3):213.

    Article  CAS  PubMed Central  Google Scholar 

  29. Hey S, Baldauf J, Opitz N, Lithio A, Pasha A, Provart N, et al. Complexity and specificity of the maize (Zea mays L.) root hair transcriptome. J Exp Bot. 2017;68(9):2175–85.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  30. Xu ZS, Yang QQ, Feng K, Yu X, Xiong AS. DcMYB113, a root-specific R2R3-MYB, conditions anthocyanin biosynthesis and modification in carrot. Plant Biotechnol J. 2020;18(7):1585–97.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  31. Hwang SG, Kim KH, Lee BM, Moon JC. Transcriptome analysis for identifying possible gene regulations during maize root emergence and formation at the initial growth stage. Genes Genomics. 2018;40(7):755–66.

    Article  CAS  PubMed  Google Scholar 

  32. Que F, Wang GL, Li T, Wang YH, Xu ZS, Xiong AS. Genome-wide identification, expansion, and evolution analysis of homeobox genes and their expression profiles during root development in carrot. Funct Integr Genomics. 2018;18(6):685–700.

    Article  CAS  PubMed  Google Scholar 

  33. Li X, Guo Z, Lv Y, Cen X, Ding X, Wu H, et al. Genetic control of the root system in rice under normal and drought stress conditions by genome-wide association study. PLoS Genet. 2017;13(7):e1006889.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  34. Chalhoub B, Denoeud F, Liu S, Parkin IA, Tang H, Wang X, et al. Plant genetics. Early allopolyploid evolution in the post-Neolithic Brassica napus oilseed genome. Science. 2014;345(6199):950–3.

    Article  CAS  PubMed  Google Scholar 

  35. Gaudinier A, Rodriguez-Medina J, Zhang L, Olson A, Liseron-Monfils C, Bagman AM, et al. Transcriptional regulation of nitrogen-associated metabolism and growth. Nature. 2018;563(7730):259–64.

    Article  CAS  PubMed  Google Scholar 

  36. Sun CH, Yu JQ, Wen LZ, Guo YH, Sun X, Hao YJ, et al. Chrysanthemum MADS-box transcription factor CmANR1 modulates lateral root development via homo-/heterodimerization to influence auxin accumulation in Arabidopsis. Plant Sci. 2018;266:27–36.

    Article  CAS  PubMed  Google Scholar 

  37. Chutia R, Abel S, Ziegler J. Iron and phosphate deficiency regulators concertedly control coumarin profiles in Arabidopsis thaliana roots during iron, phosphate, and combined deficiencies. Front Plant Sci. 2019;10:113.

    Article  PubMed  PubMed Central  Google Scholar 

  38. Zhang J, Eswaran G, Alonso-Serra J, Kucukoglu M, Xiang J, Yang W, et al. Transcriptional regulatory framework for vascular cambium development in Arabidopsis roots. Nat Plants. 2019;5(10):1033–42.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  39. Tominaga-Wada R, Wada T. CPC-ETC1 chimeric protein localization data in Arabidopsis root epidermis. Data Brief. 2018;18:1773–6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  40. Jiang K, Yung V, Chiba T, Feldman LJ. Longitudinal patterning in roots: a GATA2-auxin interaction underlies and maintains the root transition domain. Planta. 2018;247(4):831–43.

    Article  CAS  PubMed  Google Scholar 

  41. Ohtani M, Demura T, Sugiyama M. Arabidopsis root initiation defective1, a DEAH-box RNA helicase involved in pre-mRNA splicing, is essential for plant development. Plant Cell. 2013;25(6):2056–69.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  42. Mu Y, Zou M, Sun X, He B, Xu X, Liu Y, et al. BASIC PENTACYSTEINE proteins repress ABSCISIC ACID INSENSITIVE4 expression via direct recruitment of the polycomb-repressive complex 2 in Arabidopsis root development. Plant Cell Environ. 2017;58(3):607–21.

    CAS  Google Scholar 

  43. Pulwicki J, Hobill D. The dynamics of root growth: a geometric model. Bull Math Biol. 2017;79(8):1820–45.

    Article  PubMed  Google Scholar 

  44. Thomas CL, Graham NS, Hayden R, Meacham MC, Neugebauer K, Nightingale M, et al. High-throughput phenotyping (HTP) identifies seedling root traits linked to variation in seed yield and nutrient capture in field-grown oilseed rape (Brassica napus L.). Ann Bot. 2016;118(4):655–65.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  45. Thomas CL, Alcock TD, Graham NS, Hayden R, Matterson S, Wilson L, et al. Root morphology and seed and leaf ionomic traits in a Brassica napus L. diversity panel show wide phenotypic variation and are characteristic of crop habit. BMC Plant Biol. 2016;16(1):214.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  46. Liang Q, Li P, Hu C, Hua H, Li Z, Rong Y, et al. Dynamic QTL and epistasis analysis on seedling root traits in upland cotton. J Genet. 2014;93(1):63–78.

    Article  CAS  PubMed  Google Scholar 

  47. Su Y, Li M, Guo L. Different effects of phospholipase Dζ2 and non-specific phospholipase C4 on lipid remodeling and root hair growth in Arabidopsis response to phosphate deficiency. Plant J. 2018;94(2):315–26.

    Article  CAS  PubMed  Google Scholar 

  48. Mohler V, Stadlmeier M. Dynamic QTL for adult plant resistance to powdery mildew in common wheat (Triticum aestivum L.). J Appl Genet. 2019;60(3–4):291–300.

    Article  CAS  PubMed  Google Scholar 

  49. Wang X, Wang H, Long Y, Liu L, Zhao Y, Tian J, et al. Dynamic and comparative QTL analysis for plant height in different developmental stages of Brassica napus L. Theor Appl Genet. 2015;128(6):1175–92.

    Article  PubMed  Google Scholar 

  50. Su J, Yang X, Zhang F, Wu S, Xiong S, Shi L, et al. Dynamic and epistatic QTL mapping reveals the complex genetic architecture of waterlogging tolerance in chrysanthemum. Planta. 2018;247(4):899–924.

    Article  CAS  PubMed  Google Scholar 

  51. Osman KA, Tang B, Wang Y, Chen J, Yu F, Li L, et al. Dynamic QTL analysis and candidate gene mapping for waterlogging tolerance at maize seedling stage. PLoS ONE. 2013;8(11):e79305.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  52. Du B, Wang Q, Sun G, Ren X, Cheng Y, Wang Y, et al. Mapping dynamic QTL dissects the genetic architecture of grain size and grain filling rate at different grain-filling stages in barley. Sci Rep. 2019;9(1):18823.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  53. Bac-Molenaar JA, Vreugdenhil D, Granier C, Keurentjes JJ. Genome-wide association mapping of growth dynamics detects time-specific and general quantitative trait loci. J Exp Bot. 2015;66(18):5567–80.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  54. Muraya MM, Chu J, Zhao Y, Junker A, Klukas C, Reif JC, et al. Genetic variation of growth dynamics in maize (Zea mays L.) revealed through automated non-invasive phenotyping. Plant J. 2017;89(2):366–80.

    Article  CAS  PubMed  Google Scholar 

  55. Knoch D, Abbadi A, Grandke F, Meyer RC, Samans B, Werner CR, et al. Strong temporal dynamics of QTL action on plant growth progression revealed through high-throughput phenotyping in canola. Plant Biotechnol J. 2020;18(1):68–82.

    Article  CAS  PubMed  Google Scholar 

  56. Ohbayashi I, Huang S. Mitochondrial pyruvate dehydrogenase contributes to auxin-regulated organ development. Plant Physiol. 2019;180(2):896–909.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  57. Piacentini D, Falasca G, Canepari S, Massimi L. Potential of PM-selected components to induce oxidative stress and root system alteration in a plant model organism. Environ Int. 2019;132:105094.

    Article  CAS  PubMed  Google Scholar 

  58. Takei K, Yamaya T, Sakakibara H. Arabidopsis CYP735A1 and CYP735A2 encode cytokinin hydroxylases that catalyze the biosynthesis of trans-Zeatin. J Biol Chem. 2004;279(40):41866–72.

    Article  CAS  PubMed  Google Scholar 

  59. Kim MJ, Ciani S, Schachtman DP. A peroxidase contributes to ROS production during Arabidopsis root response to potassium deficiency. Mol Plant. 2010;3(2):420–7.

    Article  CAS  PubMed  Google Scholar 

  60. Renault H, El Amrani A, Berger A, Mouille G, Soubigou-Taconnat L, Bouchereau A, et al. γ-Aminobutyric acid transaminase deficiency impairs central carbon metabolism and leads to cell wall defects during salt stress in Arabidopsis roots. Plant Cell Environ. 2013;36(5):1009–18.

    Article  CAS  PubMed  Google Scholar 

  61. Estornell LH, Landberg K, Cierlik I, Sundberg E. SHI/STY genes affect pre- and post-meiotic anther processes in auxin sensing domains in Arabidopsis. Front Plant Sci. 2018;9:150.

    Article  PubMed  PubMed Central  Google Scholar 

  62. Renard J, Martínez-Almonacid I, Sonntag A, Molina I, Moya-Cuevas J, Bissoli G, et al. PRX2 and PRX25, peroxidases regulated by COG1, are involved in seed longevity in Arabidopsis. Plant Cell Environ. 2020;43(2):315–26.

    Article  CAS  PubMed  Google Scholar 

  63. Han S, Fang L, Ren X, Wang W, Jiang J. MPK6 controls H2O2-induced root elongation by mediating Ca2+ influx across the plasma membrane of root cells in Arabidopsis seedlings. New Phytol. 2015;205(2):695–706.

    Article  CAS  PubMed  Google Scholar 

  64. Bouguyon E, Brun F, Meynard D, Kubeš M, Pervent M, Leran S, et al. Multiple mechanisms of nitrate sensing by Arabidopsis nitrate transceptor NRT1.1. Nat Plants. 2015;1:15015.

    Article  CAS  PubMed  Google Scholar 

  65. Maghiaoui A, Bouguyon E, Cuesta C, Perrine-Walker F, Alcon C, Krouk G, et al. The Arabidopsis NRT1.1 transceptor coordinately controls auxin biosynthesis and transport to regulate root branching in response to nitrate. J Exp Bot. 2020;71(15):4480–94.

    Article  CAS  PubMed  Google Scholar 

  66. Jacquot A, Chaput V, Mauries A, Li Z, Tillard P, Fizames C, et al. NRT2.1 C-terminus phosphorylation prevents root high affinity nitrate uptake activity in Arabidopsis thaliana. New Phytol. 2020;228:1038–54.

    Article  CAS  PubMed  Google Scholar 

  67. Katz E, Nisani S, Yadav BS, Woldemariam MG, Shai B, Obolski U, et al. The glucosinolate breakdown product indole-3-carbinol acts as an auxin antagonist in roots of Arabidopsis thaliana. Plant J. 2015;82(4):547–55.

    Article  CAS  PubMed  Google Scholar 

  68. Francisco M, Joseph B, Caligagan H, Li B, Corwin JA, Lin C, et al. The defense metabolite, allyl glucosinolate, modulates Arabidopsis thaliana biomass dependent upon the endogenous glucosinolate pathway. Front Plant Sci. 2016;7:774.

    PubMed  PubMed Central  Google Scholar 

  69. Malinovsky FG, Thomsen MF, Nintemann SJ, Jagd LM, Bourgine B, Burow M. An evolutionarily young defense metabolite influences the root growth of plants via the ancient TOR signaling pathway. Elife Sci. 2017;6.

  70. Urbancsok J, Bones AM, Kissen R. Glucosinolate-derived isothiocyanates inhibit Arabidopsis growth and the potency depends on their side chain structure. Int J Mol Sci. 2017;18(11):2372.

    Article  PubMed Central  CAS  Google Scholar 

  71. Liao P, Lung SC, Chan WL, Bach TJ, Lo C, Chye ML. Overexpression of HMG-CoA synthase promotes Arabidopsis root growth and adversely affects glucosinolate biosynthesis. J Exp Bot. 2020;71(1):272–89.

    Article  CAS  PubMed  Google Scholar 

  72. Katz E, Bagchi R, Jeschke V, Rasmussen ARM. Diverse Allyl glucosinolate catabolites independently influence rootgrowth and development. Plant Physiol. 2020;183(3):1376–90.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  73. Fu L, Wang M, Han B, Tan D, Sun X, Zhang J. Arabidopsis myrosinase genes AtTGG4 and AtTGG5 are root-tip specific and contribute to auxin biosynthesis and root-growth regulation. Int J Mol Sci. 2016;17(6):892.

    Article  PubMed Central  CAS  Google Scholar 

  74. Ohtani M, Demura T, Sugiyama M. Particular significance of SRD2-dependent snRNA accumulation in polarized pattern generation during lateral root development of Arabidopsis. Plant Cell Physiol. 2010;51(12):2002–12.

    Article  CAS  PubMed  Google Scholar 

  75. Lan P, Li W, Wen TN, Shiau JY, Wu YC, Lin W, et al. iTRAQ protein profile analysis of Arabidopsis roots reveals new aspects critical for iron homeostasis. Plant physiol. 2011;155(2):821–34.

    Article  CAS  PubMed  Google Scholar 

  76. Lim CW, Kim JH, Baek W, Kim BS, Lee SC. Functional roles of the protein phosphatase 2C, AtAIP1, in abscisic acid signaling and sugar tolerance in Arabidopsis. Plant Sci. 2012;187:83–8.

    Article  CAS  PubMed  Google Scholar 

  77. Markakis MN, De Cnodder T, Lewandowski M, Simon D, Boron A, Balcerowicz D, et al. Identification of genes involved in the ACC-mediated control of root cell elongation in Arabidopsis thaliana. BMC Plant Biol. 2012;12:208.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  78. Riefler M, Novak O, Strnad M, Schmülling T. Arabidopsis cytokinin receptor mutants reveal functions in shoot growth, leaf senescence, seed size, germination, root development, and cytokinin metabolism. Plant Cell. 2006;18(1):40–54.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  79. Murphy E, Vu LD, Van den Broeck L, Lin Z, Ramakrishna P, van de Cotte B, et al. RALFL34 regulates formative cell divisions in Arabidopsis pericycle during lateral root initiation. J Exp Bot. 2016;67(16):4863–75.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  80. Takahashi M, Umetsu K, Oono Y, Higaki T, Blancaflor EB, Rahman A. Small acidic protein 1 and SCF(TIR)(1) ubiquitin proteasome pathway act in concert to induce 2,4-dichlorophenoxyacetic acid-mediated alteration of actin in Arabidopsis roots. Plant J. 2017;89(5):940–56.

    Article  CAS  PubMed  Google Scholar 

  81. Hoagland DR. Optimum nutrient solutions for plants. Science. 1921;52(1354):562–564

    Article  Google Scholar 

  82. Clarke WE, Higgins EE, Plieske J, Wieseke R, Sidebottom C, Khedikar Y, et al. A high-density SNP genotyping array for Brassica napus and its ancestral diploid species based on optimised selection of single-locus markers in the allotetraploid genome. Theor Appl Genet. 2016;129(10):1887–99.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  83. Pritchard JK, Stephens M, Donnelly P. Inference of population structure using multilocus genotype data. Genetics. 2000;155(2):945–59.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  84. Bradbury PJ, Zhang Z, Kroon DE, Casstevens TM, Yogesh R, Buckler ES. TASSEL: software for association mapping of complex traits in diverse samples. Bioinformatics. 2007;19:19.

    Google Scholar 

  85. Yu J, Buckler ES. Genetic association mapping and genome organization of maize. Curr Opin Biotechnol. 2006;17(2):155–60.

    Article  CAS  PubMed  Google Scholar 

  86. Wei L, Jian H, Lu K, Filardo F, Yin N, Liu L, et al. Genome-wide association analysis and differential expression analysis of resistance to Sclerotinia stem rot in Brassica napus. Plant Biotechnol J. 2016;14(6):1368–80.

    Article  CAS  PubMed  Google Scholar 

  87. Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. 2008;9:559.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

Download references

Acknowledgements

Great gratitude goes to linguistics Prof. Ping Liu from Huazhong Agriculture University, Wuhan, China for her work at English editing and language polishing.

Funding

This work was supported by the National Key Research and Development Program of China (2016YFD0100202), the Agricultural Science and Technology Innovation Project (CAAS-ZDRW202105), the Natural Foundation of Hubei Province (2019CFB617), the Technical Innovation Project of Hubei Province (2018ABA087), and China Agriculture Research System of MOF and MARA (CARS-12).

Author information

Author notes

  1. Keqi Li and Jie Wang contributed equally to this work

Authors and Affiliations

  1. Oil Crops Research Institute of the Chinese Academy of Agricultural Sciences/Key Laboratory of Biology and Genetic Improvement of Oil Crops, Ministry of Agriculture, Wuhan, 430062, China

    Keqi Li, Jie Wang, Lieqiong Kuang, Ze Tian, Xinfa Wang, Xiaoling Dun & Hanzhong Wang

  2. National Key Laboratory of Crop Genetic Improvement, Huazhong Agricultural University, Wuhan, 430062, China

    Keqi Li & Jinxing Tu

Authors
  1. Keqi Li
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Jie Wang
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Lieqiong Kuang
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Ze Tian
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Xinfa Wang
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Xiaoling Dun
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. Jinxing Tu
    View author publications

    You can also search for this author in PubMed Google Scholar

  8. Hanzhong Wang
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

All the authors planned and supervised the research. KL, JW, LK, ZT, XW, and JT performed root traits investigation and analyzed the data; KL wrote the manuscript; XD, and HW contributed to modify the manuscript. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Xiaoling Dun.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

All authors consent for publication.

Competing interests

All authors state that they have no competing interests.

Additional information

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1

: Figure S1 Phenotype of plants at different stages in B. napus. (a–e) Plants at 13 DAS, 10 DAT, 3 EL, 5 EL, and 7 EL, respectively. (f) Materials sowed on germination device. (g) Plants in germination device 6 days after sowing. (h) Plants in growth device. Scale bars = 3 cm (a–e), 4 cm (f, g), and 8 cm (h). Figure S2 Frequency distribution of root-related traits and shoot-related traits at five stages. (a–i) Frequency distribution of SFW, RFW, SDW, RDW, TNR, PRL, TRL, TSA, and TRV at the five stages (13 DAS, 10 DAT, 3 EL, 5 EL, and 7 EL), respectively. Figure S3 Analysis of population structure and kinships of 280 B. napus accessions (a) Log-likelihood data of possible clusters, K: from 1 to 10. (b) Distribution of pairwise relative kinship. (c) Population structure of 280 accessions. Figure S4 Positive correlation between RNA-seq data and qRT-PCR data. Figure S5 Principal component analysis of the transcriptome sequencing data. Figure S6 Heatmap of module eigengenes obtained by WGCNA. (a-l) Heatmaps of the expression profile of eigengenes in the purple, green, black, brown, darkorange, darkturquoise, white, darkred, red, lightyellow, saddlebrown and darkgrey modules, respectively.

Additional file 2: Table S1

Trait statistics collected at the five stages of each repetition. Table S2 Correlations in each captured trait among replication at five stages. Table S3 Summary of SNPs and LD decay on 19 chromosomes of B. napus. Table S4 Detailed information on trait-related significant SNPs identified by GWAS. Table S5 Detailed information on 48 valid QTL clusters. Table S6 Correlations among root-related traits at each examined stage. Table S7 RNA-Seq statistics of four growth types against B. napus reference genome. Table S8 FPKM of persistent DEGs. Table S9 GO enrichment results of persistent DEGs. Table S10 FPKM of stage-specific DEGs. Table S11 GO enrichment results of stage-specific expressed genes. Table S12 FPKM of genes located in haplotype blocks on 48 QTL clusters. Table S13 Information on peak SNPs overlapped with SNPs reported by Wang et al. (2017). Table S14 Primers used in this study.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Li, K., Wang, J., Kuang, L. et al. Genome-wide association study and transcriptome analysis reveal key genes affecting root growth dynamics in rapeseed. Biotechnol Biofuels 14, 178 (2021). https://doi.org/10.1186/s13068-021-02032-7

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s13068-021-02032-7

Keywords

  • Rapeseed
  • Root growth
  • Persistent
  • Stage-specific
  • GWAS
  • WGCNA

Biotechnology for Biofuels and Bioproducts

ISSN: 2731-3654

Contact us