Use of digital gene expression to discriminate gene expression differences in early generations of resynthesized Brassica napus and its diploid progenitors

Digital gene expression (DGE) profile

To investigate global transcriptome during the polyploidization of B. napus, we performed a DGE analysis on the seedling stage of resynthesized B. napus (F₁–F₄) and its diploid parents, as well as natural B. napus. Finally, DGE libraries from the leaves of four-week-old plants were generated and sequenced by Illumina technology. The sequence data are available from the GEO repository, accession number GSE43246. The statistics of the DGE tags is shown in Table 1. Approximately six million raw tags were generated for each library, and more than 97% of the raw tags were clean tags. A total of 6178564 raw tags were obtained from B. rapa (Br), 6059222 from B. oleracea (Bo), 6155227 from B. napus- F₁ (Bn-F₁), 6092805 from B. napus- F₂ (Bn-F₂), 6142098 from B. napus- F₃ (Bn-F₃), 5938583 from B. napus- F₄ (Bn-F₄), and 5964594 from natural B. napus (Bn-N). After removing the low-quality sequences and adapter sequences, 6018254, 5930726, 6022170, 5950123, 5991210, 5798939, and 5823113 clean tags were obtained with 21 nt in length in the corresponding species. Unambiguous tags were counted and normalized to the number of transcripts per million tags (TPM) to evaluate the gene expression level. The results show that the mRNA transcribed from major genes had fewer than ten copies and only a small proportion of genes were highly expressed. The distribution of clean tags in the seven libraries was determined to evaluate the normality of the dataset. The results show a consistent pattern, with most of the tags coming from highly expressed genes. In this analysis, the total number of clean tags is the sum of all clean tags and the number of distinct clean tags is the number of different clean tags (Figure 1, Additional file 1: S1 and Additional file 2: S2). We found that the percentage of distinct tags with high counts dropped dramatically, and the distinct tags with more than 100 copies accounted less than 10%. However, more than 70% of total clean tags have accounts above 100 in each library. By contrast, 50% of the distinct clean tags had copy numbers between two and five, but they only represented around 4% of the total number of clean tags, which indicates only a small number of mRNA were expressed at high abundance and the majority were expressed at very low levels [38]. The clean tags were then mapped onto the B. rapa genome [34] and the numbers of tags that could be mapped onto genes with a maximum of one base pair mismatch in Br, Bo, Bn-F₁, Bn-F₂, Bn-F₃, Bn-F_4, and Bn-N were 1964909, 1747843, 3025405, 2354059, 1652733, 2295897, and 2253347, respectively (Table 1). Statistical analysis of clean tag alignment was conducted, including analysis of total clean tags and distinct clean tags (Additional file 3: S3, Additional file 2: S2). We found that around 70% of total clean tags were mapped onto the B. rapa genome with perfect match or with 1 bp mismatch to sense or anti-sense genes, and approximately 65% of the distinct clean tags were successfully mapped. Finally, the tag mapping onto the B. rapa genome generated 19023 tag-mapped genes for Br, 18547 for Bo, 24383 for Bn-F₁, 20659 for Bn-F₂, 18881 for Bn-F₃, 20692 for Bn-F₄, and 19955 for Bn-N. In total, 25959 genes were identified from the seven libraries, which accounted for approximately 63% of genes in the annotated B. rapa genome (Additional file 4: S4).

Summary		B. rapa	B. oleracea	B. napus- F1	B. napus- F2	B. napus- F3	B. napus- F4	B. napus- N
Raw Data	Total	6178564	6059222	6155227	6092805	6142098	5938583	5964594
Raw Data	Distinct Tag	293575	243895	331406	275597	284202	272986	269285
Clean tag	Total number	6018254	5930726	6022170	5950123	5991210	5798939	5823113
Clean tag	Distinct Tag number	133499	116771	199428	134403	134721	134857	128967
Tag Mapping to Gene	Total number	1964909	1747843	3025405	2354059	1652733	2295897	2253347
Tag Mapping to Gene	Distinct Tag number	44267	36220	88270	49781	42613	51122	45358
Unambiguous Tag Mapping to Gene	Total number	1679848	1475050	2495411	2016182	1415736	1971591	1924944
Unambiguous Tag Mapping to Gene	Total % of clean tag	27.91%	24.87%	41.44%	33.88%	23.63%	34.00%	33.06%
Unambiguous Tag Mapping to Gene	Distinct Tag number	39414	31933	79648	44456	37996	45762	40561
Unambiguous Tag Mapping to Gene	Distinct Tag % of clean tag	29.52%	27.35%	39.94%	33.08%	28.20%	33.93%	31.45%
Tag-mapped Genes	number	19023	18547	24383	20659	18881	20692	19955
Tag-mapped Genes	% of ref genes	46.20%	45.05%	59.22%	50.17%	45.86%	50.26%	48.47%
Unambiguous Tag-mapped Genes	number	16574	15970	22059	18155	16479	18196	17448
Unambiguous Tag-mapped Genes	% of ref genes	40.25%	38.79%	53.58%	44.09%	40.02%	44.19%	42.38%
Mapping to Genome	Total number	2437918	2105332	1776978	2290544	2658991	2232572	2164464
Mapping to Genome	Total % of clean tag	40.51%	35.50%	29.51%	38.50%	44.38%	38.50%	37.17%
Mapping to Genome	Distinct Tag number	44076	30703	49898	42334	45680	42590	40689
Mapping to Genome	Distinct Tag % of clean tag	33.02%	26.29%	25.02%	31.50%	33.91%	31.58%	31.55%
Unknown Tag	Total number	1615427	2077551	1219787	1305520	1679486	1270470	1405302
Unknown Tag	Total % of clean tag	26.84%	35.03%	20.25%	21.94%	28.03%	21.91%	24.13%
Unknown Tag	Distinct Tag number	45156	49848	61260	42288	46428	41145	42920
Unknown Tag	Distinct Tag % of clean tag	33.82%	42.69%	30.72%	31.46%	34.46%	30.51%	33.28%

Clean tags are tags after filtering dirty tags from raw data. Distinct tags are different kinds of tags and unambiguous tags are the reminder clean tags after removing tags mapped to reference genome.

To determine whether the sequencing depth was sufficient for the transcriptome coverage, a saturation analysis was performed to check whether the number of detected genes increased with sequencing amount until the number of reads reached 2 million (Additional file 5: S5). We also analyzed the distribution of the ratio of distinct tag copy numbers in each pair of libraries and found that more than 90% of the distinct tags had ratios up to 5-fold (Additional file 6: S6). Antisense genes play an important role in gene expression and regulation, and about 36% (14775) of the genes were matched by distinct tags in antisense orientations. Up to 5726 genes were mapped by DGE tags in both sense and antisense orientations in the Br library, 5123 in the Bo library, 11903 in the Bn-F₁ library, 6568 Bn-F₂, 5692 in the Bn-F₃ library, 6887 in the Bn-F₄ library, and 6253 in the Bn-N library (Additional file 4: S4). In total, more than 56% of the genes (14775 of 25959 genes) were transcribed from both strands, which indicate the importance of RNA-mediated gene regulation in polyploidization.

Differentially expressed genes in early generations of resynthesized B. napus

To identify genes that were differentially expressed in the early generations of synthesized B. napus after polyploidization, we compared pairs of DGE profiles of the seven libraries (Br versus Bn-F₁, Bo versus Bn-F₁, Bn-F₁ versus Bn-F₂, Bn-F₁ versus Bn-F₃, Bn-F₁ versus Bn-F₄, and Bn-N versus Bn-F_1, where A was the control and B was experimental group in ‘A versus B’) to analyze gene expression variations (Figure 2 and Additional file 7: S7). The comparison of B. napus- F₁ with B. rapa revealed that 4197 DEGs were significantly upregulated and 360 DEGs were downregulated in B. napus- F₁ in comparison with B. rapa. By contrast, 4554 DEGs were downregulated and 886 DEGs were upregulated in B. napus- F₁ compared with B. oleracea. The number of upregulated DEGs in B. napus-F₁ was more than that downregulated after polyploidization, which might indicate heterosis. Apart from the DEGs among resynthesized B. napus-F₁ and its progenitors, the number of the identified DEGs differed in the early generations of synthesized B. napus, and most of the DEGs showed downregulation in the F₂–F₄ generations compared with the F₁ generation. Up to 3022 DEGs were downregulated and 507 were upregulated in B. napus- F₂ compared with B. napus- F₁; 4718 DEGs were downregulated and 502 were upregulated in B. napus- F₃ compared with B. napus- F₁; 2882 DEGs were downregulated and 545 DEGs were upregulated in B. napus- F₄ and B. napus- F₁. Comparison of resynthesized and natural B. napus showed 649 tags were downregulated and 2916 tags were upregulated in B. napus- F₁ in comparison with natural B. napus. By comparison, the lowest number of DEGs was found in the F₄ generation, which was in a relatively stable stage after polyploidization. Figure 3 shows the distribution of genes commonly expressed in B. rapa, B. oleracea, B. napus-F₁, and natural B. napus, which indicates that a number of conserved genes were identified apart from the DEGs. Distribution of genes commonly and specifically expressed in Bn-F₁, F₂, F₃ and F₄ was also analyzed (Figure 4). Among the DEGs, the 20 most abundantly expressed genes that were upregulated or downregulated during polyploidization are listed in Additional file 8: S8. We found that the genes that encode DNA binding/transcription factor (Bra039065, Bra018796) were more prominent, which were downregulated in F₁ compared with its progenitors and then upregulated in the F₂–F₃ generations compared with the F₁ generation. Genes encoding cyclin-dependent protein kinase (Bra016640), epoxycarotenoid dioxygenase (Bra020970), and glycine-rich protein (Bra040817) were upregulated in F₁ during polyploidization, and then downregulated in F₂–F₄ compared with F₁.

Functional annotation of DEGs

Annotation of the sequences using GO database yielded good results for 22059 unambiguously mapped genes identified using the Bn-F₁ DGE tags. These well-annotated sequences belonged to three main categories (cellular component, molecular function, and biological process) and distributed into 33 categories, including the most dominant pathways such as ‘binding and catalytic activity,’ ‘cell and cell part,’ ‘organelle,’ ‘cellular and metabolic processes,’ and ‘response to stimulus’ (Figure 5). A similar GO distribution was revealed in the other six libraries (Additional file 9: S9). However, we did not find any genes in the transporter cluster in Bo. Genes were not found in the extracellular region cluster in Bn-F₂. Enzyme regulators and translation regulators were not observed in Bn-F₃.

To identify the biological pathways active in the DGE libraries, we mapped all the annotated genes to terms in the KEGG database to search for significantly enriched genes involved in metabolic or signal transduction pathways (Additional file 10: S10). Among all the genes with KEGG annotation, DEGs were identified in Bn-F₁ in comparison with Br. In total, we assigned 2911 DEGs to 121 KEGG pathways. Up to 38 of these pathways were significantly enriched, with Q values ≤ 0.05 (red border region), including metabolic pathways, ribosome, and carbon fixation in photosynthetic organisms. Similar pathway enrichment was revealed in pair comparison of each libraries (Bo vs. Bn-F₁, Bn-F₁ vs. Bn-F₂, Bn-F₁ vs. Bn-F₃, Bn-F₁ vs. Bn-F₄, and Bn-N vs. Bn-F_1, where A was the control and B was experimental group in ‘A vs. B’) . The 3339 DEGs identified in Bn-F₁ in comparison with Bo were assigned to 122 KEGG pathways, 44 of which were significantly enriched. The 2293 DEGs identified in Bn-F₂ in comparison with Bn-F₁ were assigned to 122 KEGG pathways, 3322 DEGs identified in Bn-F₃ in comparison with Bn-F₁ were assigned to 120 pathways, 2229 DEGs identified in Bn-F₄ in comparison with Bn-F₁ were assigned to 119 pathways, and 2283 DEGs identified in Bn-F₁ in comparison with Bn-N were assigned to 119 pathways.

Genome-wide non-additive gene regulation in synthesized Bn-F₁

Cluster analysis of DEGs

Cluster analysis of DEGs in the early generations of synthesized B. napus was performed with the correlated expression profile (Additional file 4: S4). Generally, 7948 DEGs in all the comparisons were clustered as the union of DEGs. Up to 1358 DEGs in all the comparisons were clustered as the DEG intersections (Figure 6). Among the 12 major clusters, the upregulated transcripts were enriched in clusters A–K in the comparisons of Bn-N vs. Bn-F₁, Br vs. Bn-F₁, and Bo vs. Bn-F₁. The downregulated transcripts were enriched in clusters A–K in the comparisons of Bn-F₁ vs. Bn-F₂, Bn-F₁ vs. Bn-F₃, and Bn-F₁ vs. Bn-F₄. Only one exception indicated that cluster L showed gene upregulation in Bn-F_2, Bn-F₃ and Bn-F₄ compared with Bn-F₁, and downregulation in Bn-F₁ compared with Bn-N, Br and Bo. We also identified the enrichment of DEGs in each cluster in different functional categories and found that the genes involved in the metabolism category were enriched significantly in all clusters except cluster H. Those in the protein with binding function or cofactor requirement category were enriched significantly in clusters A–H, whereas those in the subcellular localization category were enriched in all clusters except for cluster D. All the clustered DEGs are listed in Additional file 11: S11. Multiple genes represented for proteins with binding function, such as zinc ion binding (Bra019896), ATP binding (Bra008602, Bra023518), glutathione binding (Bra022815), selenium binding (Bra032756), and FAD/NADP/NADPH binding (Bra001931), showed different expression values in different generations of resynthesized B. napus and its progenitors. Some genes that encode ribosomal proteins were also changed during polyploidization.

Differences in gene expression between synthesized B. napus-F₁ and its progenitors

Based on analysis of the leaf transcriptome data for synthesized B. napus-F₁ and its progenitors, as well as the natural B. napus generated in this study, we found that the majority of the B. napus-F₁ transcripts were conserved with the parental sequences. Up to 19023, 18547, and 24383 sequences were mapped by B. napus-F₁ and its progenitors B. rapa and B. oleracea. 3182 genes with same TPM were expressed in Bn-F₁ and Br, and 2421 genes with same TPM were expressed in Bn-F₁ and Bo, which indicate that these genes might be significantly inherited by Bn-F₁ from one of its progenitors. The genes with different TPM expressed in B. napus-F₁ and its progenitors, which indicate that the genome combination results duplicated genes in the resynthesized polyploid. However, only 14 of these commonly expressed genes showed TPM values in Bn-F₁ equal to the combination of its progenitors. Furthermore, most of the commonly expressed genes in the progenitors were non-additively expressed in Bn-F₁, which might be responsible for gene repression and activation [5, 28, 39]. The majority of the non-additively expressed genes in Bn-F₁ displayed expression values different to its progenitors, indicating that the transcriptome was reconciled during polyploidization [40]. Non-additive gene regulation should be involved in various biological pathways, which may lead to subfunctionalization of duplicated genes [41]. The gene expression in the resynthesized Bn-F₁ was more complicated than the simple combination of two genomes, which might be due to the homologous recombination between closely related genomes (A- and C-genome) [42]. Xiong et al. [30] reported the variable chromosome instability of synthesized B. napus, and chromosome loss was compensated by the gain of homologous chromosomes [30]. Szadkowski et al. [43] reported the genome blender of synthesized B. napus during first meiosis, which resulted in the rearrangement of two genomes and the restructuring in further generations [43]. Thus, resynthesized polyploids derived from closely related species result in changes in the non-additive pattern of multiple gene expression, which explains the molecular bases for hybrid vigor [44].

Besides the commonly expressed genes in polyploids and diploid progenitors, there are some species specific genes. As shown in Figure 3, 487, 800, 3325 and 549 genes were specifically expressed in Br, Bo, Bn-F₁ and Bn-N, respectively. The number of Bn-F₁ specific genes was much higher than its progenitors, indicating many new transcripts emerged during polyploidization. This was consistent with Xu et al. [45] that 20 novel transcripts were detected in leaves of resynthesized Bn-F₁[45]. Based on the clustering analysis of DEGs, we found 1333 out of the 1358 interspecific DEGs were up-regulated in Bn-F₁ compared with Br and Bo. These up-regulated genes express proteins with different biological function, including salt tolerance protein, mitochondria/chloroplast membrance binding proteins, kinds of protein kinases, etc. Confirmation of these expression differences in further study is necessary for deep into the genetic causes of interest characteristics and improved qualities. Whereas 18 out of the 1358 interspecific DEGs were down-regulated in Bn-F₁ compared with Br and Bo. Only two genes (Bra005226 and Bra026817) were up-regulated in Bn-F₁ compared with Br, whereas down-regulated in Bn-F₁ compared with Bo. Five genes (Bra040268, Bra038457, Bra033201, Bra037320 and Bra039516) were up-regulated in Bn-F₁ compared with Bo, whereas down-regulated in Bn-F₁ compared with Br (Figure 6). All these gene expression differences are of great interest for studying the genome polyploidization of Brassica species.

Gene expression differences among early generations of synthesized B. napus (F₁-F₄) and natural B. napus

Hitherto, studies on polyploidy mechanism were conducted on natural and resynthesized B. napus using proteomic, transcriptome, and cytogenetic analyses [6, 10, 29, 31, 39, 43]. However, few studies that trace genomic changes (including cytosine methylation, DNA fragment loss, genetic rearrangement) and transcriptome changes in different generations of synthesized B. napus have been reported [9, 27, 28]. As mentioned above, many novel transcripts emerged in Bn-F₁ compared its progenitors. The number of Bn-F₁ specific genes was also higher than Bn-F₁, Bn-F₂, Bn-F₃ (Figure 4). Many of these new transcripts emerged in the first generation of resynthesized B. napus disappeared after successive self hybridizations, indicating the instability of resynthesized polyploids. Based on the clustering of DEGs found in the comparisons of Br vs. Bn-F₁, Bo vs. Bn-F₁, Bn-F₁ vs. Bn-F₂, Bn-F₁ vs. Bn-F₃, Bn-F₁ vs. Bn-F₄, and Bn-N vs. Bn-F₁, we easily found that most DEGs are significantly downregulated in Bn-F₂, Bn-F₃, and Bn-F₄ compared with Bn-F₁, but these DEGs were upregulated in Bn-F₁ compared with its progenitors and natural B. napus. This phenomenon also indicated the genomic instability after polyploidization, and the discovery of multiple gene expression differences may compensate for the limitation of proteomic analysis [37]. The upregulation of DEGs in synthesized Bn-F₁ during polyploidization could explain its improved characteristics compared with its progenitors. For example, the expression pattern of Bra022585, which encodes stress-responsive protein and Bra040633, which encodes defense-related protein, were changed during polyploidization. Some of these gene expression differences might be related to the DNA methylation or gene fragment losses during the combination of two genomes [45]. For example, Bra003009, which encodes O-methyltransferase 1, and Bra006323, which encodes an N2-dimethylguanosine tRNA methyltransferase family protein, were upregulated in Bn-F₁ during polyploidization and downregulated in the F₂–F₄ generations. In the pathway analysis, we found Bra003009 participated in the flavone and flavonol biosynthesis, which functioned in the transformation of kaempferide into kaempferol. This may indicate that kaempferol level was different in the synthesized B. napus. The changes in expression include many transcription factors, binding proteins, and ribosomal proteins, must be important for the merging of the A- and C-genomes in synthesized polyploids. For the 18 DEGs, which were down-regulated in Bn-F₁ compared with Br and Bo, we found they were up-regulated in Bn-F₂, Bn-F₃ and Bn-F₄ compared with Bn-F₁. This indicated that transcriptome recovery of genome balance functioned after polyploidization. Further research is needed to verify and determine the causes of these transcriptome changes and how they correlate with the phenotypic divergence and fertility of synthesized polyploids.

Plant material

Seeds of resynthesized B. napus (F₁ generation) and its successive selfing generations (F₂-F₄), its diploid progenitors B. rapa (cv. Aikangqing) and B. oleracea (cv. Zhonghua Jielan), were generously provided by the Oil Crops Research Institute of the Chinese Academy of Agricultural Sciences. Natural B. napus cv. Yang6 was provided by the Jiangsu Institute of Agricultural Science in the Lixiahe District (China). All the plants were cultivated in climate chambers at 25°C, a 16 h light: 8 h dark photoperiod, and 70% relative humidity. The first true leaves from three plants of each genotype were pooled at the same physiologic stage (28-day-old seedlings) and frozen at 80°C for use.

RNA preparation, Illumina RNA-sequencing and data processing

Total RNA was extracted from leaves using RNAiso Plus (Takara) according to the manufacture’s protocol. RNA concentrations were measured using a Qubit Fluorometer, and integrity was confirmed via 2100 Bioanalyzer (Agilent Technologies). The DGE libraries were prepared using Illumina Gene Expression Sample Prep Kit. Single-chain molecules were fixed onto a Solexa Sequencing Chip (flowcell) and sequenced by Illumina HiSeq™ 2000 System. Millions of raw, 35 bp sequences were generated. Image analysis, base calling, generation of raw tags, and the tags were counted using the Illumina pipeline [38]. Empty tags (no tag sequence between the adaptors), adaptors, low quality tags (tags containing one or more unknown nucleotides “N”), and tags with a copy number of 1 were removed from raw sequences to obtain clean tags (21 bp) containing CATG.

Mapping of reads to the reference sequence

To identify the gene expression patterns in early generations of synthesized B. napus, all clean tags were annotated by mapping to the sequenced genome of B. rapa[34] using the SOAP2 software, with a maximum of 1 nucleotide mismatch allowed [46]. All the tags mapped to reference sequences were filtered and the remaining tags were designated as ambiguous tags. Mapping events on both sense and antisense sequences were included in the data processing. For gene expression analysis, the number of expressed tags was calculated and then normalized to TPM (number of transcripts per million tags) [47, 48]; and the DEGs were screened and used for mapping and annotation [49, 50]. Gene annotation was conducted using Blast2GO [51]. The gene ontology (GO) categorization of all DEGs was displayed as three hierarchies for cellular component, molecular function, and biological process by searching in the GO database. Web Gene Ontology Annotation Plot (WEGO) was also used for GO classification of genes mapped by each DGE library [52]. Clustering analysis of differential gene expression pattern was also conducted using hierarchical clustering explorer (HCE) [53, 54]. In this study, statistical analysis of DEGs among libraries was performed using stringent value FDR ≤ 0.001 (false discovery rate) and the absolute value of |log2Ratio| ≥ 1 as the threshold of significant difference of gene expression. The DEGs in the hierarchical clustering were grouped into functional categories based on MIPS functional catalogue using Arabidopsis orthologues [55]. Pathway enrichment analysis of differential gene expression was conducted for further understanding gene function through blasting the KEGG database. A P-value of 0.05 was selected as the threshold for considering a gene set as significantly enriched.