Crop breeding programs are expressing increased interest for root traits to improve the capture of soil resources. The optimization of the root system architecture holds the potential to sustain yield and to reduce the negative environmental impact of mineral N fertilization. To some extent, the NupE depends on the rooting pattern, which shows substantial genotypic variation in rapeseed (Wang et al. 2017; Kiran et al. 2019; Ibrahim et al. 2021; Xu et al. 2022). In modern breeding programs, diversity panels are offering valuable genetic resources. Taking advantage of such tool, we analyzed the root phenotypes of 300 rapeseed accessions grown in paper pouches with two divergent N levels (Table S1, Fig. S1). The root biomass and most root morphological traits positively correlated with shoot biomass across N levels (Fig. S3). These observations complete prior reports indicating concomitant biomass production of root and shoot organs at a seedling stage (Thomas et al. 2016b; Louvieaux et al. 2018, 2020a, Kupcsik et al. 2021). Therefore, root traits can be considered as positive indicators of biomass production. Furthermore, the two-dimensional root analysis at germination could be a reliable indicator for later developmental stages in soil (Louvieaux et al. 2020b) and for seed quality (Koscielny and Gulden 2012; Louvieaux et al. 2020a).
A substantial variability for root morphologies is leaving space for breeding targets
A considerable degree of variability in root morphologies was observed and illustrated (Fig. 1, Table S2). The winter-type (along with winter fodder and semi-winter) accessions were characterized by greater biomass production and larger root system size compared to the spring-type and the other ones (Fig. 2). The superior root characteristics of these genotypes with long life-cycle period requiring vernalization, are in line with previous observations (Rahman and McClean, 2013; Thomas et al. 2016a; Arifuzzaman and Rahman, 2017). The diversity of root morphology observed among these groups can serve as a genetic resource for breeding root morphology.
In addition, important broad sense heritability values for root morphological traits were highlighted, together with important genotype, moderate nutrition and minimal genotype-by-nutrition effects (Table S4-5). This suggests a great genetic correlation between the two N levels, opening routes for selection. In particular, varieties showing important root system size across N levels could be interesting. As such, Bellini, Conny, Kurander, and Shang-You are rather unresponsive to N level (TRL medians < 10% difference between N levels) (Table S2). These observations were consistent across experimental locations. Although precautions were taken to standardize conditions in growth chambers, differences in the panel means (Fig. S2, Table S3) between the two locations were noted for some traits (e.g., S and ΣLLR). This may be due to differences in light quality between locations.
Genetic structuration and linkage disequilibrium were investigated prior to GWAS design
The population structure and the kinship may constitute a confounding for GWAS. Hence, these were considered in the GLM to reduce the risk of false positive. The growth habit was a major determinant explaining the genetic structure of the diversity panel (Fig. S4), in accordance with previous reports (Gazave et al. 2016; Delourme et al. 2013). Investigating linkage disequilibrium enabled to evaluate the resolution of QTL mapping for every chromosome. The LD decay ranged between 0.155 Mb and 10 Mb, depending on the chromosome, with greater values for the C subgenome (Table S6), in accordance with Werner and Snowdon (2018). However, the LD substantially varied from locus to locus, as mirrored by the variation of haplotype block sizes within one chromosome. The variation in LD pattern along one chromosome is related to recombination frequency, which depends on breeding history and mating system of the species (Ecke et al. 2010). Rapeseed notably underwent diversity bottlenecks, following seed quality improvement (e.g., reducing erucic acid and glucosinolate content). This led to a decrease of genetic diversity and extension of LD throughout breeding history (Delourme et al. 2013).
Genome wide association identified 14 stable QTLs across N levels
A first mapping method (GWAS) was employed to identify QTLs associated with root traits. The genotype-to-environment interaction leads to the instability of the QTLs from one environment to another (e.g. nutrition, experimental replicates) (Corlouer et al. 2019). To strengthen QTL discovery, we applied a multi-environment design (two N levels and two experimental locations) to GWAS models and assessed QTL stability across N levels and experimental locations. The association mapping detected 161 QTLs in one single N environment. By applying a multi-environment model, the number was narrowed down to 26 main QTLs (Table 2, Table S7). Fourteen out of 26 were repeatedly detected under N + and N‒ levels, indicating a genetic stability of these QTLs which, therefore, could be promising candidates in breeding programs. These 14 QTLs were associated with traits contributing to both primary root (eg., LZ2) and lateral root phenotypes (eg., ΣLLR), which both were positively correlated (Fig. S3). Combining the corresponding favorable alleles could improve simultaneously these traits. These regions may contain either pleiotropic genes or closely linked genes showing different functions. Indeed, the genetic basis of trait correlation can be due to the pleiotropy of one single gene or to multiple genes in one LD block (Chen and Lübberstedt, 2010).
Allele frequency differences between the two bulks showed variation across the genome with a peak on A06 chromosome
In parallel, one BSA approach was conducted in one experimental F2 population after crossing two accessions with contrasting lateral root development. Recently, combining BSA and GWAS proved to be a successful strategy to identify loci regulating target traits (Takagi et al. 2013; de la Fuente Canto and Vigouroux, 2022). Nonetheless, the optimization of the BSA experimental design for QTL detection is still a matter of debate, with rather small (Branham et al. 2018) or large (Yang et al. 2013) numbers of progenies composing the bulks. Following Takagi et al. (2013) and Huang et al. (2022), we increased bulk proportion up to 15% of the population size, to increase the power of detection. The normal distribution pattern of ∑LLR measured in the F2 progenies, indicated a polygenic control of that trait (Fig. 4a). The AFD showed variation through the genome with some enrichment zones (AFD > 0.15) (Figure S5). However, only one QTL was identified on A06 (Fig. 4b). However, this QTL was located 210 kb away from main QTL 9 identified by GWAS (Fig. 3). In addition, the enrichment zone of A04 overlapped GWAS-QTLs on a region extending from ~ 21.0-21.6 Mb (Fig. 3, Table S11).
Investigation of gene content and homoeologous links revealed some duplicated QTLs
Two lists of candidate genes were generated: (i) 2,489 genes within the 26 main QTLs detected by GWAS (Table S8), (ii) 292 genes within the single QTL identified by BSA (Table S10). Given the polyploid nature of the rapeseed genome, the gene content of the main QTL regions was analyzed to seek possible homoeology between the numerous QTLs. Such homoeologous rapeseed QTLs have notably been reported for glucosinolate content (Howell et al. 2003), yield related traits (Bouchet et al. 2016b; Bilgrami et al. 2023), or the resistance to stem canker (Fopa-Fomeju et al. 2014). In the present study, some homoeologous genes were covered by QTL pairs (Table S9, Fig. 3), indicating conservation between A and C subgenomes. For example, homoeologous QTLs were found on A01 and C01 controlling TRL (regions 1–19), as well as on A04 and C04 governing LZ2, SLLR and TRL (regions 8–21). The QTL pairs could indicate candidate homoeologous genes with similar functions. Among list (i), a total of 1,933 genes were involved in QTL pairs (Table S8).
Integration of Arabidopsis orthology information provides possible indications in rapeseed root development
The functional annotations of A. thaliana orthologs were searched in relation with root development, hormonal regulation and N nutrition. Main QTL regions showing an homoeologous QTL were screened. The survey indicated B. napus gene BnaA01p14400, an ortholog of the A. thaliana At4g23100 (ROOT MERISTEMLESS 1, RML1), which encodes an enzyme in glutathione biosynthesis that is required for cell proliferation at the root tip (Vernoux et al. 2000); BnaC05p41090 (At3g22270, PAT1 HOMOLOG 1, PAT1H1) maintaining root stem niche stability (Yu et al. 2016); BnaC05p42660 (At3g21160, ALPHA-MANNOSIDASE PROTEIN 2, MNS2) playing a role in root development and cell wall biosynthesis (Liebminger et al. 2009); BnaC07p02570 (At2g19500, CYTOKININ OXIDASE 2, CKX2) encoding a dehydrogenase that catalyzes the degradation of cytokinins to influence the directional lateral root growth (Waidmann et al. 2019); BnaC07p03950 (At2g18800, XYLOGLUCAN ENDOTRANSGLUCOSYLASE/HYDROLASE 21, XTH21) involved in cell wall modification processes during primary root development (Liu et al. 2007); BnaC07p04900 (At2g18040, PEPTIDYLPROLYL CIS/TRANS ISOMERASE, PIN1AT) involved in root meristem maintenance by affecting auxin transport (Li et al. 2015b); BnaC08p45790 (At1g14350, MYB DOMAIN PROTEIN 124, MYB124) encoding a transcription factor involved in early steps of lateral root formation (Chen et al. 2015); BnaC08p45580 (At2g02020, NITRATE TRANSPORTER 1/PEPTIDE TRANSPORTER FAMILY 8.4, NPF8.4/PTR4) encoding a nitrate transporter localized at the tonoplast to regulate vacuolar storage in the vasculature (Weichert et al. 2012).
In some main QTL regions stable across N levels, BnaA06p33930 (At2g17520, INOSITOL REQUIERING ENZYME 1, IRE1) plays a role in the maintenance of root meristems under salt stress (Iwata et al. 2017); BnaA07p31500 (At1g69588, CLAVATA3/ESR-RELATED 45, CLE45) could be involved in protophloem cells fate, thereby affecting root meristem growth (Depuydt et al. 2013). In addition, some nitrate transporters were identified in main region 17. Among them, BnaA07p31940 (At1g68570, NPF3.1) is mainly expressed in root epidermis (David et al. 2016), BnaA07p31340 (At1g69870, NPF2.13) is involved in nitrate remobilization (Fan et al. 2009), and BnaA07p30370 (At1g72130, NPF5.11) is a low affinity transporter in the tonoplast (Lu et al. 2022).
In some main QTL regions showing both homoeologous QTLs and stability across N levels, BnaA02p31940 (At2g01420, PIN-FORMED 4, PIN4) codes for an auxin efflux carrier located around the quiescent center of the root meristem (Friml et al. 2002); BnaA02p31570 (At5g47750, D6 PROTEIN KINASE-LIKE 2, D6PKL2) plays a role in auxin transport, while loss-of-function results in lateral root initiation defects (Zourelidou et al. 2009); BnaA04p06590 (At3g52890, KCBP-INTERACTING PROTEIN KINASE 1, KIPK1) is a negative regulator of root growth (Humphrey et al. 2015); BnaC01p20000 (At4g24670, TRYPTOPHAN AMINOTRANSFERASE-RELATED 2, TAR2) is involved in auxin biosynthesis and modulates lateral root growth in response to N (Ma et al. 2014); BnaC03p42810 (At3g07360, PLANT U-BOX 9, PUB9) plays a role in auxin-mediated lateral root development during phosphorus starvation (Deb et al. 2014); BnaC09p55380 (At5g20820, SAUR-LIKE AUXIN-RESPONSIVE PROTEIN 76, SAUR76) regulates meristem activity in response to auxin (Markakis et al. 2013). Finally, main region19 covered BnaC01p15620 (At4g21680, NPF7.2), a nitrate transporter mediating low-affinity transport to unload nitrate from xylem vessels (Li et al. 2010).
Finally, in the QTL identified by BSA, BnaA06p02150 (At1g51950, INDOLE-3-ACETIC ACID 18, IAA18) is negatively affecting lateral root formation (Lakehal et al. 2019; Uehara et al. 2008), BnaA06p01210 (At1g5310, MITOGEN-ACTIVATED PROTEIN KINASE 18, MPK18) influences root growth dynamic through cortical microtubule function (Walia et al. 2009) and BnaA06p02570 (At1g51220, WIP DOMAIN PROTEIN 5, WIP5) is involved in cell fate determination of the quiescent center during primary root development (Crawford et al. 2015).