Phenotypic analysis of two non-heading cabbage varieties
The purple and green non-heading Chinese cabbage varieties (ZBC and LBC) were used in this study (Fig.1a, b). Total anthocyanin content of the two varieties was measured and was 0.32mg/g for ZBC, while only 0.04mg/g for LBC (Fig.1f). As the leaves developed, anthocyanins accumulated mainly in the upper leaf surface of ZBC, with almost no anthocyanin accumulation in the lower epidermis (Fig.1de). Previous studies had reported that anthocyanin accumulation was associated with low temperature [29]. To investigate whether anthocyanin accumulation in non-heading Chinese cabbage responds to low temperature, we treated the purple ZBC with 4℃ for 30 days. Compared with the control, ZBC after low temperature treatment showed internal leaf discoloration (Fig.1c). In addition, the total anthocyanin content also decreased by 43.75% compared with control, only 0.18mg/g (Fig.1f). In short, the upper epidermis of the purple ZBC accumulates a large amount of anthocyanin and is affected by low temperature.
Laser capture microdissection and RNA-Seq
To investigate genes specifically expressed on the upper epidermis of purple leaves, upper and lower epidermis tissues of leaves from the same growth period were collected from ZBC using a Laser capture microdissection frozen section method. Meanwhile, to avoid the difference other than anthocyanins between the upper and lower epidermis, the green LBC was used as a control. Hence, four sequencing libraries were constructed with three replicates. LS, LBC leaf upper epidermis (Fig.2a), LX, LBC leaf lower epidermis (Fig.2b), ZS, ZBC leaf upper epidermis (Fig.2c), ZX, ZBC leaf lower epidermis (Fig.2d). Four samples were sequenced by Illumina technology and aligned to the reference genome. After removing adapters and unknown or poor quality reads, 40,441,648-44,345,740; 41,387,390-9,937,566; 38,134,220-44,624,234 and 39,410,282-42,125,688 clean reads remained in ZS, ZX, LS and LX respectively. The Q30 score for each sample was 93.40-94.29% and mean GC content was 48.44% (Table S2).
DEG analysis was performed to identify the DEGs in each sample. The transcript levels of the genes were expressed in FPKM values. Four samples were compared each other, totally 6 groups. In group ZS_ZX, 4496 DEGs (1740 up and 2756 down) were identified. In group ZS_LS, 4461 DEGs (2546 up and 1915 down) were identified. In group ZX_LX, 5386 DEGs (3164 up and 2222 down) were identified. In group LS_LX, 4290 DEGs (1902 up and 2388 down) were identified. In group ZS_LX, 7720 DEGs (3939 up and 3781 down) were identified. In group ZX_LS, 7365 DEGs (4542 up and 2823 down) were identified (Fig.3a). The similar number of DEGs in group ZS_ZX and LS_LX indicating that the expression pattern and spatial location of these DEGs in leaves are similar. We also found that the number of DEGs in groups ZS_LX and ZX_LS was much more than that in other groups, indicating that the expression of genes in the upper and lower epidermis of different color leaves is remarkably. Anyway, the expression of genes in leaves is spatially variable, while the differences become more significant in different materials.
Functional annotation and classification of DEGs
Since anthocyanin accumulate mainly on the upper epidermis of ZBC leaves, we screened DEGs in group ZS_LS and ZS_ZX. In order to further explore DEGs expressed only on the upper epidermis of ZBC leaves. DEGs common to both group ZS_LS and group ZS_ZX were selected. Finally, we obtained 1234 genes that were differentially expressed only on the upper epidermis of ZBC leaves (Fig.3b). Then, GO functional enrichment analysis was performed on 1234 DEGs (p< 0.05). The results showed that the DEGs were mainly enriched in terms of response to UV-B (GO:0010224), biosynthesis process of anthocyanin-containing compounds (GO:0009718), response to karrikine (GO:0080167), response to sucrose (GO:0009744) and biosynthesis process of coumarin (GO:0009805) (Fig.3c; Table S3).
In addition, we also analyzed the DEGs for KEGG pathway enrichment. Metabolic pathways and biosynthesis of secondary metabolites were the most significant among all KEGG classifications, followed by flavonoid biosynthesis, phenylpropanoid biosynthesis, and glyoxylate and dicarboxylate metabolism (Fig.3d; Table S4).
Verification of RNA-Seq data by qRT-PCR
To verify the quality of the RNA sequencing results, we selected 9 genes involved from the anthocyanins synthesis for verification by quantitative reverse transcription PCR (qRT-PCR). Anthocyanin biosynthetic genes usually be divided into two subgroups: early biosynthetic genes (EBG) and late biosynthetic genes (LBG). EBGs were also flavonoid biosynthesis genes [6]. In our results, the expression levels of EBGs including CHS, CHI, LBGs including F3'H, DFR, ANS, and UGT75C1 were significantly upregulated in ZS (Fig.4). The anthocyanin-regulated genes MYB111, TT8 and EGL3 were also highly expression in ZS (Fig.4). RNA-Seq and qRT-PCR analysis showed similar expression patterns of these genes, indicating the reliability of the transcriptomic data in our study.
Expression patterns of genes involved in anthocyanin biosynthesis pathway
To enhance comprehension of the expression patterns of anthocyanin related genes in purple ZBC, we selected genes engaged in the anthocyanin synthesis pathway from 1234 DEGs and constructed an expression heat map. The biosynthetic pathway of phenylacetone involves four genes, BraC04g028280.1, BraC05g008710.1, BraC07g021110.1, BraC04g006690.1, all of which are annotated as phenylalanine ammonia lyase (PAL). There are seven genes related to the flavonoid synthesis pathway. BraC02g005310.1, BraC03g006200.1, BraC10g026540.1 are annotated as chalcone synthase (CHS), BraC07g022030.1, BraC09g053860.1 as chalcone isomerase (CHI), BraC09g049150.1 as flavanone 3-hydroxylase (F3H) and BraC10g030680.1 as flavanone 3’-hydroxylase (F3’H). These genes belong to the early biosynthetic pathway , provide precursor substrates for flavanol and anthocyanin synthesis. The anthocyanin biosynthetic pathway involves five genes, BraC09g018850.1 are annotated as dihydroflavonol 4-reductase (DFR), BraC03g052160.1, BraC01g013880.1 as anthocyanidin synthase (ANS), BraC08g010530.1, BraC10g012540.1 as UDP-glucosyltransferase (UGT). The biosynthesis and modification of anthocyanins is controlled by these late biosynthetic genes (Fig.5).
By analyzing the patterns of genes expression, we found that almost all of these anthocyanin-related genes were up-regulated in ZS, consistent with the phenotype of ZBC leaves. Notably, the expression of BcDFR (BraC09g018850) was most significantly up-regulated in ZS, consistent with the qRT-PCR results (Fig.4). These findings indicate that BcDFR may play a critical role in regulating the specific enrichment of anthocyanins on the upper epidermis of ZBC leaves.
Expression patterns of BcDFR in ZBC and LBC
To understand why BcDFR is highly expressed in ZS but barely expressed in other three parts, we cloned full-length BcDFR from ZBC and LBC. The gene region sequences of the two materials were highly conserved. Compared with ZBC, BcDFR of LBC showed an 11 base pair deletion at -585 and a 16 base pair insertion at -629 in the promoter region (Fig.6a). Two additional primers were developed based on the specific deletions and insertions sequence of BcDFR promoter, to accurately differentiate the genotypes of ZBC, LBC and F1. In contrast to ZBC, LBC showed 16bp insertion and 11bp deletion at these two loci and F1 showed heterozygous banding patterns (Fig.6bc). The insertion and deletion potentially be responsible for the low expression of BcDFR in LBC.
Analysis of regulatory elements and gene structure of DFR
To further understand the evolution of the DFR gene , we analysed the gene structure of it in Arabidopsis and six Brassica species including B. napus, B. rapa, B. carinata, B. nigra, B. juncea and B. oleracea. Eight DFR genes are highly conserved containing six exons and five introns. A domain analysis revealed 10 conserved domains within the DFR homologous gene, with only the absence of motif 7 in BjuB001305 (Fig.S1).
Cis-elements in the promoter region of the DFR gene indicate that these regulatory elements were associated with light responsiveness, MYB binding sites, the anaerobic induction, salicylic acid responsiveness, meristem expression and a variety of phytohormone responses. These findings indicate that DFR was involved in various aspects of plant growth and development and exerts a significant influence. In addition, there are more binding sites for light responsiveness and anaerobic induction elements in the DFR promoter (Fig.7), suggesting that some stress-related genes probably provide a critical role in the anthocyanin biosynthesis of plant. In addition, we also found many MYB binding sites on DFR promoter, suggesting that MYB transcription factors may act as upstream of DFR and then regulate its expression.
Prediction of BcDFR-mediated protein-protein interaction networks
To better understand BcDFR's function in anthocyanin synthesis pathway, we used STRING to predict the protein interaction network of BcDFR (Fig.S2). The results showed that BcDFR directly interacted with C4H, CHS, CHI, F3'H, FLS, UGT75C1 and TT8 proteins, while it also indirectly interacted with ANS, F3H, PAL, MYB111 and MYBL2. Interestingly, MYBL111, TT8, together with MYBL2, directly or indirectly act as catalysis proteins for the whole anthocyanin biosynthesis network. Further, RNA-Seq analysis showed that MYBL111, TT8 and MYBL2 were substantial up-regulated in ZS (Fig.S3). The specific expression of them also impacts the elevated expression of BcDFR in ZS, which may result the anthocyanins only accumulate in ZS.
Transient expression of BcDFR promotes anthocyanin accumulation
To further study BcDFR function in anthocyanin accumulation, we inject 35S:BcDFR into N. Benthamian (tobacco) leaves for the transient expression assay (Fig.8a). The total anthocyanin content of the control was 0.03 mg/g, while that of the 35S:BcDFR leaves was 0.059 mg/g (Fig.8b). This result indicates that BcDFR expression promotes anthocyanin synthesis and accumulation.