A comparative transcriptome analysis of a wild purple potato and its red mutant provides insight into the mechanism of anthocyanin transformation

In this study, a red mutant was obtained through in vitro regeneration of a wild purple potato. High-performance liquid chromatography and Mass spectrometry analysis revealed that pelargonidin-3-O-glucoside and petunidin-3-O-glucoside were main anthocyanins in the mutant and wild type tubers, respectively. In order to thoroughly understand the mechanism of anthocyanin transformation in two materials, a comparative transcriptome analysis of the mutant and wild type was carried out through high-throughput RNA sequencing, and 295 differentially expressed genes (DEGs) were obtained. Real-time qRT-PCR validation of DEGs was consistent with the transcriptome date. The DEGs mainly influenced biological and metabolic pathways, including phenylpropanoid biosynthesis and translation, and biosynthesis of flavone and flavonol. In anthocyanin biosynthetic pathway, the analysis of structural genes expressions showed that three genes, one encoding phenylalanine ammonia-lyase, one encoding 4-coumarate-CoA ligase and one encoding flavonoid 3′,5′-hydroxylasem were significantly down-regulated in the mutant; one gene encoding phenylalanine ammonia-lyase was significantly up-regulated. Moreover, the transcription factors, such as bZIP family, MYB family, LOB family, MADS family, zf-HD family and C2H2 family, were significantly regulated in anthocyanin transformation. Response proteins of hormone, such as gibberellin, abscisic acid and brassinosteroid, were also significantly regulated in anthocyanin transformation. The information contributes to discovering the candidate genes in anthocyanin transformation, which can serve as a comprehensive resource for molecular mechanism research of anthocyanin transformation in potatoes.


Introduction
Potato (Solanum tuberosum L.) has been the fifth largest crop in the world. It is nutritious with a reputation of "underground apple". Red and purple-fleshed potatoes additionally contain 2 addition, the used materials in previous researches have different genetic backgrounds, and it possibly influence the screening of differential genes and decrease the accuracy of comparative transcriptome analysis. It is urgent to use the homologous materials to study anthocyanin transformation mechanism for enriching the metabolic pathway of anthocyanin.
In this study, we fortunately obtained a red mutant SD140 from a purple wild potato SD92. Both materials were homologous, which facilitated our comparative transcriptome analysis. In order to study the mechanism of anthocyanin transformation, we identified the main types of anthocyanins and carried out a comparative transcriptome analysis of SD140 and SD92.

Plant sources
SD92, commonly known as Hei Jingang, is tetraploid with purple skin and purple flesh tuber. A red tuber mutant SD140 with red skin and flesh was obtained through in vitro regeneration of SD92 tuber slices (Fig 1). Two materials were planted in a greenhouse at 20 ± 2˚C under a 16 h light/8 h dark cycle at the Institute of Vegetables and Flowers, Shandong Academy of Agricultural Sciences at the normal sowing time (August 10, 2016) in Jinan. The fresh tubers (diameter 4-5cm) were harvested, and cleaned with sterilized water. The tubers were immediately used for pigment extraction and two biological replicates were used. The tubers for mRNA extraction were frozen in liquid nitrogen and stored at -80˚C, three biological replicates were used.

Pigment extraction and purification
Potato tubers (250 g) were washed with distilled water and smashed in 500 mL of a mixture containing 50% alcohol and 1.6% citric acid. The resulting mixture was filtrated twice and centrifuged at 4000 rpm. The supernatant was transferred into a column filled with the pretreated macroporous resin AB-8 for 6 h, and the pigments absorbed in AB-8 were eluted with 95% alcohol. The eluate was concentrated with rotary evaporator at 30˚C and dried in vacuo for 12 h. The laboratory protocol was deposited in protocols.io website (http://dx.doi.org/10.17504/ protocols.io.ma4c2gw).

Analytical high-performance liquid chromatography
The pigment extractions were analyzed by the ultra-performance liquid chromatography (UPLC) system (Waters, Milford, Massachusetts, USA) equipped with Waters acquity PDA detector and a BEH C18 column (2.1×100mm, 1.7μm particle size). Acetonitrile (mobile phase A) and 2% methanoic acid (mobile phase B) were used as the eluent at a flow rate of 0.3 mL/ min. The gradient conditions for sample elution were as follows: 5% A 95% B, 18 min; 20% A 80% B, 20 min; 100% A, 22 min; 100% A, 22.1 min; 5% A 95% B, 25 min; 5% A 95% B. The anthocyanins were detected at 530 nm using PAD. The column temperature was kept at 45˚C. The volume for each injection was 1 uL.

Mass spectrometry
Mass spectrometry (MS) analysis was performed on a quadrupole time-of-flight (Q-TOF) mass spectrometry (WATERS MALDI SYNAPT) system. The anthocyanins were analyzed in Positive-ion (PI) mode. The optimized conditions were: the capillary and sampling cone voltage of 3.5kV and 30 V, respectively; the source block temperature of 100˚C; desolvation temperature, 400˚C; desolvation gas flow, 500 L/h; cone gas flow, 50 L/h; detector voltage, 1800 V. Data were gathered in the v mode, with a scan accumulation time of 1 s and the TOF data being collected between m/z 50 and 1,500. The MS-MS experiments were performed using a collision energy of 6 eV, which was optimized for each compound.

cDNA library construction and high-throughput sequencing
According to the manufacturer's protocol, total RNA of the sample was extracted using modified Trizol reagent (Ambion, USA). And then RNA was digested with 10U DNaseI (Thermo, USA) for 1 h at 37˚C. RNA quality and purity were detected by Agilent 2100 and NanoDrop (USA). The poly (A) mRNA was enriched using Oligo (dT) magnetic beads from digested RNA and was interrupted into short fragments. Using the short mRNA fragments as templates, First-strand cDNA were synthesized using random hexamer primers. The second strand was synthesized by buffer, dNTPs, Rnase H, and DNA polymerase I. The doublestranded cDNA was obtained by PCR with illumina primers. The two-stranded cDNA was purified with magnetic beads. The cohesive ends were repaired and the base "A" was added to the 3 0 end. The fragments with suitable size were selected and were enriched by PCR amplication. The cDNA library was qualified and quantified with the Agilent 2100 Bioanalyzer (USA) and ABI StepOnePlus Real-Time PCR System (USA). The qualified library was sequenced by use of Illumina HiSeqX-Ten (USA). In total, six cDNA libraries were constructed in this study. All the raw reads were acquired by RNA-seq. In order to obtain clean reads, we filtered low quality reads, adaptor sequences and the reads in which unknown base N content was too high.
We deposited all the sequence data of the six cDNA libraries into NCBI Sequence Read Archive (http://www.ncbi.nlm.nih.gov/sra/), and the accession number is SRP125987.

Differentially expressed genes identification
All clean reads were mapped to the reference transcripts of Solanum tuberosum L. (PGSC_DM_v3.4_transcript-update, http://solanaceae.plantbiology.msu.edu/data) using the Bowtie2 (v2.2.5). Statistics and bioinformatics analyses were carried out on the number of clean reads mapped to reference genes and genome, sequencing saturation, and random sequencing distribution. And then gene expression levels were calculated based on the number of reads mapped to the reference sequences. To identify the differentially expressed genes (DEGs) in the two Samples, DEseq2 algorithm was used. For each gene, the average expression level in the SD92 group (Control-avg) and in the SD140 group (Treat-avg) were calculated, respectively. DEGs were screened according to the following criteria: |log2(Treatment-avg /Control-avg)| ! 1.00 and adjusted p-value 0.05.

Gene Ontology (GO) enrichment and KEGG pathway enrichment analysis of DEGs
The DEGs were mapped to the corresponding GO terms and KEGG Orthology (KO) pathways. We used the R software phyper perform GO functional enrichment and KEGG functional enrichment. The p-value formula (https://en.wikipedia.org/wiki/Hypergeometric_ distribution) in hypergeometric test is: Then false discovery rate (FDR) for each p-value was calculated, and the terms which FDR is not larger than 0.01 were defined as a significant enrichment.

Transcription factor (TF) coding ability prediction in DEGs
We used getorf (EMBOSS: 6.5.7.0) to detect the ORF of each DEG, then used hmmsearch (v3.0) to align the ORF to the TF protein domain (data from PlantfDB), and finally identified TF according to the regulation based on the TF family characteristics.

Real-time qRT-PCR validation of the genes expression difference
To validate the transcriptome sequencing results, the transcription level of DEGs in the transcriptome was verified by real-time qRT-PCR. RNA samples used for qRT-PCR were same as those used in high-throughput sequencing experiments. Gene specific primers were listed in S1 Table. 25 μL of 2× UltraSYBR Mixture (CWBIO, CHINA), 1 μL cDNA, and 50 μL of total reaction volume were applied to all the reactions following the manufacturer's method. The qPCR reaction was performed using an iCycler iQ (Bio-Rad, Hercules, CA, USA). The qPCR conditions were 37 cycles of 15 s at 95˚C and 1 min at 60˚C, followed by 65˚C to 95˚C melting curve detection. Three biological replicates were performed and 18S RNA was used as a reference gene.

Pigment analysis
Fig 1 showed the potato tubers appearance (A) and the cross section (B) of SD140 and SD92, respectively. SD140 was obtained by the regeneration of tuber slices of SD92. After regeneration, both skin and flesh of SD92 changed from purple to red. Pure anthocyanins were obtained through pigment extraction, concentration and purification from the fresh tubers. In order to isolate the pure anthocyanins, HPLC was performed. For the anthocyanins extracted from SD140, three major independent peaks were obtained, and their retention times were located at 9.99, 11.70 and 13.29 minutes, respectively (Fig 2A). Similarly, six major peaks were acquired from SD92 extractions, and their retention times were located at 7.13, 11.04, 11.75, 12.84, 13.78, 15.79 minutes, respectively ( Fig 2B). All these peaks were further identified by the quadrupole time-of-flight (Q-TOF) mass spectrometry, and the results were presented in Table 1. As for SD92, the MS analysis of the largest peak at the retention time of 11.70 minutes displayed an [M] + ion at m/z 933.3 and a major fragmentation in MS 2 at m/z 317.1. The peaks at 9.99 and 11.70 min revealed a major fragmentation in MS 2 at m/z 317.1 (S1 Fig). The results were confirmed on the basis of HPLC retention time, elution order, spectroscopic data, MS fragmentation pattern, MS 2 , and previous findings [26][27][28][29][30]. The MS 2 fragmentation of 317.1 at m/z corresponded to petunidin ( Table 1). The results confirmed that the main anthocyanin in  (Table 1).  Transcriptome analysis RNA sequencing and mapping of the sequence reads. Illumina sequencing technology was used to perform high-throughput RNA sequencing of SD92 and SD140, respectively. After removing the adaptor sequences and low-quality reads, 44.15-44.95 Mb clean reads were acquired from each sample. The clean reads were mapped to the genome. The mapped ratio of each sample reached 55.63%-61.49% (S2 Table). After transcripts were reconstructed, we acquired 37460 transcripts. 15232 transcripts were novel, of which 3641 were novel genes ( Table 2).
To identify each gene expression in SD92 and SD140 tubers, the clean reads were mapped to the complete reference sequences, which included the reference gene sequence and novel gene. The results demonstrated that 32.39-39.51Mb clean reads were mapped in each sample. These values accounted for 72.29%-75.42% of each library ( Table 2). Based on the mapped results, we counted the read coverage and the reads distribution of each transcript. The random reads distribution on the transcripts indicated that the read positions were evenly distributed in each gene (S3 Fig). Reads coverage on transcripts analysis showed that the average percentage of each gene covered by reads was approximately 50% (S4 Fig). Therefore, the sequencing data should accurately reflect gene expression and can be used for differential gene expression analysis.
Global gene expression pattern analysis. In SD92 and SD140, there were 22679 and 22667 genes expressed in all three biological replicates, respectively (S5 Fig). To explore the molecular mechanism involved in anthocyanin transformation, we analyzed the DEGs in three biological replicates. 267 DEGs were obtained between SD92 and SD140. Among these DEGs, 177 (60.00%) were up-regulated and 118 (40.00%) were down-regulated (Fig 3) in SD140 compared to SD92.
To facilitate the global analysis of gene expression, GO analysis of DEGs between SD92 and SD140 was performed. The results revealed 98 DEGs were classified. The DEGs were classified into the three main GO categories: biological process, cellular component and molecular function (S6 Fig). For further review of the GO classification of the DEGs, we also categorized each GO term into its sub-categories. The results indicated the sub-categories significantly enriched in the biological process category primarily contained biological regulation, localization, and response to stimulus. The DEGs in this category were associated with the regulation of catalytic activity, transcription, and macromolecule metabolic process. Cell part, cell, organelle, organelle part, membrane and membrane part were the main sub-categories belonging to the cellular component category. For the molecular function category, the main sub-categories were binding, catalytic activity, transporter activity and molecular function regulator. Catalytic activity was involved in putative RING-H2 finger protein, disulfide-isomerase protein and rasrelated RABA1f protein. Proteins in the binding sub-category included beta-galactosidase, calreticulin, luminal-binding protein, kinesin-related protein, probably inactive leucine-rich repeat receptor and heat shock 70 kDa protein.
KEGG analysis was performed to obtain the difference in the metabolic pathway in SD92 and SD140 (Fig 4). Among the DEGs, 220 genes corresponded to KO codes, including 126 significantly up-regulated genes and 94 significantly down-regulated genes in SD140. To review the function of the DEGs, we further clustered each KO term of the DEGs. The results showed that the DEGs were mostly enriched in transport and catabolism, translation, biosynthesis of other secondary metabolites and transcription. These clusters included the following pathways: protein processing in endoplasmic reticulum (13 genes), plant hormone signal transduction (12 genes), RNA transport (8 genes), flavonoid biosynthesis (3 genes), flavone and flavonol biosynthesis (2 genes), phenylpropanoid biosynthesis (7 genes) and spliceosome (7 gens).
Gene expression changes of anthocyanin biosynthesis in potato. Considering that the main type of anthocyanins in SD92 was different from SD140, and that the results of KEGG enrichment of DEGs showed differences in flavonoid metabolic pathways between SD92 and SD140, we selected 41 genes to investigate the pathway of anthocyanin biosynthesis. Of the selected genes, 8 were new genes. The criterion for selection was that the gene expression level was not less than 20 reads at least in one material (Fig 5 and S3 Table). The 41 genes mainly encoded the following enzymes: phenylalanine ammonia-lyase (PAL), 4-coumarate-CoA ligase (4CL), chalcone synthase (CHS), flavonoid 3 0 ,5 0 -hydroxylase (F3 0 5 0 H), flavonoid 3 0hydroxylase (F3 0 H), dihydroflavonol 4-reductase (DFR), leucoanthocyanidin dioxygenase (LDOX), anthocyanidin 3-O-glucosyltransferase (3GT). These enzymes had two characteristics. One was that these enzymes possessed isozymes. Phenylalanine ammonia-lyase included two isozymes: PAL and PAL1. Chalcone synthase mainly included two isozymes: CHS1B, CHS2. 4-Coumarate-CoA ligase included 4CL2 and 4CL. The other was that these enzymes belonged to a multi-gene family. PAL and 4CL families comprised 14 and 9 genes, respectively. The differences of 41 gene expression levels in anthocyanin biosynthetic pathway were investigated according to p-value. If p-value 0.05, the gene was significantly regulated. The result showed the expression of four genes significantly changed, and other genes almost didn't significantly changed. Four changed genes included two genes encoding PAL, one gene encoding 4CL, one gene encoding F3 0 5 0 H. To investigate the change range of these significantly changed genes, we calculated the rate of change by using the following formula: (Treatment-avg-Control-avg)/ Control-avg ×100% (S3 Table). The results indicated that in SD140 one gene encoding PAL was up-regulated by 144%, another gene encoding PAL was downregulated by 42.58%; one gene encoding 4CL was down-regulated by 39.83%; one gene encoding F3 0 5 0 H was down-regulated by 37.66%.

Analysis of transcription factors in DEGs.
Transcription factor is a kind of protein binding the specific nucleotide sequence of a gene and regulates the transcription of the related gene. The GO analysis indicated that there existed differences in transcription and DNA binding between SD92 and SD140. To determine the molecular mechanism of anthocyanin transformation, we examined the differential expression of transcription factors between SD92 and SD140. The results revealed the transcription factors were encoded by 18 DEGs (Table 3). These transcription factors belonged to 12 families. Alfin family, bZIP family, C3H family, FHA family, GRAS family, LOB family, MADS family were down-regulated in SD140 and   were mainly involved in histone acetyltransferase, transcription factor HBP-1, zinc finger protein ZAT11, zinc finger CCCH and uncharacterized protein. S1Fa family, Tify family, zf-HD family were up-regulated in SD140 and were mainly involved in vacuolar amino acid transporter YPQ2, plastidic glucose transporter, zinc-finger homeodomain protein 9 and protein RADIALIS. The AP2-EREBP family contained AP2 ethylene-responsive transcription factor ANT and ethylene-responsive transcription factor 4-like. ANT was up-regulated and 4-like was down-regulated. The MYB family comprised myb-related protein MYBAS1 and protein RADIALIS 3. MYBAS1 was up-regulated, and RADIALIS 3 was down-regulated in SD140.

Analysis of plant hormone signal transduction in DEGs.
Hormones are growth regulators, which affect plant germination, rooting, flowering, fruiting, sex determination, dormancy and shedding. The KEGG analysis indicated that there were differences in plant hormone signal transduction between SD92 and SD140. In order to deeply analyze the molecular mechanism of anthocyanin transformation, the genes expression differences in hormone signal transduction between SD92 and SD140 were investigated. The results indicated that 12 genes expression in hormone signal transduction exhibited significant change ( Table 4). The hormone signal transduction included GA, ABA, IAA, BR, Ethylene and Salicylic acid signal transduction. As for SD140, in GA signal transduction, one gene encoding scarecrow-like protein 8 was significantly down-regulated; in ABA signal transduction, one gene encoding abscisic acid receptor PYL4 was down-regulated; in BR signal transduction, one gene encoding D-type cyclin family 3 subgroup 3 and one gene encoding inactive leucine-rich repeat receptor-like protein kinase IMK2 were up-regulated, one gene encoding BRI1 kinase inhibitor 1 and one gene encoding LRR receptor-like serine/ threonine-protein kinase At2g24230 were down-regulated.
Validation of RNA-Seq-based gene expression by qRT-PCR. To validate the RNA-Seq results, a qRT-PCR analysis was performed on 10 genes. These genes included eight randomly selected DEGs and two genes in the anthocyanin metabolic pathway. Eight randomly selected DEGs were involved in DNA-directed RNA polymerase, ubiquitin-protein ligase, ferredoxin, translation initiation factor, xyloglucosyl transferase and ABA 8 0 -hydroxylase. Two genes in the anthocyanin biosynthetic pathway encoded F3 0 5 0 H and LDOX. The qRT-PCR result indicated that all 10 genes expression pattern was accordant with the RNA-Seq date (Fig 6).

Effect of structural genes expression differences on anthocyanin transformation in anthocyanin biosynthesis pathway
Anthocyanin biosynthesis pathway has been elucidated in many plant species. Different species have different structural genes. Many structural genes have been cloned from model species, such as maize and petunia [17]. In this study, the structural genes PAL, ACC, 4CL, C4H, CHS, CHI, F3 0 H, F3 0 5 0 H, ANS, DFR and 3GT were detected by analyzing the transcriptome of SD92 and SD140. The results indicated that the expressions levels of F3 0 H were both low in two materials. F3 0 H is responsible for biosynthesis of peonidin [14]. So the predicated peonidin levels were low in the two materials, which were consistent with our HPLC-MS analysis. These results further confirmed that F3 0 H was responsible for the biosynthesis of peonidin in potato. Additionally, the differential expression analysis of structural genes in anthocyanin biosynthesis revealed that one PAL, one 4CL and one F3 0 5 0 H were significantly down-regulated in SD140; one PAL was significantly up-regulated. PAL and 4CL belong to anthocyanin biosynthetic upstream genes, and they are involved in multiple metabolic pathways [31,32]. These genes expression changes may not cause changes in anthocyanin biosynthesis. We paid more attention to the gene change in downstream. F3 0 5 0 H belongs to anthocyanin biosynthetic downstream gene. F3 0 5 0 H plays a key role in regulating different types of anthocyanin biosynthesis [14]. We found that it was significantly down-regulated in anthocyanin transformation.

Effect of F3 0 5 0 H on anthocyanin transformation
In this work, petunidin-3-O-glucoside and pelairgonidin-3-O-glucoside were found to be responsible for purple tuber and red tuber, respectively, and these findings were consistent with previous results [9,15]. F3 0 H and F3 0 5 0 H, belonging to the cytochrome P-450 total branch, are key enzymes controlling different branches of anthocyanin biosynthesis in some species [14]. When both F3 0 H and F3 0 5 0 H are present, F3 0 5 0 H is predominant. In other words, F3 0 5 0 H not only promotes anthocyanins synthesis pathway to 3 0 5 0 -hydroxylated anthocyanins but also stops the biosynthetic pathway to 3 0 -hydroxylated anthocyanins [33,34]. Besides, F3 0 5 0 H expression was much higher and F3 0 H expression was much lower, so F3 0 5 0 H had a main effect on the color of potato tuber. F3 0 5 0 H controls the biosynthesis of petunidin-based anthocyanins [17]. Silenced F3 0 5 0 H, inactivation and degeneration of F3 0 5 0 H can decrease petunidin-based anthocyanins [35,36]. The decreased petunidin in SD140 was caused by the decreased expression of F3 0 5 0 H, which was verified by qRT-PCR. Moreover, the main type of anthocyanin in the purple tubers and the red tubers were petunidin and pelargonidin, respectively. Therefore, the transformation from petunidin to pelargonidin in tuber was related to the decrease of F3 0 5 0 H expression level and dihydrokaempferol directly synthesize pelargonidin [14]. So it was possible that the transcript of F3 0 5 0 H was influenced.

Effect of transcription factors on anthocyanin transformation
Transcription factors are adaptor molecules. It can be used to detect regulatory sequences in the DNA and target the assembly of protein complexes controlling gene expression [37]. Transcription factors play important roles on controlling plant development and differentiation. In this study, the transcription factor families, including AP2-EREBP family, C2H2 family, MYB family, S1Fa family, Tify family, GRAS family, Alfin family, C3H family, FHA family, LOB family, MADS family and zf-HD family, significantly changed during anthocyanin transformation.
In plants, R2R3MYB family, a Myb-type transcription factor, is closely related to flavonoid metabolic pathways. The R2R3MYB family includes a lot of subfamilies in different species. The subfamilies, such as Stan2, VvMYBPA1, VvMYB5b and VvMYB5a, can regulate anthocyanin biosynthesis by activating structural genes. All the activated genes include F3 0 5 0 H [11,38,39]. There are also inhibitory factors in the MYB family, such as SlTRY (CPC) and MYB182. They inhibit the anthocyanin biosynthesis [40,41]. In our study, the expressions of the genes encoding the MYB family members were significantly different in potatoes with different types of anthocyanins. It was possible that these MYB family members regulated anthocyanin transformation from petunidin to pelargonidin by regulating F3 0 5 0 H, which switched from 3 0hydroxylation to 3 0 5 0 -hydroxylation anthocyanin synthesis pathway. In MYB family, MYBAS1 was up-regulated, which suggested it promoted the transformation from petunidin to pelargonidin by repressing F3 0 5 0 H expression. RADIALIS 3 was down-regulated, which suggested it inhibited the transformation from petunidin to pelargonidin.
Besides, bZIP family also played an important role in anthocyanin transformation. bZIP, an alkaline leucine zipper transcription factor, comprised a basic part and a leucine zipper part [42]. bZIP transcription factors G/HBF-1, LONG HYPOCOTYL 5 (HY5) and HY5 HOMO LOG (HYH) can increase the structural genes expression of CHS, CHI, F3H and DFR and promote anthocyanin biosynthesis by binding the promoter of the structural genes [43,44,45]. bZIP transcription factors also improve anthocyanin biosynthesis by binding the R2R3MYB transcription factor PAP1 promoter region [46,47]. R2R3MYB can regulate the expression of F3 0 5 0 H, which regulates transformation from petunidin to pelargonidin [36]. Therefore, bZIP can regulate anthocyanin transformation from petunidin to pelargonidin by F3 0 5 0 H. In this study, a gene encoding bZIP transcription factor HBP-1b (c38) -like isoform X1, X3, X6, X2 was down-regulated. It suggested that the HBP-1b (c38) inhibited the transformation from petunidin to pelargonidin by MYB.
MADS, a transcription factor, consists of the MADS-box, the I, the K and the C regions [48]. MADS family is closely related to anthocyanin biosynthesis [49][50][51]. MADS transcription factors increase anthocyanin biosynthesis by activating the genes encoding CHS, CHI, F3H, DFR, ANS and UFTG [49,52]. MADS transcription factor VmTDR4 plays an important role in accumulation of anthocyanins through the up-regulation of R2R3MYB family [53]. R2R3MYB can regulate anthocyanin transformation through F3 0 5 0 H [36]. So MADS can regulate anthocyanin transformation through F3 0 5 0 H. In this study, a gene encoding MADS was significantly down-regulated in SD140. It was possible that the MADS suppressed the transformation from petunidin to pelargonidin through regulating F3 0 5 0 H.
In this study, the transcription factors were significantly regulated, such as AP2-EREBP family, C2H2 family, S1Fa family, Tify family, GRAS family, Alfin family, C3H family, FHA family, LOB family, and zf-HD family. The functions of these transcription factors were worth further research.

Effect of plant hormone on anthocyanin transformation
Plant hormone can influence anthocyanin biosynthesis. GA signal is involved in anthocyanin biosynthesis pathway. GA can decrease anthocyanin accumulation under low temperature or phosphate starvation [54,55]. The depression effect under low temperature was determined by two bZIP transcription factors (HY5 and HYH) [54]. In this study, one DEG encoding scarecrow-like protein 8 in GA signal pathway, and one DEG encoding bZIP transcription factor HBP-1b(c38) were both down-regulated. The results suggested that scarecrow-like protein 8 inhibited the transformation through HBP-1b(c38). ABA can independently improve anthocyanin accumulation and also has a synergic effect with sugar [56,57]. ABA regulates anthocyanin biosynthesis via MYB, bHLH and bZIP [58][59][60]. In this work, one DEG with the function of ABA receptor PYL4 in ABA signal pathway was down-regulated. Additonallly, the aforementioned analysis indicated that MYB and bZIP regulated anthocyanin transformation. Therefore, it could be concluded that ABA receptor PYL4 inhibited the transformation from petunidin to pelargonidin by MYB and bZIP. BR also affects JA-induced anthocyanin accumulation by the WD-repeat/MYB/bHLH transcriptional complexes [61]. In current work, D-type cyclin family 3 and inactive leucine-rich repeat receptor-like protein kinase IMK2 were up-regulated, BRI1 kinase inhibitor 1 and LRR receptor serine/threonine-protein kinase At2g24230 were down-regulated. Besides, the above-mentioned findings confirmed that MYB and bZIP regulated anthocyanin transformation. As a result, it could be infered that with the aid of the MYB and bZIP, the two up-regulated BR response proteins inhibited anthocyanin transformation and the two down-regulated BR response proteins promoted anthocyanin transformation.
Moreover, the response proteins in signal transduction of IAA and ethylene were significantly regulated in anthocyanin transformation in this work, The functions of these hormone response proteins were worth further research.

Conclusion
In this study, both pigments and the comparative transcriptome analysis of purple potato tuber and its red mutant were carried out. It revealed the decreased expression level of F3 0 5 0 H was mainly responsible for anthocyanin transformation from petunidin to pelargonidin in potato tubers. Transcription factors bZIP family, MSD family regulated the transformation through MYB family, which regulated F3 0 5 0 H expression. Response proteins of ABA, GA, BR also affected the transformation through transcription factors.The structural genes, transcription factors and their related factors formed a complex that affected the transformation from petunidin to pelargonidin in potato tuber (Fig 7). The transcriptome analysis provided valuable information for anthocyanin transformation.