- Open Access
Comparative transcriptomic analysis reveals the molecular mechanism underlying seedling biomass heterosis in Brassica napus
BMC Plant Biology volume 22, Article number: 283 (2022)
Heterosis is an important biological phenomenon in which the hybrids exceed the parents in many traits. However, the molecular mechanism underlying seedling heterosis remains unclear.
In the present study, we analyzed the leaf transcriptomes of strong hybrids (AM, HM) and weak hybrids (CM, HW) and their parents (A, C, H, M, and W) at two periods. Phenotypically, hybrids had obvious biomass heterosis at the seedling stage, with statistically significant differences between the strong and weak hybrids. The transcriptomic analysis demonstrated that the number of differentially expressed genes (DEGs) between parents was the highest. Further analysis showed that most DEGs were biased toward parental expression. The biological processes of the two periods were significantly enriched in the plant hormone signal transduction and photosynthetic pathways. In the plant hormone signaling pathway, DEG expression was high in hybrids, with expression differences between strong and weak hybrids. In addition, DEGs related to cell size were identified. Similar changes were observed during photosynthesis. The enhanced leaf area of hybrids generated an increase in photosynthetic products, which was consistent with the phenotype of the biomass. Weighted gene co-expression network analysis of different hybrids and parents revealed that hub genes in vigorous hybrid were mainly enriched in the plant hormone signal transduction and regulation of plant hormones.
Plant hormone signaling and photosynthesis pathways, as well as differential expression of plant cell size-related genes, jointly regulate the dynamic changes between strong and weak hybrids and the generation of seedling-stage heterosis. This study may elucidate the molecular mechanism underlying early biomass heterosis and help enhance canola yield.
Heterosis refers to the phenomenon that hybrids perform better than their parents in various traits, such as biomass and resistance , and it has been successfully applied in field crops. It is a complicated phenomenon that may be coordinated by one or numerous genes and has always attracted the attention of breeders and molecular biologists. In the past century, various hypotheses have been proposed, including the dominance hypothesis , overdominance hypothesis , and epistasis hypothesis . These hypotheses are limited to explaining the genetic basis of heterosis at the DNA level but do not reveal the biological basis for its occurrence. Furthermore, the mechanism underlying phenotypic heterosis is not completely understood.
The rapid development of sequencing technology has made it possible to explore heterosis at the molecular level , and it has been used in many species, such as Arabidopsis [6–8], rice [9–11], maize [12, 13], and wheat [14, 15]. These studies combined morphological and omics techniques to analyze the changes in the vegetative and reproductive growth phases of hybrids relative to their parents. Nonetheless, these studies did not determine the changes in important genes that participated in the adjustment of the biological pathways of biomass heterosis.
The seedling stage is a critical period for plants to accumulate nutrients and grow, and the coordinated source-sink-flow relationship contributes to the formation of yield at the mature stage . Researchers have focused on the molecular mechanisms underlying early heterosis [17–20]. Leaves are important photosynthetic organs and are critical for biomass heterosis . Transcriptomic analyses of Arabidopsis hybrids and their parents showed that the increase in the expression of indole-3-acetic acid (IAA)-targeted genes and the decrease in the expression of salicylic acid biosynthetic pathway-related genes in the hybrids were consistent with the phenotype of increased leaf cell size . Transcriptomic analyses of canola hybrids and their parents four and eight days after sowing (DAS) revealed that the early biomass heterosis of seedlings in hybrids was related to the photosynthetic pathway and auxin biosynthetic pathway genes . Liu et al. studied the temporal dynamic transcriptome of Col-0 and Per-1 and their hybrid seedlings during early development and found that the dominant expression complementation of the basic biological pathways helps the formation of heterosis . However, most studies only focus on the heterosis of one or two hybrids and their parents, ignoring the comparative study of multiple pairs of hybrids and parents at different time points.
The utilization of canola (Brassica napus) heterosis effectively enhances the yield of edible vegetable oil. Canola and the model plant genus Arabidopsis belong to the family Brassicaceae. B. napus is developed through an interspecific cross between B. rapa (AA) and B. oleracea (CC). It is a young crop that originated approximately 7,500 years ago . Moreover, expression level dominance (ELD) has been widely reported in allopolyploid plants [14, 20]. Transcriptomic, small RNA, and methylation sequencing analyses have systematically revealed the epigenetic mechanism underlying the heterosis of yield traits in canola , research on canola (B. napus) heterosis is still limited [24, 24, 25].
In the present study, we analyzed the leaf transcriptomes of AM, CM, HM, and HW hybrids and their parents (A, C, H, M, and W) at 21 and 24 DAS. This study demonstrates that affinities and distinctions in gene expression changes occur between strong and weak hybrids and their parents in the early stages of seedling development. A comparative analysis of the whole genome transcriptome comparison will help us to better understand the mechanism underlying gene expression changes in the hybrid and provide new ideas for hybrid breeding in canola.
Phenotypic analysis of hybrids and parents
To better compare the biomass heterosis between the different hybrids, we calculated the mid-parent heterosis (MPH) and high-parent heterosis (HPH) of traits and found that the MPH (83.77%) and HPH (68.75%) of the AM combination (dry weight) were significantly higher than those of the parents at 21 DAS (Table S1). The MPH and HPH ranges in the dry weight were 18.64 – 83.77% and 6.98 – 68.75%, respectively, at 21 DAS (Table S1).
We observed that the total leaf area of all hybrids differed significantly from the corresponding mid-parent value (MPV) at 21 and 24 DAS (Fig. 1A, B; Fig. S1). Furthermore, there were significant differences in the total leaf area of strong hybrid AM relative to CM and the advantageous hybrid HM relative to HW in the two periods (Fig. 1A, B). Interestingly, when we analyzed the fresh weights of the two periods, we observed the similar changes mentioned above (Fig. 1C, D; Fig. S1; Tables S1, S2). The results indicated that AM and HM were vigorous hybrids. The results also showed that the two pairs of hybrids and parents with obvious differences were suitable for the comparative transcriptomic analysis and regulation mechanism study on biomass heterosis at the seedling stage of B. napus.
Identification and analysis of differentially expressed genes (DEGs)
The range of clean reads obtained from all sequencing samples was 33,395,822 to 65,055,276. The clean data of each sample reached 6.52 Gb on average. The range of clean bases for all sequencing samples was 5,009,373,300 to 9,758,291,400. The average Q20 value of our samples was 97.29% and the average GC content was 47.82%. Among the clean reads, the percentage of reads uniquely mapped to the reference genome averaged 79.26% (Table S3). The Pearson correlation coefficient between biological replicates was high, with an average of 98.69% (Fig. S2). This indicates that our transcriptome data are credible.
Herein, we analyzed the dynamic changes in the transcriptome of two pairs of strong and weak hybrids (AM-CM and HM-HW) and their parents. Venn diagram analysis conducted on the DEGs of 21 DAS and 24 DAS hybrids apropos of their parents (Fig. 2A, B). For example, 21 DAS AMvsA, AMvsM, and AvsM had 2111, 3,285, and 4702 unique DEGs, respectively (Fig. 2A). In addition, 24 DAS AMvsA, AMvsM, and AvsM had 668, 3105, and 4706 unique DEGs, respectively (Fig. 2B). The numbers of upregulated and downregulated genes between F1 hybrids and their parents at two periods were further analyzed (Fig. 2C, Tables S4, S5). The results showed that at 21 DAS, the upregulated and downregulated genes in HMvsH, HMvsM, and HvsM were 3492 and 1715, 6247 and 3726, and 9941 and 8438, respectively (Fig. 2C, Tables S4, S5). The total DEGs of HMvsH and HMvsM were 3437 and 13,811 at 24 DAS, respectively (Fig. 2C, Tables S4, S5).
Interestingly, we discovered that the number of upregulated genes in most hybrids was markedly higher than that of downregulated genes (Fig. 2C, Tables S4, S5). The number of DEGs between parents was the highest. The number of DEGs also increased as the plants grew, and more drastic transcriptomic changes were observed at 24 than at 21 DAS. The number of DEGs between all hybrids and the male parents was less than that between the hybrids and female parents (Fig. 2C).
Analysis of the expression level of hybrids and their parents
To explore the possible changes and directions, we divided the DEGs into 12 possible groups based on previous studies (Fig. 3A) [26, 27]. Many ELD genes were identified at 21 and 24 DAS. Among them, additivity, transgressive downregulation, and transgressive upregulation genes accounted for only a small part of the total (Fig. 3B, C). Most of the gene expression in the hybrids was similar to that of the parents (Fig. 3B, C, Fig. S3, Tables S6, S7). In the two periods examined, ELD male (ELD-M) and ELD female (ELD-F) were classified as having high expression in all hybrids. Among the expressed genes of ELD classification, the parental-ELD genes at 21 and 24 DAS accounted for an average of up to 85% and 90%, respectively (Fig. 3B, C). The genes classified by high parental-ELD in the two periods had a high proportion of 58% and 62%, respectively (Fig. S3). In conclusion, the analysis of expression patterns showed that ELD was ubiquitous. Parental ELD had a strong positive correlation with seedling biomass. In particular, the dominant effect was the main contributor to the leaf growth of hybrids at the seedling stage.
Gene ontology (GO) enrichment analysis of parental-ELD (P-ELD) gene
To explore the biological function of the P-ELD gene, the four types of genes (II + XI + IV + IX) in each F1 hybrid in the two periods were analyzed by GO enrichment analysis (Tables S8, S9). At 21 DAS, the molecular functions of the P-ELD gene of the hybrids were significantly enriched in the structural components of RNA binding and ribosomes when compared to those of the parents (Fig. S4A, Table S8). Through cell component enrichment analysis, we found that they were significantly enriched in chloroplasts and the chloroplast matrix. (Fig. S4B, Table S8). At 24 DAS, the molecular functions of the hybrids relative to the P-ELD gene in both parents were significantly enriched in NAD, NADP binding, and oxidoreductase (Fig. S5A, Table S9). The cell components were significantly enriched in the chloroplast, chloroplast matrix, and photosynthetic membrane (Fig. S5B, Table S9).
We focused on the biological process of the significant enrichment of the P-ELD genes of the four hybrids in the two periods (Tables S8, S9). The Venn diagram showed that the two periods had a common and unique biological pathway (BP) term. The number of significant BP terms in the T2 period (24 DAS) was much greater than that in the T1 period (21 DAS) (Fig. S6A, B). Interestingly, the significantly enriched biological pathways were mainly enriched in photosynthesis (GO:0,015,979), carbon metabolism (GO:0,005,975), plant hormone response (GO:0,009,725), and other pathways at 21 and 24 DAS (Fig. S6C, D).
Comparative analysis of strong and weak hybrids
To explore the changes in strong and weak hybrids, we analyzed the differential expression of two pairs of strong and weak hybrids (AM-CM and HM-HW) in two periods. The two pairs of strong and weak hybrids had 875 DEGs in common at 21 DAS (Fig. 4A, Table S10). The GO enrichment inquiry the common DEGs disclosed that they primarily benefited from biological processes, such as response to hormones, regulation of cell size, and carbohydrate metabolism (Fig. 4B, Table S11). The Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis of shared DEGs showed they were mainly enriched in the carbohydrate metabolism and signal transduction pathways (Fig. 4C, Table S12).
We performed the same analysis at 24 DAS. A total of 563 DEGs were shared by the two pairs of strong and weak hybrids (Fig. S7A, Table S13). The GO enrichment analysis of these shared DEGs demonstrated that they were largely enriched in the biological processes of hormone metabolism, circadian rhythm, and carbohydrate metabolism (Fig. S7B, Table S14). The KEGG enrichment analysis performed on the shared DEGs showed that they were mainly enriched in pathways such as carbohydrate metabolism and starch and sucrose metabolism (Fig. S7C, Table S15). These outcomes suggest that the contrasting expression of these pathways can lead to differences between strong and weak hybrids.
Analysis of plant hormone transduction between hybrids and parents
According to the enriched biological process, the plant hormone signal transduction pathway was significantly enriched (Fig. S6C, D). We focused on genes related to the auxin and cytokinin transduction pathways (Tables S8, S9). The regulatory pathways are shown in Fig. 5a and b. Significant differences in small auxin upregulated RNA (SAUR) were observed in the auxin signaling pathway, including SAUR4, SAUR14, SAUR15, SAUR19, SAUR20, SAUR21, SAUR23, SAUR33, and SAUR51 (Fig. 5A). Interestingly, based on our analysis, up to 10 SAUR21 genes (BnaA03g06810D, BnaA10g16790D, BnaA10g16800D, BnaC02g07800D, BnaC03g08730D, BnaC03g08740D, BnaC03g08750D, BnaC03g08780D, BnaC09g39830D, and BnaC09g52070D) in the pathway were significantly expressed in the hybrids, and they were expressed in high amounts relative to the MPV.
We compared the alterations in the gene expression of the four hybrids (AM, CM, HM, and HW) relative to MPV in two periods (T1 and T2) and observed significant differential expression of the genes of B-Arabidopsis Response Regulators (ARR) and A-ARR (Fig. 5B). Notably, the B-ARR gene was predominantly expressed in AM and HM hybrids relative to MPV; the same gene was more highly expressed in AM and HM hybrids than in CM and HW hybrids. Similar changes were observed in A-ARR genes. These results indicate that genes in the signal transduction pathway of auxin and cytokinin may change leaf expansion and eventually evoke changes in hybrid rape biomass vigor.
Changes in cell size-related gene expression and leaf phenotypes in hybrids
Many of the DEGs of hybrids promoted an increase in the cell size (Fig. 6A). Compared with those in the CM and HW hybrids, the genes that promoted cell size in the strong hybrids AM and HM were mainly upregulated and more highly expressed (Fig. 6A). Hybrids with heterosis at the seedling stage were likely to have either or both increased leaf length and leaf width. Consistent with the results of biomass, the AM and HM hybrids had significant heterosis in the second true leaf area relative to the MPV and had significant differences between AMvsCM and HMvsHW (Fig. 6B). The AM, CM, and HM hybrids showed an increase in leaf length and width correlated with the MPV (Fig. 6C, D). We found that the contribution of leaf width to the leaf area of the enlarged second true leaf was greater than that of the leaf length (Fig. 6E). By performing the same analysis as above on the third true leaf at 24 DAS (Fig. S8, S9), we discovered that genes related to leaf cell size, such as BXL1 (BnaA02g30740D and BnaC02g39030D), LNG2 (BnaC03g32810D), and PIF4 (BnaA03g19970D, BnaC04g48630D, BnaA04g24760D, and BnaC01g40470D) were highly expressed in strong dominant hybrids (AM and HM) and significantly different in the CM and HW hybrids (Fig. S8). Similar results were obtained regarding the phenotypic data of the third true leaf area, leaf length, and breadth (Fig. S9A, B, C). Except for that to the HW hybrid, the contribution of leaf width to the leaf area of other hybrids was greater than that of leaf length (Fig. S9D).
Comparison among genes in the photosynthetic pathways of hybrids
According to previous reports, photosynthesis is very important for plant heterosis [7, 12]. The significant enrichment process of many P-ELD genes was related to photosynthesis (Fig. 7A), and several DEGs constituted photosystem I, photosystem II, and light-harvesting complex (LHC) (Fig. 7B, C, D; Tables S8, S9). We analyzed the changes in DEGs related to the photosynthesis of hybrids and their corresponding MPVs in the T1 and T2 periods (Fig. 7B, C, D). For LHC, we observed that the genes of BnaA09g02120D (LHCB2), BnaC09g01520D (LHCB2), BnaA05g05540D (LHCB4.3), BnaC04g05260D (LHCB4.3), and BnaA08g22200D (LHCA6) showed upregulated expression relative to the MPV at 21 and 24 DAS. AM had higher expression levels of these genes than CM, and HM had higher expression levels of these genes than HW (Fig. 7B). Moreover, the PSII and PSI genes showed similar expression trends (Fig. 7C, D). The results showed that the upregulated expression of genes associated with photosynthesis contributes to the production of biomass heterosis.
Correlation network analysis of leaf development in hybrids and parents at the seedling stage
Results suggested that the DEGs shared by all hybrids and parents are closely related to the biomass phenotype in the two periods. Consequently, weighted gene co-expression network analysis (WGCNA) was executed for related genes in all samples to identify the key genes and underlying mechanism highly related to the biomass phenotype. Ultimately, 17 modules were identified, as shown in the dendrogram (Fig. S10). The eigengenes of the 17 modules were associated with different samples (Fig. S11).
This study focused specifically on the module closely related to leaf development in the strong heterosis hybrid. Therefore, we associated modules with different samples. We discovered that the yellow module was markedly correlated with biomass (Fig. S11). The yellow module contained 487 genes and highly expressed transcripts in the vigorous hybrid (AM) (Fig. S12, Table S16). KEGG enrichment analysis was conducted on the genes of this module, which were primarily enriched in plant hormone signal transduction and various metabolic processes (Fig. S13, Table S17).
To further analyze the yellow module, we identified 53 hub genes (membership degree value [KME] > 0.90, P < 10–6) (Table S18). GO enrichment analysis indicated that these hub genes were overexpressed in functional categories related to the regulation of plant growth and development, such as hormone-mediated signal transduction, the response to auxin, and regulation of the hormone level (Fig. S14, Table S19). These results support the role of the hub genes in early leaf growth and development.
qRT-PCR verification of the transcriptome
We randomly selected nine candidate genes in the transcriptomic data at 21 and 24 DAS to perform qRT-PCR experiments to certify the reliability of the transcriptomic outcomes. Nine genes, including BnaC07g24660D, BnaA09g49810D, BnaC06g28840D, and BnaA09g48920D, were used for quantitative analysis. Detailed information on the primers used is shown in Table S20. The quantitative results were highly consistent with the trends of the transcriptomic data, further confirming the accuracy of transcriptomic outcomes (Fig. 8 and Fig. S15).
Heterosis is a widely known event wherein hybrids display vigor greater than that of their parents [11, 28]. The vegetative growth of plant seedlings is critical to the increase in yield at a later time, and the advantage of the seedling stage is to eliminate environmental interference and better reflect the mechanism underlying heterosis. In this study, we used comparative transcriptomic analysis of the leaves of two pairs of hybrids and their parents (AM-CM and HM-HW) at two stages and explored the molecular mechanism underlying biomass heterosis at the seedling stage.
Interestingly, we measured the leaf length and width of the second and third true leaves at the same time point at 20–30 DAS and found that the days with the fastest daily growth rate of the second and third true leaves were 21 DAS and 24 DAS, respectively (Fig. S16). Therefore, we analyzed the transcriptomic data of the second and third true leaves in two sampling periods from four hybrids (AM, CM, HM, and HW) and their parents (A, C, H, M, and W). Herein, we propose a model in which we observe the differential expression of genes in the photosynthetic pathway and the plant hormone signal transduction pathway between hybrids and parents. The changes in the cell size-related genes led to changes in leaf length and leaf width and ultimately affected the differences in biomass heterosis between strong and weak hybrids, as well as between the hybrids and their parents (Fig. 9).
Possible roles of hormone transduction pathway genes in seedling heterosis
In the present study, two strong rape hybrid combinations (AM and HM) in the auxin and cytokinin signal transduction pathways (especially in the former) exhibited many genes with upregulated expression relative to their parents. Recent studies have provided evidence for a correlation between auxin and cytokinin and heterosis. The high expression of auxin-related genes has played a positive role in promoting heterosis in crops such as Arabidopsis , maize , and cotton . These reports further support our conclusions. We detected the differential expression of the SAUR gene in the auxin signal transduction pathway. Consistent with our results, Huang et al. who analyzed the size of Euryale ferox Salisb hybrid seeds, revealed that the SAUR gene can affect the distribution of auxin and play a potential regulatory role in hybrid seed size .
Cytokinins play a crucial role in regulating plant growth and development. They can promote the proliferation and expansion of leaf cell development and coordinate with auxins to regulate plant growth . In our study, two genes related to the response regulation of Arabidopsis were upregulated in strongly dominant hybrids. These genes are thought to be related to an increase in biomass. However, the factors that upregulate auxin and cytokinin signal transduction pathway genes in hybrids remain unclear. The differential expression of different genes in plant IAA and cytokinin signal transduction pathways may lead to different levels of hybrid seedling vigor. This indicates the high expression of genes related to auxin and cytokinin signaling pathways in hybrids contributes to the occurrence of heterosis in canola at the seedling stage.
Role of cell size-related genes in plant growth and development
Reports that focus on the relationship between cell size-related genes and heterosis are of particular interest. In poplars, the advantage of leaf size in triploid hybrids is positively correlated with an increase in the expression of genes regulating cell size . In Arabidopsis, changes related to genes of leaf size during hybrid rosette development have been observed  as well as the upregulated expression of xyloglucan endotransglucosylase/hydrolase-related genes in the three hybrid systems . However, the factors that drive the upregulation of cell size-related genes in hybrids remain unclear. In our seedling transcriptomic data, genes that promote cell size were highly expressed in rape hybrids relative to their parents. In addition, we provided further physiological evidence that hybrids have more advantages in terms of leaf area than their parents. The leaf length and leaf width of the hybrids and parents were further measured. The difference between the strong and weak hybrids of the same male parent and different female parents may be caused by maternal effects in the early stages of development. This effect on heterosis may indicate there are various cytoplasmic interactions between different parents.
Photosynthesis plays a crucial role in early biomass hybrid vigor
Photosynthesis is an essential basic approach, which supplies plants with energy and helps balance the energy distribution . In Brassicaceae hybrids, the early activation of related genes in the photosynthetic pathway is universal. In Arabidopsis, photosynthesis significantly promotes an increase in the heterosis of hybrid biomass [7, 21, 37]. In Chinese cabbage hybrid seedlings, high parental expression of related genes in hybrids during photosynthesis, such as in photosystem I and photosystem II, has been detected [7, 38]. Our results show that in rape hybrids, the photosynthesis genes were significantly expressed in strong dominant biomass hybrids at 21 and 24 DAS relative to their parents. In addition, we compared the gene expression changes in the photosynthetic pathway between strong and weak hybrids. Based on our results, it is not yet possible to explain whether the contribution of the female or male parent to the hybrid is superior, and further research is required. Subsequent methylation sequencing can be performed on hybrids and parents to explore whether the promoter region of the core gene of interest is regulated by methylation or demethylation. This affects the differential expression of core genes, which may lead to the dominant phenotype of the hybrid.
We analyzed the transcriptomic data of different strong and weak hybrids and their parents at multiple time points to illustrate the important role of dominant genes with high parental expression levels, which are mainly involved in the formation of biomass heterosis in B. napus seedlings through plant hormone signal transduction and the photosynthetic pathway. However, our results do not negate other possible situations that affect seedling heterosis. Future studies should explore the mechanisms through which the genes of related pathways promote heterosis production. This study provides novel insights into biomass heterosis during the vegetative stage and may help provide theoretical support for the cross-breeding of canola.
Plant samples and growth situations
The seeds used in this experiment were provided by the rape Laboratory of Huazhong Agricultural University and created by Dr. Zhao [39, 1]. Through phenotypic analysis of 30 DAS biomass heterosis , we screened strong and weak hybrid combinations with the same parents. Finally, two pairs of strong and weak hybrids (AMvsCM and HMvsHW) and their five inbred parents (A, C, H, M, and W) were used to study biomass heterosis. M and W were female parents; A, C, and H were male parents; and AM, CM, HM, and HW were F40 hybrids. Among them, the female parent was sterile and male parent was fertile. All materials germinated in a greenhouse in ordered surroundings (temperature: 21–23 °C, humidity: 30–60%, 16 h light/8 h darkness) under the support of sterile gauze. Seven days later, seedlings with consistent growth were selected for hydroponic growth . During this period, the nutrient solution was changed every seven days. After two weeks, the second and third true leaves of all experimental materials were collected at 21 and 24 DAS, respectively. After the samples were collected, they were instantly placed in liquid nitrogen and transferred to -80 °C for storage. Three biological replicates were obtained for each sampling, and each biological replicate contained at least 8–10 samples of the same material. F42 hybrids and the parents used in the experiment were placed in the same pot.
Measurement and analysis of phenotypic traits
The fresh weight, total leaf area, and dry weight of the parents and hybrids were measured at 21 and 24 DAS. Three biological replicates were used, each containing at least 8–10 seedlings. At 21 and 24 DAS, the leaf length, leaf width, and leaf area of the second (third) true leaves of the parents and hybrids were measured. The dry weight was measured by initially measuring the fresh weight; next, the samples were placed in a craft paper bag, heated at 105 °C for 30 min, and then dried at 80 °C to an invariant weight. The fresh and dry weights were measured using an electronic balance. Leaf area was evaluated using ImageJ (http://rsb.info.nih.gov/ij/).
Total RNA was obtained from the materials using the TRIzol method. After total RNA collection and DNase I therapy, RNA degradation and contamination were detected on a 1% agarose gel solution. NanoDrop (Thermo Fisher Scientific, Inc., Waltham, MA, USA) and Agilent 2100 (Agilent Technologies, Santa Clara, CA, USA) were used to verify the purity and completeness of the RNA. Library building and RNA-seq were conducted by the BGI Biotechnology Company (Wuhan, China) under the manufacturer’s guidance. All libraries were sequenced to achieve 150 nucleotide paired-end reads on a HiSeq 4000 platform (Illumina, San Diego, CA, USA).
Analysis of DEGs
Raw data from the Illumina HiSeq 4000 sequencer was collated to obtain clean data by eliminating reads comprising adapters (N > 10%) or reads with low-quality nucleotides (> 50%) using fastp v0.19.4 software . The clean data were aligned to the “Darmor-bzh” reference genomes  using the HISAT2 v2.1.0 software . Finally, the uniquely mapped reads were obtained. Subsequently, the reads mapped to genes or exons were applied to count the gene expression level using FeatureCounts v1.5.0 software and illustrated with transcripts per kilobase of exon model per million mapped reads (TPM) . Genes with > 30 counts from the sum of the three biological replicates were considered expressed genes. The DESeq2 R package was used to analyze DEGs . The DEGs were selected with a false discovery rate (FDR) < 0.05 and |log2(fold change)|> 1. We also calculated Pearson correlation coefficients between pairwise samples with the TPM of all expressed genes.
Weighted gene co-expression network analysis
WGCNA was performed on differentially expressed genes to explore co-expression modules that were highly correlated with biomass. Using the correlation network of the WGCNA R package, correlation analysis between highly correlated DEGs and samples was carried out to identify modules with high specific expression . The first principal component of the expression profile of each module was calculated and expressed based on the module eigengenes, which were used to estimate the association between the module and the biomass phenotype of different samples. Moreover, the KME of each gene in the module was used to mine the hub gene.
GO and KEGG analysis
The complete GO annotation information of B. napus was annotated according to the description by Wu et al. . GO and KEGG  enrichment investigation of DEGs was executed by the Tbtools v1.098689 software . The functional classes with a corrected p-value < 0.05 were subsequently analyzed.
Quantitative analysis using qRT-PCR
To assess the reliability of the transcriptomic data, we randomly selected nine genes in two periods for quantification. PerlPrimer software was deployed by constructing gene-specific quantitative primers. The original sequencing sample was used for quantitative analysis. First-strand cDNA synthesis was performed using DNase I management reagent (DP441; Tiangen, Beijing, China) and a first-strand cDNA synthesis kit (K1622; Thermo Fisher Scientific, Inc.). qRT-PCR was performed using SYBR Green Real-Time PCR Master Mix (Toyobo, Osaka, Japan) and a CFX96 real-time system (Bio-Rad, Hercules, CA, USA). UBC10 was used as a housekeeping gene. The 2−ΔΔCT method was used to analyze the data . All reactions were performed for three biological replicates, and every biological replicate contained three technical replicates.
MPH and HPH were estimated using the following formula: MPH = (F1-MP)/MP × 100%, HPH = (F1-HP)/HP × 100%, where MP is defined as the average performance of two parents, and HP is defined as the performance of the best parent. GraphPad Prism 8 (GraphPad Software, Inc.) was used to analyze the significant differences between hybrid and mid-parent values and strong and weak hybrids using unpaired t-tests.
Availability of data and materials
These sequence data reported in this paper have been deposited in the Genome Sequence Archive (https://ngdc.cncb.ac.cn/gsa/) with accession number CRA005745.
Differentially expressed genes
Days after sowing
Transcripts per million
False discovery rate
Kyoto Encyclopedia of Genes and Genomes
Small auxin-upregulated RNA
Type B-Arabidopsis response regulator
Type A-Arabidopsis response regulator
Weighted gene co-expression network analysis
Quantitative real-time polymerase chain reaction
Shull GH. What is heterosis. Genetics. 1948;33(5):439–46.
Davenport CB. Degeneration, albinism and inbreeding. Science. 1908;28:454–5.
East EM. Heterosis. Genetics. 1936;21(4):375–97.
Powers L. An expansion of Jones’s theory for the explanation of heterosis. Am Nat. 1944;78:275–80.
Metzker ML. Applications of next-generation sequencing sequencing technologies - the next generation. Nat Rev Genet. 2010;11(1):31–46.
Botet R, Keurentjes JJB. The role of transcriptional regulation in hybrid vigor. Front Plant Sci 2020; 11(410).
Liu W, He G, Deng XW. Biological pathway expression complementation contributes to biomass heterosis in Arabidopsis. P Natl Acad Sci Usa 2021; 118(e202327811816).
Wang L, Wu LM, Greaves IK, Zhu A, Dennis ES, Peacock WJ. PIF4-controlled auxin pathway contributes to hybrid vigor in Arabidopsis thaliana. P Natl Acad Sci Usa. 2017;114(17):E3555–62.
Huang X, Yang S, Gong J, et al. Genomic architecture of heterosis for yield traits in rice. Nature. 2016;537(7622):629.
Fujimoto R, Uezone K, Ishikura S, Osabe K, Peacoc WJ, Dennis ES. Recent research on the mechanism of heterosis is important for crop and vegetable breeding systems. Breeding Sci. 2018;68(2):145–58.
Liu J, Li M, Zhang Q, Wei X, Huang X. Exploring the molecular basis of heterosis for plant breeding. J Integr Plant Biol. 2020;62(3):287–98.
Li Z, Zhu A, Song Q, Chen HY, Harmon FG, Chen ZJ. Temporal regulation of the metabolome and proteome in photosynthetic and photorespiratory pathways contributes to maize heterosis. Plant Cell. 2020;32(12):3706–22.
Luo J, Wang M, Jia G, He Y. Transcriptome-wide analysis of epitranscriptome and translational efficiency associated with heterosis in maize. J Exp Bot. 2021;72(8):2933–46.
Li A, Liu D, Wu J, et al. mRNA and small RNA transcriptomes reveal insights into dynamic homoeolog regulation of allopolyploid heterosis in nascent hexaploid wheat. Plant Cell. 2014;26(5):1878–900.
Liu YJ, Gao SQ, Tang YM, et al. Transcriptome analysis of wheat seedling and spike tissues in the hybrid Jingmai 8 uncovered genes involved in heterosis. Planta. 2018;247(6):1307–21.
Ahammed GJ, Gantait S, Mitra M, Yang Y, Li X. Role of ethylene crosstalk in seed germination and early seedling development: A review. Plant Physiol Bioch. 2020;151:124–31.
Meyer RC, Witucka-Wall H, Becher M, et al. Heterosis manifestation during early Arabidopsis seedling development is characterized by intermediate gene expression and enhanced metabolic activity in the hybrids. Plant J. 2012;71(4):669–83.
Ko DK, Rohozinski D, Song Q, et al. Temporal shift of circadian-mediated gene expression and carbon fixation contributes to biomass heterosis in maize hybrids. Plos Genet 2016; 12(e10061977).
Zhu A, Greaves IK, Liu P, Wu L, Dennis ES, Peacock WJ. Early changes of gene activity in developing seedlings of Arabidopsis hybrids relative to parents may contribute to hybrid vigour. Plant J. 2016;88(4):597–607.
Zhu A, Wang A, Zhang Y, Dennis ES, Peacock WJ, Greaves IK. Early establishment of photosynthesis and auxin biosynthesis plays a key role in early biomass heterosis in Brassica napus (canola) hybrids. Plant Cell Physiol. 2020;61(6):1134–43.
Liu PC, Peacock WJ, Wang L, Furbank R, Larkum A, Dennis ES. Leaf growth in early development is key to biomass heterosis in Arabidopsis. J Exp Bot. 2020;71(8):2439–50.
Groszmann M, Gonzalez-Bayon R, Lyons RL, et al. Hormone-regulated defense and stress response networks contribute to heterosis in Arabidopsis F1 hybrids. P Natl Acad Sci Usa. 2015;112(46):E6397–406.
Chalhoub B, Denoeud F, Liu S, et al. Plant genetics. Early allopolyploid evolution in the post-Neolithic Brassica napus oilseed genome. Science 2014; 345(6199): 950–3.
Shen Y, Sun S, Hua S, et al. Analysis of transcriptional and epigenetic changes in hybrid vigor of allopolyploid Brassica napus uncovers key roles for small RNAs. Plant J. 2017;91(5):874–93.
Shalby N, Mohamed I, Xiong J, et al. Overdominance at the gene expression level plays a critical role in the hybrid root growth of Brassica napus. Int J Mol Sci 2021; 22(17).
Rapp RA, Udall JA, Wendel JF. Genomic expression dominance in allopolyploids. Bmc Biol. 2009;7:18.
Yoo MJ, Szadkowski E, Wendel JF. Homoeolog expression bias and expression level dominance in allopolyploid cotton. Heredity (Edinb). 2013;110(2):171–80.
Hochholdinger F, Baldauf JA. Heterosis in plants. Curr Biol. 2018;28(18):R1089–92.
Chen Y, Zhou Q, Tian R, et al. Proteomic analysis reveals that auxin homeostasis influences the eighth internode length heterosis in maize (Zea mays). Sci Rep. 2018;8(1):7159.
Shahzad K, Zhang X, Guo L, et al. Comparative transcriptome analysis of inbred lines and contrasting hybrids reveals overdominance mediate early biomass vigor in hybrid cotton. BMC Genomics. 2020;21(1):140.
Huang Z, Bao K, Jing Z, et al. Small Auxin Up RNAs influence the distribution of indole-3-acetic acid and play a potential role in increasing seed size in Euryale ferox Salisb. Bmc Plant Biol. 2020;20(1):311.
Wu W, Du K, Kang X, Wei H. The diverse roles of cytokinins in regulating leaf development. Hortic Res. 2021;8(1):118.
Zhang Y, Wang B, Qi S, et al. Ploidy and hybridity effects on leaf size, cell size and related genes expression in triploids, diploids and their parents in Populus. Planta. 2019;249(3):635–46.
Groszmann M, Gonzalez-Bayon R, Greaves IK, et al. Intraspecific Arabidopsis hybrids show different patterns of heterosis despite the close relatedness of the parental genomes. Plant Physiol. 2014;166(1):265–80.
Wang L, Wu LM, Greaves IK, Dennis ES, Peacock WJ. In Arabidopsis hybrids and hybrid mimics, up-regulation of cell wall biogenesis is associated with the increased plant size. Plant Direct. 2019;3(11): e174.
Matuszyńska A, Saadat NP, Ebenhöh O. Balancing energy supply during photosynthesis - a theoretical perspective. Physiol Plant. 2019;166(1):392–402.
Fujimoto R, Taylor JM, Shirasawa S, Peacock WJ, Dennis ES. Heterosis of Arabidopsis hybrids between C24 and Col is associated with increased photosynthesis capacity. Proc Natl Acad Sci U S A. 2012;109(18):7109–14.
Kong X, Chen L, Wei T, et al. Transcriptome analysis of biological pathways associated with heterosis in Chinese cabbage. Genomics. 2020;112(6):4732–41.
Bu SH, Zhao X, Yi C, Wen J, Tu J, Zhang YM. Interacted QTL mapping in partial NCII design provides evidences for breeding by design. PLoS ONE. 2015;10(3): e121034.
Zhao X, Li B, Zhang K, et al. Breeding signature of combining ability improvement revealed by a genomic variation map from recurrent selection population in Brassica napus. Sci Rep. 2016;6:29553.
Hu Y, Xiong J, Shalby N, et al. Comparison of dynamic 3D chromatin architecture uncovers heterosis for leaf size in Brassica napus. J Adv Res 2022.
Tocquin P, Corbesier L, Havelange A, et al. A novel high efficiency, low maintenance, hydroponic system for synchronous growth and flowering of Arabidopsis thaliana. Bmc Plant Biol. 2003;3:2.
Chen S, Zhou Y, Chen Y, Gu J. Fastp: an ultra-fast all-in-one FASTQ preprocessor. Bioinformatics. 2018;34(17):i884–90.
Kim D, Langmead B, Salzberg SL. HISAT: a fast spliced aligner with low memory requirements. Nat Methods. 2015;12(4):357–60.
Liao Y, Smyth GK, Shi W. FeatureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics. 2014;30(7):923–30.
Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(12):550.
Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. 2008;9:559.
Wu J, Zhao Q, Yang Q, et al. Comparative transcriptomic analysis uncovers the complex genetic network for resistance to Sclerotinia sclerotiorum in Brassica napus. Sci Rep. 2016;6:19007.
Kanehisa M. KEGG: Kyoto Encyclopedia of Genes and Genomes. Nucleic Acids Res. 28(1):27–30. https://doi.org/10.1093/nar/28.1.27.
Chen C, Chen H, Zhang Y, et al. Tbtools: an integrative toolkit developed for interactive analyses of big biological data. Mol Plant. 2020;13(8):1194–202.
Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2(T)(-Delta Delta C) method. Methods. 2001;25(4):402–8.
This work was funded by the National Key Research and Development Program of China (2021YFF1000100) and the Hubei Provincial Science and Technology Major Project (2021ABA011).
Ethics approval and consent to participate
This study complies with relevant institutional, national, and international guidelines and legislation in ethics section.
Consent for publication
The authors declare that they have no conflict of interest.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Figure S1. Photos showing the seedling phenotypes of the parents and hybrids in the T1 and T2 periods.T1: 21DAS; T2: 24 DAS; Male parent: A, C, H; Female parent: M, W; Hybrid F1:AM, CM, HM, HW. Figure S2. Pearson correlation heat map between all RNA sequencing samples. The red color indicates a higher correlation between samples. -1, -2, and -3 represent three biological replicates at each time point. T1: 21 DAS; T2: 24DAS; -l: leaf. Figure S3. Number of the expression level of dominant(ELD) genes in hybrid canola seedlings in each of the 12 DEG types in the T1 and T2 periods. Genes with an expression degree in the F1 hybrids analogous to that of the female parent are defined as ELD-F; genes with an expression degree in the F1 hybrids analogous to that of the male parent are defined as ELD-M. T1: 21 DAS;T2: 24 DAS. Figure S4. Venn diagram analysis of the significant molecular function (MF) and cellular composition (CC) of the parental-ELD gene at 21 DAS. A and B display Venn diagrams of MF and CC, respectively. MF: molecular function; CC: cell component Figure S5. Venn diagram analysis of the significant molecular function (MF) andcellular composition (CC) of the parental-ELD gene at 24 DAS. A and B show Venn diagrams of MF and CC, respectively. MF: molecular function; CC: cell component. Figure S6. Remarkable biological processes of the parental-ELD gene at T1 and T2 in different hybrids.A and B exhibit Venn diagrams of the parental-ELD with significantly enriched BP in theT1 and T2 periods, respectively. C and D show heat maps of the BP terms where the parental-ELD gene is significantly enriched in different hybrids at 21 DAS and 24 DAS, respectively. T1: 21 DAS; T2: 24 DAS. Figure S7. Analysis of differentially expressed genes between strong and weak hybrids at 24 DAS.A Venn diagram of the number of unique and shared DEGs between the two groups of strong and weak hybrids at 24 DAS. B GO terms with significant enrichment of DEGs shared between the two strong and weak hybrids at 24 DAS. C The significantly enriched KEGG terms of DEGs shared between the two strong and weak hybrids at 24 DAS. Figure S8. Heat map of genes that promoted cell size in F1 hybrids relative to their parents at 24 DAS. T2: 24 DAS. Figure S9. Phenotypic changes of the leaf area of the third true leaf of all hybrids relative to their parents at 24 DAS. The third true leaf area A, leaf length B, and leaf width C of all F1 hybrids and their MPV at 24 DAS. D,The histogram shows the percentage of the contribution of leaf length and leaf width to the third true leaf area in four hybrids at 24 DAS. The data are expressed as mean ± SD, derived from the results of three biological replicates. AM, CM, HM, and HW are F1 hybrids; MPV: mid-parent value; *P<0.05; **P<0.01; ***P<0.001. Figure S10. Hierarchical clustering tree of 17 modules obtained by WGCNA. Figure S11. Heatmap of the correlation analysis of modular traits. Each row represents a module and each column represents a different sample. The correlation coefficients and P-values are shown in the figure. Red and green represent positive and negative correlations, respectively. Figure S12. Heatmap showing the eigengene expression profile of the yellow module in different samples. Figure S13. Significantly enriched KEGG terms in yellow modules. Figure S14. Circle chart of GO terms showing significant enrichment of hub genes in the yellow module. The first circle is the enriched classification, and the outer circle is the coordinate ruler of the number of genes. The second circle is the number of the category in the background gene. The third circle is the number enriched in the category. The fourth circle is the value of the rich factor for each category. Figure S15. qRT-PCR validation of the expression degrees of 10 genes in the photosynthetic and plant hormone signal transduction pathways. The blue histogram shows the relative expression levels (mean ± SD) of genes obtained through quantitative verification. The red line graph depicts the gene expression level value obtained by the transcriptome. Male: A, C, H; Female: M, W; F1 hybrid: AM, CM, HM, HW. The first four genes were expressed at 21 DAS, whereas the other genes were expressed at 24 DAS. Figure S16. Determination of the leaf sampling period. Figures A and B represent the line graphs of the measured daily growth rate of the leaf length and breadth of the second and third true leaves, respectively, at 20–30 DAS.
Table S1. Phenotypic data of T1 are expressed as mean ± SD, T1 refers to 21 DAS. The data represents the results of three biological repetitions. Table S2. Phenotypic data of T2 are expressed as mean ± SD, T2 refers to 24 DAS. The data represents the results of three biological repetitions. Table S3. Quality evaluation of RNA-seq reads and mapped to reference genome. Table S4. Information on upregulated and downregulated DEGs of 21 DAS hybrid F1 (AM, CM, HM, HW) and their parents (A, C, H, M, W). Table S5. Information on upregulated and downregulated DEGs of 24 DAS hybrid F1 (AM, CM, HM, HW) and their parents (A, C, H, M, W). Table S6. Specific gene information for the classification of expression patterns of 12 types of ELD genes at 21 DAS. Table S7. Specific gene information for the classification of expression patterns of 12 types of ELD genes at 24 DAS. Table S8. Results of significant GO enrichment of F1 relative to the parent of the parental-ELD gene at 21 DAS. Table S9. Results of significant GO enrichment of F1 relative to the parent of the parental-ELD gene at 24 DAS. Table S10. Information on DEGs of strong hybrids and weak hybrids at 21 DAS. Table S11. The results of the significantly enriched GO terms of the DEGs shared by the two pairs of strong and weak hybrids at 21 DAS. Table S12. The results of the significantly enriched KEGG pathways of the DEGs shared by the two pairs of strong and weak hybrids at 21 DAS. Table S13. Information on DEGs of strong hybrids and weak hybrids at 24 DAS. Table S14. The results of the significantly enriched GO terms of the DEGs shared by the two pairs of strong and weak hybrids at 24 DAS. Table S15. The results of the significantly enriched KEGG pathways of the DEGs shared by the two pairs of strong and weak hybrids at 24 DAS. Table S16. Total genes in the yellow module. Table S17. KEGG results of significant enrichment of genes in the yellow module. Table S18. Hub genes in the yellow module. Table S19. GO terms that were significantly enriched with respect to hub genes in the yellow module. Table S20. Specific information about the sequence of the quantitative primers of randomly selected genes.
About this article
Cite this article
Xiong, J., Hu, K., Shalby, N. et al. Comparative transcriptomic analysis reveals the molecular mechanism underlying seedling biomass heterosis in Brassica napus. BMC Plant Biol 22, 283 (2022). https://doi.org/10.1186/s12870-022-03671-0
- Brassica napus
- Plant hormones
- Cell size