Understanding the genetic and molecular architecture underlying complex quantitative traits or phenotypic variation is the central goal of current biology (Kliebenstein, 2009). Changing transcript abundance arising from sequence polymorphism is one possible underlying molecular mechanism for complex traits (Rockman and Kruglyak, 2006). A wide range of complex traits within multiple species have been shown to be controlled by their underlying gene expression (Kliebenstein, 2009). A large-scale analysis of the relationship between complex diseases and known non-synonymous single nucleotide polymorphisms (SNPs) found that an increased number of SNPs did not find a significant number of associations, further indicating the potential role of variation in gene expression in controlling phenotypic traits (Burton et al., 2007). Together, these findings support the hypothesis that gene expression acts as the intermediate molecular phenotype, connecting genotypic variation to phenotypic variation (Rockman and Kruglyak, 2006).

The traditional quantitative trait loci (QTLs) mapping approach has been widely used to dissect the molecular mechanism underlying complex traits. However, most complex traits are controlled by multiple genes (Kliebenstein, 2009), and multiple QTLs with marginal effects (Mackay et al., 2009), which can be difficult to dissect with traditional QTL analysis. With the advent of advanced high-throughput technologies, including expression arrays and genome sequencing technologies, a new approach emerged to meet the challenge of complex traits, gene expression quantitative trait locus (eQTL) analyses. This approach is characterized by its combination of population genomic and quantitative genetics (Jansen and Nap, 2001). Briefly, QTL analysis is applied to the transcript abundances measured from a well-genotyped segregating population, allowing genetic loci (eQTL) controlling variation of each transcript to be identified on a genome-wide level (Keurentjes et al., 2007). A locus of sequence polymorphism affecting gene expression that is tightly linked to the physical position of the gene itself is defined as a cis-eQTL; all other eQTLs act in trans (trans-eQTLs) (Brem et al., 2002). The heritability of transcript variation and the genetic effect on this variation can also be inferred from eQTL analysis. Studies from Arabidopsis (West et al., 2007; Keurentjes et al., 2007), yeast (Brem and Kruglyak, 2005), maize (Schadt et al., 2003), and barley (Potokina et al., 2007), revealed an uneven distribution of trans-eQTLs with some genomic regions identified as eQTL hotspots (Schadt et al., 2003), enriched with numerous trans-eQTLs. Identification and dissection of eQTL hotspots may help to identify underlying master regulators or pleiotropic genes, and concomitantly the affected downstream genes or biological pathways. Thus, global eQTL analysis has the potential to yield great insight into the complex genetic architecture of transcript variation, and when linked to phenotypic variation (pQTLs), molecular mechanisms underlying complex traits can be directly inferred.

Several in silico large-scale data analyses have demonstrated the power of gene co-expression networks in dissecting a range of biological problems such as secondary wall biosynthesis (Persson et al., 2005), glucosinolate accumulation (Hirai et al., 2007), and seed germination and dormancy (Bassel et al., 2011) at the systems level. Gene regulatory networks were also hypothesized to constitute important information that can be inferred from eQTL data (Jansen and Nap, 2001); yet, this component has often been overlooked. Studies on seed development of Arabidopsis demonstrated that seed filling was controlled by a dynamic and distinct developmental transcriptional program, and genes involved in central carbon transport and multiple other metabolic pathways were implicated (Ruuska et al., 2002). However, most eQTL analyses have been limited to single gene-level analysis.

Fibre in seed meal negatively affects digestion and processibility for animal feed, and reducing the seed fibre content in oilseed rape could result in increased seed oil and protein content (Wittkop et al., 2009; Snowdon et al., 2010; Liu et al., 2012, 2013; Yu et al., 2016; Behnke et al., 2018; Miao et al., 2019; Chao et al., 2022; Yusuf and Möllers, 2024). Therefore, the low fibre trait was identified as a high-value target for canola (Brassica napus) breeding programs. Crude fibre is mainly comprised of lignin, including neutral detergent fibre (NDF), acid detergent fibre (ADF), acid detergent lignin (ADL), and polysaccharides including cellulose and hemicellulose (Van Soest et al., 1991; Wittkop et al., 2012). Genetics studies revealed a few major QTLs for lignin traits concentrated on chromosomes A09, C09, and C05 (Badani et al., 2006; Liu et al., 2012, 2013; Stein et al., 2013; Yu et al., 2016; Behnke et al., 2018; Miao et al., 2019; Chao et al., 2022). Lignin biosynthetic pathway(s) are complex. Forward and reverse genetics and biochemical analyses have identified more than 600 involved genes, with a subset of key enzymes that include phenylalanine ammonia-lyase1 (PAL1) and PAL2; cinnamate 4-hydroxylase (C4H); 4-coumarate: CoA ligase1 (4CL1) and 4CL2; caffeoyl-CoA O-methyltransferase1 (CCoAOMT1); cinnamoyl-CoA reductase1 (CCR1); ferulate 5-hydroxylase (F5H1); caffeic acid O-methyltransferase (COMT); and cinnamyl alcohol dehydrogenase6 (CAD6) (Vanholme et al., 2010). Previous studies also reported that QTLs for seed colour colocalize with lignin content in B. napus seeds (Badani et al., 2006; Liu et al., 2012, 2013; Yu et al., 2016; Chao et al., 2022). In Arabidopsis, genes in the flavonoid pathway, including those from the TRANSPARENT TESTA (TT) family (eg. TT1-TT19) and transcription factor complex MYB-bHLH-WD40, are key regulators controlling seed colour (Lepiniec et al., 2006; Xu et al., 2014, 2015). The biosynthetic pathways of lignin and pigments share the same precursors (Vanholme et al., 2012), indicating that the two complex biosynthetic pathways may act together to control both lignin and seed colour, and possibly other quality traits in B. napus seeds.

As Arabidopsis is a very close relative of B. napus within the Brassicaceae family, it is reasonable to infer that seed quality traits of B. napus, e.g., lignin, are also controlled by interactions of multiple metabolic and biological pathways. Therefore, a systems genetic approach that combines quantitative genetics, regulatory network analysis, and eQTL analysis would be more appropriate to investigate seed quality traits within B. napus ( Kliebenstein, 2009 and Marshall-Colón and Kliebenstein, 2019; Kliebenstein, 2020). In the present study, we aimed to fully exploit the power of this approach to decipher the regulatory network underlying lignin traits in B. napus seeds.

Complex traits or quantitative traits are typically controlled by complex genetic architecture consisting of hundreds to thousands of genes (Kliebenstein, 2009). These genes form a genetic regulatory network that controls the trait, which challenges the single gene or traditional QTL analysis. Therefore, we utilized a systems genetic approach that integrates the quantitative genetic, gene co-regulatory network, network-QTL, and eQTL analysis to uncover the causal regulatory network underlying lignin content in B. napus seeds.

Three major QTLs were identified on A09 spanning a physical interval of approximately 13 Mb. All three QTLs conferred stable and relatively large effects on lignin content. QTLs for NDF, ADF, and ADL colocalized at the same genetic locus, indicating these three components of lignin shared the same genetic architecture. Numerous QTLs for seed quality traits including seed colour, oil, protein, and fibre were previously identified on the A09 chromosome (Wittkop et al., 2009; Liu et al., 2012, 2013; Yu et al., 2016; Behnke et al., 2018; Miao et al., 2019; Chao et al., 2022; Yusuf and Möllers, 2024), indicating that A09 is a hotspot that harbours several genes for these seed quality traits. However, this study and all previous research have used a relatively small population size, from ~100 to 300 lines, making it challenging to precisely identify QTLs within this hotspot. Future, larger and advanced-structured genetic populations could be used to dissect this hotspot. In addition, a stable and large effect homologous QTL to the qLignin_A09_2 was identified on C08, indicating both the A and C genome contributed to the variation of lignin content in our population. Currently, limited studies, including Badani et al. (2006) and Chao et al. (2022), have reported the QTL on C08 for seed colour and lignin content. The relative simplicity of the C08 QTL, compared to A09, makes it a better target to uncover genes underlying lignin QTL in the future. Furthermore, since the C08 QTL is homologous to the A09 lignin QTL, uncovering its genetic architecture may also shed light on the A09 hotspot.

The genome-wide transcriptome datasets for developed seeds of the SG population, allowed a seed development biological network to be successfully built ( Figure 2 ). Biological network construction assumes that it will follow the topological scale-free property, in which the distribution of the nodes and edges follows a power-law distribution (Barabási and Oltvai, 2004). Although definitions of the properties of biological networks is still under discussion, the seed developmental regulatory network did exhibit the expected scale-free topology, described as a small number of nodes with a large number of connections, while most nodes have few connections ( Figure 2 ). The distribution of the nodes and edges in our network exhibited a power-law distribution, supporting its topological free property, and was in agreement with the network property discovered in Arabidopsis (Atias et al., 2009; Mao et al., 2009), rice (Ficklin et al., 2010; Ficklin and Feltus, 2011) and maize (Ficklin and Feltus, 2011).

The structure of genome-wide co-expression networks are hierarchical with modular units (Aoki et al., 2007). Given that genes with similar biological functions are typically co-expressed together and formed functional modules, disassembling biological networks into functional modules or clusters is more informative and practical for biologists to interpret the network (Mutwil et al., 2010). Here, the seed network was divided into 83 computational biological modules. If phenotype variation is attributed to the variation of genes within the module, a correlation should be observed between modules and phenotypes. This was the case for two modules (modules 33 and 35) that were significantly positively correlated with lignin, and one (module 67) that was negatively coordinated with lignin. With network QTL analysis, it was found that lignin QTLs (pQTLs) and module QTLs shared the same genetic loci on A09, supporting the suggestion that lignin variation was determined by the interactions of these three functional modules. These findings also support the hypothesis that entire plant networks can be inherited as a phenotypic trait (Marshall-Colón et al., 2019). Large-scale evaluation of 1,671 genes on agronomic performance in maize using a single gene engineering approach found that only one gene encoding a transcription factor, zmm28, showed promise, indicating the marginal effect of a single gene for the improvement of complex traits (Simmons et al., 2021). The defined lignin network could enable breeders to target a specific regulatory network and potentially use an engineering approach to rewire and fine-tune the network to improve crop performance.

Since modules 33, 35, and 67 were altogether associated with seed lignin, these three modules were assembled into one large lignin-related functional module, called “lignin causative network”, to perform further analysis. This network was enriched with secondary metabolism-related genes. One group of these genes includes PAL2, C4H, and CCR1, encoding three key enzymes controlling critical steps of the lignin pathway (Vogt, 2010). The other group, which mainly consists of flavonoid pathway genes, is involved in the accumulation of flavonoids in the seed coat, such as TT3, TT7, TT6, TT8, TT18, TT19, and BAN. All these genes are significantly and positively associated with seed lignin variation. Enrichment of these secondary metabolite genes indicated the high reliability of this causative regulatory network underlying lignin content.

eQTL analysis was performed on genes within this network to identify the causal genes underlying the network for lignin content. The lignin gene, CCR1, was detected as a cis-eQTL, while both PAL and C4H were found as trans-eQTLs. Flavonoid genes including TT4, TT6, and TT8 were identified as cis-eQTLs, while the remaining transparent testa genes were trans-eQTLs. In addition, all these cis-eQTLs colocalize with the fibre QTL from chromosome A09, indicating they are prioritized as causal genes underlying this network (Kliebenstein, 2009). Liu et al. (2012) pinpointed CCR1 from chromosome A09 as the candidate gene for ADF in B. napus seeds. Zhang et al. (2022) identified a CCR-like gene, CCRL, from C07 as one candidate gene for seed content in B. napus. CCR1 has been found to encode the key enzyme for lignin accumulation in stems (Jones et al., 2001). The cis-controlled eQTL impacting the variation of CCR1, the co-localization with phenotype QTL, and the unique known gene from the monolignol pathway located on our lignin causative network as well as previous findings, indicated that CCR1 is one of the causative genes underlying the seed lignin trait. Of the three cis-regulated flavonoid genes, TT8 was found to be a key regulator for seed oil and colour in Arabidopsis and B. napus (Chen et al., 2014; Xu et al., 2014, 2015; Zhai et al., 2020; Zhang et al., 2022) and lignin in B. napus (Zhang et al., 2022). In addition, Zhai et al. (2020) found TT8, as a BHLH transcription factor, regulated both the flavonoid and lignin pathways. All these findings pinpoint TT8 as another key regulator in this network. Li et al. (2023) reported that TT4, which encodes an enzyme acting at an earlier step of the flavonoid pathway and located on chromosomes A03 and C02, is key regulator for oil and protein content in B. napus seeds, while we identified TT4 from A09 as playing a key role for lignin content. Therefore, it is hypothesized that multi-homologous and homoeologous TT4 genes in different genotypes control seed quality traits. In addition, another earlier flavonoid biosynthetic gene TT6 was identified as a potential new key regulator for lignin content in B. napus. Thus, the lignin traits associated with the A09 loci are controlled by a complex regulatory network with multiple genes.

Notably, we did not find any cis-eQTL from the key lignin and flavonoid genes associated with the homologous C08 lignin QTL. Network QTL analysis only found that the QTL of module 67 colocalized with the qLignin_C08 QTL. Genes for the development of B. napus seed are spatially and temporally regulated (Zhai et al., 2020). It is possible that the A and C genome coordinately regulate seed lignin content in a temporal manner, while the regulatory network built here, was only from one key time point, 21 days after flowering. We observed, under greenhouse conditions, that the A donor line of PSA12 flowered and matured earliest, followed by DH12075 and PSA12, while the C genome donor line of PSA12 flowered and matured very late (unpublished data), supporting this hypothesis, since positive alleles at the C08 locus were contributed by PSA12. A further spatial-temporal regulatory network built on seed development will be more powerful for dissecting the genetic architecture controlling seed quality traits and should be considered for future studies.

Another interesting finding of the lignin network was that it was enriched with biotic stress genes, including key transcript factors, WRKY25, 33 and SIGMA FACTOR BINDING PROTEIN 1 (SIB1), which are involved in triggering resistance to the necrotrophic pathogen Botrytis cinerea in Arabidopsis (Andreasson et al., 2005; Zheng et al., 2006; Eulgem and Somssich, 2007; Pandey and Somssich, 2009; Lai et al., 2011). However, all these genes were significantly negatively associated with lignin accumulation, which is contradictory to the positive role of lignin on enhancing the resistance response (Vanholme et al., 2010). This paradoxical role of lignin in the plant defence mechanism compelled us to hypothesize that cell wall integrity plays a key role in this process. Reduced lignin affects the integrity of the wall which mimics disease infection, leading to the release of defence-related elicitors. When a normal plant is infected by disease, cell wall integrity is damaged and in response to this, the cell wall induces lignin biosynthesis to maintain its integrity. This theory is further supported by studies which have shown that cell wall damage directly triggers lignin biosynthesis (Denness et al., 2011).

It is worth noting a few limitations of the current study. First, the population size was small, and may not provide enough statistical power for identifying minor QTL. Second due to the limitation of array technology, many unknown genes within our causative network may be missing from our research. Third, this regulatory network was built only from one stage of seed development, and could not capture a complete temporal picture of development. New advanced structured genetic populations, and genome sequencing tools including long read and pangenome sequence analysis will empower these systems genetic approaches to uncover the genetic mechanism of complex traits.

An example of the power of systems genetics that combines quantitative genetic, regulatory network analysis, network QTL, and eQTL analysis to identify causal regulatory networks and causal genes for lignin traits was provided. Further insights were revealed into the complex genetic architecture that regulates lignin content, specifically that controlled by the locus on A09 of B. napus. The analyses provide more prioritized genes to target to improve seed quality trait in B. napus. Moreover, the research here may also provide a new avenue for crop communities to use regulatory networks to select or design crops with better performance and quality in a changing climate.

The data presented in the study are deposited in the NCBI data repository, accession number GSE268725.

WZ: Data curation, Formal analysis, Investigation, Methodology, Validation, Visualization, Writing – original draft, Writing – review & editing. EH: Data curation, Formal analysis, Investigation, Methodology, Validation, Visualization, Writing – review & editing. SR: Data curation, Methodology, Writing – review & editing. WC: Data curation, Methodology, Writing – review & editing. KB: Formal analysis, Visualization, Writing – review & editing. AS: Conceptualization, Funding acquisition, Writing – review & editing. PF: Conceptualization, Funding acquisition, Investigation, Project administration, Supervision, Writing – review & editing. IP: Conceptualization, Funding acquisition, Methodology, Project administration, Resources, Supervision, Validation, Writing – original draft, Writing – review & editing.