Abstract

Background

   Brassica includes many successfully cultivated crop species of
   polyploid origin, either by ancestral genome triplication or by
   hybridization between two diploid progenitors, displaying complex
   repetitive sequences and transposons. The U’s triangle, which consists
   of three diploids and three amphidiploids, is optimal for the analysis
   of complicated genomes after polyploidization. Next-generation
   sequencing enables the transcriptome profiling of polyploids on a
   global scale.

Results

   We examined the gene expression patterns of three diploids (Brassica
   rapa, B. nigra, and B. oleracea) and three amphidiploids (B. napus, B.
   juncea, and B. carinata) via digital gene expression analysis. In
   total, the libraries generated between 5.7 and 6.1 million raw reads,
   and the clean tags of each library were mapped to 18547–21995 genes of
   B. rapa genome. The unambiguous tag-mapped genes in the libraries were
   compared. Moreover, the majority of differentially expressed genes
   (DEGs) were explored among diploids as well as between diploids and
   amphidiploids. Gene ontological analysis was performed to functionally
   categorize these DEGs into different classes. The Kyoto Encyclopedia of
   Genes and Genomes analysis was performed to assign these DEGs into
   approximately 120 pathways, among which the metabolic pathway,
   biosynthesis of secondary metabolites, and peroxisomal pathway were
   enriched. The non-additive genes in Brassica amphidiploids were
   analyzed, and the results indicated that orthologous genes in
   polyploids are frequently expressed in a non-additive pattern.
   Methyltransferase genes showed differential expression pattern in
   Brassica species.

Conclusion

   Our results provided an understanding of the transcriptome complexity
   of natural Brassica species. The gene expression changes in diploids
   and allopolyploids may help elucidate the morphological and
   physiological differences among Brassica species.

Electronic supplementary material

   The online version of this article (doi:10.1186/s12870-015-0417-5)
   contains supplementary material, which is available to authorized
   users.

   Keywords: Brassica spp, Polyploidization, Sequencing, Digital gene
   expression (DGE)

Background

   Polyploidy is an important factor in the evolution of many plants and
   has attracted considerable scientific attention for a long period of
   time. Many important economical crops are of polyploid origin,
   including wheat, cotton, and rapeseed [[35]1]. Cruciferae includes the
   model species Arabidopsis thaliana and the economically important
   Brassica crops. These important crops include three diploid Brassica
   species, namely, B. rapa (AA, 2n = 20; Chinese cabbage, turnip, turnip
   rape), B. nigra (BB, 2n = 16; black mustard), and B. oleracea (CC,
   2n = 18; cauliflower, broccoli, kale), and three allopolyploids
   spontaneously derived from pairwise hybridization of the diploids,
   which are B. napus (AACC, 2n = 38; oilseed rape, swede), B. juncea
   (AABB, 2n = 36; abyssinian or Ethiopian mustard), and B. carinata
   (BBCC, 2n = 34; Indian or brown mustard) [[36]2]. Lysak et al. (2005)
   confirmed the chromosome triplication history of Brassica that
   corresponds to that of Arabidopsis [[37]3]. Cheung et al. (2009) found
   that the divergence of Arabidopsis and Brassica lineage was
   approximately 17 Mya [[38]4], and the replicated Brassica subgenomes
   (probably the divergence of A/C from B genome) was diverged by 14.3 Mya
   [[39]4]. In addition, the A and C genomes were estimated with 3.7 Mya
   divergence. Woodhouse et al. (2014) stated that the B. rapa lineage
   underwent a whole-genome triplication of 5–9 Mya [[40]5]. For the
   allopolyploids, B. napus probably arose from the natural hybridization
   of A and C genomes around 10,000 years ago. However, when the
   hybridization between A and B genomes and between B and C genomes
   happened is still unclear. The precise ancestors of B. napus, B.
   juncea, and B. carinata are not yet identified [[41]6]. The duplication
   of gene segments reported on Brassica is explained as frequent loss,
   remote genome duplication, or unbalanced homologous recombination
   [[42]7]. During the divergence of Brassica species, the
   sub-functionalization and/or neo-functionalization of these paralogs
   coupled with novel gene interactions contribute significantly to genome
   evolution [[43]8]. Moreover, genetic mapping and sequencing analysis
   substantiate the triplication hypothesis of diploid Brassica genomes
   [[44]9-[45]12]. The comparative mapping of Brassica by using genetic
   markers has successfully revealed homologous rearrangements,
   translocations, and fusions that are crucial to the diversification of
   the A, B, and C genomes from A. thaliana [[46]13-[47]15].

   Many linkage maps and karyotype analysis have identified extensive
   collinearity and genomic polymorphisms among Brassica genomes. Given
   the complexity of the gene copy number and syntenic conservation caused
   by polyploidization, Brassica genomes are difficult to study
   [[48]16,[49]17]. Identifying the genes related to specific traits based
   on the linkage maps is also challenging because of the complexity of
   the homologs and paralogs in polyploidy genomes [[50]15,[51]18].
   Profiling arrays of A. thaliana are useful in the transcriptome
   analysis of Brassica [[52]6]. However, A. thaliana-based microarrays
   lack the resolution of Brassica specific genes and paralogous genes.
   Furthermore, Brassica microarrays were developed to confirm
   Brassica-specific expressed genes [[53]19]. Identifying different
   homologous copies of Brassica sequences is challenging for microarray
   expression analysis [[54]20]. Next-generation sequencing is an optimal
   method for genomic and transcriptomic studies and provides
   opportunities for polyploidy studies and enables the extensive genome
   profiling of crops with complex genomes, such as soybean, potato,
   tomato, cotton, maize, and common bean [[55]21-[56]26]. This technology
   also promotes sequencing analysis in Brassica; the genome sequence of
   B. rapa has already been released and annotated [[57]12]. The genome
   sequencing of B. oleracea, B. nigra, and B. napus is still in progress.
   However, the genome sequences of B. oleracea are available in the Basic
   Local Alignment Search Tool in Brassica database
   ([58]www.brassica.info). The transcriptome profiling of B. napus has
   been analyzed via RNA sequencing [[59]27-[60]29]. This information is
   valuable for the investigation of Brassica genome evolution. Many
   technologies have been applied to quantify transcript abundance,
   including microarray, serial analysis of gene expression, digital gene
   expression (DGE), and RNA-seq. DGE and RNA-seq have been widely used to
   identify the molecular information of plant transcriptome and gene
   expression variation between comparable samples. DGE, as a well-known
   technique suitable to directly quantify transcript abundance counts, is
   optimized over RNA-seq because of its cost efficiency. RNA-seq is a
   flexible approach that can detect full-transcript sequence, alternative
   splicing, exon boundaries, and transcript abundance. However, each
   transcript in RNA-seq can be mapped multiple times, and the sequencing
   depth of RNA-seq is correlated with but is not equal to transcript
   abundance. Each read in DGE is expected with a sole hit on an RNA
   molecule. Therefore, DGE is better to represent rare transcripts or
   exclude transcripts with less interest than RNA-seq [[61]30].

   Many studies have analyzed the genomic and phenotypic changes in
   synthesized Brassica, particularly B. napus and hexaploid Brassica
   [[62]31-[63]33]. However, limited information is available for the
   natural species of Brassica. In the present research, we performed DGE
   analysis on three diploid Brassica species (B. rapa, B. nigra, and B.
   oleracea) and three allopolyploids (B. napus, B. juncea, and B.
   carinata) to determine the transcriptome changes after natural
   polyploidization. The expression profile of the genes in the six
   Brassica species was reported, and the multiple gene expression
   differences were observed. Differentially expressed genes (DEGs) are
   involved in a wide range of stress resistance and development
   processes. This study is the first transcriptomic research that
   identifies DEGs and the pathways involved in the natural
   polyploidization of the six Brassica species.

Results

DGE profile

   This research investigates the transcriptome profiling of diploids and
   spontaneous allopolyploids in Brassica by performing DGE analysis on
   the seeding stage of the six Brassica species, namely, B. rapa (Br), B.
   nigra (Bg), B. oleracea (Bo), B. napus (Bn), B. juncea (Bj), and B.
   carinata (Bc). DGE libraries from the leaves of four-week-old plants
   were generated and sequenced by an Illumina technology. The sequence
   data are available from the GEO repository with an accession number of
   [64]GSE43246. The statistics of the DGE tags are shown in Table [65]1.
   Approximately six million raw tags were generated for each library.
   Clean tags were obtained after removing the low-quality sequences and
   adaptor sequences from the raw data. 6178564, 5881618, 6059222,
   5964594, 6076830, and 5795234 clean tags were obtained in Br, Bg, Bo,
   Bn, Bj, and Bc, respectively. Unambiguous tags (tags that were uniquely
   matched to one gene of reference genome with no more than one mismatch)
   were counted and normalized to TPM to evaluate the gene expression
   level. To evaluate the normality of the DGE data, the distribution of
   the total tags and distinct clean tags (tags with specific nucleotide
   sequence) over different tag copy numbers was analyzed. The
   distribution of the tag expression was similar for each library.
   Moreover, the distribution of clean tags in the six libraries showed
   that most of the tags are from highly expressed genes (Figure [66]1,
   Additional files [67]1 and [68]2). The percentage of distinct tags with
   high counts dropped dramatically, and the distinct tags with more than
   100 copies accounted for less than 8%. However, more than 67% of the
   total clean tags accounted for more than 100 copies in each library. By
   contrast, more than 43% of the distinct clean tags had copy numbers
   between two and five, which represented approximately 96% of the total
   number of clean tags. Generally, a small number of categories of mRNA
   showed high abundance, whereas the other majority had a quite low
   expression level. This finding indicates that only a small number of
   mRNAs are expressed at high abundance, and majority of them are
   expressed at very low levels [[69]34].

Table 1.

   Statistics of categorization and abundance of DGE tags
   Summary B. rapa B. nigra B. oleracea B. napus B. juncea B. carinata
   Raw Tag Total 6178564 5881618 6059222 5964594 6076830 5795234
   Raw Tag Distinct Tag 293575 214427 243895 269285 400134 278768
   Clean Tag Total number 6018254 5772449 5930726 5823113 5858527 5657697
   Clean Tag Distinct Tag number 133499 106552 116771 128967 181965 142281
   Tag Mapping to Gene Total number 1964909 1990442 1747843 2253347
   1857572 1915305
   Tag Mapping to Gene Distinct Tag number 44267 30413 36220 45358 56289
   40425
   Unambiguous Tag Mapping to Gene Total number 1679848 1635594 1475050
   1924944 1531974 1594991
   Unambiguous Tag Mapping to Gene Total% of clean tag 27.91% 28.33%
   24.87% 33.06% 26.15% 28.19%
   Unambiguous Tag Mapping to Gene Distinct Tag number 39414 26114 31933
   40561 49892 35285
   Unambiguous Tag Mapping to Gene Distinct Tag% of clean tag 29.52%
   24.51% 27.35% 31.45% 27.42% 24.80%
   Tag-mapped Genes number 19023 16687 18547 19955 21995 19436
   Tag-mapped Genes % of ref genes 46.20% 40.53% 45.05% 48.47% 53.42%
   47.20%
   Unambiguous Tag-mapped Genes number 16574 13867 15970 17448 19424 16645
   Unambiguous Tag-mapped Genes % of ref genes 40.25% 33.68% 38.79% 42.38%
   47.18% 40.43%
   Mapping to Genome Total number 2437918 1147106 2105332 2164464 2047451
   1462061
   Mapping to Genome Total% of clean tag 40.51% 19.87% 35.50% 37.17%
   34.95% 25.84%
   Mapping to Genome Distinct Tag number 44076 15159 30703 40689 50304
   29547
   Mapping to Genome Distinct Tag% of clean tag 33.02% 14.23% 26.29%
   31.55% 27.64% 20.77%
   Unknown Tag Total number 1615427 2634901 2077551 1405302 1953504
   2280331
   Unknown Tag Total% of clean tag 26.84% 45.65% 35.03% 24.13% 33.34%
   40.30%
   Unknown Tag Distinct Tag number 45156 60980 49848 42920 75372 72309
   Unknown Tag Distinct Tag% of clean tag 33.82% 57.23% 42.69% 33.28%
   41.42% 50.82%
   [70]Open in a new tab

   Clean tags are tags after filtering low-quality tags from raw data.
   Distinct tags are different tags and unambiguous tags are the remaining
   clean tags after removing tags mapped to more than one locus of
   reference genome.

Figure 1.

   Figure 1
   [71]Open in a new tab

   Distribution of total tag and distinct tag counts over different tag
   abundance categories from the six libraries.

   The clean tags were then mapped onto the B. rapa genome with a maximum
   of one base-pair mismatch [[72]12]. Table [73]1 shows that the 1964909,
   1990442, 1747843, 2253347, 1857572, and 1915305 tags in Br, Bg, Bo, Bn,
   Bj, and Bc were mapped to B. rapa genome, respectively. Statistical
   analysis of clean tag alignment was conducted, including the analysis
   of total clean tags and distinct clean tags (Additional files [74]2 and
   [75]3). More than 54% of the total clean tags and 42% of the distinct
   clean tags in each library were mapped onto the B. rapa genome.
   However, the tags mapped in the DGE library of Bg and Bc were lower
   than those in the other four libraries, which might be due to the
   divergence of the B genome to the A/C genome. Moreover, the tag mapping
   onto the B. rapa genome generated 19023 tag-mapped genes for Br, 16687
   for Bg, 18547 for Bo, 19955 for Bn, 21995 for Bj, and 19436 for Bc. In
   total, approximately 61% of the genes in the B. rapa genome (25298
   genes) could be mapped with unique tags (Additional file [76]4).
   Furthermore, we mapped all the clean tags of each DGE library to the
   genome of A. thaliana, and the summary and details of the mapping
   result are listed in Additional file [77]5. Only approximately 47% of
   A. thaliana genes (19557 genes) were successfully mapped, and the
   percent of unambiguous tag-mapped genes in A. thaliana is much lower
   than in B. rapa. The number of DGE tags in each library that well
   matched with Arabidopsis genome is also lower than that mapped to B.
   rapa. The difference in mapping rate is in accordance with the
   prediction that the A, B, and C genomes of Brassica diverged after the
   divergence of Arabidopsis and Brassica lineages [[78]6]. Thus, we chose
   the mapping information that used B. rapa as reference for further
   analysis. Saturation analysis was performed to check if the number of
   detected genes increased with sequencing amount. The result showed that
   the number of detected genes stopped increasing when the number of
   reads reached 2 million (Additional file [79]6). The distribution of
   the ratio of distinct tag copy numbers in each pair of libraries was
   analyzed. More than 90% of the distinct tags had ratios up to five
   folds (Additional file [80]7).

DEGs in Brassica diploids

   The DEGs in Brassica diploids (Br, Bg, and Bo) were compared (Br vs.
   Bo, Bg vs. Br, Bg vs. Bo, where A was the control group and B was the
   experimental group in “A vs. B”) to analyze the gene expression
   variations (Figure [81]2 and Additional file [82]8). A comparison of Br
   and Bo showed that 1352 and 1282 DEGs were significantly up-regulated
   and down-regulated, respectively. By contrast, 2278 DEGs were
   down-regulated and 2391 DEGs were up-regulated in Br compared with Bg.

Figure 2.

   Figure 2
   [83]Open in a new tab

   Number of differentially expressed genes in each comparison of Brassica
   diploids. The numbers of up-regulated (in red) and down-regulated genes
   (in green) are presented. Br, Bg and Bo are abbreviations of B. rapa,
   B. nigra and B. oleracea, respectively.

   Moreover, 2140 DEGs were down-regulated and 2437 DEGs were up-regulated
   in Bo compared with Bg. The number of DEGs in Bg compared with Br/Bo
   was more than Br vs. Bo, which indicates that the A and C genomes of
   Brassica were closer than the B genome. Among the 20 most abundantly
   expressed genes that were up-regulated or down-regulated in the pair
   comparison of the three diploids (Additional file [84]8), Bra015187,
   Bra026992, Bra017452, Bra029372, Bra028406, Bra017112, Bra036352,
   Bra000377, and Bra016934 were up-regulated in Bg compared with Br/Bo.
   Moreover, Bra023103, Bra011285, Bra014371, Bra031070, Bra028805, and
   Bra006083 were down-regulated in Bg compared with Br/Bo. Most DEGs
   between Br and Bo differed from those between Br and Bg as well as
   between Bo and Bg. Figure [85]3A shows the distribution of the genes
   commonly expressed in Br, Bg, and Bo, and 8932 genes were co-expressed
   in the three diploid libraries apart from the DEGs. A second filter of
   expression differences (at least twofold or greater) was performed in
   the pairwise comparisons to reduce the total number of significant
   changes. This analysis revealed 6401 significantly expressed genes,
   such as Br vs. Bg = 4669, Br vs. Bo = 2634, and Bg vs. Bo = 4577
   (Figure [86]3B). The functional significance of the genes expressed in
   each library was examined, and the GO analysis results are shown in
   Figure [87]3C. The well-annotated gene sequences were assigned to 33
   functional groups of the three main GO categories (cellular component,
   molecular function, and biological process). The results showed that
   DGEs were primarily involved in the cell and organelle, in the binding,
   catalytic, cellular, and metabolic process, as well as in response to
   stimulus. Two specific processes, the extracellular region part and the
   molecular transducer, were unique to Bo. However, Bo lacked a
   transporter, and Bg lacked anatomical structure formation.

Figure 3.

   Figure 3
   [88]Open in a new tab

   Distribution of expressed mRNAs in Brassica diploids. A. Venn diagram
   of genes expressed in Br, Bg and Bo. B. Venn diagram of unique
   expressed genes between pairwise of Br, Bg and Bo. C. GO classification
   of genes in Br, Bg and Bo.

DEGs among allopolyploids and ancestral diploid progenitors

   To identify the DEGs in allopolyploids and ancient diploid progenitors,
   the DGE profiles of Br vs. Bn, Bo vs. Bn, Br vs. Bj, Bg vs. Bj, Bg vs.
   Bc, and Bo vs. Bc were compared to analyze the gene expression
   variations after natural polyploidization (Figure [89]4 and Additional
   file [90]8). The results showed that 1230 DEGs were up-regulated and
   324 DEGs were down-regulated in Bn compared with Br, whereas 1872 DEGs
   were up-regulated and 797 DEGs were down-regulated in Bn compared with
   Bo. After natural polyploidization, 1519 DEGs were induced in Bj
   compared with Br, whereas 508 DEGs were down-regulated. Moreover, 2692
   DEGs were induced in Bj compared with Bg, whereas 1393 DEGs were
   down-regulated. With regard to Bc, 2099 DEGs were up-regulated and 1344
   were down-regulated compared with Bg, and 1691 DEGs were up-regulated
   and 1070 were down-regulated compared with Bg. The variations in the
   gene expression among the diploids and amphidiploids are essential to
   the diversity of phenotype, growth, and production. The 20 most
   abundantly expressed genes that were up-regulated or down-regulated in
   the pair comparison of amphidiploids and diploids are listed in
   Additional file [91]8. The distribution of the genes that were commonly
   and uniquely expressed in amphidiploid and its ancestral diploids is
   shown in Figure [92]5. The results further show that 11810 genes were
   conserved in Br, Bo, and Bn, whereas 1362, 1666, and 1824 genes were
   specifically expressed in Br, Bo, and Bn, respectively (Figure [93]5A).
   A similar pattern was observed when Bj was compared with Br/Bg
   (Figure [94]5B) and Bc with Bg/Bo (Figure [95]5C). The expressional
   differences of species-specific genes might be caused by the genome
   interaction during natural polyploidization. The GO pattern of the
   genes in amphidiploid and ancestral diploids is shown in Figure [96]6.
   Based on Figure [97]6A, the numbers of DGEs in the envelope,
   extracellular region, macromolecular complex, biological regulation,
   cellular component biogenesis, death, multicellular organism process,
   and pigmentation were different in Br, Bo, and Bn. GOs of molecular
   transducer was found in Bo only. Apparent GO difference was observed
   among Br, Bg, and Bj (Figure [98]6C). As shown in Figure [99]6C, GOs of
   transporter were found in Bg only, and anatomical structure formation
   was not present in Bg.

Figure 4.

   Figure 4
   [100]Open in a new tab

   Number of differentially expressed genes in comparison of diploids and
   amphidiploids. The numbers of up-regulated (in red) and down-regulated
   genes (in green) are presented. Br, Bg, Bo, Bn, Bj and Bc are
   abbreviations of B. rapa, B. nigra, B. oleracea, B. napus, B. juncea
   and B. carinata, respectively.

Figure 5.

   Figure 5
   [101]Open in a new tab

   Distribution of the genes commonly and specifically expressed in
   diploids and amphidiploids. A. Venn diagram of genes expressed in Br,
   Bo and Bn. B. Venn diagram of genes expressed in Br, Bg and Bj. C. Venn
   diagram of genes expressed in Bg, Bo and Bc.

Figure 6.

   Figure 6
   [102]Open in a new tab

   GO classification of genes in diploids and amphidiploids. A. GO
   classification of genes expressed in Br, Bo and Bn. B. GO
   classification of genes expressed in Br, Bg and Bj. C. GO
   classification of genes expressed in Bg, Bo and Bc.

Functional annotation of DEGs

   Pathway enrichment analysis was performed on the expressed transcripts
   of the six DGE libraries. This analysis was performed by mapping all
   annotated genes in the KEGG database to search for significantly
   enriched genes involved in the metabolic or signal transduction
   pathways (Additional file [103]9). Among the genes with KEGG
   annotation, DEGs were identified in Bn compared with Br. In total, 894
   DEGs were assigned to 109 KEGG pathways, and 13 of these pathways were
   significantly enriched with Q values ≤ 0.05 (red border region). The
   enriched pathways include metabolic, biosynthesis of secondary
   metabolites, and peroxisome. A similiar pathway enrichment was
   discovered in pair comparison of other libraries (Bo vs. Bn, Br vs. Bj,
   Bg vs. Bj, Bg vs. Bc, and Bo vs. Bc). The 1562 DEGs identified in Bn
   vs. Bo were assigned to 122 KEGG pathways, 15 of which were
   significantly enriched. The 1171 DEGs identified in Bj vs. Br were
   assigned to 116 KEGG pathways, the 2373 DEGs identified in Bj vs. Bg
   were assigned to 121 pathways, the 1975 DEGs identified in Bc vs. Bg
   were assigned to 120 pathways, and the 1639 DEGs identified in Bc vs.
   Bo were assigned to 117 pathways. All these pathways are involved in
   the process of plant growth and stress reaction, which are important
   for the morphological and physiological differences among the Brassica
   species. The biosynthesis of unsaturated fatty acids, which was
   significantly enriched in Bo vs. Bn and Bg vs. Bc, explains the
   different fatty acid compositions in Brassica species
   [[104]35,[105]36]. The DEGs identified in the peroxisome pathway were
   related to the improved stress ability of Bn and Bj.

Clustering of DEGs

   Hierarchical cluster analysis of the DEGs in Br, Bg, Bo, Bn, Bj, and Bc
   were performed to examine the similarity and diversity of gene
   expression (Additional file [106]4). All results were displayed by Java
   Treeview, where red and green represent the positive and negative
   values of gene expression, respectively. Generally, 651 genes with
   differential expression in Br, Bg, and Bo were clustered as DEG
   intersections (Figure [107]7A, Additional file [108]10). The comparison
   of Br, Bg, and Bo showed that 5417 DEGs were clustered as the union of
   DEGs (Additional file [109]11). Moreover, 285 DEGs in Bn, Br, and Bo
   were also clustered (Figure [110]7B and Additional file [111]9), and
   3786 DEGs differentially expressed in Bn and Br/Bo are listed in
   Additional file [112]11. The 630 DEGs in Bj, Br, and Bg were also
   clustered (Figure [113]7C and Additional file [114]9), and 5590 DEGs
   differentially expressed in Bj and Br/Bg are listed in Additional file
   [115]11. In addition, 726 DEGs in Bc, Bg, and Bo were clustered
   (Figure [116]7D and Additional file [117]9), and 5264 DEGs
   differentially expressed in Bc and Bg/Bo are listed in Additional file
   [118]11. The functional categories of these DEGs showed similar
   enrichment patterns for each comparison, including categories of
   metabolism, development, cellular transport, and interaction with the
   environment (data not shown). Genes with binding function were enriched
   in each comparison, which is different from previous reports
   [[119]32,[120]33].

Figure 7.

   Figure 7
   [121]Open in a new tab

   Hierarchical cluster analysis of differentially expressed transcripts.
   A. Clustering of genes expressed in diploids of Brassica. B. Clustering
   of genes expressed in Br, Bo and Bn. C. Clustering of genes expressed
   in Br, Bg and Bj. D. Clustering of genes expressed in Bg, Bo and Bc.

Analysis of methyltransferase genes differentially expressed in Brassica

   Epigenetic variation has an important function in the evolution of
   plants. DNA methylation is included in this variation and has received
   much attention in previous years. Proteins such as methyltransferase
   are considered important for plant methylation [[122]37,[123]38]. Thus,
   the putative methyltransferase and methylase genes from all DEGs in
   each comparison were filtered to understand the mechanism of the
   changes in DNA methylation in Brassica (Additional file [124]12). Two
   methyltransferase genes (Bra003928 and Bra020452) were differentially
   expressed in Br, Bg, and Bo, and 30 genes exhibited differential
   expression in Br vs. Bo/Bg vs. Bo/Bg vs. Br. One methyltransferase gene
   (Bra008507) was differentially expressed in Bn, Br, and Bo, and 23
   genes exhibited differential expression in Br vs. Bn/Bo vs. Bn/Br vs.
   Bo. Five methyltransferase genes (Bra003396, Bra004391, Bra010977,
   Bra022603, and Bra024271) were differentially expressed in Bj, Br, and
   Bg, and 36 genes exhibited differential expression in Br vs. Bj/Bg vs.
   Bj/Bg vs. Br. Three methyltransferase genes (Bra003928, Bra004391, and
   Bra012494) were differentially expressed in Bc, Bg, and Bo, and 33
   genes exhibited differential expression in Bg vs. Bc/Bo vs. Bc/Bg vs.
   Bo. The results showed that Bra003928 was significantly down-regulated
   in Br compared with Bg/Bo, which was up-regulated in Bn compared with
   Br and down-regulated in Bn compared with Bo. The expression of
   Bra003928 in Bj was higher than in Br and lower than in Bg. The
   expression of this methyltransferase gene in Bc was significantly
   up-regulated than in Bg and Bo. Moreover, Bra020452 was obviously
   down-regulated in Bo compared with Br/Bg. Different expression values
   were also examined in Brassica amphidiploids compared with their
   ancestral diploid parents. The methyltransferase gene was up-regulated
   in Bn compared with Br and Bo, which was also up-regulated in Bc
   compared with Bg and Bo. However, the expression value of Bra020452 in
   Bj was similar to that of Br and Bg.

Non-additive genes expressed in Brassica amphidiploids

   The expression value of genes in Brassica amphidiploids (Bn, Bj, and
   Bc) were compared with the relative mid-parent value (MPV) to identify
   the genes with differential expression pattern. Up to 19844 genes in Bn
   showed differences in expression from MPV, 9605 (48.4%) of these genes
   showed higher expression levels, whereas 10239 (51.6%) showed lower
   expressions than MPV. Among the non-additively expressed genes, 9519
   (48%) genes were expressed at higher levels, whereas 10325 (52%) genes
   were expressed at lower levels in Br than in Bo (Table [125]2). This
   finding is similar to the data reported by Jiang et al. (2013)
   [[126]32]. With regard to Bj, 20317 genes showed differences in
   expression from MPV, 11173 (55%) of these genes were expressed higher
   in Br than in Bg, and 9144 (45%) genes were expressed at lower levels.
   Moreover, 19921 genes in Bc showed differences in expression from MPV,
   8189 (46.1%) of them were expressed higher in Bg than in Bo, whereas
   10732 (53.9%) genes were expressed lower. Significantly more
   non-additive genes than additive genes in amphidiploids implied the
   complicated evolution history of Brassica. In this study, no standard
   RNA sample was included in library construction. Given that two samples
   often differ in the total number of transcripts per cell, the
   differences in transcriptome size, not just the differences in RNA
   yields during extraction, can introduce biases [[127]39-[128]41]. In
   addition, polyploidization of plant species is a complicated process
   that is unequal to simple duplication or combination of ancient
   genomes; fractionation of duplicated genes would result in both gene
   and genome preferences in stabilized Brassica polyploids [[129]5]. The