Plant materials were collected from yellow-seeded and black-seeded B. napus in different genetic backgrounds (Table 1, Figure 1). The yellow seed degree are calculated using the published reports (Li et al., 2012). The B₁ and B₂ were inbred lines ZY821 and ZS11, and Y₁ and Y₂ were GH06 and 05E258-1, which were used to represent typical black-seeded and yellow-seeded genetic types of B. napus, respectively. B₃ and Y₃ were represent two near isogenic lines of black- and yellow-seeded B. napus, selected from the successive backcross of the sixth generation with parent ZS11 and recurrent parent 06E241 (Qu et al., 2015). The plants were grown under normal field conditions at Chongqing Rapeseed Technology Research Center (CRTRC) in 2014. Field management essentially followed normal agronomic procedures. Various tissue organs, including stems (St), leaves (Le), flowers (Fl), buds (Bu), silique pericarps (SP), and seeds harvested at 10, 20, 30, 40, and 50 days after pollination (DAP) were sampled and stored at −80°C until used.

The full genome sequences of Brassicaceae were downloaded from multiple databases, including the BRAD database (http://brassicadb.org/brad/ftpTrans.php; Cheng et al., 2011), the Phytozome database (http://phytozome.jgi.doe.gov/pz/portal.html#!search; Goodstein et al., 2012), PGDD (PLANT GENOME DUPLICATION DATABASE; http://chibba.agtec.uga.edu/duplication/index/files; Lee et al., 2013), the Raphanus sativus Genome DataBase (http://radish.kazusa.or.jp/; Kitashiba et al., 2014), and the B. napus database (http://www.genoscope.cns.fr/brassicanapus/; Chalhoub et al., 2014). From amongst the species with full genome sequences, the following 17 species were selected: Aethionema arabicum (Aa), Arabidopsis halleri (Ah), Arabidopsis lyrata (Al), A. thaliana (At), Brassica napus L. (Bn), B. oleracea (Bo), B. rapa (Br), Boechera stricta (Bs), Capsella grandiflora (Cg), Capsella rubella (Cr), Camelina sativa (Cs), Leavenworthia alabamica (La), Raphanus sativus L. (Rs), Sisymbrium irio (Si), Schrenkiella parvula (Sp), Thellungiella halophile (Th), and Thellungiella salsuginea (Ts). All coding sequences (CDSs) and amino acid sequences were stored in a local Brassicaceae database using Geneious Pro 4.8.5 software (http://www.geneious.com/; Biomatters, Auckland, New Zealand). To identify the flavonoid biosynthesis pathway genes and their homologous genes, all amino acid sequences of flavonoid biosynthesis pathway genes in A. thaliana retrieved from TAIR 10 (http://www.arabidopsis.org/; Lamesch et al., 2012) were used as queries to search against the Brassicaceae protein models with HMMER3 (version 3.1b2 with Pfam HMM library Pfam 28.0; Finn et al., 2011). To identify the flavonoid biosynthesis genes, unique protein sequences of these genes from the A. thaliana genome were used as query using BLASTP program (Altschul et al., 1997) in the local Brassicaceae database developed in this research. All taxa were named using two-letter acronyms and gene type was used as the species gene name. Briefly, the first uppercase letter represents the genus, the second the species, and the following the gene name. A number at the end indicates the copy number. For example, A. thaliana phenylalanine ammonia lyase 1, which has only one copy in A. thaliana, is indicated by AtPAL1.

To assign the location of flavonoid biosynthesis pathway genes in the Brassicaceae species genomes, the GFF genome files were downloaded from the aforementioned databases. Then, MapChart 2.0 was used to draw graphic representations of their corresponding physical position on pseudo-molecular chromosomes of Brassica crops.

Based on previously described methods, all sequence alignments for each flavonoid biosynthesis pathway gene superfamily were performed using ClustalW2 software (Larkin et al., 2007), and phylogenetic analysis was carried out using Molecular Evolutionary Genetics Analysis (MEGA) 6.0 (Tamura et al., 2013) with a maximum likelihood (ML). In the ML method, phylogenetic trees were constructed using the JTT+I+G substitution model in PhyML version 3.0.1. To ensure the accuracy of the phylogenetic tree, each tree was subjected to bootstrap analysis with 1000 replicates (Guindon et al., 2010). Finally, all the phylogenetic trees were visualized using FigTree v1.4.2 (http://tree.bio.ed.ac.uk/software/figtree/).

To characterize differences in expression of the 18 genes associated with the flavonoid biosynthesis pathway between the yellow- and black-seeded varieties, we designed the primers in consensus region based on alignment the gene sequence (Table 3). Then the total RNA was extracted from various tissues using the RNAprep Pure Plant Kit (TIANGEN BIOTECH, Beijing, China) according to manufacturer's instructions. Then, 1 μg RNA sample was reverse transcribed with the Oligo dT-Adaptor Primer using the RNA PCR Kit (AMV) Ver. 3.0 (TaKaRa, http://www.takara.com.cn). To monitor sample uniformity of initial RNA input and RT efficiency, Brassica napus 26S rRNA was used as the internal control according to the previously described method (Qu et al., 2013).

Real-time PCR was performed using SYBR® Premix Ex Taq™ II (Perfect Real Time) (TaKaRa, China) in a 20 μl volume that included 10 μl of SYBR® Premix Ex Taq™ II, 2 μl (100 ng) of template cDNA, and 0.4 μM of each PCR primer. All primer sequences used for the qRT-PCR are listed in Table 3, designed according to the methods described in our previous research (Qu et al., 2013). Then the specific primers used in this study and cycling conditions were 95°C for 2 min, followed by 40 cycles at 95°C for 10 s (denaturation) and 60°C for 20 s (annealing and extension). The melting curves of each PCR application were obtained using the following cycling conditions: 95°C for 10 s followed by a constant increase in temperature between 65 and 95°C at an increment of 0.5°C/cycle, and samples were run on the Bio-Rad CFX96 Real Time System (USA). The relative expression of the target genes was analyzed using the 2^−ΔΔCt method with BnACTIN7 (EV116054) and BnUBC21 (EV086936) as internal controls (Wu et al., 2010). Three biological replicates for each sample were used for real-time PCR analysis and three technical replicates were analyzed for each biological replicate. Then the values represent the average ± SD of three biological replicates with three technical replicates of each tissue and organ. Relative gene expression levels were normalized according to the expression values in black-seeded at 10 DAP.

To identify members of the Brassicaceae gene subfamily that are involved in the flavonoid biosynthesis pathway, the whole genome sequences of 17 species were downloaded from multiple public databases (Materials and Methods). In total, 649 nucleotide sequences of 21 genes and their associated amino acid sequences were respectively were identified using the HMM profile (Finn et al., 2011). In addition, we used the nucleotide and encoded amino acid sequences of four genes from B. juncea, Bj_PAL1 (ACX31148.1), Bj_TT3 (ADB45307.1), Bj_TT5 (ADB45305.1), and Bj_TT18 (ACH58397.1), that we identified and aligned previously (Qu et al., 2013). The copy numbers of each gene varied from 20 to 91, and the genes were widely distributed in the 17 Brassicaceae species examined. The encoded amino acid sequences ranged from 84 aa (Bn_UGT2d) to 1676 aa (Sp_TT10; Supplementary Table S2). In addition, the number of gene families also varied by species; for example, only 10 copies of PAL3 were found in 6 species, but 43 C4H gene copies occurred in 17 species with as many as 10 orthologs (BnC4H) in B. napus (Supplementary Table S2). However, TT4, TT6, and AHA10 were not identified in all Brassicaceae species. For example, TT4 was not found in A. halleri and C. rubella, TT6 was absent from A. halleri, and AHA10 was not present in A. lyrata (Supplementary Table S2). We used the sequences of flavonoid biosynthesis genes from A. thaliana as query to search genome databases, and identified highly conserved sequences, including LESS ADHESIVE POLLEN5/6 (LAP5/6) and TT4 in A. halleri and C. rubella, DMR6-LIKE OXYGENASE1/2 (DLO1/2) and TT6 in A. halleri, and members of the Autoinhibited H(+)-ATPase (AHA) superfamily, which might be homologs of TT4, TT6, and AHA10 and function in flavonoid biosynthesis in A. lyrata. Here, 36 copies of 21 genes of the flavonoid biosynthesis pathway were identified in A. thaliana, but 53 were identified in B. rapa, 50 in B. oleracea, and 95 in B. napus, in accordance with the fact that sequences present as a single copy in the A. thaliana were present in 2–8 copies in B. napus (Cavell et al., 1998), and with the observation that excessive gene loss is typical after polyploidization in eukaryotes (Sankoff et al., 2010; Wang et al., 2011).

All of the flavonoid biosynthesis pathway genes in 17 Brassicaceae species were mapped onto pseudo-molecules or chromosomes using GFF files of their nucleotide sequences (Supplementary Table S3). Brassica crops are the ideal model for studying genome evolution (Wang et al., 2011). To intuitively assign the physical position to the chromosomes of B. rapa, B. oleracea, and B. napus, all gene members of flavonoid biosynthesis pathway were mapped to their chromosomes [52 (98.1%) gene copies in B. rapa, 39 (78.0%) gene copies in B. oleracea, and 84 (88.4%) gene copies in B. napus] and 22 (1 in B. rapa, 11 in B. oleracea, and 11 in B. napus) were distributed on the unanchored scaffolds (Figure 2, Supplementary Table S3), which showed strong collinearity between the A subgenomes from B. rapa and B. napus and the C subgenomes from B. oleracea and B. napus (Figure 2). This result indicates that the gene copies are distributed in orthologous blocks in each genome, and that substantial genome reshuffling had occurred. For example, copies of TT10 and TT16 were not identified on B. napus chromosome A02, and the differential gene copies located in orthologous blocks on chromosome C06 differed greatly between B. oleracea and B. napus (Figure 2). These findings are consistent with the fact that Brassicaceae genomes underwent Brassicaceae-lineage-specific whole genome triplication, followed by diploidization after divergence from their common ancestor (Lysak et al., 2005; Town et al., 2006; Mun et al., 2009; Wang et al., 2011; Cheng et al., 2013).

Despite different flavonoid biosynthesis pathway genes having similar functions in Brassicaceae species, the phylogenetic relationships of members of each superfamily have not been comprehensively studied. To gain insight into the functional roles of flavonoid biosynthesis pathway genes that may have arisen during evolution, we performed a phylogenetic analysis using the putative amino acid sequence based on each member of the flavonoid biosynthesis pathway superfamily in Brassicaceae species.

Initially, we conducted a BLASTP analysis against the local database, which included 17 Brassicaceae species, followed by HMM searching (Finn et al., 2011). We identified 649 sequences of 21 flavonoid biosynthesis pathway gene family members, and found that the number of gene copies varied from 1 to 10 in each family (Supplementary Table S2), which were mapped onto pseudo-molecules or chromosomes using GFF files of their nucleotide sequences (Supplementary Table S3). To intuitively assign the physical position to the chromosomes of B. rapa, B. oleracea, and B. napus, all gene members of flavonoid biosynthesis pathway were mapped to their chromosomes [52 (98.1%) gene copies in B. rapa, 39 (78.0%) gene copies in B. oleracea, and 84 (88.4%) gene copies in B. napus] and 22 (1 in B. rapa, 11 in B. oleracea, and 11 in B. napus) were distributed on the unanchored scaffolds (Figure 2, Supplementary Table S3), which showed strong collinearity between the A subgenomes from B. rapa and B. napus and the C subgenomes from B. oleracea and B. napus (Figure 2). This result indicates that the gene copies are distributed in orthologous blocks in each genome, and that substantial genome reshuffling had occurred. For example, copies of TT10 and TT16 were not identified on B. napus chromosome A02, and the differential gene copies located in orthologous blocks on chromosome C06 differed greatly between B. oleracea and B. napus (Figure 2). These findings are consistent with the fact that Brassicaceae genomes underwent Brassicaceae-lineage-specific whole genome triplication, followed by diploidization after divergence from their common ancestor (Lysak et al., 2005; Town et al., 2006; Mun et al., 2009; Wang et al., 2011; Cheng et al., 2013).

In addition, we performed a phylogenetic analysis to identify each flavonoid biosynthesis pathway gene type, and constructed their corresponding phylogenetic trees using MEGA 6.0 with the ML method and modified the tree using FigTree v1.4.2. According to the phylogeny generated using the ML method, all representative sequences for each gene member from the Brassicaceae species formed a well-supported clade, which was classified into the Arabidopsis or Brassicaceae subgroups, but each gene was assigned to monophyletic clades using other substitution models (Figures 3–5, and Supplementary Figures S1–S10). For example, we identified 4 members of PAL, which encodes key enzymes of the phenylpropanoid pathway, that were clearly grouped into four subclades (Figure 3). The gene members, Bn_PAL1b and Si_PAL3 belonged to the same PAL2 subcategories, and Aa_PAL4 and Cs_PAL4c clustered alone or in the vicinity of the PAL3 and PAL4 group (Figure 3; PAL1, Gray; PAL2, Light blue; PAL3, Pink; PAL4, Green). We could thus predict the functions of specific flavonoid biosynthesis pathway genes in Brassicaceae species, because genes with similar functions tend to be retained with orthologous genes (Koonin, 2005). In addition, Brassica species are an ideal model for systematically studying polyploidy genome evolution. We identified one clade of genes of the flavonoid biosynthesis pathway, and the copies in B. rapa, B. oleracea, and B. napus (e.g., C4H, TT4, TT6, TT18, TT19, and UGT2) were divided into different subclades or monophyletic subclades by phylogenetic analysis (Figures 4A–F). However, multiple gene copies of TT12, TTG2, and BAN were present in the subclades in B. rapa, B. oleracea, and B. napus, which have high levels of divergence from their common ancestor with A. thaliana (Figures 5A–C).

Based on the collinearity of orthologs, we identified syntenic orthologous genes of the flavonoid biosynthesis pathway among the genomes of the Brassicaceae species (Supplementary Table S4). In the genomes of B. rapa, B. oleracea, and B. napus, the copies of C4H were distributed in the LF, MF1, and MF2 subgenomes. Besides, the copies of C4H were also identified outside the genome triplication segments in the genomes of Brassica crops. However, 10 of 25 gene members had a single copy in the B. rapa and B. oleracea genomes, but multiple copies in B. napus. Moreover, we identified 10 pseudo-copies of five genes in seven species, i.e., 3 of which in C. rubella (C4H, TT10, and UGT2), 2 copies of TT1 in L. alabamica and S. irio, and 5 copies of TT6 in A. arabicum, B. rapa, B. oleracea, and B. napus (Supplementary Table S4). These findings can be used to systematically study gene retention in the triplicated genomes of Brassicaceae species, as well as to understand the evolutionary history of these orthologous genes among Brassicaceae species.

In this study, we analyzed the expression patterns of 18 TT-type genes involved in flavonoid biosynthesis by qRT-PCR in various tissues of black- and yellow-seeded lines of B. napus with different genetic backgrounds. Except for BnTT4, BnTT5, BnTT6, and BnTT7, the expression level of 11 structural genes was much higher in developing seeds than in other tissues in B. napus (Figure 6). The expression level of BnBAN, BnTT12, BnTT19, and BnAHA10 peaked during the early stages of seed development, which were the EBGs (Figures 6G–I,K, Table 4). BnTT3, and BnTT18 were expressed in a similar pattern as BnBAN, BnTT12, BnTT19, and BnAHA10, but peaked later in development, which were the LBGs (Figures 6E,F,J, Table 4). Thus, these genes may be regulated by the same upstream gene or they may have a synergistic effect on the flavonoid biosynthetic pathway during seed development. In addition, the expression levels of BnTT5 and BnTT7 did not differ significantly (Student's t-test, P > 0.05) among the rapeseed tissues, whereas BnTT4, BnTT5, BnTT6, and BnTT7 expression was higher in the buds and flowers of B. napus (Figures 6A–D, Table 4). Furthermore, the expression level of five regulatory genes, BnTT2, BnTT8, BnTT16, BnTTG2, and BnTT1, was much higher in developing B. napus seeds, and peaked during the early and middle stages of seed development, which may be the EBGs (Figures 6L,M,O,P,Q, Table 4). However, BnTT15 was expressed at higher levels in flowers than in other organs, indicating that this gene may play an important role in flower formation (Figure 6R).

In this study, most of the key genes involved in flavonoid biosynthesis were found to exhibit similar expression patterns, and the expression levels of genes varied greatly between the yellow- and black-seeded lines of B. napus (Figure 7). The structural genes involved in flavonoid biosynthesis, i.e., BnTT3, BnTT18, BnBAN, BnTT12, BnTT19, and BnAHA10, had similar expression patterns in the developing seeds of different genetic backgrounds, but had significant differences (Student's t-test, P < 0.01) in expression level between the black- and yellow-seeded lines (Figures 7C,E–I,K, Table 4). BnTT4, BnBAN, BnTT12, BnTT19, and BnAHA10 expression peaked at 20 DAP, which were EBGs (Figures 7A,G–I,K, Table 4), and BnTT6, BnTT3, and BnTT18 expression peaked at 40 DAP, were belong the LBGs (Figures 7E,F, Table 4). These genes were expressed at higher levels in the black-seeded than in the yellow-seeded lines, especially at 20 DAP (Figures 7A,H,K, Table 4), 30 DAP (Figures 7G,I), and 40 DAP (Figures 7B,E,F, Table 4), but the expression patterns of these genes were hardly affected by genetic background. Additionally, BnTT7, BnTT5, and BnTT10 expression peaked at different development stages both in the black- and yellow-seeded lines (Figures 7B,D,J, Table 4, and Supplementary Table S5), indicating that the expression of these genes was affected by the genetic background. As for the expression patterns of structural genes, regulatory genes, such as BnTT2, BnTT8, BnTT16, BnTTG2, BnTT1, and BnTT15, had similar expression patterns in the black- and yellow-seeded B. napus, with the genetic background having a minimal effect, and the expression of these genes peaked during the early and middle developmental stages (20 and 30 DAP, EBGs; Figures 7L,M,O–R, Table 4). In addition, the expression patterns and levels of BnTTG1 showed obvious variations during seed development of different lines of B. napus, suggesting that they were significantly affected by genetic background (Figure 7N).

Although the expression patterns of most of the genes involved in flavonoid biosynthesis did not significantly differ (Student's t-test, P > 0.05) between the developing seeds of the black- and yellow-seeded B. napus, the expression levels of these genes did differ and peaked at different developmental stages (Figure 7), indicating that these genes can also be classified as Early Biosynthetic Genes and Late Biosynthetic Genes, as in Arabidopsis (Lepiniec et al., 2006). The expression of genes that function upstream in the flavonoid biosynthesis pathway peaked before those that functioned later, suggesting that the downstream genes were influenced by the upstream genes. For example, the expression of BnTT4, which encodes a protein that catalyzes the first committed step of flavonoid biosynthesis (Albert et al., 1997; Tang et al., 1997), peaked at 20 DAP, which was before the expression of the downstream genes BnTT5, BnTT6, and BnTT7 peaked (Figure 7). All flavonoids and isoflavonoids are derived from the nargingenin chalcone generated by this first enzyme. Moreover, the expression of most of the regulatory genes also peaked at 20 DAP (Figure 7), suggesting that they may be essential for regulating the expression of genes involved in flavonoid biosynthesis during the early stages of seed development. Although, the expression of BnTTG2 peaked later than did that of BnTT2 and BnTT8, the expression of BnTT2 was greater in yellow- than in black-seeded B. napus, and BnTTG1 expression was also largely affected by the genetic background (Figures 7L–N,P), indicating that the mechanism of BnTTG2 is different from that of its orthologs in A. thaliana (Lepiniec et al., 2006). In addition, the expression of BnTT6, BnTT18, and BnTT2 peaked sooner in the yellow-seeded than in the black-seeded lines (Figures 7C,F,L).

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.