Uncovering the mechanism of anthocyanin accumulation in a purple-leaved variety of foxtail millet (Setaria italica) by transcriptome analysis

Plant material and growth conditions

The seeds of foxtail millet varieties B100 (purple-leaved variety) and Yugu 1 (green-leaved variety) were sown in an experimental plot located at 37°N, 112°E. The leaves were collected for analysis when plants reached a mature four-leaf stage.

Morphological observations

The leaves were sliced using two, side-by-side blades for morphoanatomical evaluation. The slices were embedded in water and were then observed under a microscope.

Measurement of plant pigments

Anthocyanin was extracted from leaves with 0.1% HCl-methanol at 4 °C for 12 h in dark. The absorption values of the extraction were measured at 530 and 700 nm wavelengths using the mixture of supernatant and reagent 1 (50 mmol L⁻¹ KCl, 150 mmol L⁻¹ HCl, pH = 1.0), followed by the mixture of supernatant and reagent 2 (100 mmol L⁻¹ NaAC, 240 mmol L⁻¹ HCl, pH = 4.5). The method of calculation was as follows (Luo et al., 2017):

Anthocyanin (mg/g FW) = [(A₅₂₀ − A₇₀₀)_reagent 1 − (A₅₂₀ − A₇₀₀)_{reagent 2}] × M × V × n/(ɛ × m).

M: relative molecular weight of anthocyanins, 449.2 g/mol; V: extraction volume; ɛ: molar extinction coefficient of anthocyanins, 2.69 × 104 L/mol/cm; n: dilution times; m: sample weight.

Chlorophyll and carotenoids were extracted from leaves using 80% acetone for 24 h. Then the absorption values were measured from the extraction at 663, 646 and 470 nm wavelengths (Porra, 2002). The method of calculation was as follows (V: extraction volume, W: sample weight):

Chlorophyll a (mg/g Fw) = (12.21A₆₆₃ − 2.81A₆₄₆) × V/W × 1,000

Chlorophyll b (mg/g Fw) = (20.13A₆₄₆ − 5.03A₆₆₃) × V/W × 1,000

Chlorophyll (mg/g Fw) = chlorophyll a (mg/g Fw) + chlorophyll b (mg/g Fw)

Carotenoid (mg/g Fw) = [1,000A₄₇₀ − 3.27 × (12.21A₆₆₃ − 2.81A₆₄₆) − 104 × (20.13A₆₄₆ − 5.03A₆₆₃)]/229 × V/W × 1,000.

The raw data of the pigment measurement is shown in Table S5.

Total antioxidant capacity assay

Antioxidants were extracted from leaves with 80% methanol. The antioxidant capacity was measured using the ferric-reducing antioxidant power (FRAP) method (Benzie & Strain, 1996).

Trolox solutions were prepared using methanol as the solvent of gradient concentrations. The Trolox solutions and FRAP solution were mixed for reaction for 1 h in the dark. Then we measured the absorption value of the mixture at 593 nm wavelength. The standard curve was drawn by taking the Trolox concentration as the x-axis and the absorbance as the y-axis.

The leaf extracts and FRAP solution were mixed for reaction for 1 h in the dark. Then we measured the absorption value of the mixture at 593 nm wavelength. The antioxidant capacity was assayed according to the standard curve and absorbance. The raw data of this assay is shown in Table S6.

Measurement of soluble sugar content

The leaves were homogenized with 80% ethanol and were extracted at 80 °C for 30 min. The supernatant was transferred to a fresh, clean tube after centrifugation at 3,500 g for 10 min. The sediment was then extracted twice using the above procedure with 80% ethanol. The supernatant was added into an anthrone-sulfuric acid reagent and was boiled for 10 min. After cooling to room temperature, the absorption value of the extraction was measured at 620 nm (Leng et al., 2016). The raw data of this assay is shown in Table S6.

Transcriptome sequencing and analyses

The total RNA was extracted from leaves using RNAiso Plus (9109; Takara). The RNA preparations were evaluated with agarose gel electrophoresis and were quantified using a microultraviolet–visible fluorescence spectrophotometer.

Three biological replicates of both Yugu1 and B100 were used in RNA-sequencing. The statistical power of this study was calculated using RNASeqPower (https://bioconductor.org/packages/release/bioc/html/RNASeqPower.html) as P < 0.05. The quality control of raw data was performed as in a previous study (Wang et al., 2017). We retrieved 144,677,546 clean reads, including 43.4 G total bases. The sequencing depth was 100x. The clean data were aligned to the foxtail millet genome (Yugu1, version 2.2: https://data.jgi.doe.gov/refine-download/phytozome?organism=Sitalica&expanded=312) using HISAT2 software (Kim, Langmead & Salzberg, 2015), and the mapped reads were analyzed by String Tie software (Pertea et al., 2015). The abundance of transcripts was normalized using the fragments per kilobase of transcript per million mapped reads (FPKM) method (Florea, Song & Salzberg, 2013). Prior to differential gene expression analysis, for each sequenced library, the read counts were adjusted using the edgeR program package with one scaling normalized factor. Differential expression analysis of the two samples was performed using the EBSeq R package. The resulting false discovery rate (FDR) was adjusted using the posterior probability of being DE (PPDE). The FDR < 0.05 & |log2 (fold change)|≥1 was set as the threshold for significantly differential expression. Differential expression analysis of the two conditions/groups was performed using the DESeq R package (1.10.1). DESeq provided statistical routines for determining the differential expression in digital gene expression data using a model based on the negative binomial distribution. The resulting P values were adjusted using the Benjamini and Hochberg’s approach for controlling the FDR. Genes with an adjusted P-value < 0.05 by DESeq were assigned as differentially expressed.

We performed Gene Ontology (GO) enrichment analysis by the GOseq R package method (Young et al., 2010). The KEGG database (http://www.genome.jp/kegg/) was used to discover the regulatory pathways of genes (Kanehisa et al., 2008). We used KOBAS (Mao et al., 2005) to determine the enrichment of DEGs in KEGG pathways.

The raw data of the RNA sequencing was submitted to an SRA database (BioProject accession number: PRJNA777600).

Quantitative real-time PCR analysis

The total RNA was isolated as mentioned above. The first-strand cDNA was synthesized using the PrimeScript RT reagent kit with gDNA Eraser (RR047A; Takara). qRT-PCR was performed using the Bio-Rad CFX96TM Touch real-time PCR detection system with SYBR Premix Ex TaqTM II (Tli RNaseH Plus) (RR820A; Takara). Transcript levels were normalized to the SiActin (Seita.5G464000). The raw data of qRT-PCR is shown in Table S7. The $2\stackrel{ˆ}{}$-ΔΔCt method was used to calculate the expression level of the genes and the one-way ANOVA was used to identify significant differences (http://vassarstats.net/anova1u.html).

Purple-leaved variety B100 accumulated more anthocyanin than YG1 in epidermal cells

In a previous study, we obtained a few foxtail millet varieties with different leaf colors, stems and panicles. We selected two representative varieties B100 (purple-leaved variety) and Yugu1 (YG1, green-leaved variety) to analyze the variation of anthocyanin distribution in foxtail millet. The morphological characteristics showed that B100 presented the visible purple color of the leaves at the seedling stage but the leaves of YG1 were green (Fig. 1A). Furthermore, the leaves of B100 and YG1 showed a different color at the maturation stages (Fig. 1B). The purple color of the purple leaves was deepest in the vein, which implied that anthocyanin was transported through the vascular bundle. Anatomical observations of the leaf revealed that in the leaves of B100, purple pigmentation was present in the adaxial and abaxial epidermal cells as well as adjacent mesophyll cells, whereas the purple pigmentation was not significantly present in the leaves of YG1 (Figs. 1C–1D).

The content of anthocyanin, not other pigments, was higher in B100 than YG1

We quantitatively analyzed the content of some pigments in B100 and YG1 to examine whether the visible purple was due to an increase in anthocyanin accumulation. The results showed that the anthocyanin content in B100 was significantly higher than YG1 at the seedling stage and the maturation stage (Fig. 2A), while the chlorophyll and carotenoid contents were higher in YG1 than in B100 (Figs. 2B–2C). Furthermore, all three pigments contents were higher at the maturation stage than at the seedling stage and the anthocyanin content was similar at the two stages in YG1.

B100 displayed enhanced stress resistance

Since anthocyanin synthesis and accumulation are closely related to the stress responses in plants, we explored the functions of anthocyanin in regulating stress resistance in foxtail millet. The total antioxidant capacity analysis was performed using the ferric-reducing antioxidant power (FRAP) method and we established the standard curve (Fig. S1A). B100 was shown to accumulate more anthocyanin and was better able to scavenge free radical and reductive ferrous ions than YG1 (Fig. 3A). Meanwhile, B100 could accumulate more soluble sugar in both the seedling and maturation stages compared to YG1 (Fig. 3B). The standard curve of the quantitative determination of soluble sugar is shown in Fig. S1B.

Differentially expressed genes were identified through the analysis of transcriptome characteristics in purple- and green-leaved foxtail millet varieties

We conducted RNA sequencing to elucidate the molecular mechanism underlying the phenotypic and functional differences between the purple- and green-leaved varieties. The total RNA of B100 and YG1 was extracted from the top second leaf at the maturation stage. The quality of the RNA was detected by agarose gel electrophoresis (Fig. S2A), and all RNA of corresponding samples that fulfilled the requirements were used for RNA-sequencing. A total of 144,677,546 clean reads were obtained, with a total base number of 43.4 G and GC contents were between 58.24% and 59.01%. The quality values of the sequencing data were statistically assessed for each sample of three biological replicates. The Q30 values were all above 92.84%, indicating that the transcriptome data were of high quality and could be further analyzed (Table S1).

The clean reads were sequentially compared with the reference genome (Setaria italica v2.2, Yugu1). Approximately 90.38%-91.50% of the clean reads in YG1 were sequentially compared to the reference genome. In the purple-leaved variety B100, 88.05%-89.97% of reads were compared to the reference genome. Among them, the proportion of clean reads with multiple locations was between 1.42% and 1.72%, and the number of reads compared to the positive and negative strands of the genome was similar. Therefore, the selected reference millet genome assembly may meet the needs for analysis (Table S2). We then obtained 1,184 significant differentially expressed genes (DEGs) of YG1 and B100, in which 598 genes were up-regulated and 588 genes were down-regulated in B100 compared with YG1 (Fig. S2B). A cluster analysis of these genes was conducted by hierarchical clustering. According to the different expression patterns, these DEGs were clustered into three groups. The genes in cluster I were down-regulated in B100, and the genes in cluster III were up-regulated in B100 (Fig. 4).

Subsequently, enrichment analysis of DEGs was performed based on Gene Ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis. GO analysis showed that DEGs were significantly enriched in the metabolic process, cell and cell part, binding, and catalytic activity (Fig. 5, Table S8). KEGG pathway analysis revealed that the DEGs were mainly enriched in the metabolic pathways. Among these DEGs, we identified some genes encoding synthetase of flavonoid biosynthesis and phenylpropane biosynthesis which may be related to anthocyanin metabolism (Fig. 6).

Nine structural genes related to anthocyanin biosynthesis may function in the accumulation of anthocyanin

We screened nine structural genes involved in anthocyanin biosynthesis from the results of the significant differentially expressed genes analysis that were up-regulated in B100 compared to YG1. We performed quantitative real-time PCR to verify the results of RNA sequencing in foxtail millet to determine the major role of these genes. As shown in Fig. 7, the relative expression levels of all nine structural genes in B100 were significantly higher than YG1 at the maturation stage. Among these up-regulated genes, DFR (Seita.5G237900), LDOX-2 (Seita.1G000700) and AT (Seita.9G002300) had higher expression levels in B100 than YG1 at both the seedling and maturation stages. With the development of foxtail millet, most of these structural genes displayed an increasing expression level. However, the plant accumulated more transcripts of PAL (Seita.1G240200) and UFGT (Seita.3G190000) at the seedling stage. All of the nine DEGs in the B100 purple-leaved variety were up-regulated, which was consistent with the results of RNA-sequencing. Moreover, the correlation between qRT-PCR and RNA-sequencing was 0.9625 (Fig. S3). These results showed that the sequencing results had high accuracy and reliability, and may be used for further investigations.