Integrating genome-wide association study with transcriptomic data to predict candidate genes influencing Brassica napus root and biomass-related traits under low phosphorus conditions

Abstract

Background: Rapeseed (Brassica napus L.) is an essential source of edible oil and livestock feed, as well as a promising source of biofuel. Breeding crops with an ideal root system architecture (RSA) for high phosphorus use efficiency (PUE) is an effective way to reduce the use of phosphate fertilizers. However, the genetic mechanisms that underpin PUE in rapeseed remain elusive. To address this, we conducted a genome-wide association study (GWAS) in 327 rapeseed accessions to elucidate the genetic variability of 13 root and biomass traits under low phosphorus (LP; 0.01 mM P +). Furthermore, RNA-sequencing was performed in root among high/low phosphorus efficient groups (HP1/LP1) and high/low phosphorus stress tolerance groups (HP2/LP2) at two-time points under control and P-stress conditions.

Results: Significant variations were observed in all measured traits, with heritabilities ranging from 0.47 to 0.72, and significant correlations were found between most of the traits. There were 39 significant trait-SNP associations and 31 suggestive associations, which integrated into 11 valid quantitative trait loci (QTL) clusters, explaining 4.24-24.43% of the phenotypic variance observed. In total, RNA-seq identified 692, 1076, 648, and 934 differentially expressed genes (DEGs) specific to HP1/LP1 and HP2/LP2 under P-stress and control conditions, respectively, while 761 and 860 DEGs common for HP1/LP1 and HP2/LP2 under both conditions. An integrated approach of GWAS, weighted co-expression network, and differential expression analysis identified 12 genes associated with root growth and development under LP stress. In this study, six genes (BnaA04g23490D, BnaA09g08440D, BnaA09g04320D, BnaA09g04350D, BnaA09g04930D, BnaA09g09290D) that showed differential expression were identified as promising candidate genes for the target traits.

Conclusion: 11 QTL clusters and 12 candidate genes associated with root and development under LP stress were identified in this study. Our study's phenotypic and genetic information may be exploited for genetic improvement of root traits to increase PUE in rapeseed.

Keywords: Biomass traits; GWAS; PUE; Phosphorus; QTL; Rapeseed.

Fig. 1

Correlation analysis of the investigated traits. a Correlations of studied traits under low phosphorus stress. The frequency distribution for each trait was displayed on the diagonal. The upper and lower parts represent the correlation coefficient and scatter plots between two diagonal traits, respectively. b Correlations of each investigated trait between control and low phosphorus stress. Red and blue indicate positive and negative correlations, respectively. ***, ** and * denote significance at the 0.1%, 1% and 5% levels of probability, respectively

Fig. 2

a–c Manhattan plots of the phenotype–genotype association analysis for 13 root and shoot biomass traits of B. napus by MLM with BLUE values. The x-axis displays the chromosome label, and the y-axis displays −log10 (p-value). The solid gray line shows significant associations between SNPs and phenotype value with a threshold level of the p-value (−log10 1/20, 131 = 4.30 × 10−5). The colour dots above the threshold values indicate the significant SNPs for root and shoot biomass traits. d–f QQ plots represent MLM analysis of the investigated traits

Fig. 3

Differential gene expression analysis. a Phenotypic performance of high/low phosphorus efficient groups and high/low phosphorus stress tolerance groups. b Venn diagram of the DEGs in the selected groups. c Upand downregulated DEGs in different groups. d Correlation between qRT-PCR and RNA-seq data. ** and * denote significance at the 1% and 5% levels of probability, respectively. ns, not significant

Fig. 4

Gene ontology (GO) analysis of differentially expressed genes. a–c GO terms correspond to HP1/LP1-specific, HP1CK/LP1CK-specific, HP1/LP1/HP1CK/LP1CK-common, HP2/LP2-specific, HP2CK/LP2CK-specific, and HP2/LP2/HP2CK/LP2CK-common, respectively. d–f GO terms correspond to HP2/LP2-specific, HP2CK/LP2CK-specific, and HP2/LP2/HP2CK/LP2CK-common, respectively. The Y-axis is −log10 (p-value). The relevant p-value decreases as the bar chart height increases. Red, blue, and green colours correspond to molecular function, cellular components, and biological processes

Fig. 5

WGCNA of gene expression matrix. a Gene-based co-expression network analysis dendrogram. b Module–sample association; each row represents a module labelled with the same colour as in (a), and each column represents a sample. c Overview of identified genes corresponds to each module

Fig. 6

Correlation of networks in MEbrown module. The yellow colour in the network indicates the candidate genes overlapped by GWAS and WGCNA