Transcriptome and methylome profiling reveals relics of genome dominance in the mesopolyploid Brassica oleracea
- Isobel AP Parkin1Email author,
- Chushin Koh2,
- Haibao Tang3,
- Stephen J Robinson1,
- Sateesh Kagale1, 2,
- Wayne E Clarke1,
- Chris D Town3,
- John Nixon1,
- Vivek Krishnakumar3,
- Shelby L Bidwell3,
- France Denoeud4, 5, 6,
- Harry Belcram7,
- Matthew G Links1,
- Jérémy Just7,
- Carling Clarke2,
- Tricia Bender1,
- Terry Huebert1,
- Annaliese S Mason8,
- J Chris Pires9,
- Guy Barker10,
- Jonathan Moore10,
- Peter G Walley10,
- Sahana Manoli8,
- Jacqueline Batley8,
- David Edwards8,
- Matthew N Nelson11,
- Xiyin Wang12,
- Andrew H Paterson12,
- Graham King13,
- Ian Bancroft14,
- Boulos Chalhoub7Email author and
- Andrew G Sharpe2Email author
© Parkin et al.; licensee BioMed Central Ltd. 2014
Received: 27 January 2014
Accepted: 10 June 2014
Published: 10 June 2014
Brassica oleracea is a valuable vegetable species that has contributed to human health and nutrition for hundreds of years and comprises multiple distinct cultivar groups with diverse morphological and phytochemical attributes. In addition to this phenotypic wealth, B. oleracea offers unique insights into polyploid evolution, as it results from multiple ancestral polyploidy events and a final Brassiceae-specific triplication event. Further, B. oleracea represents one of the diploid genomes that formed the economically important allopolyploid oilseed, Brassica napus. A deeper understanding of B. oleracea genome architecture provides a foundation for crop improvement strategies throughout the Brassica genus.
We generate an assembly representing 75% of the predicted B. oleracea genome using a hybrid Illumina/Roche 454 approach. Two dense genetic maps are generated to anchor almost 92% of the assembled scaffolds to nine pseudo-chromosomes. Over 50,000 genes are annotated and 40% of the genome predicted to be repetitive, thus contributing to the increased genome size of B. oleracea compared to its close relative B. rapa. A snapshot of both the leaf transcriptome and methylome allows comparisons to be made across the triplicated sub-genomes, which resulted from the most recent Brassiceae-specific polyploidy event.
Differential expression of the triplicated syntelogs and cytosine methylation levels across the sub-genomes suggest residual marks of the genome dominance that led to the current genome architecture. Although cytosine methylation does not correlate with individual gene dominance, the independent methylation patterns of triplicated copies suggest epigenetic mechanisms play a role in the functional diversification of duplicate genes.
The Brassicaceae family encompasses extensive species diversity with a wide range of intra- and inter-specific morphological and phytochemical profiles, which has contributed to their global importance in crop production. The most agriculturally significant species are found within the Brassiceae tribe where selection for targeted traits in both seed and vegetative tissues has generated high-value oilseed, vegetable and condiment crops. Brassica oleracea represents an important species integral to human diets that encompasses multiple cultivar groups that are classified based on the specialized morphology of their edible structures, namely kales, cabbages, broccoli, cauliflower, Brussels sprouts and kohl rabi. In addition to this radical morphological variation, the cultivars and their related wild species accumulate a range of secondary metabolites that have been implicated in promoting human health, such as compounds functioning as anti-carcinogens and anti-oxidants. B. oleracea (n = 9, CC genome) along with its close relative Brassica rapa (n = 10, AA genome) are the extant descendants of the two diploid species, which fused to form the amphidiploid Brassica napus (canola or oilseed rape; n = 19, AACC genome), arguably the most economically valuable species among the Brassiceae.
Brassica research has been enhanced through the recent publication of the B. rapa genome sequence . In the current report we present the genome of the close relative B. oleracea, which diverged from B. rapa some 4 million years ago (Mya) . The completion of the genome sequences of these diploid species and the anticipated assembly of the amphidiploid B. napus genome will provide foundational resources for not only crop improvement but also studying polyploid evolution. The importance of polyploidy as a driving force in plant speciation, diversification and evolution has been widely acknowledged . The Brassica species carry remnants of three whole genome duplication (WGD) events, the α, β and γ events , the most recent of which, the α event, dates from approximately 47 Mya. In addition, the Brassica diploid species are mesopolyploids, having experienced an hexaploidization event approximately 23 Mya . As such, Brassica species offer a unique opportunity to study the impacts of WGD on genome evolution, and in particular the mechanisms employed to mitigate or exploit the impact of maintaining multiple copies of each gene. Analysis of the B. rapa genome suggested that the hexaploidy event occurred in two discrete steps (from diploid to tetraploid, then from tetraploid to hexaploid) with a reduction in gene number at each stage, resulting in one genome (least-fractionated) copy maintaining a higher percentage of genes than the other two . It has also been suggested that genome dominance has led to the inequity in the maintenance of gene copies after the WGD events, such that the least-fractionated genome of B. rapa exhibits higher levels of gene expression than its nuclear counterparts . However, the mechanisms for maintaining such dominance are unknown, although epigenetic phenomena have been suggested as possible determinants.
The current report describes the generation of a reference genome for the mesopolyploid B. oleracea along with transcriptome and methylome characterization. The genome was compared with its close relative B. rapa to provide insights into the observed difference in genome size, with B. oleracea being approximately 20% larger than B. rapa. In addition, the legacy of the mesopolyploidy event relative to maintenance of gene copy number, relative gene expression and cytosine methylation levels was assessed.
Assembly statistics for Brassica oleracea
Number of sequences a
Number of bases
Max length (bp)
Genetic anchoring of the genome assembly
Genetic anchoring of TO1000 assembly to B. oleracea pseudochromosomes (C1 to C9)
Number of linked scaffolds
Number of base pairs anchored
Percentage of base pairs anchored
Number of anchored scaffolds
Number of base pairs (%) in pseudochromosomes
Genome annotation of Brassica oleracea
A combination of homology-based and de novo gene prediction was used to annotate 59,225 gene models in the B. oleracea genome . Complete protein coding genes were predicted for 54,475 models. The gene models were annotated based on homology to proteins in public plant databases (Table S3 in Additional file 2). Orthologues for 94.6% (56,082) of all gene models were identified within the various UniProt databases, of which 52,774 (94%) were found in species from the Brassicaceae (Table S3 in Additional file 2). RNASeq data generated from four tissue types provided evidence of expression from 62.7% of the gene models (Table S4 in Additional file 2).
The number of B. oleracea gene models was 1.4-fold greater than for its close relative B. rapa (41,174). However, different annotation pipelines were utilized, and although repeat masking was employed, it is possible that some of the gene expansion may be due to the presence of uncharacterized transposon-related sequences . Reciprocal alignments between the annotated genes of the two Brassica diploids using BLAST found 32,061 orthologous pairs. A further 20,343 B. oleracea gene models showed similarity to an annotated B. rapa gene, suggesting expansion of gene families or tandem duplication, and 3,507 were homologous to B. rapa intergenic sequence, suggesting the presence of potentially un-annotated genes in the B. rapa genome . The remaining 3,314 genes could represent B. oleracea-specific genes or are possible artifacts of avaricious annotation.
Summary of repeat elements annotated in B. rapa and B. oleracea
Number of copies
Number of base pairs
Percentage genome (not N)
Number of copies
Number of base pairs
Percentage genome (not N)
Comparative genome organization of B. oleraceawith close relatives
Single base resolution of the B. oleraceamethylome
Sequence data from bisulfite-treated DNA extracted from TO1000 leaf tissue was analyzed to provide a picture of the cytosine DNA methylation landscape across the genome (Figure 3). From the almost 69 million trimmed sequence reads, 47.7% were uniquely mapped to the B. oleracea genome, providing an average read depth of 7.29 (per strand) at 93.5% of the cytosines in the assembled genome. Methylated cytosine bases were determined by applying a binomial test to data extracted from reads at each cytosine, using unmethylated lambda DNA as a control.
Cytosine methylation can be divided into three types based on sequence context, where mCG and mCHG (where H is A, C or T) possess symmetry and mCHH is asymmetric. The level of DNA methylation observed throughout the B. oleracea genome was 54.9% of all CGs, 9.4% of all CHGs and 2.4% of all CHHs. Varying reports of such levels have been found for plants ranging, for mCG, from 22% in A. thaliana[20, 21], to 51% in soybean , and 59% in rice . The B. oleracea genome was enriched for mCG when compared to levels more commonly observed for A. thaliana but demonstrated similar levels of cytosine methylation in the mCHG and mCHH contexts. The pattern observed in B. oleracea was in contrast to that noted recently for the paleopolyploid soybean, where higher levels of both mCG and mCHG were reported . The distribution of methylated cytosines across the genome mirrored the distribution of repetitive elements, with methylation enrichment observed in particular towards the centromeric regions (Figure 3). This biased distribution is a reflection of the higher levels of cytosine methylation found in all contexts at the individual repetitive element level (Table S8 in Additional file 2). The average level of mCG methylation for repeat elements ranged from 88% to 100%, with SINEs showing the highest levels. mCHG and mCHH are commonly associated with silencing of repetitive elements, and although showing more variation compared to mCG, they were consistently two to seven times the observed genome-wide levels across all repeat types (Table S8 in Additional file 2).
Correlation of cytosine methylation and expression levels
Impact of the Brassiceae-specific hexaploidy event
The most recent WGD event(s) in related B. rapa resulted in a hexaploid structure composed of three sub-genomes differentiated not only by extensive chromosomal restructuring but also varying levels of gene fractionation (or loss). This WGD preceded the divergence of B. oleracea from B. rapa and except for a limited number of additional gross chromosomal rearrangements both species have maintained a similar genome structure. In addition, comparable levels of biased gene fractionation have resulted in one sub-genome maintaining higher numbers of orthologous genes relative to the other two. It has been suggested that genome dominance post-polyploidy determines the fate of duplicated genes and gene expression has been positively correlated with maintenance of gene copies in other species and may play a similar role in Brassica species .
The B. oleracea species include a wide range of important vegetable crops, with diverse morphological variation and an eclectic mix of phytochemicals and secondary metabolites, many with health promoting properties. Of these morphotypes, TO1000, a kale-like plant, was selected as an excellent experimental model since it is rapid cycling, self-compatible, uncommon within B. oleracea, and has associated genomic and genetic resources, including a BAC-based physical map , a doubled haploid mapping population  and a mutagenized population . Thus, the reference genome sequence of TO1000 will provide an excellent tool for dissecting the molecular basis of the remarkable variation found within the species. In addition, B. oleracea is a close relative of the C genome species that hybridized with B. rapa to form the economically important oilseed crop B. napus; as such, the genome sequence will provide a foundation for further research in the allopolyploid relative.
The assembled nuclear genome sequence of TO1000 indicates that the increased genome size of B. oleracea compared to B. rapa results from multiplication of both repetitive and genic sequences. However, the two genomes evolved from a common ancestor and, as such, are remarkably similar with minimal large-scale chromosomal rearrangements having occurred since their divergence approximately 4 Mya . The two genomes share the Brassiceae-specific hexaploidization event and 77.8% of the annotated genes from B. rapa have a retained B. oleracea orthologue. Thus, even with on-going independent fractionation post-polyploidization and species divergence there is significant conservation of syntenic gene content between the A and C genomes, which is reflected in the composition of ancestral blocks along with their fractionation status in the two genomes (Figure 5). The maintenance of the syntenic block structure reflects the recent shared ancestry of the two genomes. The mechanism by which the non-syntenic genes have evolved is unclear, although they now show evidence of pseudogenization. It has been suggested that such changes could result from transposon activity causing transduplication , and there was a slight bias in prevalence of transposable elements flanking non-syntenic genes; however, 30% of syntenic genes also had transposable elements annotated within 500 bp of the coding region, which largely reflects the higher density of such structures in the B. oleracea genome (Table S9 in Additional file 2). It appears that the majority of the gene loss observed in the two Brassica diploids occurred prior to species divergence. Although the chromosomal rearrangements and independent repetitive element evolution may be sufficient to ensure inter-specific incompatibility, the majority of the species morphological characteristics presumably stem from allelic diversification and expansion of particular gene families.
Polyploidy is a prevalent evolutionary mechanism and it has been suggested that most angiosperms have experienced one or more polyploidy events in their history . Brassica species are excellent models for the study of this mode of genome evolution since at least four WGD events underlie their current genome architecture, three shared with the Brassicaceae experimental model A. thaliana and the most recent (estimated at approximately 23 Mya) specific to Brassiceae. The genome sequence of B. rapa allowed a detailed analysis of the three ancestral genomes that had fused to form the Brassiceae tribe, revealing an asymmetrical pattern of gene loss, suggesting a two-step polyploidy event with the last genome to join retaining the highest number of ancestral genes . The genome of B. oleracea mirrors this observation. It has been proposed that the genome dominance inferred to regulate this biased gene loss stems from imbalances in gene expression  and our analyses of transcriptome data from four B. oleracea tissue types supports this hypothesis, since higher expression of genes from the LF genome was observed for 43% of the fully retained triplets showing a significant sub-genome by tissue type interaction. However, study of the retained triplicated genes provides more compelling evidence for functional diversification of gene duplicates, with 83% of the triplicated genes showing independent expression patterns across the four tissue types. The fate of duplicate gene copies following WGD has been the subject of much debate and numerous models have been proposed to explain their retention or loss with discussion seemingly favoring dosage balance, neo-functionalization and sub-functionalization [29, 30]. Perhaps at one time controlled through dosage balance, which would require coordinated expression, the maintenance of duplicate genes post-hexaploidy in B. oleracea appears now to result from neo- or sub-functionalization.
The mechanisms that respond to WGD to provide meiotic stability post-polyploidy are unknown but the rapidity of change offered by epigenetic effects has suggested their potential role . Indeed, methylation re-patterning post-hybridization has been found in specific plant lineages, including Brassica species . Cytosine methylation has been shown to have an impact on gene expression and its role in controlling the activity of transposable elements could influence genome stability. In addition, cytosine methylation itself can be a source of further molecular variation through spontaneous deamination of methylated cytosines. Through single base resolution we have provided the first cytosine methylation map for the reference B. oleracea genome, which parallels findings in other eukaryotic genomes. Biased distribution reflecting the role of cytosine methylation in repeat element suppression was observed and patterns of methylation across gene structures were similar to those observed for multiple species . Although methylation across the three ancestral sub-genomes was imbalanced, with the least fractionated genome experiencing lower levels, this did not appear to be reflected at the level of the retained triplicate genes. Since the triplicate genes appear to have undergone a significant level of functional diversification, the influences of residual genome dominance from cytosine methylation seem minimal at the gene level. If cytosine methylation was an influencing factor in restoring balance in the nucleus of the ancestral polyploid B. oleracea, its effects are no longer apparent. It has been shown that methylation patterns can be conserved across orthologous genes within lineages  yet the patterns observed for the fully retained triplicated homologues in B. oleracea appeared largely independent. Thus, it is possible that cytosine methylation played a significant role early in the ancestry of B. oleracea to quickly establish relative expression differences but may now be involved in diversification of gene expression among the three sub-genomes.
We present a reference genome for the mesopolyploid B. oleracea, a resource for further study of polyploid evolution and a platform for genetic improvement of an array of Brassica vegetable crop types. Associated transcriptome data suggest the functional diversification of duplicate gene copies within the genome could be a source of the rich diversity within this species. In addition, the cytosine methylation landscape for this mesopolyploid provides insights into the role of this genome-wide mechanism in the evolution of the B. oleracea genome. Our findings show empirical links between gene fractionation, expression and epigenetic phenomena as major factors shaping the evolution of a polyploid genome.
Materials and methods
Library construction and sequencing
A homozygous doubled haploid line TO1000DH3 derived from B. oleracea cultivar TO1434 was chosen for sequencing. High quality nuclear DNA was extracted from young leaves using a megabase-sized DNA isolation protocol as described in . Briefly, approximately 40 g of fresh leaf tissue was homogenized in 200 ml buffer (HB; 0.01 M Trizma base, 0.08 M KCL, 0.01 M EDTA, 1 mM spermidine, 1 mM spermine, 0.5 M sucrose plus 0.15% β-mercaptoethanol, pH 9.4 to 9.5). The homogenate was filtered and the nuclei pelleted by centrifugation (1,800 g at 4°C for 20 minutes). The pellet was resuspensed (1 × HB plus 0.5% Triton-×100) and centrifuged three times. Finally, the nuclei were resuspended in 10 ml lysis buffer (100 mM TrisCl, 100 mM NaCl, 50 mM EDTA, 2% SDS). High molecular weight genomic DNA was extracted by traditional proteinase K (0.05 mg/ml; 65°C for 2 h) digestion followed by RNAase A treatment, two cycles of phenol/chloroform extraction and ethanol precipitation. Quantification of genomic DNA was performed using PicoGreen dsDNA kit (Molecular Probes, Life Technologies Inc., Burlington, ON, Canada). Genomic DNA (5 to 40 μg) was randomly sheared using one of: Covaris S2 ultrasonicator (Covaris Inc., Woburn, MA, US); Hydroshear (Genomic Solutions Inc., Ann Arbor, MI, US); or gas-driven nebulizers. For Illumina sequencing, four paired-end (PE) libraries (with median insert sizes of 273, 335, 418 and 532 bp; Table S1 in Additional file 2) and five short-span mate-paired (MP) libraries (from 2.5 to 8.5 kb) were constructed following the manufacturer’s instructions (TruSeq DNA sample preparation and MP library preparation kit v2 (Illumina, San Diego, CA, US), respectively). Libraries were size selected using the Pippin prep automated gel electrophoresis system (Sage Science, Beverly, MA, US), quantified using a BioAnalyzer (Agilent Technologies, Mississauga, ON, Canada) and KAPA library quantification kit for Illumina (KAPA Biosystems, Wilmington, MA, US), and sequenced from both ends (paired-end) for 100 cycles on an Illumina HiSeq 2000 instrument. For 454 pyrosequencing, two medium span MP libraries with median insert sizes of 8 and 17 kb (Table S1 in Additional file 2) were constructed following the method described in the GS FLX Titanium 20 kb span PE library preparation manual from Roche and sequenced using a Roche 454 FLX Titanium sequencer (454 Life Sciences, Branford, CT, US).
Initially, all Illumina and 454 reads were filtered for adapter contamination, PCR duplicates, ambiguous residues (N residues) and low quality regions. The initial backbone of the draft genome was assembled with Illumina reads using De Bruijn graph-based SOAPdenovo (version 1.05) assembler [8, 35], run with a k-mer parameter of 47 and each library ranked according to insert size from smallest to largest. The gaps within assembled scaffolds were filled with the short insert PE reads using GapCloser (version 1.12). The resulting assembly consisted of a total of 35,436 contigs and short scaffolds, with a sequence span of 488 Mb and an N50 size of 265 kb. BAC end sequences for TO1434 were downloaded from NCBI (LIBGSS_011756) and trimmed for quality, ambiguous bases and adapter sequences. Bambus  was used to overlay all the 454 MP information and the BAC end sequence data onto SOAPdenovo scaffolds to improve scaffold lengths as described in . In short, all 454 MP reads and BAC end sequence reads (Table S1 in Additional file 2) were aligned to the scaffolds using a genomic mapping and alignment program (GMAP) . The output from GMAP was used to create a Bambus-compatible GDE formatted contig file that indicated scaffold links. Redundant or multi-mapped mates, mates where only one read mapped, and those where both mates mapped to a single scaffold were considered invalid links. Each link was considered in Bambus in ascending order of their length, with scaffolding parameters including a redundancy level of 3 and link size error of 5%. Any potentially ambiguous scaffolds were resolved using the 'untangle' utility of Bambus. Bambus was able to order, orient and merge 2,623 of these pre-assembled SOAPdenovo scaffolds into 646 superscaffolds, resulting in a greatly improved assembly with an N50 size of 850 kb.
Construction of a high density genetic map and anchoring of the genome
A high density genetic map representing nine linkage groups was constructed using a mapping population of 94 doubled haploid lines (DH) derived previously from a cross between TO1000 and Early Big . A total of 2,299 polymorphic loci (SNPs, simple sequence repeats and insertion/deletion polymorphisms) identified using the RAD approach  were used to integrate assembled scaffolds with the genetic map. Collinearity between B. oleracea provisional pseudomolecules and A. thaliana and/or B. rapa was used to further assist with ordering and orientation of scaffolds for which there was paucity of adequate genetic recombination and markers. A total of 66 instances of false joins or insertions within Bambus superscaffolds were identified based on marker discontiguity and collinearity information. These scaffolds were split and the correct position of each of the fragments was determined based on marker and collinearity information. Final scaffolds were renamed as ‘Scaffold’ and numbered sequentially based on their length from longest to shortest. The order and orientation of scaffolds within each pseudomolecule were determined based on marker order within each scaffold, and marker contiguity pattern between adjoining scaffolds. Scaffolds with too few markers were ordered and oriented using collinearity information. The final version of the draft genome representing nine pseudochromosomes and 32,919 unanchored scaffolds was collated using a custom Perl script (Additional file 3), and the ordering and orientation information of scaffolds within each pseudochromosome was compiled in AGP files. The quality of the assembled genome was ascertained by performing several independent tests (Table S3 in Additional file 2; Figure S2 in Additional file 1).
Half-tetrad genetic mapping of centromeres
Centromere positions were mapped using half-tetrad analysis in a population of 49 plants derived from first-division restitution unreduced gametes of several B. napus × B. carinata hybrids . A total of 13,098 SNPs polymorphic between the B. napus and B. carinata parent genotypes for at least two progeny sets were selected from the recently released Illumina Brassica 60 K SNP chip and used to assess offspring heterozygosity at each C-genome locus. Increasing heterozygosity towards each centromere was used to predict centromere locations (see  for details).
The whole plant transcriptome of B. oleracea based on Illumina RNAseq data was characterized for four different tissue samples (leaf, root, flower, and pod; Table S4 in Additional file 2). Total RNA was isolated using the RNeasy plant mini kit (QIAGEN, Toronto, ON, Canada), including on-column DNase digestion, according to the manufacturer’s instructions. The integrity and quantity of total RNA was assessed using RNA 6000 Nano labchip on the BioAnalyzer (Agilent). Sequencing libraries were constructed following the standard TruSeq RNA sample preparation guide (Illumina), and PE sequencing was performed using the Illumina Hiseq 2000 platform. A total of 12.5 Gb raw RNAseq data were generated. Prior to read mapping, all reads were filtered for adapter contamination, ambiguous residues (N residues) and low quality regions. B. oleracea gene expression was assessed using a previously described TopHat and Cufflinks-based method . For both TopHat and Cufflink analysis, maximum and minimum intron lengths were set at 2,500 and 20 bp, respectively. Transcript abundance was measured as FPKM values.
For accurate annotation of gene models, an integrated computational approach based on two major genome annotation pipelines, Maker  and PASA , was adopted. Maker provides a simplified process for aligning expressed sequence tags (ESTs) and proteins to the genome, and integrates this external homology evidence with ab initio gene predictions to produce polished gene annotations with evidence-based quality statistics. Inputs for Maker included the repeat-masked B. oleracea genome assembly (masked against ‘te_proteins.fasta’ in the Maker package, which contains a generic list of common transposable elements (TEs)), PlantGDB ESTs from B. oleracea, B. rapa and B. napus, and Uniprot (SwissProt + TrEMBL) plant protein database. Ab initio gene predictions were made by Fgenesh  and Augustus . Maker gene structure annotations were further updated by PASA using evidence from Sanger ESTs and multiple de novo RNA-seq assemblies. Post-Maker processing included splitting potentially fused genes, extending genes, resolving internal rearrangements, trimming overlapping genes, and removing proteins of less than 50 amino acids and with no BLAST match to A. thaliana. The final annotation set contained a total of 59,225 gene models. Protein names are assigned using AHRD pipeline  with names extracted from BLASTP hits in TAIR v10, Swissprot and TrEMBL databases (queried in March 2012).
A TE database constituted from de novo analysis of B. napus (Chalhoub et al., in preparation) was merged with databases of TEs previously constructed from analysis of B. rapa and B. oleracea. The TEs were classified into major subclasses and superfamilies (Table 3) based on their structural features. Inside each superfamily, elements sharing more than 80% sequence identity over more than 80% of their length, and at least 80 bp, are considered as belonging to the same family . The merged TE database was used for comparative repeat masking  of the B. oleracea and B. rapa genomes.
Sequence homology was detected by BLASTP of the predicted proteins against the A. thaliana proteome. BLAST hits with E-value 1E-20 or better and within the top 40% from the best bit score were kept for further analysis. The chains of syntenic B. oleracea-A. thaliana gene pairs were computed by DAGChainer  using the default parameters. In the case of a B. oleracea gene participating in more than one syntenic chain due to duplication in the A. thaliana genome, the B. oleracea-A. thaliana pair in the weaker scoring chain was removed from the analysis. The syntelog table was generated by placing the syntenic chains onto the B. oleracea chromosomes.
Analysis of expression divergence
All calculations were completed using R . For the 2,242 triplets of (fully retained) B. oleracea genes the raw expression values (FPKM) from each of the four tissue types (T) and three sub-genomes (G) were transformed by adding 1 and taking the natural logarithm. All subsequent calculations were done with the transformed values. The means for each sub-genome were taken over the four tissue types, and the genome with the highest expression was recorded. With no replication interaction (G×T) cannot be separated from random error, so the gene triplets were sorted by the P-value for an extra interaction term (Tukey)  that is the product of the main effects terms in the ANOVA model for the expression data for a gene triplet. The cumulative frequency of each genome with the highest expression was plotted against the Tukey P-values.
Hierarchical clustering of gene expression was carried out using 1 minus the Pearson sample correlation coefficient as the distance measure, which requires the variances of the expression values for each gene to be non-zero. This clusters correlated sets of expression values regardless of the absolute magnitude of the expression levels. The average linkage method was used to form the clusters.
Protein family classification
An all-against-all BLASTP (E-value cutoff of 1E-10)  similarity search was performed among predicted protein sequences of four completely sequenced Brassicaceae species, A. thaliana, A. lyrata, B. rapa and B. oleracea. The pairwise similarities generated by this analysis were parsed and stored in a matrix that served as input for clustering by the TRIBE-MCL approach . The groups of homologous proteins were clustered using Markov cluster (MCL) algorithm at an inflation rate (I) of 3.0 and other default parameters. A χ2 test was performed to determine if significant expansion of gene families among the Brassicaceae species occurred. χ2 probability values were adjusted by Bonferroni correction and the post-hoc test was conducted only if the P-value after Bonferroni correction was <0.05.
Single base resolution methylome sequencing
Genomic DNA was isolated from B. oleracea TO1000 nuclei prepared according to the method found in Zhang et al.  from leaf tissue as above. TO1000 genomic DNA (5 μg), along with 0.5% w/w of unmethylated lambda DNA (Promega, Madison, WI, US) included to assess bisulfite conversion efficiency, was sheared using the Diagenode Bioruptor Plus to produce an average fragment size of 300 bp. Library construction was performed using the Illumina Paired-End Sample Prep Kit and QIAGEN PCR Purification kits according to the Illumina Whole-genome Bisulfite Sequencing for Methylation Analysis protocol. Methylated adapters from the Illumina TruSeq DNA Sample Prep Kit v2 were substituted for the unmethylated adapters within the Illumina Paired-End Sample Prep Kit. Ligation products equivalent to 450 bp were purified using the BluePippin DNA size selection system (Sage Science). Bisulfite conversion of this DNA was performed using the QIAGEN EpiTect Bisulfite kit. The libraries were quantified using a BioAnalyzer (Agilent) and KAPA library quantification kit for Illumina (KAPA Biosystems), and sequenced from both ends (paired-end) for 100 cycles on either an Illumina HiSeq 2000 or Illumina HiScanSQ.
BSMAP (v2.74) and an associated python script (methratio.py)  were used to analyze cytosine methylation ratios. All reads were trimmed for quality and adapter sequence, potential PCR duplicates were removed and only those reads mapping uniquely to the genome were retained. Sequencing of cytosines in lambda DNA included in each library provided an error rate for non-conversion; this rate was used in a custom Perl script to identify methylated cytosines using the binomial probability distribution at a false discovery rate of <5% (Additional file 4).
Perl was used for text manipulation to assign pseudochromosome coordinates to gene features, 5′ promoter, 5′ UTR, gene, exon, intron, 3′ UTR and 3′ regions, and incorporate these annotations in a custom .gff file. A combination of Perl and R was used to perform statistical analyses on these data. Graphics summarizing the observed statistical relationships were developed using ggplot2, RGL and corrgram in R . All relevant Perl and R scripts are open source and freely available upon request.
All the raw data used in the assembly of the genome, the RNASeq data used in the expression analysis, and the bisulfite sequence data have been deposited at DDBJ/EMBL/GenBank short read archives under accession numbers SRP040796, SRS586190, SRS586272, SRS586273. The whole genome shotgun project has been deposited at DDBJ/EMBL/GenBank under the accession JJMF00000000. The version described in this paper is version JJMF01000000.
bacterial artificial chromosome
expressed sequence tag
fragments per kilobase of exon per million fragments mapped
million years ago
restriction site-associated DNA
short interspersed nuclear element
single nucleotide polymorphism
whole genome duplication.
We thank Jacek Nowak and Faouzi Bekkaoui for project management support. Work in Canada was supported by the Canadian Canola Genome Sequencing Initiative (http://canseq.ca), Genome Alberta, industry partners and the Agricultural Flexibility Fund (Agri-Flex, AAFC). Work at the J. Craig Venter Institute, the University of Missouri and the University of Georgia was supported by grant number 0638536 from the US National Science Foundation. Work in Australia was supported by the Australian Research Council (projects LP0882095, LP0883462 and LP110100200). Support is also acknowledged from the Australian Genome Research Facility (AGRF) and the Queensland Cyber Infrastructure Foundation (QCIF). Work in the UK was supported by the UK Biotechnology and Biological Sciences Research Council (BBSRC BB/E017363/1) and the DEFRA grant IF0157 as part of the Vegetable Genetic Improvement Network (VeGIN).
- Wang X, Wang H, Wang J, Sun R, Wu J, Liu S, Bai Y, Mun JH, Bancroft I, Cheng F, Huang S, Li X, Hua W, Wang J, Wang X, Freeling M, Pires JC, Paterson AH, Chalhoub B, Wang B, Hayward A, Sharpe AG, Park BS, Weisshaar B, Liu B, Li B, Liu B, Tong C, Song C, Duran C, et al: The genome of the mesopolyploid crop species Brassica rapa. Nat Gen. 2011, 43: 1035-1039.View ArticleGoogle Scholar
- Navabi ZK, Huebert T, Sharpe AG, O’Neill CM, Bancroft I, Parkin IA: Conserved microstructure of the Brassica B Genome of Brassica nigra in relation to homologous regions of Arabidopsis thaliana, B. rapa and B. oleracea. BMC Genomics. 2013, 14: 250-PubMedPubMed CentralView ArticleGoogle Scholar
- Leitch AR, Leitch IJ: Genomic plasticity and the diversity of polyploid plants. Science. 2008, 320: 481-483.PubMedView ArticleGoogle Scholar
- Bowers JE, Chapman BA, Rong J, Paterson AH: Unravelling angiosperm genome evolution by phylogenetic analysis of chromosomal duplication events. Nature. 2003, 422: 433-438.PubMedView ArticleGoogle Scholar
- Beilstein MA, Nagalingum NS, Clements MD, Manchester SR, Mathews S: Dated molecular phylogenies indicate a Miocene origin for Arabidopsis thaliana. Proc Natl Acad Sci U S A. 2010, 107: 18724-18728.PubMedPubMed CentralView ArticleGoogle Scholar
- Tang H, Woodhouse MR, Cheng F, Schnable JC, Pedersen BS, Conant G, Wang X, Freeling M, Pires JC: Altered patterns of fractionation and exon deletions in Brassica rapa support a two-step model of paleohexaploidy. Genetics. 2012, 190: 1563-1574.PubMedPubMed CentralView ArticleGoogle Scholar
- Cheng F, Wu J, Fang L, Sun S, Liu B, Lin K, Bonnema G, Wang X: Biased gene fractionation and dominant gene expression among the subgenomes of Brassica rapa. PLoS One. 2012, 7: e36442-PubMedPubMed CentralView ArticleGoogle Scholar
- Luo R, Liu B, Xie Y, Li Z, Huang W, Yuan J, He G, Chen Y, Pan Q, Liu Y, Tang J, Wu G, Zhang H, Shi Y, Liu Y, Yu C, Wang B, Lu Y, Han C, Cheung DW, Yiu SM, Peng S, Xiaoqian Z, Liu G, Liao X, Li Y, Yang H, Wang J, Lam TW, Wang J: SOAPdenovo2: an empirically improved memory-efficient short-read de novo assembler. Gigascience. 2012, 1: 18-PubMedPubMed CentralView ArticleGoogle Scholar
- Pop M, Kosack DS, Salzberg SL: Hierarchical scaffolding with Bambus. Genome Res. 2004, 14: 149-159.PubMedPubMed CentralView ArticleGoogle Scholar
- Parra G, Bradnam K, Korf I: CEGMA: a pipeline to accurately annotate core genes in eukaryotic genomes. Bioinformatics (Oxford, England). 2007, 23: 1061-1067.View ArticleGoogle Scholar
- Baird NA, Etter PD, Atwood TS, Currey MC, Shiver AL, Lewis ZA, Selker EU, Cresko WA, Johnson EA: Rapid SNP discovery and genetic mapping using sequenced RAD markers. PLoS One. 2008, 3: e3376-PubMedPubMed CentralView ArticleGoogle Scholar
- Poland JA, Brown PJ, Sorrells ME, Jannink JL: Development of high-density genetic maps for barley and wheat using a novel two-enzyme genotyping-by-sequencing approach. PLoS One. 2012, 7: e32253-PubMedPubMed CentralView ArticleGoogle Scholar
- Cantarel BL, Korf I, Robb SM, Parra G, Ross E, Moore B, Holt C, Sanchez Alvarado A, Yandell M: MAKER: an easy-to-use annotation pipeline designed for emerging model organism genomes. Genome Res. 2008, 18: 188-196.PubMedPubMed CentralView ArticleGoogle Scholar
- Town CD, Cheung F, Maiti R, Crabtree J, Haas BJ, Wortman JR, Hine EE, Althoff R, Arbogast TS, Tallon LJ, Vigouroux M, Trick M, Bancroft I: Comparative genomics of Brassica oleracea and Arabidopsis thaliana reveal gene loss, fragmentation, and dispersal after polyploidy. Plant Cell. 2006, 18: 1348-1359.PubMedPubMed CentralView ArticleGoogle Scholar
- Lespinet O, Wolf YI, Koonin EV, Aravind L: The role of lineage-specific gene family expansion in the evolution of eukaryotes. Genome Res. 2002, 12: 1048-1059.PubMedPubMed CentralView ArticleGoogle Scholar
- Kawabe A, Hansson B, Hagenblad J, Forrest A, Charlesworth D: Centromere locations and associated chromosome rearrangements in Arabidopsis lyrata and A. thaliana. Genetics. 2006, 173: 1613-1619.PubMedPubMed CentralView ArticleGoogle Scholar
- Parkin IA, Gulden SM, Sharpe AG, Lukens L, Trick M, Osborn TC, Lydiate DJ: Segmental structure of the Brassica napus genome based on comparative analysis with Arabidopsis thaliana. Genetics. 2005, 171: 765-781.PubMedPubMed CentralView ArticleGoogle Scholar
- Schranz ME, Lysak MA, Mitchell-Olds T: The ABC’s of comparative genomics in the Brassicaceae: building blocks of crucifer genomes. Trends Plant Sci. 2006, 11: 535-542.PubMedView ArticleGoogle Scholar
- Cheng F, Mandakova T, Wu J, Xie Q, Lysak MA, Wang X: Deciphering the diploid ancestral genome of the mesohexaploid Brassica rapa. Plant Cell. 2013, 25: 1541-1554.PubMedPubMed CentralView ArticleGoogle Scholar
- Feng S, Cokus SJ, Zhang X, Chen PY, Bostick M, Goll MG, Hetzel J, Jain J, Strauss SH, Halpern ME, Ukomadu C, Sadler KC, Pradhan S, Pellegrini M, Jacobsen SE: Conservation and divergence of methylation patterning in plants and animals. Proc Natl Acad Sci U S A. 2010, 107: 8689-8694.PubMedPubMed CentralView ArticleGoogle Scholar
- Cokus SJ, Feng S, Zhang X, Chen Z, Merriman B, Haudenschild CD, Pradhan S, Nelson SF, Pellegrini M, Jacobsen SE: Shotgun bisulphite sequencing of the Arabidopsis genome reveals DNA methylation patterning. Nature. 2008, 452: 215-219.PubMedPubMed CentralView ArticleGoogle Scholar
- Schmitz RJ, He Y, Valdes-Lopez O, Khan SM, Joshi T, Urich MA, Nery JR, Diers B, Xu D, Stacey G, Ecker JR: Epigenome-wide inheritance of cytosine methylation variants in a recombinant inbred population. Genome Res. 2013, 23: 1663-1674.PubMedPubMed CentralView ArticleGoogle Scholar
- Schnable JC, Springer NM, Freeling M: Differentiation of the maize subgenomes by genome dominance and both ancient and ongoing gene loss. Proc Natl Acad Sci U S A. 2011, 108: 4069-4074.PubMedPubMed CentralView ArticleGoogle Scholar
- Wang X, Torres MJ, Pierce G, Lemke C, Nelson LK, Yuksel B, Bowers JE, Marler B, Xiao Y, Lin L, Epps E, Sarazen H, Rogers C, Karunakaran S, Ingles J, Giattina E, Mun JH, Seol YJ, Park BS, Amasino RM, Quiros CF, Osborn TC, Pires JC, Town C, Paterson AH: A physical map of Brassica oleracea shows complexity of chromosomal changes following recursive paleopolyploidizations. BMC Genomics. 2011, 12: 470-PubMedPubMed CentralView ArticleGoogle Scholar
- Iniguez-Luy FL, Lukens L, Farnham MW, Amasino RM, Osborn TC: Development of public immortal mapping populations, molecular markers and linkage maps for rapid cycling Brassica rapa and B. oleracea. Theor Appl Genet. 2009, 120: 31-43.PubMedView ArticleGoogle Scholar
- Himelblau E, Gilchrist EJ, Buono K, Bizzell C, Mentzer L, Vogelzang R, Osborn T, Amasino RM, Parkin IA, Haughn GW: Forward and reverse genetics of rapid-cycling Brassica oleracea. Theor Appl Genet. 2009, 118: 953-961.PubMedView ArticleGoogle Scholar
- Cheung F, Trick M, Drou N, Lim YP, Park JY, Kwon SJ, Kim JA, Scott R, Pires JC, Paterson AH, Town C, Bancroft I: Comparative analysis between homoeologous genome segments of Brassica napus and its progenitor species reveals extensive sequence-level divergence. Plant Cell. 2009, 21: 1912-1928.PubMedPubMed CentralView ArticleGoogle Scholar
- Jiao Y, Wickett NJ, Ayyampalayam S, Chanderbali AS, Landherr L, Ralph PE, Tomsho LP, Hu Y, Liang H, Soltis PS, Soltis DE, Clifton SW, Schlarbaum SE, Schuster SC, Ma H, Leebens-Mack J, de Pamphilis CW: Ancestral polyploidy in seed plants and angiosperms. Nature. 2011, 473: 97-100.PubMedView ArticleGoogle Scholar
- Innan H, Kondrashov F: The evolution of gene duplications: classifying and distinguishing between models. Nat Rev Genet. 2010, 11: 97-108.PubMedView ArticleGoogle Scholar
- Conant GC, Wolfe KH: Turning a hobby into a job: how duplicated genes find new functions. Nat Rev Genet. 2008, 9: 938-950.PubMedView ArticleGoogle Scholar
- Madlung A, Wendel JF: Genetic and epigenetic aspects of polyploid evolution in plants. Cytogenet Genome Res. 2013, 140: 270-285.PubMedView ArticleGoogle Scholar
- Lukens LN, Pires JC, Leon E, Vogelzang R, Oslach L, Osborn T: Patterns of sequence loss and cytosine methylation within a population of newly resynthesized Brassica napus allopolyploids. Plant Physiol. 2006, 140: 336-348.PubMedPubMed CentralView ArticleGoogle Scholar
- Takuno S, Gaut BS: Gene body methylation is conserved between plant orthologs and is of evolutionary consequence. Proc Natl Acad Sci U S A. 2013, 110: 1797-1802.PubMedPubMed CentralView ArticleGoogle Scholar
- Kagale S, Koh C, Nixon J, Bollina V, Clarke WE, Tuteja R, Spillane C, Robinson SJ, Links MG, Clarke C, Higgins E, Huebert T, Sharpe AG, Parkin IAP: The emerging biofuel crop Camelina sativa retains a highly undifferentiated hexaploid genome structure. Nat Comm. 2014, 5: 3706-View ArticleGoogle Scholar
- Short Oligonucleotide Analysis Package. [http://soap.genomics.org.cn/soapdenovo.html]
- Wu TD, Watanabe CK: GMAP: a genomic mapping and alignment program for mRNA and EST sequences. Bioinformatics. 2005, 21: 1859-1875.PubMedView ArticleGoogle Scholar
- Mason AS, Nelson MN, Castello M-C, Yan G, Cowling WA: Genotypic effects on the frequency of homoeologous and homologous recombination in Brassica napus × B. carinata hybrids. Theor Appl Genetics. 2011, 122: 543-553.View ArticleGoogle Scholar
- Park TH, Kim JB, Hutten RC, van Eck HJ, Jacobsen E, Visser RG: Genetic positioning of centromeres using half-tetrad analysis in a 4×-2× cross population of potato. Genetics. 2007, 176: 85-94.PubMedPubMed CentralView ArticleGoogle Scholar
- Trapnell C, Roberts A, Goff L, Pertea G, Kim D, Kelley DR, Pimentel H, Salzberg SL, Rinn JL, Pachter L: Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks. Nat Protoc. 2012, 7: 562-578.PubMedPubMed CentralView ArticleGoogle Scholar
- Haas BJ, Delcher AL, Mount SM, Wortman JR, Smith RK, Hannick LI, Maiti R, Ronning CM, Rusch DB, Town CD, Salzberg SL, White O: Improving the Arabidopsis genome annotation using maximal transcript alignment assemblies. Nucleic Acids Res. 2003, 31: 5654-5666.PubMedPubMed CentralView ArticleGoogle Scholar
- Stanke M, Tzvetkova A, Morgenstern B: AUGUSTUS at EGASP: using EST, protein and genomic alignments for improved gene prediction in the human genome. Genome Biol. 2006, 7: S11 11-18View ArticleGoogle Scholar
- Automated Assignment of Human Readable Descriptions (AHRD). [https://github.com/groupschoof/AHRD/]
- Zhao M, Du J, Lin F, Tong C, Yu J, Huang S, Wang X, Liu S, Ma J: Shifts in the evolutionary rate and intensity of purifying selection between two Brassica genomes revealed by analyses of orthologous transposons and relics of a whole genome triplication. Plant J. 2013, 76: 211-222.PubMedGoogle Scholar
- Wicker T, Sabot F, Hua-Van A, Bennetzen JL, Capy P, Chalhoub B, Flavell A, Leroy P, Morgante M, Panaud O, Paux E, SanMiguel P, Schulman AH: A unified classification system for eukaryotic transposable elements. Nat Rev Genet. 2007, 8: 973-982.PubMedView ArticleGoogle Scholar
- RepeatMasker. [http://www.repeatmasker.org]
- Haas BJ, Delcher AL, Wortman JR, Salzberg SL: DAGchainer: a tool for mining segmental genome duplications and synteny. Bioinformatics. 2004, 20: 3643-3646.PubMedView ArticleGoogle Scholar
- R Development Core Team: R: A Language and Environment for Statistical Computing. 2010, Vienna, Austria: R Foundation for Statistical Computing, [http://www.R-project.org/]Google Scholar
- Tukey JW: One degree of freedom for non-additivity. Biometrics. 1949, 5: 232-242.View ArticleGoogle Scholar
- Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ: Basic local alignment search tool. J Mole Biol. 1990, 215: 403-410.View ArticleGoogle Scholar
- Enright AJ, Van Dongen S, Ouzounis CA: An efficient algorithm for large-scale detection of protein families. Nucleic Acids Res. 2002, 30: 1575-1584.PubMedPubMed CentralView ArticleGoogle Scholar
- Zhang H-B, Zhao X, Ding X, Paterson AH, Wing RA: Preparation of megabase-size DNA from plant nuclei. Plant J. 1995, 7: 175-184.View ArticleGoogle Scholar
- Xi Y, Li W: BSMAP: whole genome bisulfite sequence MAPping program. BMC Bioinformatics. 2009, 10: 232-PubMedPubMed CentralView ArticleGoogle Scholar
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.