Transcriptome Analysis of Interspecific Hybrid between Brassica napus and B. rapa Reveals Heterosis for Oil Rape Improvement

Abstract

The hybrid between Brassica napus and B. rapa displays obvious heterosis in both growth performance and stress tolerances. A comparative transcriptome analysis for B. napus ([superscript]An[/superscript] [superscript]An[/superscript] CC genome), B. rapa ([superscript]Ar[/superscript] [superscript]Ar[/superscript] genome), and its hybrid F1 ([superscript]An[/superscript] [superscript]Ar[/superscript] C genome) was carried out to reveal the possible molecular mechanisms of heterosis at the gene expression level. A total of 40,320 nonredundant unigenes were identified using B. rapa (AA genome) and B. oleracea (CC genome) as reference genomes. A total of 6,816 differentially expressed genes (DEGs) were mapped in the A and C genomes with 4,946 DEGs displayed nonadditively by comparing the gene expression patterns among the three samples. The coexistence of nonadditive DEGs including high-parent dominance, low-parent dominance, overdominance, and underdominance was observed in the gene action modes of F1 hybrid, which were potentially related to the heterosis. The coexistence of multiple gene actions in the hybrid was observed and provided a list of candidate genes and pathways for heterosis. The expression bias of transposable element-associated genes was also observed in the hybrid compared to their parents. The present study could be helpful for the better understanding of the determination and regulation of mechanisms of heterosis to aid Brassica improvement.

Jinfang Zhang 1 and Guangrong Li 2 and Haojie Li 1 and Xiaobin Pu 1 and Jun Jiang 1 and Liang Chai 1 and Benchuan Zheng 1 and Cheng Cui 1 and Zujun Yang 2 and Yongqing Zhu 3 and Liangcai Jiang 1

Academic Editor:Graziano Pesole

1, Crop Science Institute, Sichuan Academy of Agricultural Sciences, Chengdu, Sichuan 6110066, China
2, School of Life Science and Technology, University of Electronic Science and Technology of China, Chengdu, Sichuan 610054, China
3, Institute of Agro-Products Processing Science and Technology, Sichuan Academy of Agricultural Sciences, Chengdu, Sichuan 6110066, China

Received 3 May 2015; Revised 3 August 2015; Accepted 5 August 2015; 10 September 2015

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

Brassica napus (rapeseed, AACC, [figure omitted; refer to PDF] ) is one of the main oil crops used for human consumption and is widely grown in China, Canada, Europe, and Australia and increasingly grown in South America [1]. B. napus is likely derived from the natural hybridization and genome doubling from its two diploid parents, B. rapa (genome AA, [figure omitted; refer to PDF] ) and B. oleracea (genome CC, [figure omitted; refer to PDF] ) along the Mediterranean coastline in Southern Europe approximately 10,000 years ago [2]. However, there is evidence that the genetic diversity in rapeseed has been continuously reduced by extensive breeding efforts [3]. Approaches to increase the genetic variation of rapeseed have been made by introducing the genetic components from related species as B. napus can easily cross with one of its ancestral parents, B. rapa, to produce viable interspecies hybrids [4]. Due to the high crossability between B. napus and B. rapa , and low aneuploidy of their interspecific hybrids, different B. rapa accessions have been widely hybridized to rapeseed breeding program and novel agronomic traits from B. rapa have been successfully transmitted into commercial B. napus varieties [5]. Strong heteroses affecting both biomass and seed yield have been observed in hybrids derived from these interspecific crosses [6]. However, the genetic and molecular mechanism of heterosis in the interspecific hybrids has not been investigated.

Heterosis is a prevalent phenomenon in evolution and breeding process of plants [7]. The hybrid offspring obtains the advantages in many agricultural and developmental traits including biomass yield, plant height, vigor, and stress tolerances from their parents, thereby improving the adaptation of the crop and subsequently increasing the global area sown [8]. Although heterosis has been widely utilized, the genetic and molecular basis of heterotic improvements remains unclear. Various genetic models have been proposed for explaining heterosis and include dominance, overdominance, and pseudo overdominance [9]. Heterotic interactions in the tetraploid B. napus can be derived from two types of allelic interactions, namely, those between the genomes of the two diploid progenitors as well as those between the two parental genomes making up a hybrid. Gene expression is potentially altered not only due to a dosage effect but also with gene functional divergence [8]. Epigenetic modification in hybrids genomes can also account for heterosis [10]. Emerging new technologies make it possible to investigate the DNA sequence, RNA transcripts, proteins, and metabolism products genome-wide, allowing for new insights into benefits of heterotic recombinations [7, 9]. Next generation sequencing technology utilizing RNA-seq has emerged as a powerful tool for investigating gene expression data at the whole genome level. Both the gene transcript and expression levels can be detected in RNA-seq which could provide new insights into the molecular mechanisms of heterosis [11].

Chinese cabbage, Brassica rapa (syn. B. campestris ), contains two subspecies named pekinensis and chinensis. These two distinct groups are used as leafy vegetables widely across eastern Asia and have much genetic variation. In particular, B. rapa ssp. chinensis has wide morphological variation and is found in China, Korea, and Japan [12]. Recently, a number of genes responsible for developmental and stress tolerances were isolated [13-15], showing that B. rapa ssp. chinensis could be a novel resource for B. napus improvement.

In this study, the transcriptomes of B. rapa ssp. chinensis Makino, B. napus, and their interspecific hybrid were sequenced to investigate the molecular basis of heterosis. Gene expression profiles were used to identify the number of differentially expressed genes, and gene expression levels were analyzed and compared among two parents and their hybrid. This data was used to give a better understanding of additive and nonadditive effects of genes and how these contribute to heterosis, thus providing new insights into the genetic and epigenetic mechanism of subgenomic heterosis in B. napus .

2. Materials and Methods

2.1. Sample Preparation and RNA Isolation

Three lines were grown in a green house at 22°C under 12 hours of light and 12 hours of darkness. Leaves of the seedling were harvested after 4 weeks by cutting off the sixth leaf of each seedling. The leaves of four plants were combined for RNA extraction. Samples were immediately frozen in liquid nitrogen and stored at -80°C until RNA was extracted. Total RNA of each sample was isolated using RNAprep Pure Plant Kit (Tiangen, Shanghai, China). RNA quality was characterized initially by NanoDrop ND1000 spectrophotometer (NanoDrop Technologies, Wilmington, DE, USA) and then further assessed by RIN (RNA Integrity Number) value (>9.5) using Agilent 2100 Bioanalyzer (Santa Clara, CA, USA). Equal quantities of high-quality RNA from each sample were pooled for cDNA synthesis. Cytological observation of the F1 hybrid was performed as described by Zhan et al. [16].

2.2. cDNA Library Construction for Illumina Sequencing

The cDNA library was constructed following the manufacturer's instructions of mRNA-Seq Sample Preparation Kit (cat. number RS-930-1001, Illumina, Inc., San Diego, CA). Briefly, the poly(A) mRNA was isolated from total RNA samples with Magnetic Oligo (dT) Beads. The mRNA was then fragmented into small pieces using RNA fragmentation kit (Ambion). Using these short fragments as the templates, the first cDNA strand was synthesized using random hexamer primers and reverse transcriptase (Invitrogen), and the second-strand cDNA was synthesized using DNA polymerase I and RNase H. The cDNA fragments were purified using the QiaQuick PCR extraction kit (Qiagen) and resolved with EB buffer for end reparation and poly(A) addition. The short fragments were then connected with sequencing adapters, and the products were subsequently purified and amplified via PCR. Libraries were prepared from a 400-500 bp size-selected fraction following adapter ligation and agarose gel separation. The quality control analysis on the sample library was performed to quantify the DNA concentration and validate the library. After validation with an Eppendorf Mastercycler Real-Time PCR System, the cDNA libraries were sequenced on the Illumina Hiseq 2000 platform with read length of 2 × 100 bp. The sequencing-derived raw image data were transformed by base calling into sequence data using Illumina Pipeline Software v1.6. The raw sequencing reads have been submitted to NCBI Short Read Archive under the accession number of SRR2131203, SRR2134440, and SRR2134444.

2.3. Sequence Data Analysis and Expression Analysis

The raw reads were cleaned by removing adapter sequences, low-quality sequences (reads with ambiguous bases " [figure omitted; refer to PDF] "), and reads with more than 10% [figure omitted; refer to PDF] bases. To analyze transcript abundance levels, the uniquely mapped reads for a specific gene were counted by mapping reads to de novo assembled distinct sequences using SOAP2 software [17], and the RPKM (Reads Per Kb per Million reads) values were computed as proposed by Mortazavi et al. [18]. Gene transcript abundance differences were obtained from RPKM values using a method modified from Audic's proposal [19]. Gene expression levels were calculated with the RPKM method which is able to eliminate the influence of different gene lengths and sequencing discrepancy within the calculation of gene expression.

2.4. Statistics

Statistical analysis of the phenotypic data was performed with a Welch two-sample [figure omitted; refer to PDF] -test in R (R Development Core Team 2005) [20]. Midparent heterosis (MPH) values were calculated using the formula [figure omitted; refer to PDF] (whereas [figure omitted; refer to PDF] is the value of the hybrid and P1 and P2 are the values of the one and the other parental line). To test for significance of MPH values the contrast [figure omitted; refer to PDF] was used [21].

2.5. Quantitative RT-PCR Validation of Differentiation Gene Expression

RNA-seq data were further validated using quantitative real-time PCR analysis (qRT-PCR) for a selected number of genes using gene-specific primer sets. The plant leaves and growth condition for total RNA extraction used for qRT-PCR were the same as these for RNA-seq experiments. Three independent RNA extractions per sample were performed for biological replicates. The reverse transcribed into cDNA using PrimeScritH RT reagent kit with gDNA Eraser (Takara, Dalian, China). Primer pairs of 10 unigenes were designed using Invitrogen's (Carlsbad, CA, USA) OligoPerfect Designer software [22]. Specificity of the primer sets and their product length was verified by agarose gel electrophoresis. The qRT-PCR reaction mixture consisted of the SYBR Premix EX-Taq II Kit (TakaRa, Japan) on iCycler iQ (BIO-RAD) for three repeats of each sample. The thermal cycling conditions were as follows: 95 uC 2 min and 40 cycles at 95°C for 10 s for denaturation and 65°C 20 s and 72°C 30 s for annealing and extension. The expression of TIPS-41 with primer sequences 5[variant prime]-TGAAGAGCAGATTGATTTGGCT-3[variant prime] and 5[variant prime]-ACACTCCATTGTC AGCCAGTT-3[variant prime] was used as an internal control for normalization to compare the gene expression level between the accessions. The relative levels of gene expression were calculated using the [figure omitted; refer to PDF] method.

3. Results

3.1. Sample Preparation, Sequencing, and Data Filtering

The hybridization of B. napus cv. CWH-2 (P1) and B. rapa ssp. chinensis cv. Qianjin (P2) was made in 2010-2011 oilseed growing season. Root-tip chromosome assessments of F1 showed that all F1 plants had a chromosome number of 29 (Figure 1), suggesting that they were hybrids. The seeds' sizes of F1 hybrids were intermediate to those of the parents, but the seed color resembled that of B. napus (Figure 2(a)). Root growth of the F1 over the first 24 h was similar to P1. Leaf morphology of the F1 was more similar to B. rapa than to B. napus (Figure 2). In 4-week-old seedlings, the F1 plant displayed stronger roots than either of the parents (Figure 2(b)) while adult F1 plants displayed the similar heading time to B. napus .

Figure 1: Root tip squash showing chromosomes from the F1 hybrid (chromosome number of 29).

[figure omitted; refer to PDF]

Figure 2: Seed germination (a) and 2-week-old seedlings (b) of P1, P2, and F1 germplasm.

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

To identify transcripts of P1, P2, and F1, three cDNA libraries were constructed from their leaves and sequenced by Illumina paired-end sequencing. Sequence data from each of the libraries consisted of 6.0 Gb from B. rapa (AA genome), 5.2 Gb from B. napus (AACC genome), and 6.7 Gb from the hybrid (AAC genome). High-quality cleaned raw reads from all of those data were aligned to reference sequences. In total, 67,951,208 reads and 46,207,184 reads were mapped in the B. rapa and B. oleracea reference genomes, respectively. In the B. rapa genome, 38,522 unigenes were assembled with a minimum scaffold size of 200 bp and a total length of 51,995,312 bp, with an average length of 1,349 bp (Figure 3(a)). The B. oleracea genome had a total of 42,009 unigenes with a minimum scaffold size of 300 bp, a total length of 51,691,329 bp, and an average length of 1,230 bp (Figure 3(b)).

Figure 3: The size distribution of transcriptomic unigenes assembled by the B. rapa (a) and B. oleracea (b).

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

3.2. Heterosis of Differential Gene Expression in B. rapa , B. napus , and Their Hybrid

The numbers of expressed gene in each genome varied slightly. There were 29,579 (71.8%), 27,939 (67.9%), and 29,447 (71.7%) from B. rapa, B. napus, and the hybrid genome, respectively. The unique mapping rates were similar among three Brassica samples, suggesting that their genomes were quite alike. Some genes had expression that was germplasm specific. The B. rapa genome had approximately about 6% (1,639) of genes that were unique, whilst B. napus and the hybrid had only 2% (420) and 3% (951) unique sequences, respectively. There were 25,985 genes expressed in all three lines. A total of 489 genes were expressed in both B. rapa and B. napus but silenced in the hybrid. Gene expression profiles of the three lines revealed that the number of expressed genes did not show obvious difference in the three sets of material. The mode of gene action for the differentially expressed genes (DEGs) was analyzed. About 27.00% genes (1840 of 6816) exhibited an expression pattern that was not distinguishable from additivity, while the other 72.56% (4946 of 6816) of genes showed nonadditive expression patterns (Table 1). The nonadditive number of DEGs from the cross was further classified into five distinct classes: high-parent dominance (HPD), low-parent dominance (LPD), overdominance (ODO), underdominance (UDO), and partial-dominance. There were 7% of transcripts (545) that were absent in their parents, indicating that there was transcriptional activation of new genes in hybrid F1. There was also transcriptional silencing with 3.2% (218) of parentally expressed gene not being expressed in the hybrid.

Table 1: Summary of dominance patterns of differentially expressed genes.

DGs	Number	Percentage
Additivity	1840	27
Nonadditivity	4946	72.56
Others	30	0.44
High-parent dominance (HPD)	353	7.14
Low-parent dominance (LPD)	191	3.86
Overdominance (ODO)	593	11.99
Underdominance (UDO)	193	3.9
Positive partial-dominance (PPD)	1635	33.06
Negative partial-dominance (NPD)	1981	40.5
Total	6816

Additivity: F1 [approximate] 1/2(P1 + P2); nonadditivity: F1 > 1/2(P1 + P2) or F1 < 1/2(P1 + P2). High-parent dominance (HPD): F1 [approximate] P1 > P2 or F1 [approximate] P2 < P1; low-parent dominance (LPD): F1 [approximate] P1 < P2 or F1 [approximate] P2 < P1; overdominance (ODO): F1 > P1 and F1 > P2; underdominance (UDO): F1 < P1 and F1 < P2.

There were 2894 unigenes that mapped to the A genome, 2107 that mapped to the C genome (Figure 4), and 593 that were not assigned to either genome. Chromosomal distribution of DEG numbers in hybrids showed that low-parent domiance (LPD) genes were significantly accumulated on chromosomes A06 and C07 ( [figure omitted; refer to PDF] ), while high-parent dominance (HPD) genes were on chromosomes A09 and C08 ( [figure omitted; refer to PDF] ). Gene expression variation caused by interspecific hybridization was not randomly distributed along the chromosomes. Therefore, different B. napus -B. rapa chromosome introgression lines can be developed to the further study of heterosis pertaining to specific chromosomes.

Figure 4: Chromosomal distribution of nonadditive differentially expressed genes with high-parent dominance (HPD), low-parent dominance (LPD), overdominance (ODO), and underdominance (UDO) on the A genome (a) and C genome (b).

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

3.3. Contribution of A and C Subgenomes to Gene Differential Expression in the Hybrid

There were 1386 significant differentially expressed genes between P1 and F1 samples, 858 genes had increased expression in the F1, and 528 genes had decreased expression. There were 5077 significant differentially expressed genes between P2 and F1 samples, with 3241 genes having increased expression in the F1 and 1836 genes being decreased in expression (Figure 5). Changes to expression from P2 were approximately three times the number of genes compared to from P1. There were more genes with increased expression in the F1 (4098) than with decreased expression (2464) when compared to expression levels in the respective parents. We performed gene ontology (GO) enrichment analysis to test the functional categories for DEGs among P1, P2, and F1 (Figure 6) by AgriGO online tool [23] and R statistical software [20]. There were 10, 12, and 18 functional categories of these transcripts that belonged to the cell component (CC), molecular function (MF), and biological process (BP), respectively. The GO analysis showed that the composition of functional pathways associated with genes showing differential expression was similar when comparing changes from either P1 or P2 with F1.

Figure 5: Total numbers of differentially expressed genes between P1, P2, and F1 by Venn diagram analyses (a) and statistics of up- or downregulated genes (b).

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

Figure 6: The heat map from GO enrichment analysis of DEGs between the P1, P2, and F1 combinations. Color scale represents [figure omitted; refer to PDF] values of enrichment test.

[figure omitted; refer to PDF]

The biological process category was further subdivided into functional classes in those DEGs with high-parent dominance (HPD, [figure omitted; refer to PDF] ) based on gene ontology. Table 2 highlights the 18 gene ontology classes that were differentially expressed in the [figure omitted; refer to PDF] . Combined with the results in Figure 6, we found that the GOs such as regulation of transcription, metabolic process, defense response, multicellular organismal development, and transporters were enriched in A genomes transcripts. The genes responsible for stress were striking enriched in the differential expressed genes from C genome. It possibly implied that the hybrid expressed protein-coding genes from A genomes of both B. rapa and B. napus potentially participated in metabolism and development, while those from C genome of B. napus were largely involved in stress resistance.

Table 2: Significant gene ontology terms of high-parent dominance (HPD) in the biological process category from A and C genomes.

Gene ID	Gene O terms
Bra028791	Regulation of transcription (GO:0006355)
Bra005111	Metabolic process (GO:0008152)
Bra032185	Defense response (GO:0002679)
Bra010724	Multicellular organismal development (GO:0007275)
Bra007111	Metabolic process (GO:0010408)
Bra031515	Single-organism transport (GO:0044765)
Bra031678	Transport (GO:0006810)
Bra027171	Transmembrane transport (GO:0006855)
Bra024033	Peroxidase activity (GO:0004601)
Bra002550	Chromosome segregation (GO:0007059)
Bol015697	Response to molecule of bacterial origin (GO:0002237)
Bol032712	Defense response to insect (GO:0002213)
Bol014351	Defense response to bacterium (GO:0042742)
Bol045822	Defense response (GO:0002679)
Bol034224	Response to wounding (GO:0009611)
Bol038369	Single-organism process (GO:0044699)
Bol045822	Regulation of transcription (GO:0006355)
Bol015697	Response to molecule of bacterial origin (GO:0002237)

3.4. Transposon Active or Inactive in Hybrid

Compared to the B. rapa genome, there were 25 transposable elements (TE) like genes expressed in P1 and P2, while the 11 (44%) genes displayed differential expression in [figure omitted; refer to PDF] (Table 3). The En/Spm type of transposon, and the retrotransposon-like transcript were significantly enriched. It is likely that the interspecific hybrids significantly modify the transcription of transposon-like genes, which may directly or indirectly influence gene expression of other genes in the heterotic material.

Table 3: The differential expressed transposon-like genes in hybrid [figure omitted; refer to PDF] .

Gene ID	Regulation	Chromosome	Similar to transposon
B.rapa.newgene.1068	Inactive	A09	Retrotransposon-like protein (Arabidopsis thaliana )
B.rapa.newgene.871	Inactive	A07	Putative transposon protein (Arabidopsis thaliana )
B.rapa.newgene.1142	Inactive	A09	Retrotransposon-like protein (Arabidopsis thaliana )
B.rapa.newgene.957	Inactive	A08	Retrotransposon (Arabidopsis thaliana )
Bra040302	Up	A08	Putative transposon protein (Arabidopsis thaliana )
B.rapa.newgene.307	Up	A03	En/Spm transposon protein (Arabidopsis thaliana )
B.rapa.newgene.726	Up	A06	En/Spm-like transposon protein (Arabidopsis thaliana )
Bra034711	Active	A05	hAT transposon superfamily protein (Arabidopsis thaliana )
Bra027867	Active	A09	Pol polyprotein transposon element Bs1 (Zea mays )
B.rapa.newgene.8	Active	A01	Tam1 transposon protein (Arabidopsis thaliana )
Bra035011	Active	A07	Putative transposon protein (Arabidopsis thaliana )
Bol019342	Inactive	C09	Retrotransposon Tto1 DNA (Nicotiana tabacum )
Bol041429	Up	C07	Putative Tam3-like transposon protein (Zea mays )
Bol035479	Down	C03	En/Spm-related transposon protein (Brassica oleracea )
Bol005562	Down	C03	En/Spm-related transposon protein (Brassica oleracea )
Bol019765	Down	C09	En/Spm-related transposon protein (Brassica oleracea )
Bol006373	Down	C07	Transposon protein-like (Arabidopsis thaliana )
Bol035958	Up	C02	Transposon-like ORF (Brassica oleracea )
Bol030703	Active	C03	Retrotransposon Tto1 DNA (Nicotiana tabacum )
Bol030374	Active	C09	En/Spm-related transposon protein (Brassica oleracea )
Bol004880	Active	C04	Retrotransposon Tto1 DNA (Nicotiana tabacum )
Bol043208	Active	C07	Similar to retrotransposon (Arabidopsis thaliana )

3.5. Validation of Transcriptome Data by qRT-PCR

To assess the accuracy of RNA-seq data, ten differentially expressed unigenes including stress responsive genes, secondary metabolism biosynthesis genes, and epigenetic modifying genes were selected. Three genes belong to underdominance (UDO) and seven genes are belonging to overdominance (ODO) types, which are represented for heterosis analysis. Primer pairs of 10 unigenes for qRT-PCR were designed and listed in Table 4. We tested the similarity between differential gene expression identified by transcriptome and those identified by qRT-PCR. As shown in Figure 7, the qRT-PCR revealed that 8 of 10 genes (except Bol022348 and Bol026880) that showed the differential gene expression level agreed well with the expression patterns of DGE data. Hence, the qRT-PCR results showed general agreement with their transcript abundance changes determined by RNA-seq, which suggested the reliability of the transcriptomic profiling data among the three samples.

Table 4: The qPCR primers for selected genes.

Unigenes	Primers	Biological process	Nonadditivity
Bra007087 F	TGCAGCGCTTGATTTACCT	Oxidation-reduction process	HPD
Bra007087 R	GCAAACTCCAGAGCTATGT
Bra020376 F	AGGTCATTCTGGTGAGCCACA	Chaperones	HPD
Bra020376 R	TGGAGACTTTGGAAGGATACT
Bra031210 F	ACGAGGCTCAGTCTCGTGGT	Response to water deprivation	ODO
Bra031210 R	TCCGCTGCGGTATCCACCA
Bra031666 F	TGCCAAGGACAACAACTTGGACT	Lignin biosynthetic process	HPD
Bra031666 R	AGTTGGTTGTAGGACTGGTCCA
Bol016209 F	TAACCTACCAGAAGCACGGT	Response to abscisic acid	ODO
Bol016209 R	ACAACTTCAACGGTGCACGACT
Bol022348 F	AGGTACCTTACGAGTCTCGT	Response to karrikin	HPD
Bol022348 R	CACGAACCTGATAGAAGCTCGT
Bol026880 F	TGTAAGGCTACGAAGGGACAT	Defense response to fungus	HPD
Bol026880 R	TCCAACATATCCATATGTCCGT
Bol026959 F	TAGCTCTTCCTCTTCAAGCGAT	Defense response to fungus	ODO
Bol026959 R	TCCTTCTTCCTCTTCTCACCA
Bol041012 F	TGACACATTGTGGATGGAACT	Cellular hyperosmotic response	ODO
Bol041012 R	GGCATACCAATCATTGGAACT
Bol044348 F	TACTACTGGACCTTTGGTGCT	Calcium iontransport	HPD
Bol044348 R	GATCTTCATCTGAAGGTCACT

Figure 7: Genes with expression levels validated by qRT-PCR.

[figure omitted; refer to PDF]

[figure omitted; refer to PDF]

4. Discussion

Interspecific hybridization is an important approach to obtaining and utilizing novel agronomic traits from related species for crop improvement. Previous researches revealed strong heterosis between B. napus and B. rapa on biomass and seed yield [5-7]. Recently, Li et al. [24] reported that B. napus -like individuals in the F3 and F4 generations, from interspecific hybridization between B. napus and B. oleracea , showed diverse genetic variation relative to current B. napus and strong heterotic potential. We selected a number of backcrossed offspring derived from the F1 by being crossed with P1 (B. napus ) and new traits such as wrinkle leaves, six sepals which different from their parents were observed in the progenies. In particular, we identified several lines which displayed high yield up to 4,500 kg per hectare in BC1F2 population. It suggested that the novel lines from interspecific hybridization between B. napus and B. rapa with high yield good quality can be potentially useful for oilseeds production in the near future.

Intraspecific hybrids of Arabidopsis thaliana [25, 26], rice [27-29], and maize [30, 31] revealed that there are both additive and nonadditive effects of gene expression in hybrids, and heterosis is strongly influenced by genetic differences between two parents. A transcriptome approach has been utilized to investigate the heterosis of intraspecific hybrids of A. thaliana , rice, maize [29, 30]. Our research of the transcriptomes from B. napus and B. rapa demonstrated the genetic effect of interspecies heterosis at the gene expression level. More genes specific to certain pathways in further experiments are needed to fully validate the heterotic association between the gene expression patterns and the target agronomic traits. In the present study, we found that both additive and dominance effects contributed to interspecies heterosis between B. napus and B. rapa . Analysis of the number of DEGs in hybrid and parental germplasm revealed that the nonadditive effect contributed more than the additive effect. Previous research reveals interspecies heterosis between B. napus and B. rapa through genetic analysis but the main focus was on phenotypic variation [6, 32]. Our study provided evidence of differential gene expression and shed new insights on the molecular mechanism for this heterosis in the interspecific hybrids.

Brassica species constitute an elite system for investigating the changes in genomic structure and functional divergence of duplicate genes in the process of hybridization and polyploidization among the diploid species and their amphidiploids [32, 33]. Subgenomes in Brassica species are designated with superscripts with [figure omitted; refer to PDF] and [figure omitted; refer to PDF] representing A genomes from B. napus and B. rapa , respectively [33]. We have shown intersubgenomic heterosis between [figure omitted; refer to PDF] and [figure omitted; refer to PDF] genomes, supporting work by others that show that such heterosis can be useful in creating genetic diversity and be used in breeding. Studies indicated that the midparent heterosis contributes intersubgenomic heterosis [8, 34, 35]. Therefore, the utilization of intersubgenomic heterosis can be creating genetic diversity and breeding. In the present study, we found that the nonadditive DEGs from the A genome in the hybrids mainly participated in metabolism and development, while those from the C genome were largely involved in stress resistance. It is likely that the neofunctionalization of subgenome contributes the heterosis or the genomic adaptation in hybrid. Recently, Li et al. [36] investigated the total homoeolog expression level of hexaploid wheat (AABBDD) and compared it to its ancestral parents T. turgidum (AABB) and Ae. tauschii (DD), finding that genes expressed from the AABB potentially participated in development and those from Ae. tauschii (DD) were more likely involved in adaptation. Differential contributions from the parental sources may also be occurring in the Brassica hybrids, and this may hint at the molecular mechanisms behind heterosis. Transposable elements (TEs) have major effects on different plant genomes and can significantly contribute differences in genome size of over 1,000-fold [37]. Recent sequencing of the B. rapa and B. oleracea genomes revealed that the amplification of TEs is one of the main factors inducing the difference in genome size [37]. The overall TEs in either Brassica genomes are expressed at very low levels, and the expression levels of different TE categories and families vary among different organs [38]. Recently, the variations in TE-associated sequences during the process of allopolyploidisation were detected between Brassica rapa (AA) and B. oleracea (CC), as well as in successive generations of self-pollinated progeny [39]. Our study revealed that the expression activities of TE categories or families were clearly modified in the hybrids between the two Brassica species, and different types of TEs showed different patterns of variation during the process of hybridization. Therefore, the association of activation and inactivation of TEs and of TE related genes expression in relation to heterotic gene expression needs to be further investigated.

5. Conclusions

Development of interspecific hybrids has been widely exploited for the heterosis breeding of Brassica crops. The interspecific hybrid between Brassica napus and B. rapa displayed obvious heterosis. In the present study, the comparative transcriptome analysis for the parents B. napus and B. rapa and its hybrid F1 was sequenced using the platform of Illumina/Solexa with de novo assemblage. The nonredundant unigenes were identified using B. rapa (A genome) and B. oleracea (C genome) as reference genomes. The nonadditively differentially expressed genes were potentially related to the heterosis. The results indicated that the differentially expressed genes in hybrids from the A genome mainly participated in metabolism and development, while those from the C genome were largely involved in stress resistance. The coexistence of multiple gene actions and transposable element-associated genes regulations in the hybrid were observed. The present study could be helpful for the better understanding of the mechanisms of heterosis, and a list of candidate genes will be used for future Brassica breeding.

Acknowledgments

The authors thank the Earmarked Fund for Modern Agro-industry Technology Research System of China, the National High Technology Research and Development Program "863" of China (2011AA10A104), and National Science and Technology Supporting Plan of China (2011BAD35B04), The Ministry of agriculture experimental observation of the upper reaches of the Yangtze River Oil Crop Science Station (09203020), and Sichuan Crop Breeding Community (2011NZ0098-5, 2011NZ0098-16, and 2011JYGC04-013) for the financial support.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

References

[1] J. A. Foley, N. Ramankutty, K. A. Brauman, E. S. Cassidy, J. S. Gerber, M. Johnston, N. D. Mueller, C. O'Connell, D. K. Ray, P. C. West, C. Balzer, E. M. Bennett, S. R. Carpenter, J. Hill, C. Monfreda, S. Polasky, J. Rockström, J. Sheehan, S. Siebert, D. Tilman, D. P. M. Zaks, "Solutions for a cultivated planet," Nature , vol. 478, no. 7369, pp. 337-342, 2011.

[2] U. Nagaharu, "Genome analysis in Brassica with special reference to the experimental formation of B. napus and peculiar mode of fertilization," Japanese Journal of Botany , vol. 7, pp. 389-452, 1935.

[3] A. Abbadi, G. Leckband, "Rapeseed breeding for oil content, quality, and sustainability," European Journal of Lipid Science and Technology , vol. 113, no. 10, pp. 1198-1206, 2011.

[4] T. Jesske, B. Olberg, A. Schierholt, H. C. Becker, "Resynthesized lines from domesticated and wild Brassica taxa and their hybrids with B. napus L.: genetic diversity and hybrid yield," Theoretical and Applied Genetics , vol. 126, no. 4, pp. 1053-1065, 2013.

[5] H. L. Liu Genetics and Breeding in Rapeseed , Chinese Agricultural Universities Press, Beijing, China, 2000.

[6] R. Liu, W. Qian, J. Meng, "Association of RFLP markers and biomass heterosis in trigenomic hybrids of oilseed rape (Brassica napus × B. campestris )," Theoretical and Applied Genetics , vol. 105, no. 6-7, pp. 1050-1057, 2002.

[7] D. Fu, M. Xiao, A. Hayward, G. Jiang, L. Zhu, Q. Zhou, J. Li, M. Zhang, "What is crop heterosis: new insights into an old topic," Journal of Applied Genetics , vol. 56, no. 1, pp. 1-13, 2014.

[8] J. A. Birchler, D. L. Auger, N. C. Riddle, "In search of the molecular basis of heterosis," Plant Cell , vol. 15, no. 10, pp. 2236-2239, 2003.

[9] Z. J. Chen, "Genomic and epigenetic insights into the molecular bases of heterosis," Nature Reviews Genetics , vol. 14, no. 7, pp. 471-482, 2013.

[10] N. R. Hofmann, "A global view of hybrid vigor: DNA methylation, small RNAs, and gene expression," Plant Cell , vol. 24, no. 3, article 841, 2012.

[11] R. Zhai, Y. Feng, H. Wang, X. Zhan, X. Shen, W. Wu, Y. Zhang, D. Chen, G. Dai, Z. Yang, L. Cao, S. Cheng, "Transcriptome analysis of rice root heterosis by RNA-Seq," BMC Genomics , vol. 14, article 19, 2013.

[12] S. Yu, F. Zhang, X. Wang, X. Zhao, D. Zhang, Y. Yu, J. Xu, "Genetic diversity and marker-trait associations in a collection of Pak-choi (Brassica rapa L. ssp. chinensis Makino) Accessions," Genes and Genomics , vol. 32, no. 5, pp. 419-428, 2014.

[13] X. L. Zhang, L. G. Zhang, "Molecular cloning and expression of the male sterility-related CtYABBY1 gene in flowering Chinese cabbage (Brassica campestris L. ssp chinensis var. parachinensis )," Genetics and Molecular Research , vol. 13, no. 2, pp. 4336-4347, 2014.

[14] N. Wang, F. Li, B. Chen, K. Xu, G. Yan, J. Qiao, J. Li, G. Gao, I. Bancroft, J. Meng, G. J. King, X. Wu, "Genome-wide investigation of genetic changes during modern breeding of Brassica napus ," Theoretical and Applied Genetics , vol. 127, no. 8, pp. 1817-1829, 2014.

[15] S. Lin, H. Dong, F. Zhang, L. Qiu, F. Wang, J. Cao, L. Huang, "BcMF8, a putative arabinogalactan protein-encoding gene, contributes to pollen wall development, aperture formation and pollen tube growth in Brassica campestris ," Annals of Botany , vol. 113, no. 5, pp. 777-788, 2014.

[16] H. Zhan, X. Zhang, G. Li, Z. Pan, J. Hu, X. Li, L. Qiao, J. Jia, H. Guo, Z. Chang, Z. Yang, "Molecular characterization of a new wheat-Thinopyrum intermedium translocation line with resistance to powdery mildew and stripe rust," International Journal of Molecular Sciences , vol. 16, no. 1, pp. 2162-2173, 2015.

[17] R. Li, C. Yu, Y. Li, T.-W. Lam, S.-M. Yiu, K. Kristiansen, J. Wang, "SOAP2: an improved ultrafast tool for short read alignment," Bioinformatics , vol. 25, no. 15, pp. 1966-1967, 2009.

[18] A. Mortazavi, B. A. Williams, K. McCue, L. Schaeffer, B. Wold, "Mapping and quantifying mammalian transcriptomes by RNA-Seq," Nature Methods , vol. 5, no. 7, pp. 621-628, 2008.

[19] S. Audic, J.-M. Claverie, "The significance of digital gene expression profiles," Genome Research , vol. 7, no. 10, pp. 986-995, 1997.

[20] RC Team R: A Language and Environment for Statistical Computing , R Foundation for Statistical Computing, Vienna, Austria, 2013.

[21] H. Ding, C. Qin, X. Luo, L. Li, Z. Chen, H. Liu, J. Gao, H. Lin, Y. Shen, M. Zhao, T. Lübberstedt, Z. Zhang, G. Pan, "Heterosis in early maize ear inflorescence development: a genome-wide transcription analysis for two maize inbred lines and their hybrid," International Journal of Molecular Sciences , vol. 15, no. 8, pp. 13892-13915, 2014.

[22] Custom Primers-OligoPerfect Designer , Invitrogen, Carlsbad, Calif, USA, 2006.

[23] Z. Du, X. Zhou, Y. Ling, Z. Zhang, Z. Su, "AgriGO: a GO analysis toolkit for the agricultural community," Nucleic Acids Research , vol. 38, no. 2, pp. W64-W70, 2010.

[24] Q. Li, Q. Zhou, J. Mei, Y. Zhang, J. Li, Z. Li, X. Ge, Z. Xiong, Y. Huang, W. Qian, "Improvement of Brassica napus via interspecific hybridization between B. napus and B. oleracea ," Molecular Breeding , vol. 34, no. 4, pp. 1955-1963, 2014.

[25] J. Wang, L. Tian, H.-S. Lee, N. E. Wei, H. Jiang, B. Watson, A. Madlung, T. C. Osborn, R. W. Doerge, L. Comai, Z. J. Chen, "Genome-wide non-additive gene regulation in Arabidopsis allotetraploids," Genetics , vol. 172, no. 1, pp. 507-517, 2006.

[26] R. Fujimoto, J. M. Taylor, S. Shirasawa, W. J. Peacock, E. S. Dennis, "Heterosis of Arabidopsis hybrids between C24 and Col is associated with increased photosynthesis capacity," Proceedings of the National Academy of Sciences of the United States of America , vol. 109, no. 18, pp. 7109-7114, 2012.

[27] X. Ge, W. Chen, S. Song, W. Wang, S. Hu, J. Yu, "Transcriptomic profiling of mature embryo from an elite super-hybrid rice LYP9 and its parental lines," BMC Plant Biology , vol. 8, article 114, 2008.

[28] G. Wei, Y. Tao, G. Liu, C. Chen, R. Luo, H. Xia, Q. Gan, H. Zeng, Z. Lu, Y. Han, X. Li, G. Song, H. Zhai, Y. Peng, D. Li, H. Xu, X. Wei, M. Cao, H. Deng, Y. Xin, X. Fu, L. Yuan, J. Yu, Z. Zhu, L. Zhu, "A transcriptomic analysis of superhybrid rice LYP9 and its parents," Proceedings of the National Academy of Sciences of the United States of America , vol. 106, no. 19, pp. 7695-7701, 2009.

[29] G.-S. Song, H.-L. Zhai, Y.-G. Peng, L. Zhang, G. Wei, X.-Y. Chen, Y.-G. Xiao, L. Wang, Y.-J. Chen, B. Wu, B. Chen, Y. Zhang, H. Chen, X.-J. Feng, W.-K. Gong, Y. Liu, Z.-J. Yin, F. Wang, G.-Z. Liu, H.-L. Xu, X.-L. Wei, X.-L. Zhao, P. B. F. Ouwerkerk, T. Hankemeier, T. Reijmers, R. Van Der Heijden, C.-M. Lu, M. Wang, J. Van Der Greef, Z. Zhu, "Comparative transcriptional profiling and preliminary study on heterosis mechanism of super-hybrid rice," Molecular Plant , vol. 3, no. 6, pp. 1012-1025, 2010.

[30] N. Hoecker, B. Keller, H.-P. Piepho, F. Hochholdinger, "Manifestation of heterosis during early maize (Zea mays L.) root development," Theoretical and Applied Genetics , vol. 112, no. 3, pp. 421-429, 2006.

[31] A. Thiemann, J. Fu, T. A. Schrag, A. E. Melchinger, M. Frisch, S. Scholten, "Correlation between parental transcriptome and field data for the characterization of heterosis in Zea mays L," Theoretical and Applied Genetics , vol. 120, no. 2, pp. 401-413, 2010.

[32] W. Qian, R. Liu, J. Meng, "Genetic effects on biomass yield in interspecific hybrids between Brassica napus and B. rapa ," Euphytica , vol. 134, no. 1, pp. 9-15, 2003.

[33] F. Cheng, T. Mandakova, J. Wu, Q. Xie, M. A. Lysak, X. Wang, "Deciphering the diploid ancestral genome of the mesohexaploid Brassica rapa ," Plant Cell , vol. 25, no. 5, pp. 1541-1554, 2013.

[34] J. Zou, J. Zhu, S. Huang, E. Tian, Y. Xiao, D. Fu, J. Tu, T. Fu, J. Meng, "Broadening the avenue of intersubgenomic heterosis in oilseed Brassica ," Theoretical and Applied Genetics , vol. 120, no. 2, pp. 283-290, 2010.

[35] W. Qian, X. Chen, D. Fu, J. Zou, J. Meng, "Intersubgenomic heterosis in seed yield potential observed in a new type of Brassica napus introgressed with partial Brassica rapa genome," Theoretical and Applied Genetics , vol. 110, no. 7, pp. 1187-1194, 2005.

[36] A. Li, D. Liu, J. Wu, X. Zhao, M. Hao, S. Geng, J. Yan, X. Jiang, L. Zhang, J. Wu, L. Yin, R. Zhang, L. Wu, Y. Zheng, L. Mao, "mRNA and small RNA transcriptomes reveal insights into dynamic homoeolog regulation of allopolyploid heterosis in nascent hexaploid wheat," Plant Cell , vol. 26, no. 5, pp. 1878-1900, 2014.

[37] X. Wang, H. Wang, J. Wang, "The genome of the mesopolyploid crop species Brassica rapa ," Nature Genetics , vol. 43, no. 10, pp. 1035-1039, 2011.

[38] M. X. Zhao, B. Zhang, S. L. Liu, J. Ma, "Transposon expression and potential effects on gene regulation of Brassica rapa and B. oleracea genomes," Hereditas , vol. 35, no. 8, pp. 1014-1022, 2013.

[39] Z. An, Z. Tang, B. Ma, A. S. Mason, Y. Guo, J. Yin, C. Gao, L. Wei, J. Li, D. Fu, "Transposon variation by order during allopolyploidisation between Brassica oleracea and Brassica rapa ," Plant Biology , vol. 16, no. 4, pp. 825-835, 2014.