One B. napus inbred lines which exhibited purple leaves (abbreviated as PL) at the seedling stage and ZS11 (Zhongshuang 11) which showed green leaves (abbreviated as GL), as well as the F₁ progeny by a cross between these two lines were used in this study. Seeds of PL, GL, and F₁ were sown directly on September 30th, 2020, and three seedlings of each were maintained in the same pot (43 cm length × 40 cm wide × 25 cm high) to minimize the developmental or environmental effects. All the plants were grown under natural lighting in a greenhouse located at the Hunan University of Science and Technology (27°91′05″N, 112°92′60′'E, China). Considering that the leaf color of PL showed a dynamic change from light purple to dark purple and then to green, the fresh leaves at the seedling stage (<101 days after sowing) and bolting stage (150 days after sowing) were collected and stored at −80°C.

The total anthocyanin content was measured using the method described previously (Li et al., 2016; Zhang et al., 2020). In brief, the young leaves PL, GL, and F₁ with three replicates were ground into powder after low-temperature freeze drying (SJIA-10N, Ningbo, China). Then 0.1 g powder was mixed with 1 ml of 95% ethanol: 1.5 M HCl (85:15, v/v) for 20 h at 4°C in the dark with moderate shaking. The aqueous phase was collected after centrifuging at 12,000 ×g for 10 min, and the absorbance was subsequently measured at 530 nm using a Cary 60 UV-Vis (Agilent Technologies, Palo Alto, United States).

Previously, the ABGs in B. rapa were well-characterized based on the syntenic analysis results between B. rapa and A. thaliana (Guo et al., 2014). Following the method provided by Cheng et al. (2012) and Guo et al. (2014 and 2019), syntenic analysis across B. napus, A. thaliana, B. rapa, and B. oleracea was performed, and ABGs were identified using SynOrths (http://brassicadb.cn/#/Download/) with the same parameters (NumQ = 20, NumR = 100, RatioQR = 0.2). Gene pairs that are the best hits or with Blastp e-values <1E⁻²⁰ were selected for further analysis. The reference genome of A. thaliana was downloaded from BRAD (http://brassicadb.cn), as described before. The B. oleracea (genotype HDEM) and B. rapa (genotype Z1) reference genomes were obtained from https://www.genoscope.cns.fr/externe/plants (Rousseau-Gueutin et al., 2020). The B. napus (genotype ZS11) reference genome was obtained from http://cbi.hzau.edu.cn/bnapus/(Song et al., 2020), respectively. Integrative software TBtools was used for determining the chromosomal locations of ABGs of B. napus (Chen et al., 2020).

The entire of fully expanded young leaves from PL, GL, and F₁ with three biological replicates at 41, 91, and 101 days after sowing were collected for further RNA-seq, respectively. The total RNA of samples was extracted using the RNA Easy Fast Plant Tissue Kit (TIANGEN, Beijing, China) according to a standard protocol. After that, 1.5% agarose gel, a NanoPhotometer® spectrophotometer (IMPLEN, CA, United States), and a Bioanalyzer 2,100 system (Agilent Technologies, CA, United States) were also used to check the purity and integrity of RNA. Samples with RIN scores≥7.5 and OD260/280 = 1.8–2.2 were utilized for further sequencing. Sequencing libraries were generated using the NEBNext® UltraTM RNA Library Prep Kit for Illumina® (NEB, United States), following the manufacturer’s recommendations. The libraries were sequenced on an Illumina NovaSeq 6,000 platform and 150-bp paired-end reads were generated by Shanghai Majorbio Bio-pharm Technology Co., Ltd.

After the high-throughput sequencing, the low-quality reads which contained adapters and ploy-N were trimmed using Trimmomatic. Then, the clean reads were mapped to the high-quality B. napus (genotype ZS11; http://cbi.hzau.edu.cn/bnapus/) reference genomes using TopHat2 (http://tophat.cbcb.umd.edu/) with the default parameters. Software RSEM was used for read counting, and DESeq2 was used for differential expression analysis. Genes with adjusted p-values < 0.05 were classified as differentially expressed. The expression levels of individual genes were also quantified using the transcripts per million (TPM) method. The KEGG enrichment analysis of differentially expressed genes was implemented using the free online platform of Majorbio Cloud Platform (www.majorbio.com). KEGG terms with corrected p-values less than 0.05 were considered to be significantly enriched.

The extracted RNA was also used for qRT–PCR. First-strand cDNA synthesis was performed using a Tsingke Goldenstar RT6 cDNA Synthesis Kit ver.2 (Tsingke, Beijing, China). Differentially expressed ABGs were selected, and specific primers were designed for qRT–PCR using Oligo 7 and NCBI Primer-BLAST (https://www.ncbi.nlm. nih. gov/tools/primer-blast/). The β-actin gene (GenBank accession No. AF111812) was used as an internal control, and SYBR Premix Ex TaqII with a Bio-Rad CFX96 Real-Time Detection System was used as described previously (Zhang et al., 2018).

During the seedling stage, the PL displayed purple pigments in the surface and thin vein of leaves, while the F₁ only showed faint purple in the leaf surface (Figures 1A,B). The microscopic observation of sections revealed that the purple pigments were mainly accumulated in the adaxial epidermis and subepidermal cells of PL leaves, while it decreased in the top cell layer of F₁ (Figures 1C,D). In contrast, GL with common green color only exhibited green chlorophyll layers in leaves (Figure 1E).

Although the average anthocyanin contents in PL (1.09 mg/kg DW) were about 3.9-fold higher than those in GL (0.28 mg/kg DW, p < 0.05, paired Students’ t-test), they varied at different developmental stages in the leaves of PL (Figure 1F). During the seedling stage (<101 days), the total anthocyanin contents of PL were stable at initial weeks (0.88 mg/kg DW, 28–41 days after sowing) and increased after that (76–91 days) but dropped dramatically from 1.83 mg/kg (91 days) to 1.10 mg/kg (101 days). No significant difference (p > 0.05, t-test) was observed between PL (0.39 mg/kg) and GL (0.37 mg/kg) at the bolting stage (150 days), which was correlated with the color change from purple to green in the leaves of PL. The average anthocyanin contents in F₁ (0.64 mg/kg DW) were higher than those in GL but significantly lower than those in PL (p < 0.05, paired t-test), suggesting that incompletely dominant gene or genes might control the purple leaf trait in B. napus.

Based on the method from a previous study (Guo et al., 2019; Zhang et al., 2020), a total of 74, 93, and 152 genes were identified to be ABGs in B. rapa, B. oleracea, and B. napus (Supplementary Table S1). As compared with 41 ABGs in A. thaliana, FLS2, FLS4, FLS5, FLS6, CPC, MYB111, MYB113, and MYB114 have not been identified in B. rapa, B. oleracea, and B. napus; these observations suggested that these genes were lost following the whole-genome triplication. Among the 152 genes identified in B. napus, 98 were classified into structural genes, 50 were regulatory genes, and 4 were transport genes. The number of these three types of genes was quadrupled as compared with A. thaliana, which may be caused by whole-genome triplication and subsequent allopolyploidization. Using the total A. thaliana and B. napus genes as references, the 152 ABGs represent 0.15% of the 100,919 annotated B. napus genes, which was similar to that in A. thaliana (0.15%). Meanwhile, the overall expansion levels of structural, regulatory, and transport genes also exhibited no significant difference to the whole genome (Table 1, chi-squared test, p > 0.05).

Among these 152 ABGs, 150 were anchored to the 19 chromosomes of B. napus, in which 66 genes were distributed on the A subgenome and other 84 genes were located on the C subgenome, ranging from 1 to 16 on each chromosome (Figure 2A). Among the 150 anchored ABGs, 33.3% (50/150) were originated from tandem duplication events. The 50 tandemly duplicated genes could be divided into 22 groups, with 9 groups on the A subgenome (19/50) and 13 groups on the C subgenome (31/50), of which 19 groups contained two genes each, two groups contained 3 genes each, and 1 group contained 6 genes each. Except for one set of tandemly duplicated gene on the chromosome A02 and two set of tandemly duplicated genes on chromosome C02 and chromosome C05, respectively, all the tandemly duplicated genes had homologous copies both in A and C subgenomes. These indicated that most tandem duplication events occurred before polyploidization. The syntenic analysis of ABGs among A. thaliana, B. rapa, B. oleracea, and B. napus showed that the whole-genome triplication of Brassica and subsequent allopolyploidization contributed to the expansion of ABGs in B. napus (Figure 2B). Although some genes were lost as a result of gene fractionation that occurred following the whole-genome triplication, a large proportion of ABGs exhibited syntenic between B. napus and B. rapa/B. oleracea.

To investigate the molecular mechanisms of anthocyanin accumulation, young leaves with three biological replicates at three different developmental stages (41, 91, 101 days) from GL, F₁, and PL were collected for library construction (3 × 3 × 3 = 27 libraries) and high-throughput RNA sequencing (Supplementary Table S2). In summary, about 127.0 G clean reads were obtained after filtering, and more than 90% of them can map on the B. napus reference genome released recently (Song et al., 2020). Gene expression levels were calculated using the TPM methods, and gene expression correlation showed good consistence between replications (Supplementary Figure S1).

The comparative transcriptome analysis of B. napus among GL, F₁, and PL showed that the number of DEGs was varied in different comparisons at different developmental stages (Figure 3A). The highest number of DEGs was found between PL and GL at 91 days when PL showed highest anthocyanin content, for 6,263 genes were upregulated and 3,795 genes were downregulated. While the number of DEGs were decreased between F₁ and the parental lines (F₁ vs GL and PL vs F₁) with the development of plants, and the lowest number of DEGs was found between F₁ and GL at 101 days when the anthocyanin content showed no significant difference. About 8,167, 10,058, and 8,291 genes were differentially expressed between PL and GL at 41, 91, and 101 days after sowing, and 4,054 of them were shared in the aforementioned comparisons (Figure 3B). The KEGG analysis of those DEGs between PL and GL suggested that the genes involved in anthocyanin biosynthesis were significantly enriched at 41, 91, and 101 days (Figure 3C).

To identify key genes responsible for leaf color variation in B. napus, the expression patterns of ABGs were investigated. Among the 152 ABGs, the expression levels of 57 genes were extremely low (average TPM value < 1) or undetected (Supplementary Table S3). Of the other 95 expressed ABGs, 22 genes were differentially expressed between PL and GL. These 22 DEGs are equally distributed in the A subgenome (11/22) and C subgenome (11/22), while the average expression level of the 11 DEGs in C subgenome was higher than that in the A subgenome (p < 0.05, two-tailed Student’s t-test; Figure 4B). Based on these results, the expressed levels of DEGs in the anthocyanin biosynthetic pathway from different leaf color samples at different stages were demonstrated (Figure 4). The majority of upstream genes in the anthocyanin biosynthetic pathway (45/46 in the phenylpropanoid pathway and 32/35 in EBGs) were not differentially expressed, implying that these genes might not play important roles in anthocyanin accumulation differences. While these genes which are involved in LBGs, such as DFR (3/3), ANS (3/3), and UFGT (6/9), as well as TT19 (3/4) and TT8 (2/3) were significantly upregulated in PL, as compared with GL at 41 and 91 days after sowing (Supplementary Table S3, Figure 4A). As to F₁, the overall expression level of 22 ABGs was significantly higher than that in GL at 41, 91, and 101 days but lower than that in PL at 41 and 91 days (p < 0.05, two-tailed Student’s t-test; Figure 4B). Except for one C4H gene (BnaA03G0147400ZS), all the other genes (21/22) involved in early or late anthocyanin biosynthetic exhibited the similar trend of expression level PL > F₁>GL at 41 and 91 days (Figure 4A). Quantitative RT-PCR results also confirmed the overexpression of these genes in PL and F₁, suggesting that these genes might be the essential factors for anthocyanin accumulation at the early stage (Supplementary Figure S2, Supplementary Table S4). It is also worthy to note that the expression levels of most upregulated (p < 0.01, paired t-test) genes in PL at 41 and 91 days were decreased due to the similar levels of GL (p = 0.1175, paired t-test) at 101 days, which was consistent with the decrease in anthocyanin contents and the leaf color phenotype in PL (Figure 4B and Figure 1F).

To further understand the relationship between expression patterns of the 22 differentially expressed ABGs and anthocyanin accumulation, correlation coefficients were calculated (Supplementary Figure S3). Significant correlations (p < 0.05, Pearson correlation coefficient) between expression and anthocyanin contents were observed in UFGT (BnaA08G0082300ZS, r = 0.60), TT8 (BnaA09G0260600ZS, r = 0.62), TT19 (BnaC09G0492500ZS, r = 0.76; BnaA10G0196600ZS, r = 0.82; BnaA02G0068400ZS r = 0.72), F3H (BnaC08G0314600ZS, r = 0.79), and F3 'H (BnaA10G0256900ZS, r = 0.64).