Performance of Eight Lines under P-Concentration Gradients
Eight lines, H1, H2, H3, H4, H5, H6, H7, and H8, were randomly chosen from the B. napus association panel and cultivated hydroponically in two independent experiments with eight P levels ranging from P1 to P8 with 1, 0.1, 0.05, 0.025, 0.01, 0.007, 0.005 and 0.003 mM P+, respectively. To determine the effect of different levels of phosphorus deficiency on rapeseed growth, the shoot fresh weight (SFW) of the eight lines was examined (Additional file 1; Fig. S1a). As shown (Additional file 1; Fig. S1a-b), SFW decreased with decreasing treatment concentrations, demonstrating that P stress profoundly affects plant growth. The SFW of all eight lines declined significantly at 0.01 mM P+ (P5) compared with one mM P+ control, with a ratio of less than 50% (Additional file 1; Fig. S1b). Furthermore, there were significant differences between the genotypes regarding their SFW ratio, with values ranging from 14.3–40.9% (Additional file 1; Fig. S1b), suggesting that phosphorus uptake efficiency differs significantly among genotypes of B. napus.
Phenotypic variations of root and biomass traits under LP stress in the association population
Thirteen root and biomass-related traits (Table 1) were investigated in 327 B. napus association mapping panel under control (1 mM P+) and low phosphorus stress (0.01 mM P+) treatments. Under control and LP conditions, the coefficient variation of the 13 traits ranged from 14.5–45.6% and 10.0–32.0%, respectively. For all the studied traits, moderate to high heritability of 0.47–0.72 was observed in both conditions (Table 2). Furthermore, our findings demonstrated that phosphorus stress had a detrimental effect on normal root and shoot growth, resulting in a decrease in shoot fresh weight (SFW) by 131.8%, shoot dry weight (SDW) by 76.8%, total fresh weight (TFW) by 89.2%, and total dry weight (TDW) by 60.7% (Table 2). Some of the studied traits were increased significantly under LP stress; primary root length (PRL) by 8.5%, root fresh weight (RFW) by 13.4%, total root length (TRL) by 14.9%, total root surface area (TSA) by 17.1%, total root volume (TRV) by 15.4%, total root numbers (TRN) by 10.3%, root dry weight (RDW) by 16.1%, the root-shoot ratio in fresh weight (RSRF) by 62.9% and root-shoot ratio in dry weight (RSRD) by 53.8% (Table 2).
The Pearson’s correlation coefficients between the measured traits studied under LP and CK treatments are presented in Fig. 1a-b. Obviously, SFW was found to be positively and significantly correlated with RFW under LP stress (r = 0.62, P < 0.01), which was consistent with the correlation under the control condition [27]. Root morphology contributed significantly to the formation of biomass-related traits, as demonstrated by the significant and strong correlation between SFW and TRL (r = 0.47, P < 0.01), TSA (r = 0.33, P < 0.01), TRV (r = 0.47, P < 0.01) and TRN (r = 0.39, P < 0.01). Moreover, negative correlations were also detected with the values of -0.37 between SFW and RSRF and − 0.43 between SDW and RSRD, respectively (Fig. 1a). In addition, the Pearson correlation for root and biomass traits between control and LP-stress treatments showed the highest correlation for PRL (r = 0.53, P < 0.01) followed by RFW (r = 0.50, P < 0.01) and TRL (r = 0.48, P < 0.01). At the same time, TRN exhibited the lowest correlation (r = 0.15, P < 0.01) (Fig. 1b). Further, the study illustrated that multiple traits interrelated to PUE should be considered for a comprehensive assessment of RSA traits. Furthermore, normal and continuous frequency distributions were observed under LP stress for the studied traits, indicating that the investigated accessions would be appropriate for subsequent association studies (Fig. 1a).
Marker trait-association analysis for root and biomass traits under LP-stress
Brassica 50K Illumina Infinium SNP array containing 45,708 SNPs was used to genotype the association panel. Following SNP filtering, 20,131 SNP markers were used further to find an association between traits and SNPs [27]. Since the results of association analysis under normal conditions were already reported [28], this study only conducted GWAS with BLUE values from three trials under LP stress (Fig. 2a-f). SNPs with close proximity (within 1 MB) and an LD r2 > 0.2 were grouped together since they belong to the same QTL [40].
In total, 39 marker-trait associations were found with a significant threshold of -log10P-value > 4.30 distributed over seven linkage groups, varying from 1 on chromosome C02 to 12 on chromosome A09 (Additional file 2; Table S1). The important SNPs related to examined traits were displayed on Manhattan and QQ plots (Fig. 2a-f). We determined suggestive trait-SNP associations (3.50 < − log10 P ≤ 4.30) to avoid missing SNPs due to the complex nature of RSA traits and the strict MLM criteria [28]. As a result, 39 significant SNP associations, with 15 significant SNP markers and 31 suggestive trait-SNP associated, were integrated into 11 valid QTL clusters (Tables 3, Additional file 2; Table S1), which included at least two investigated root and biomass-related traits. These SNP markers were detected except PRL, TRL, and RSRD for ten investigated traits. These QTLs were distributed on seven chromosomes; A04, A09, A10, C02, C03, C07, and C09. The highest number of loci were identified on A09 and C03, containing 28 and 19 loci (Additional file 2; Table S1).
Because root and biomass traits demonstrated significant and strong relationships, numerous pleiotropic genetic loci were revealed, particularly QTL clusters RT-A09-1, RT-A09-2, RT-A10-1, RT-C03-2, RT-C03-3, RT-C07-1, and RT-C09-1, which influenced both root development and aboveground biomass formation. Two pleiotropic SNP markers (seq-new-rs41996 and seq-new-rs46512) in the QTL clusters RT-A09-2 and RT-C07-1 revealed high phenotypic variance of 24.43% and 24.04%, respectively, for RFW and were associated with both RMT and BM. These loci that simultaneously influence root and shoot biomass traits could be promising loci for marker-assisted breeding upon validation.
Differential gene expression analysis between high/low phosphorus efficient groups and high/low phosphorus stress tolerance groups
Based on the marked differences in SFW at 7 (T1) and 14 (T2) days after transplanting, five accessions (S57, S178, S203, S233, S283) with higher SFW were classified into the high phosphorus efficient (HP1). Five accessions (S55, S159, S202, S205, S269) with lower SFW were classified as low phosphorus efficient (LP1). In addition, other ten accessions were evaluated based on the phosphorus stress tolerance index (STI) at 14 days under control and LP treatments and were classified into two groups; high phosphorus stress tolerance (HP2, including S126, S178, S203, S219, S262) and low phosphorus stress tolerance (LP2, including S33, S55, S59, S192, S311) (Fig. 3a). As a result, 48 RNA-seq libraries, including three biological replicates of the HP1, LP1, HP2, and LP2 groups under the low phosphorus stress, HP1CK, LP1CK, HP2CK, and LP2CK groups under the control condition at T1 and T2 time points were generated. The total, mapped, and unique reads to the reference B. napus genome are shown (Additional file 2; Table S2). After filtering the adaptor sequences and those with a low base quality, the Illumina sequencing produced 2280.3 million clean reads. At the same time, the average guanine-cytosine (GC) content was 47.4%, and the Phred quality score (Q30) was 88.8%. Using principal component analysis (PCA) and correlation analysis of gene expression levels, it was demonstrated that the correlation between individuals within the same group was much higher than the correlation between individuals within different groups (Additional file 1; Fig. S2a-b), indicating that the three biological repeats used were sufficiently accurate for the experiment.
A pairwise approach was used to determine the differentially expressed genes (DEGs) of the group-I (HP1 and LP1) under control and LP stress, respectively, including HP1 vs. LP1 and HP1CK vs. LP1CK under T1 and T2 points, respectively. Then the common DEGs of Group-I under both two-time points were found using a Venn diagram between HP1/LP1-T1 vs. HP1/LP1-T2 and HP1CK/LP1CK-T1 vs. HP1CK/LP1CK-T2. The same approach was used to screen the common DEGs of group-II (HP2/LP2) that occur under both control and LP stress conditions at T1 and T2 time points (Fig. 3b). The DEGs between these groups were determined using a false discovery rate (FDR) ≤ 0.05, and an absolute value of |log2 (fold change) | was used as the threshold. Under the LP stress, 1453 common DEGs of group-I from the pairwise comparisons of HP1/LP1-T1 vs. HP1/LP1-T2 were identified, including 842 upregulated and 611 downregulated DEGs (Fig. 3c). Under the control condition, 1837 common DEGs of group-I from the pairwise comparisons of HP1CK/LP1CK-T1 vs. HP1CK/LP1CK-T2 were identified (1107 upregulated and 730 downregulated). Furthermore, the Venn diagram showed 761 DEGs (502 upregulated and 259 downregulated) overlapped between HP1/LP1/HP1CK/LP1CK-T1T2 (Fig. 3b-c). This means that there were 692 DEGs specific to the HP1/LP1 group under the LP stress condition, while 1076 DEGs were specific under the control condition (Additional file 2; Table S3).
Similarly, 1508 common DEGs were identified from group-II between HP2/LP2-T1 vs. HP2/LP2-T2 under the LP stress, including 825 upregulated and 683 downregulated DEGs and 1794 DEGs under the control condition (1018 upregulated and 776 downregulated) were identified between HP2CK/LP2CK-T1 vs. HP2CK/LP2CK-T2. Furthermore, 860 DEGs (507 upregulated and 353 downregulated) were regarded as common DEGs for HP2/LP2/HP2CK/LP2CK-T1/T2 (Fig. 3b-c). Considering this, HP2/LP2 showed 648 specific DEGs under the LP stress condition compared to 934 DEGs under the control condition (Additional file 2; Table S3).
DEGs' complete names, FPKM values, and accompanying descriptive information are presented in (Additional file 2; Table S3). Normalized FPKM values ranging from − 1 to 1 were used to generate a heatmap based on expression profile similarity and diversity to classify high/low P-specific DEGs (Additional file 1; Fig. S3a-f). According to the heatmap, the expression patterns of DEGs were clearly shown as clusters of upregulated and downregulated genes. The qRT-PCR results for seven DEGs in all samples were strikingly similar to those from RNA-Seq, demonstrating the accuracy of the RNA-Seq data (Fig. 3d).
Functional classification of DEGs involved in high/low phosphorus efficiency and high/low phosphorus stress tolerance index
To determine the functional significance of the DEGs in each group at two time points (7 and 14DAT), Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) classifications were identified for the DEGs in HP1/LP1-specific, HP1CK/LP1CK-specific, HP1/LP1/HP1CK/LP1CK-common, HP2/LP2-specific, HP2CK/LP2CK-specific, and HP2/LP2/HP2CK/LP2CK-common. Annotated genes were further divided into three major functional classifications; molecular functions (MF), cellular components (CC), and biological processes (BP).
Among the HP1/LP1 DEGs specific under LP stress, 543 (78.46%) of 692 received GO functional annotation. The most over-represented GO terms were displayed in Fig. 4a (according to their p-values). In the MF category, the most abundant GO terms were; benzoic acid carboxyl methyl transferase activity, methyltransferase activity, and carboxylate reductase activity. In the CC category, endoplasmic reticulum membrane and cell, followed by THO complex, were the most significant GO terms. Under the BP category, the DEGs were significantly enriched in; rRNA 5'-end processing, negative regulation of reductive pentose-phosphate cycle, and negative regulation of posttranscriptional gene silencing (Fig. 4a, Additional file 2; Table S4). Similarly, GO function annotations were obtained for 984 of the 1076 (91.44%) DEGs unique to HP1CK/LP1CK. Most transcripts enriched in the MF category were; sulfate transmembrane transporter activity, sinapine esterase activity, glutathione binding, and nutrient reservoir activity. Enrichment for DEGs associated with the CC category was; ribosome, RNA nuclear export complex, followed by proton-transporting two-sector ATPase complex. Maintenance of seed dormancy, sulfate transport, intracellular iron ion sequestering, root meristem identity maintenance, and photosynthesis was highly enriched in the BP category (Fig. 4b, Additional file 2; Table S4). The GO-term enrichment analysis for HP1/LP1/HP1CK/LP1CK-common DEGs also revealed important GO terms (Fig. 4c, Additional file 2; Table S4). The common GO terms were related to; microtubule-severing ATPase activity, transporter activity, phosphatide phosphatase activity, membrane, photosystem I reaction center, cellular chloride ion homeostasis, positive regulation of response to oxidative stress, photosynthesis, hormone-mediated signaling pathway, and carbohydrate transmembrane transport, sucrose synthase activity, activation of MAPK activity and glutathione metabolic process.
Similarly, GO enrichment analysis was performed for HP2 and LP2 at T1 and T2 under both LP stress and control treatments. A total of 648 HP2/LP2-specific, 934 HP2CK/LP2CK-specific, and 860 HP2/LP2/HP2CK/LP2CK-common DEGs were assigned to 203, 186, and 194 significant GO terms, respectively (Additional file 2; Table S4). For HP2/LP2 specific transcripts, the dominant terms were; methyltransferase activity, chitinase activity, L-glutamate transmembrane transporter activity, NADP binding, chloroplast, and RNA polymerase III complex, polysaccharide catabolic process, chitin catabolic process, cellular response to reactive oxygen species, regulation of phosphate transport, lateral root branching, and positive regulation of abscisic acid-activated signaling pathway (Fig. 4d, Additional file 2; Table S4). In group HP2CK/LP2CK-specific, the most enriched GO terms were; chlorophyll-binding, nutrient reservoir activity, photosystem I, photosystem II, photosynthesis, photosystem II assembly, and root hair cell tip growth (Fig. 4e, Additional file 2; Table S4). Moreover, 860 DEGs (Fig. 3b) that overlapped between HP2/LP2 and HP2CK/LP2CK were significantly enriched in different GO terms such as protein methyltransferase activity, NADP binding, proton-transporting ATPase activity, SUMO binding, DNA replication factor A complex, regulation of chlorophyll biosynthetic process, cellular response to glucose stimulus, cellular response to phosphate starvation, and carbohydrate metabolic process (Fig. 4f, Additional file 2; Table S4).
Based on the GO and KEGG results, sucrose synthase activity, DNA-binding transcription factor activity, auxin transport, lipid biosynthetic process, fatty acid degradation, phenylpropanoid biosynthesis, and tryptophan metabolism were unique pathways for low phosphorus uptake. In contrast, glutathione synthase activity, nutrient reservoir activity, sucrose-phosphate synthase activity, cellular response to oxidative stress, and glutathione metabolism were the important specific pathways for high phosphorus uptake. ATPase activity, phosphate binding, photosynthesis, response to oxidative stress, starch biosynthetic process, and selenocompound metabolism were common pathways for low/high phosphorus uptake. Furthermore, on the other hand, amino acid binding, electron transfer activity, NADP binding, lateral root branching, propanoate metabolism, and cellular response to sucrose starvation were the critical pathways for low/high phosphorus stress tolerance. In contrast, chlorophyll-binding, potassium ion binding, translation elongation factor activity, entrainment of the circadian clock, carbohydrate metabolic process, and amino sugar and nucleotide sugar metabolism were the most crucial pathways for low/high phosphors stress tolerance under control conditions. The most important pathways that were common in low/high phosphorus stress tolerance; transcription coactivator activity, SUMO binding, regulation of chlorophyll biosynthetic process, carbohydrate metabolic process, negative regulation of iron ion transport, cyanoamino acid metabolism, and glycerolipid metabolism. Therefore, it can be speculated that the genes associated with these pathways play an important role in regulating the low phosphorus responsiveness of rapeseed.
Candidate Genes Prediction and Prioritization by Integrating GWAS, DEGs, and WGCNA
An analysis of co-expression gene modules from 83,088 annotated genes with p > 0.05 was conducted to investigate the gene regulatory network under low phosphorus stress using WGCNA. In Fig. 5a, the dendrogram identified 21 modules based on gene correlation, and Fig. 5b depicts the relationships between modules and samples. A total of 48,454 genes were determined to be involved in these 21 modules, ranging from 292 in the ‘MEgrey’ module to 10006 in the ‘MEturquoise’ module (Fig. 5c), indicating the complexity of gene regulation during phosphorus stress treatment. The MEtan, MEpurple, MEmidnightblue, and MEcyan modules were found to be highly correlated with HP1/CK-T1T2, LP1/CK-T1T2, HP2/CK-T1T2, and LP2/CK-T1T2, respectively (Fig. 5b). Based on the heatmaps, genes within one module were significantly expressed in samples that were strongly correlated with that module (Additional file 1; Fig. S4a-d).
The GO and KEGG analysis suggested that the significant GO terms of genes in the MEbrown module were related to glutamate-cysteine ligase activity, phosphatidate phosphatase activity, auxin binding, ATPase activator activity, response to ethylene, regulation of the salicylic acid metabolic process, cytokinin-activated signaling pathway, lateral root branching and formation and post-embryonic root development (Additional file 2; Table S5). Pyruvate metabolism, riboflavin metabolism, propanoate metabolism, N-Glycan biosynthesis, and inositol phosphate metabolism were significantly enriched in the “MEbrown” module by KEGG pathway enrichment analysis (Additional file 2; Table S6). These critical pathways were crucial in the phosphorus metabolism and assimilation process.
The LD approach was used to identify candidate genes within 300 kb upstream and downstream of the significant peak SNPs associated with the traits studied [41, 42]. Based on GWAS results, 1060 genes were identified around each peak SNP in 11 QTL clusters located within the 300 kb region (Additional file 2; Table S7). We explored five potential genes that overlapped with GWAS and WGCNA genes based on the substantial and consistent correlation between WGCNA genes and each module (Table 4, Fig. 7). These five hub genes were found within MEbrown having high and consistent correlation across the samples under both phosphorus stress and control conditions. In addition, GWAS and DEG integration allowed us to identify seven genes simultaneously as common candidate genes. These candidate genes were detected within seven QTL clusters; BnaA04g23490D and BnaA04g23610D (RT-A04-1), BnaA09g04320D, BnaA09g04350D and BnaA09g04930D (RT-A09-1), BnaA09g08440D, BnaA09g08940D and BnaA09g09290D (RT-A09-2), BnaC02g23620D (RT-C02-1), BnaC03g41900D (RT-C03-1), BnaC07g30970D (RT-C07-1) and BnaC09g18800D (RT-C09-1). Within these 12 genes, six genes (BnaA04g23490D, BnaA09g08440D, BnaA09g04320D, BnaA09g04350D, BnaA09g04930D, BnaA09g09290D) that showed differential expression were identified previously related to phosphorus use efficiency, phosphorus utilization, assimilation, and root growth and development. These findings evaluated the efficacy of an integrated GWAS, WGCNA, and differential expression analysis strategy for screening potential genes.