Plant materials sequenced in this study were provided by the Chinese Cabbage Research Group of the Institute of Vegetables and Flowers, Chinese Academy of Agricultural Sciences (CAAS-IVF), from the Dutch Crop Genetic Resources Center (CGN) in Wageningen UR and from Bejo and Rijk Zwaan breeding companies, comprising a total of 60 Brassica accessions (Supplementary Table S1). For the accessions that were inbred or DH lines, DNA was isolated from 50 to 100 two-week-old seedlings (hypocotyls and cotyledons) or from leaves from a single plant. For GeneBank accessions, DNA was isolated from single plants, as accessions are generally heterogeneous. In the paper, we use the term accession for both the single plants representing GeneBank accessions and the DH or inbred lines from breeding companies. Total DNA was extracted from seedlings (hypocotyls and cotyledons) or fresh leaves of the 60 Brassica accessions by the conventional cetyltrimethylammonium bromide (CTAB) method (Chen and Ronald, 1999) and stored at -80°C. Chloroplast genome sequences of A. thaliana (NC_000932), A. lyrata (LN877383), B. rapa (DQ231548), B. juncea (KT581449), B. oleracea (KR233156), B. nigra (KT878383) and B. napus (KP161617, GQ861354, KM454973 and KJ872515) were obtained from the original publication.

Whole genome sequences were generated using Illumina HiSeq2000 and HiSeq3000. In brief, short-insert (400-bp) libraries were constructed using the TruSeq^TM DNA sample preparation kit (Illumina, USA) for Illumina HiSeq sequencing and following the manufacturer’s protocol (Illumina). DNA from the different samples was indexed by tags and pooled together in one lane of Illumina’s Genome Analyzer for sequencing. Raw reads were first filtered to obtain the high-quality clean data by removing adaptor sequences and low-quality reads with Q-value ≤20.

The chloroplast genomes were assembled using a highly automated five stage pipeline (manuscript in preparation). Briefly this pipeline consists of the following stages: (1) Extraction of a k-mer frequency table from raw data and detection of two linked peaks at frequency N and 2^∗N, representing single copy and IR regions of the chloroplast genome, respectively, and extraction of underlying k-mers. (2) Extraction of reads containing these k-mers and initial assembly using a parameter spread. (3) Refinement of initial assembly by iterative selection of additional reads from raw data based on previous assembly and subsequent reassembly. (4) Local selection and reassembly of additional reads near scaffold ends to alleviate problems caused by variable representation of the chloroplast genome in the sequencing library. (5) Gap filling by local selection and de novo assembly of these gap-filling reads and insertion of gap-filling contigs into the scaffolds. This pipeline takes into account the circular and quadripartite [(LSC, IR (twice) and SSC] nature of the chloroplast genome, and may, if data supports this, produce finished circular chloroplast genome assemblies. The assembled chloroplast genomes were deposited to BRAD¹.

The chloroplast genes were annotated using an online DOGMA tool (Wyman et al., 2004), using default parameters to predict protein-coding genes, transfer RNA (tRNA) genes, and ribosome RNA (rRNA) genes. Start and stop codons of protein-coding genes were searched and determined by BLASTX against the NCBI protein database, with B. rapa as a guide. Genome maps were drawn with OGDraw (Lohse et al., 2007). Multiple sequence alignments were performed in MAFFT version 5 (Katoh et al., 2005) under standard parameters, and adjusted manually where necessary. Full alignments with annotations were visualized using the VISTA viewer (Frazer et al., 2004). Pairwise assembly alignments were made between the reference genome and the de novo assembled genomes using MUMmer (Kurtz et al., 2004).

Sample sequences were mapped against the B. rapa chloroplast genome (Chiifu). The mapping procedure applied in this study is summarized as follows: Burrows–Wheeler Aligner (BWA 0.6.2) software was used to map paired-end reads in this study (Li and Durbin, 2010), and reads with low-quality alignments were then filtered out from mapping using SAMtools (Li et al., 2009) with the default parameters. The final set of mapping files (BAM) was used to perform a multisample variant calling by using BCFtools. The calling parameters were adjusted to obtain both SNVs and InDels.

To evaluate consistency of phylogenetic trees produced from regions with different molecular evolutionary rates (Palmer, 1985), we extracted three subsets (LSC, SSC, and IRs) from the complete chloroplast data set, and combined these to produce three types of trees.

The program JModeltest 2 was used to find the optimal substitution model for each subset (Posada and Crandall, 1998; Darriba et al., 2012), using both the Bayesian information criterion (Felsenstein, 1988) and the Akaike information criterion (Akaike, 1974). Models for three datasets were explored: LSC+SSC dataset (GTR+G model), IR dataset (GTR+I model) and whole chloroplast genome (GTR+I+G model). Assembled chloroplast sequences were aligned with the help of MAFFT software (Katoh et al., 2005). The maximum likelihood (ML) tree was computed in PhyML (Guindon et al., 2010) using settings of the best-fitting model tested by JModeltest program (Darriba et al., 2012). The method used was a maximum-likelihood iterative model and a bootstrap of 1,000 repetitions was used to assess the reliability of the phylogeny reconstructed.

Divergence times of Brassica species were estimated using the BEAST program with the Bayesian method (Drummond et al., 2012). Uncorrelated relaxed lognormal clock was used in our analysis that allows different rates to be optimized independently on each branch of the tree. The Yule process described the likelihood of speciation for this analysis (Yule, 1925), and branching rates were determined under an optimal GTR+I+G nucleotide model. Two independent Markov chain Monte Carlo runs in BEAST were conducted for 500,000,000 generations sampling every 50,000 generations for the chloroplast genome analysis. To calibrate the phylogenetic tree of the chloroplast genomes of the accessions belonging to the six Brassica species, estimations of dates of speciation events calculated by valid fossil evidence were used (Hohmann et al., 2015). In our experiment, Arabidopsis was the outgroup, and the root age was constrained by a normal distribution with a mean of 23.35 Mya. In addition, the divergence age of A. thaliana and A. lyrata was set at 5.97 Mya. We then used Tracer v. 1.6 to explore the output of BEAST (Rambaut et al., 2014). TreeAnnotator v. 1.8 (part of the BEAST package) was used to produce maximum clade credibility trees from the postburn-in trees and to determine the 95% posterior density of ages for all nodes in the tree.

Sixty sequenced chloroplast genomes of Brassica genotypes together with seven published chloroplast genomes were used in the present study, including chloroplast genomes of 30 B. rapa accessions representing nine different morphotypes, of 17 B. oleracea accessions representing seven different morphotypes, of eight B. napus accessions representing oilseed and swede morphotypes, of eight B. juncea accessions representing five morphotypes, of one B. nigra and of one B. carinata; as outgroup one chloroplast genome each of A. thaliana and A. lyrata was included (Supplementary Table S1). The raw reads were mapped onto the B. rapa subsp. pekinensis (Chiifu, DQ231548.1) chloroplast genome as described in the Section “Materials and Methods.” Mapped reads of each genotype for the sequences of 60 whole chloroplast genomes are listed in Supplementary Table S1. The chloroplast genome coverage of each sample was very high (Supplementary Table S1) because of the numerous copies of chloroplasts in leaf cells (Zhang et al., 2011).

In this study, we used k-mer frequency histograms to identify and extract the chloroplast reads from total DNA sequencing data (Supplementary Figure S1) and assembled chloroplast genomes for these accessions using a highly integrated and automated pipeline. Dot plot analyses using MUMmer software was performed to compare de novo assembled genomes with the reference genome (Figure 1). Assembled genomes of all six Brassica species are collinear to previously published chloroplast genomes of Brassica (Wu et al., 2012; Hao et al., 2014), as no rearrangements were identified. Their genome sizes are very similar, ranging from 152860 to 153700 bp (Supplementary Table S1).

A total of 2646 SNVs and 464 insertions and deletions (InDels) were detected by Burrows–Wheeler Aligner software (BWA) and SAMtools, summing up to a total of 3110 high-quality variant positions (Supplementary Tables S2 and S3). Manual curation and the high coverage used in this study allowed the detection of discrepancies with the reference genome that can most probably be attributed to Sanger sequencing errors in the original sequence. For example, position 116457 (A) of the B. rapa subsp. pekinensis (Chiifu) chloroplast genome that we re-sequenced was different from the reference (G, DQ231548.1). The number of changes observed for all the four nucleotides among the 2646 SNVs were 52.4% A or T, 47.6% G or C. In total, 56.3% SNVs were transversions and 43.7% SNVs were transitions. Of the transitions, 48 and 52% were A/G and T/C transitions, respectively. While A/C, A/T, C/G, and G/T transversions constituted 34.09, 16.31, 8.26, and, 41.34% of the cases, respectively. Also 25 triallelic positions were observed in these Brassica species.

The average SNV density in the whole chloroplast genome was 17.2 SNVs per kb, while the average SNV density for each region (LSC, SSC, and IRs) was 21.9, 33.0 and 4.5 SNVs per kb, respectively. The IR regions showed a much lower SNV density in comparison with the chloroplast average. The average SNV density for exons (total size 87356 bp), introns (total size 12574 bp) and intergenic regions (total size 53552 bp) was 12.5, 14.0 and 25.7 SNVs per kb, respectively. About 52.0% of variant positions (1376) were found in intergenic regions. Furthermore, some genes such as matK, rpoC2, rpoB, ycf2.1, ndhF, and ycf1.2, displayed a remarkably high number of SNVs in their exons (46, 71, 39, 40, 49, and 237, respectively; Supplementary Table S4). A 500 bp window size was used to calculate and draw the SNVs density distribution along the chloroplast genome, to identify regions with increased SNV densities that we consider hotspots of variation of chloroplast genomes (Figure 2). In addition, circular maps were drawn to show the hotspots of variation of chloroplast genomes in the three diploid Brassica species separately (Supplementary Figure S2). In B. rapa, based on 30 chloroplast genomes, the average SNV density in the chloroplast genome was 2.2 SNVs per kb, while the average SNV density for each region (LSC, SSC, and IRs) was 3.0, 3.9, and 0.5 SNVs per kb, respectively. In B. oleracea, based on 17 chloroplast genomes and the reference B. oleracea chloroplast genome, the average SNV density in the chloroplast genome was much lower with 0.1 SNVs per kb, while the average SNV density for each region (LSC, SSC, and IRs) was 0.2, 0.2, and 0 SNVs per kb, respectively. In B. nigra, based on two chloroplast genomes only the average SNV density was 0.5 SNVs per kb, while the average SNV density for each region (LSC, SSC, and IRs) was 0.6, 1.3, and 0.1 SNVs per kb, respectively. In total, 343 SNVs and 125 InDels were identified in the chloroplast genomes of the 30 B. rapa accessions. Interestingly, the chloroplast genomes of four B. rapa broccoletto accessions contributed many more SNVs and InDels than chloroplast genomes of other B. rapa accessions (Supplementary Table S1).

The average InDel density in the chloroplast genome, when chloroplast genomes of all six Brassica species were included was 3.0 InDels per kb, while the average InDel density for each region (LSC, SSC, and IRs) was 4.6, 3.3, and 0.4 SNVs per kb, respectively. The average InDel density for exons (total size 87356 bp), introns (total size 12574 bp) and intergenic regions (total size 53552 bp) was 0.1, 4.9, and 7.3 InDel per kb, respectively. In B. rapa, the average InDel density in the chloroplast genome was 0.8 InDel per kb, while the average InDel density for each region (LSC, SSC, and IRs) was 1.2, 1.1, and 0.1 InDel per kb, respectively. In B. oleracea, the average InDel density in the chloroplast genome was 0.04 InDel per kb, while the average InDel density for each region (LSC, SSC, and IRs) was 0.06, 0.06, and 0 InDel per kb, respectively. In B. nigra, the average InDel density in the chloroplast genome was 0.16 InDel per kb, while the average InDel density for each region (LSC, SSC, and IRs) was 0.28, 0.06, and 0.02 InDel per kb, respectively. Analysis of the distribution of the 464 InDels revealed that 392 (84.5%) were situated in intergenic regions whereas only 72 (15.5%) were located in genes (Supplementary Table S4; Supplementary Figure S3). Sixty one out of the 72 InDels were in introns whereas the remaining 11 InDels affected the coding regions of trnK-UUU, trnQ-UUG, rpoC2, trnS-UGA, and accD genes. A total of 464 manually curated high-quality InDels, ranging from 1 to 38 bp (including 41 triallelic positions), were detected. Only 6 out of 464 InDels were longer than 10 bp. Furthermore, 237 InDels were single-base InDels.

As the single copy regions (LSC and SSC) and IRs regions have very different molecular evolutionary rates (Palmer, 1985), three subsets of the data were produced, containing either the combination of LSC and SSC, only the IR or the whole chloroplast genome sequence.

The phylogenetic relationships derived from the LSC+SSC dataset under GTR+G model, IR dataset under GTR+I model and whole chloroplast genome under GTR+I+G model are shown in Figure 3. The topologies estimated from ML analyses of different data subsets were similar to each other, with four clades clearly distinguished. The only differences were the poorly supported branches within lowly divergent clades (Figure 3). Overall, phylogenetic analyses of the three different datasets did not generate strongly supported topological conflict. The phylogenetic resolution of the LSC+SSC subset (Figure 3a) and the whole chloroplast genome sequences subset (Figure 3c) was higher than that of the IR subset (Figure 3b). Clade I contained the chloroplast genomes of mustards (B. juncea), of most B. rapa accessions and of two B. napus accessions. Clade I based on the LSC+SSC subset was divided into five subclades (subclade a–e), while divided in three subclades (subclade a–c) when based on the whole chloroplast genome. These subclade groupings do not fully agree with each other. The subclade ‘a’ in the LSC+SSC analysis, containing chloroplast genomes of all mustard accessions, and of several B. rapa accessions (Japanese leafy types and turnips, one Italian broccoletto, and a European turnip), generally corresponds to the subclade ‘a’ of the whole chloroplast analysis, except for the chloroplast genome of one European turnip (5). The subclade b of the whole chloroplast analysis, containing chloroplast genomes of several B. rapa morphotypes, corresponds to subclade b, c, and d of the LSC+SSC based analysis. Subclade ‘e’ with chloroplast genomes of sarsons (B. rapa) and B. napus is identical for the two analyses. The fact that the results of different datasets were almost identical, in addition to their corresponding internal bootstrap support, shows the reliability of the phylogenies. Moreover, ML analyses of the LSC+SSC chloroplast sequences demonstrated increased resolution in the tree topology with higher support values (Figure 3a) compared to the other two datasets. Thus, we mainly focus on results from the LSC+SSC region of chloroplast genome analyses.

We generated a ML phylogeny from the 67 single copy LSC+SSC regions of chloroplast genomes including 60 de novo assembled sequences and seven previously published chloroplast genomes. The phylogenetic tree was rooted with A. thaliana and A. lyrata as outgroup and showed a clustering topology in general well supported with high bootstrap values (Figure 4). Our phylogenetic analyses of chloroplasts sequences indicated a clear division of the genus Brassica into two separate evolutionary linages, including four different clades (Figure 4), which was quite similar to the phylogeny using whole chloroplast sequences (Supplementary Figure S4).

One of these two separate linages contained three clades formed by the chloroplast genomes of B. rapa, B. oleracea, B. napus and B. juncea accessions (Figure 4). The clade I (strongly supported with 100% BP), contained the chloroplast genomes of 26 B. rapa accessions representing all different B. rapa morphotypes, together with chloroplast genomes of two B. napus accessions and of all the eight B. juncea accessions. Furthermore, the chloroplast genomes of these two B. napus accessions were grouped with those of three B. rapa (sarsons) in this clade. The chloroplast genomes of all the B. juncea accessions were grouped with those of the B. rapa accessions in clade I, supporting that B. rapa is the maternal ancestor of B. juncea. In clade II, the chloroplast genomes of all the 17 B. oleracea accessions we examined were placed together. Interestingly, clade III in the phylogenetic tree is composed of the chloroplast genomes of six B. napus and four B. rapa broccoletto accessions, significantly supported by bootstrap values. The chloroplast genomes of the six B. napus accessions that include both rapeseed and swede accessions, were closely related to those of the four broccoletto accessions and less related to those of the other two rapeseed B. napus accessions that were placed in clade I, close to the chloroplast genomes of the B. rapa sarsons. Clade IV represents the other lineage, and corresponds to B. nigra (black mustard) and B. carinata (Ethiopian mustard) chloroplast genomes, supporting that B. nigra is the maternal ancestor of B. carinata. The branch length of the tree indicated that B. nigra and B. carinata chloroplast genomes are more distant from B. oleracea, B. napus, B. juncea and B. rapa chloroplast genomes than the chloroplast genomes of these four species are from each other, suggesting much earlier divergence of B. nigra from the other diploid Brassica’s.

To estimate when separation of these chloroplast lineages in Brassica occurred, we used two calibration points of speciation events calibrated by valid fossil evidence (Hohmann et al., 2015). The results of our dating experiments are presented in Figure 5.

Our results indicate that, the ancestors of the current major Brassica groups were generated in several speciation events that started with a first radiation approximately 13.7 Mya separating clade IV (B. nigra and B. carinata) from the rest of Brassica species. B. nigra and B. carinata species diverged from each other 1.07 Mya. In a second period of speciation, the clade III (broccoletto and B. napus) separated from clade I (B. rapa, B. juncea and B. napus) and clade II (B. oleracea) at 2.98 Mya. In a third period, approximately 2.18 Mya, the chloroplast line shared by clade I (B. rapa, B. juncea and B. napus) diverged from clade II (B. oleracea). Later on, the Sarson chloroplast lines diverged from the B. rapa cluster approximately 0.86 Mya.