Comparative Transcriptome Analysis of Different Actinidia arguta Fruit Parts Reveals Difference of Light Response during Fruit Coloration

Simple Summary Kiwifruit (A. arguta) color is one of the most important quality characters. Exploring the coloration mechanism is significant for the genetic improvement of color quality and the breeding of new germplasms. As a critical environmental factor, light plays a key role in fruit coloration. However, the effecting mechanism of light on A. arguta coloration remains unclear. In current research, different A. arguta parts with different treatments were performed high throughput RNA sequencing, based on which candidate genes and corresponding annotations were obtained. Finally, AaMYB308like was screened as an R2R3-MYB typed TF involved in light-inducible fruit coloration through the result analysis of bioinformatics and molecular biology experiments. Our study provides insights into the photoreponse mechanisms in A. arguta coloration. Abstract Kiwifruit coloration is an important agronomic trait used to determine fruit quality, and light plays a vital role in the coloration process. The effect of light on fruit coloration has been studied in many species, but differences in the photoresponse of different fruit parts during fruit coloration is unclear in kiwifruit (Actinidia arguta). In this study, peel and core with bagging and non-bagging treatment at two stages were selected to perform high throughput RNA sequencing. A total of 100,417 unigenes (25,186 unigenes with length beyond 1000 bp) were obtained, of which 37,519 unigenes were annotated in functional databases. GO and KEGG enrichment results showed that ‘plant hormone signal transduction’ and ‘carbon metabolism’ were the key pathways in peel and core coloration, respectively. A total of 27 MYB-related TFs (transcription factors) were differentially expressed in peel and core. An R2R3-MYB typed TF, AaMYB308like, possibly served as a candidate objective, which played a vital role in light-inducible fruit coloration based on bioinformatics analysis. Transient overexpression of AaMYB308like suggested overexpression of AaMYB308like elevated transcription level of NtCHI in Nicotiana tabacum leaves. Integration of all these results imply that AaMYB308like might be served as a light-responsive transcription factor to regulate anthocyanin biosynthesis in A. arguta. Moreover, our study provided important insights into photoreponse mechanisms in A. arguta coloration.


Introduction
Kiwifruit belongs to Actinidiaceae, genus Actinidia, which possess rich germplasm resources, including 54 species and 21 varieties [1]. Some kiwifruit cultivars are red due to their high concentration of anthocyanins, metabolites with high antioxidant properties that contribute to improved vascular elasticity and protect against liver injury [2]. However, most kiwifruit cultivars are green, so there is a need to develop new red varieties given their potential benefit to human health. The whole-red Actinidia arguta (Sieb. et Zucc.) Planch. et Miq. is one of the new cultivated species in recent years. Due to the characteristics in each biological replicate. Fruits at 30 days after full bloom (DAFB) were bagged using two-layer light-impermeable bags to ensure that the fruit in bags receives no light at all ( Figure 1A). Untreated fruits receiving normal light were sampled at the same stages as the control. Sampling was carried out at two stages. The first sampling was done at 70 DAFB (S1) when the fruits were still green, and the second sampling was done at 110 DAFB (S2) when the fruit color was obviously different between bagged and unbagged fruits ( Figure 1B). The peel and core were separated from fruits using a lab-used blade and were immediately frozen in liquid nitrogen and stored at −80 °C for subsequent use. Figure 1. Sample treatment and preparation for RNA-seq. (A) Fruits were bagged using two-layer light-impermeable paper bags at one month after full bloom. (B) Peel and core were collected from bagged and unbagged fruit making a total of eight typed samples used for RNA-seq. (C) The anthocyanin content of these eight typed samples. Note: S1-fruit stage at 70 days after full bloom, S2-fruit stage at 110 days after full bloom, TP70-peel from bagged fruit at S1, TX70-core from bagged fruit at S1, TP110-peel from bagged fruit at S2, TX110-core from bagged fruit at S2, WP70-peel from unbagged fruit at S1, WX70-core from unbagged fruit at S1, WP110-peel from unbagged fruit at S2, WX110-core from unbagged fruit at S2. Data were analyzed with Student's t-test (** p < 0.01).
The measurement of anthocyanin content was conducted using previous methods [2], and all operations were performed for three replicates.

RNA Extraction, cDNA Synthesis, Library Construction, and Sequencing
RNA was extracted from sample using RNAprep Pure Plant Kit (DP441 TIANGEN, Beijing, China), after which the purity and concentration of RNA were tested by a NanoDrop 2000c spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA), and the integrity of RNA was tested using agarose gel electrophoresis, ensuring high qualification for the next procedure. The first cDNA strand was synthesized using random hexamer primer and M-MuLV Reverse Transcriptase (NEB, Rowley, NE, USA), and subsequently the second strand cDNA was obtained by DNA Polymerase I and RNase H Figure 1. Sample treatment and preparation for RNA-seq. (A) Fruits were bagged using two-layer light-impermeable paper bags at one month after full bloom. (B) Peel and core were collected from bagged and unbagged fruit making a total of eight typed samples used for RNA-seq. (C) The anthocyanin content of these eight typed samples. Note: S1-fruit stage at 70 days after full bloom, S2-fruit stage at 110 days after full bloom, TP70-peel from bagged fruit at S1, TX70-core from bagged fruit at S1, TP110-peel from bagged fruit at S2, TX110-core from bagged fruit at S2, WP70-peel from unbagged fruit at S1, WX70-core from unbagged fruit at S1, WP110-peel from unbagged fruit at S2, WX110-core from unbagged fruit at S2. Data were analyzed with Student's t-test (** p < 0.01).
The measurement of anthocyanin content was conducted using previous methods [2], and all operations were performed for three replicates.

RNA Extraction, cDNA Synthesis, Library Construction, and Sequencing
RNA was extracted from sample using RNAprep Pure Plant Kit (DP441 TIANGEN, Beijing, China), after which the purity and concentration of RNA were tested by a Nan-oDrop 2000c spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA), and the integrity of RNA was tested using agarose gel electrophoresis, ensuring high qualification for the next procedure. The first cDNA strand was synthesized using random hexamer primer and M-MuLV Reverse Transcriptase (NEB, Rowley, NE, USA), and subsequently the second strand cDNA was obtained by DNA Polymerase I and RNase H (Thermo Fisher Scientific, Waltham, MA, USA). To get cDNA fragments of preferentially 240 bp in length, the library fragments were purified with AMPure XP system (Beckman Coulter, Beverly, MA, USA), and then 3 µL USER Enzyme (NEB, Rowley, NE, USA) was used with size-selected, adaptor-ligated cDNA at 37 • C for 15 min followed by 5 min at 95 • C before PCR. Afterwards, PCR was performed with Phusion High-Fidelity DNA polymerase (Roche, Basel, Switzerland), Universal PCR primers and Index (X) Primer. Finally, the PCR products were purified by AMPure XP system (A63880, Beckman Coulter, CA, USA) and library quality was assessed on the Agilent Bioanalyzer 2100 system. The prepared library was sequenced on an Illumina platform and paired-end reads were generated. There were three biological replicates performed for the RNA sequencing during whole RNA-seq analysis. In addition, all raw data has already been submitted to the public database NCBI Sequence Read Archive (SRA accession: PRJNA680442).

Quality Control and Gene Functional Annotation
Raw reads of fastq format were firstly processed through fastp software 0.21.0 (HaploX Biotechnology, Shenzhen, China) with specific parameters setting (-q 10 -u 50 -y -g -Y 10 -e 20 -l 100 -b 150 -B 150), during which the clean reads were obtained by removing adapter, poly-N, and low-quality reads. Meanwhile, Q20, Q30, GC-content, and sequence duplication level of clean reads were calculated to ensure downstream analyses based on high quality clean reads. The left files (read one file) from all libraries were pooled into one big left.fq file, and right files (read two files) into one big right.fq file. Transcriptome assembly was accomplished based on the left.fq and right.fq using Trinity with min_kmer_cov set to 2 by default and all other parameters set to default. Six databases, Nr (NCBI non-redundant protein sequences), Nt (NCBI non-redundant nucleotide sequences), Pfam (Protein family), KOG/COG (Clusters of Orthologous Group of proteins), Swiss-Prot (A manually annotated and reviewed protein sequence database), KO (KEGG Ortholog database), and GO (Gene Ontology), were used for gene functional annotation.

Quantification of Gene Expression and Differential Expression Analysis
FPKM (fragments per kilobase of transcript per million fragments mapped) was selected for estimation of gene expression quantification. DEseq was used for differential expression analysis [16,17]. DEseq provide statistical routines for determining differential expression in gene expression data using a model based on the negative binomial distribution. The resulting p value were adjusted using the Benjamini and Hochberg's approach for controlling the false discovery rate. Genes with an adjusted p-value < 0.01 found by DEseq were assigned as differentially expressed.

GO and KEGG Enrichment Analysis
GO (gene ontology) enrichment analysis of the DEGs (differentially expressed genes) was implemented by the GOseq R packages based Wallenius non-central hyper-genometric distribution [18], which can adjust for gene length bias in DEGs. KEGG (Kyoto Encyclopedia of Genes and Genomes) is a database resource for understanding high-level functions and utilities of the biological system, such as the cell, the organism, and the ecosystem, from molecular-level information, especially large-scale molecular datasets generated by genome sequencing and other high throughput experimental technologies (http://www.genome.jp/kegg/, accessed on 1 January 2000). KOBAS software (3.0 version, Peking University, Beijing, China) was used for testing the statistical enrichment of differential expression genes in KEGG pathways [19].

Vector Construction and Transient Overexpression of AaMYB308like in N. tabacum
The CDS of AaMYB308like was amplified from A. arguta by specific primers 5 -GCTCTAGAATGGGGCGATCACCATGTTG-3 (forward) and 5 -CGGGATCCTCTACACA ACCCTTGATTGC-3 (reverse) containing Xba I and BamH I restriction enzymatic sites. The PCR product was recombined with the plant binary expression vector pBI121 to form CaMV 35S:AaMYB308like-GFP. Empty pBI121 vector only with the GFP gene (35S:GFP) was used as control. Two constructs were introduced into A. tumefaciens strain EHA105 using a freeze-thaw method. A. tumefaciens strains were kept at 28 • C in LB medium with kanamycin antibiotic, re-suspended in infiltration buffer containing 10 mM MgCl 2 , 10 mM MES and 200 µM acetosyringone to OD 600 of 0.6-1.0, and placed at room temperature for 2 h before infiltration.
The overexpression vector AaMYB308like-pBI121 constructed above was used for A. tumefaciens-mediated infiltration in tobacco leaves using a 1 mL needleless syringe. The empty vector pBI121-injected and non-treated leaves were both served as the control. Phenotypic observation and gene expression were assayed 7 days after transformation. Experimental operation was carried out on three biological repeats.

qRT-PCR Analysis
A total of 20 µL PCR mixture was used, including 5 µL ddH 2 O, 10 µL SYBR Green I master mix (Roche, Basel, Switzerland), 1 µL forward primer, 1 µL reverse primer, and 3 µL cDNA template. Real time PCR reaction was run on the LightCycler ® 480 system with a 96-well plate accompanied by the PCR procedure as follows: 95 • C for 5 min, followed by 45 cycles of 10 s at 95 • C, 20 s at 60 • C, and 20 s at 72 • C. The Nicotiana tabacum alpha-tubulin (NtTubA1 gene) was selected as control for the transient overexpression qRT-PCR. The relative expression level of genes was estimated by the 2 −∆∆Ct method [20].

Changes in Phenotype of Peel and Core between Bagging and Non-Bagging Treatment
The bagged and unbagged fruits were sampled at 70 DAFB and 110 DAFB. A total of eight types of samples from peel and core of fruits, namely TP70 (peel from bagged fruit at S1), TX70 (core from bagged fruit at S1), TP110 (peel from bagged fruit at S2), TX110 (core from bagged fruit at S2), WP70 (peel from unbagged fruit at S1), WX70 (core from unbagged fruit at S1), WP110 (peel from unbagged fruit at S2), and WX110 (core from unbagged fruit at S2) were used for RNA-seq analysis. The color of peel in TP70 and WP70 was still green ( Figure 1B). However, the peel color of WP110 turned to red, while TP110 still remained green, which indicated that bagging treatment significantly suppressed peel coloration. TX110 core accumulated light red pigments compared to the WX110 core, which indicated that bagging failed to completely restrain the coloration process in core. The anthocyanin content of WP110 was significantly higher than that of TP110, meanwhile, the anthocyanin content of WX110 was significantly higher than that of TX110 ( Figure 1C), which was consistent with the phenotypic presentation.

Transcriptome Assembly
Transcriptome analysis was performed on a total of 24 'TY' samples, based on which 229,397 transcripts with N50 of 1856 and 100,417 unigenes with N50 of 1600 bp were obtained. The length distribution showed that with an increase of length of transcript and unigenes, their number and proportion gradually decreased. The specific distribution data for transcripts showed that 133,270 (58.09%) were shorter than 1000 bp, 55,806 (24.33%) were ranged between 1000 and 2000 bp, 40,321 (17.58%) were longer than 2000 bp, while data for unigenes presented that 75,231 (74.92%) were shorter than 1000 bp, 13,838 (13.78%) were ranged between 1000 and 2000 bp, 11,348 (11.30%) were longer than 2000 bp ( Figure S1, Table 1). All in all, we obtained 100,417 unigenes, of which 25,186 (longer than 1000 bp) were used for further analysis.

Functional Annotation
To gain relevant functional information, all the unigenes obtained from RNA-seq were mapped to various functional databases. Among 100,417 unigenes, a total of 37,519 unigenes were annotated in databases including COG, GO, KEGG, KOG, Pfam, Swissprot, eggNOG, and NR (Table S1). Most of the unigenes could be annotated in NR and eggNOG databases, followed by Swissprot, Pfam, KOG, and GO. The distribution of unigenes among different databases was as follow: 34,561 were annotated in NR, 34,480 were annotated in eggNOG, 22,300 were annotated in Pfam, 22,123 were annotated in Swissprot, 21,560 were annotated in GO, 21,038 were annotated in KOG, 13,391 were annotated in KEGG and 10,167 were annotated in COG. The NR and eggNOG were top two among the databases, which include most of the annotated unigenes ranging from 300 to 1000 bp. The unigenes composed of more than 1000 bp followed a similar pattern ( Table 2). The abovementioned results about annotated unigenes inferred that most of the unigenes could be successfully annotated in different databases, which could provide invaluable functional basic for further exploration of genes.

Expression Analysis of DEGs among Different Pairwise
The present study aimed to assess the differences in coloration process among different fruit parts. Therefore, a co-analysis was performed to find out common DEGs in six different comparisons for peel including WP110 vs.   (Table S3). Co-analysis results showed that 205 DEGs were commonly expressed among six comparisons in peel ( Figure S2a), while only 5 common DEGs were found in the core ( Figure S2b). The presence of more common DEGs in the peels than that in the core indicated its significantly stronger response to light during fruit coloration, which was consistent with the phenotypic results ( Figure 1B). The largest phenotypic difference occurred in 110 DAFB based on Figure 1B,C, so the DEGs were analyzed using volcano plot presence. A total of 2424 and 1928 transcripts were up-regulated and down-regulated expression in WP110 vs. TP110, respectively ( Figure S2c, Table S4). A total of 3177 and 2771 transcripts were up-regulated and down-regulated expression in WX110 vs. TX110, respectively ( Figure S2d, Table S5).

GO and KEGG Analysis of DEGs
To find out the possible response of genes to light during peel and core coloration, the DEGs in WP110 vs. TP110 and WX110 vs. TX110 were selected to conduct GO and KEGG analysis. A total of 2849 DEGs were assigned to three GO terms including biological process, cellular component, and molecular function in WP110 vs. TP110. The 'metabolic process', 'cellular process', and 'single-organism process' were the top 3 leveltwo terms in biological process, followed by 'biological regulation', 'localization', and 'response to stimulus'. The 'cell', 'cell part', and 'membrane' were the top 3 level-two terms in cellular component, followed by 'organelle', 'membrane part', and 'organelle part'. The 'catalytic activity', 'binding', and 'transporter activity' were the top 3 leveltwo terms in molecular function, followed by 'structural molecule activity', 'nucleic acid binding transcription factor activity', and 'molecular function regulation' (Figure 2A). In order to obtain candidate genes involved in biological pathway, all unigenes were used to conduct KEGG analysis. A total of 883 DEGs were involved in 121 KEGG pathways in WP110 vs. TP110. Among these, 72, 71, and 70 DEGs were classified into 'plant hormone signal transduction', 'ribosome', and 'biosynthesis of amino acids', respectively, followed by 52 DEGs in 'starch and sucrose metabolism', 48 DEGs in 'protein processing in endoplasmic reticulum', and 48 DEGs in 'carbon metabolism' ( Figure S3a, Table S6). These results suggested the possible involvement of the abovementioned pathways during peel coloration in response to light. For WX110 vs. TX110, there were 3886 DEGs assigned to three GO terms including biological process, cellular component, and molecular function. The specific level-two terms showed similar a distributed rule with that in WP110 vs. TP110 ( Figure 2B). A total of 1225 DEGs were involved in 124 KEGG pathways in WX110 vs. TX110, among which 89, 83, and 83 DEGs were classified into 'biosynthesis of amino acids', 'carbon metabolism', and 'ribosome', respectively, followed by 80 DEGs in 'starch and sucrose metabolism', 76 DEGs in 'protein processing in endoplasmic reticulum', and 75 DEGs in 'plant hormone signal transduction', indicating that these pathways played a key role in core coloration ( Figure S3b, Table S6). The genes and main pathways were different in peel and core during fruit coloration, which indicated the presence of different response mechanism in peel and core to light. In addition, to find out the interesting genes involved in anthocyanin biosynthesis, hormone signal transduction and carbon metabolism of A. arguta fruit subjected to light treatment, 12 possible genes were screened from transcriptome data (Table S7, Figure S4), which would provide important reference for future study.

The Possible MYB-Related Genes Controlling Coloration of Different Fruit Parts
As one of the largest families in plant transcription factors, MYB typed TFs play a vital role in plant growth and development as well as fruit coloration. To find out if the candidate MYB-typed TFs, WP110 vs. TP110 and WX110 vs. TX110 continued to be served as the excavating point for further exploration. There was a total of 45 and 60 MYB-related DEGs in WP110 vs. TP110, and WX110 vs. TX110, respectively. Besides, a total of 27 common DEGs were found between each other ( Figure 3A, Table S8). There were two kinds of different expression patterns for these 27 DEGs, 18 of which were up-regulated after bagging treatment, while the other 9 DEGs presented decreasing trend ( Figure 3B), indicating that these MYB-related DEGs could respond to light. To find out which gene was involved in anthocyanin biosynthesis, the analysis of functional protein association networks was performed in online STRING database (https://string-db.org/cgi, accessed on 2 November 2010). Interestingly, among 27 DEGs, only one gene, c72412.graph c1 (the closest ortholog gene of MYB4 in Arabidopsis), interacted with anthocyanin-related genes ( Figure 3C). To explore the specific information about c72412.graph c1, NCBI-blast was conducted to search the top 20 homologous genes. The phylogenetic tree showed that c72412.graph c1 was very likely a AaMYB308like gene ( Figure 3D), based on which we named c72412.graph c1 as AaMYB308like.

The Possible MYB-Related Genes Controlling Coloration of Different Fruit Parts
As one of the largest families in plant transcription factors, MYB typed TFs play a vital role in plant growth and development as well as fruit coloration. To find out if the candidate MYB-typed TFs, WP110 vs. TP110 and WX110 vs. TX110 continued to be served as the excavating point for further exploration. There was a total of 45 and 60 MYB-related DEGs in WP110 vs. TP110, and WX110 vs. TX110, respectively. Besides, a total of 27 common DEGs were found between each other ( Figure 3A, Table S8). There were two kinds of different expression patterns for these 27 DEGs, 18 of which were up-regulated after bagging treatment, while the other 9 DEGs presented decreasing trend ( Figure 3B), indicating that these MYB-related DEGs could respond to light. To find out which gene was involved in anthocyanin biosynthesis, the analysis of functional protein association networks was performed in online STRING database (https://string-db.org/cgi, accessed on 2 November 2010). Interestingly, among 27 DEGs, only one gene, c72412.graph c1 (the closest ortholog gene of MYB4 in Arabidopsis), interacted with anthocyanin-related genes ( Figure 3C). To explore the specific information about c72412.graph c1, NCBI-blast was conducted to search the top 20 homologous genes. The phylogenetic tree showed that c72412.graph c1 was very likely a AaMYB308like gene ( Figure 3D), based on which we named c72412.graph c1 as AaMYB308like. Biology 2021, 10, x 9 of 14

Transient Overexpression of AaMYB308like in Nicotiana tabacum
In order to further investigate the specific role of AaMYB308like in anthocyanin regulation, transient overexpression of AaMYB308like was carried out in N. tabacum leaves by Agrobacterium-mediated transformation. Phenotypic observation of tobacco after infiltration showed that the leaves of transiently over-expressed AaMYB308like showed no obvious changes compared with the control with the empty vector or no treatment ( Figure  4A). However, the transcription levels of structural genes including NtPAL, NtCHI, NtDFR, NtLDOX, and NtUFGT investigated in corresponding leaves showed a different

Transient Overexpression of AaMYB308like in Nicotiana tabacum
In order to further investigate the specific role of AaMYB308like in anthocyanin regulation, transient overexpression of AaMYB308like was carried out in N. tabacum leaves by Agrobacterium-mediated transformation. Phenotypic observation of tobacco after infiltration showed that the leaves of transiently over-expressed AaMYB308like showed no obvious changes compared with the control with the empty vector or no treatment ( Figure 4A). However, the transcription levels of structural genes including NtPAL, NtCHI, NtDFR, NtLDOX, and NtUFGT investigated in corresponding leaves showed a different expression level. In over-expressed tobacco leaves, the over-expression of AaMYB308like significantly up-regulated expression levels of NtCHI, but did not change the expression of other genes ( Figure 4B), which suggested that overexpression of AaMYB308like could activate NtCHI expression in tobacco.
Biology 2021, 10, x 10 of 14 expression level. In over-expressed tobacco leaves, the over-expression of AaMYB308like significantly up-regulated expression levels of NtCHI, but did not change the expression of other genes ( Figure 4B), which suggested that overexpression of AaMYB308like could activate NtCHI expression in tobacco.

Discussion
Fruit color is an important agronomic trait and is responsible for market value. As ready-to-eat fruit, A. arguta coloration plays a key role in ensuring fruit quality and attracting consumers. Light is a key environmental factor that mediates fruit coloration by regulating related genes [21,22]. To gain molecular information of the photoresponse mechanism of different parts of A. arguta, we performed high throughput RNA-seq for a total of 24 samples. Finally, we obtained 100,417 unigenes with an average length of 863.77 bp and N50 of 1600 bp, of which 37,519 unigenes were annotated in functional databases (Tabel 2). Co-analysis presented 205 DEGs in peel and 5 DEGs in core ( Figure S2), which indicated a significantly stronger response of peel than that of core to light during fruit coloration. To find out the key differences between peel and core coloration process, WP110 vs. TP110 and WX110 vs. TX110 were selected as two main comparisons for further analysis. The 2849 and 3886 DEGs were classified into GO terms in WP110 vs. TP110 and

Discussion
Fruit color is an important agronomic trait and is responsible for market value. As ready-to-eat fruit, A. arguta coloration plays a key role in ensuring fruit quality and attracting consumers. Light is a key environmental factor that mediates fruit coloration by regulating related genes [21,22]. To gain molecular information of the photoresponse mechanism of different parts of A. arguta, we performed high throughput RNA-seq for a total of 24 samples. Finally, we obtained 100,417 unigenes with an average length of 863.77 bp and N50 of 1600 bp, of which 37,519 unigenes were annotated in functional databases (Tabel 2). Co-analysis presented 205 DEGs in peel and 5 DEGs in core ( Figure S2), which indicated a significantly stronger response of peel than that of core to light during fruit coloration. To find out the key differences between peel and core coloration process, WP110 vs. TP110 and WX110 vs. TX110 were selected as two main comparisons for further analysis. The 2849 and 3886 DEGs were classified into GO terms in WP110 vs. TP110 and WX110 vs. TX110, respectively. Similarly, 883 and 1225 DEGs were assigned into KEGG pathways in WP110 vs. TP110 and WX110 vs. TX110, respectively (Figure 2, Figure S3). Different DEGs involved in GO and KEGG suggested that peel and core had a different photoresponse mechanism during fruit coloration. Based on the KEGG enrichment results, we can conclude that 'plant hormone signal transduction' is the most important pathway in peel coloration, while 'carbon metabolism' is the most important pathway in core coloration. As the outermost part of A. arguta, the peel coloration is an important aspect for fruit appearance and market value.
To further explore the specific genes involved in peel coloration, a light responsive transcription factor AaMYB308like was cloned from A. arguta peel. In order to further confirm its function during the coloring process, the AaMYB308like-PBI121 expression vector was constructed and transiently transformed into N. tabacum leaves. The results showed that transient overexpression of AaMYB308like in N. tabacum leaves significantly up-regulated NtCHI (EBG) expression, which suggests the possible involvement of AaMYB308like in anthocyanin biosynthesis. However, there were no changes of color and anthocyanin accumulation, which indicated the inadequacy of AaMYB308like in activating the whole anthocyanin biosynthetic pathway ( Figure 4). Our results were at part with a previous study, where transient transformation of AcMYB5-1/5-2/A1-1 cloned from A. chinenses in N. benthamiana leaves showed no change in leaf color and accumulation of anthocyanin [23].
To perform the homology analysis between AaMYB308like identified above and other anthocyanin-related MYB typed transcription factors, the MYB TFs sequences in model plants, such as Arabidopsis thaliana and Solanum lycopersicum, and several fruit plants including Malus domestica, Citrus sinensis, Prunus persica, Pyrus pyrifolia, Vitis vinifera, Fragaria ananassa, Fragaria vesca, Actinidia chinensis, and Actinidia arguta were collected from NCBI database. A total of 18 MYB TFs were selected for the sequence and homology analysis. The homology of AaMYB308like and these 18 MYB TFs were maintained at a middle level with the identity ranging from 45.68% to 69.83%, which indicates that AaMYB308like was indeed a MYB typed TF and played a regulating role in anthocyanin biosynthesis (Table 3). MdMYB1 was a key TF regulating anthocyanin biosynthesis in apple (Malus domestica) by not only activating expression of downstream structural genes, but also interacting with other TFs including MdERF3 (ETHYLENE RESPONSE FACTOR 3) and MdEIL1 (ETHYLENE INSENSITIVE 3-LIKE 1) to mediate ethylene responding fruit coloration. Additionally, as a light response factor, MdMYB1 was also served as a light response factor participating in the light-mediated anthocyanin regulation [24,25]. In our study, AaMYB308like was screened as a candidate light response factor to regulate anthocyanin biosynthesis, showing the possible consistency with MdMYB1. However, whether AaMYB308like works in the same way with MdMYB1 involved in fruit coloration still needs further confirmation. Besides in other species, MYB transcription factors in kiwifruit have also been found to be involved in anthocyanin regulation. AcMYB75 was identified to be an R2R3-MYB TF to regulate anthocyanin accumulation by activating the promoter of AcANS in A. chinensis [26]. Whether AaMYB308like regulates anthocyanin biosynthesis in a similar way needs to be confirmed. Overall, the candidate AaMYB308like TF could be used as the key factor to participate in light-mediated anthocyanin regulation, while the specific regulatory mechanism needs to be further explored.

Conclusions
In the present study, we revealed the molecular differences in the coloration process between A. arguta peel and core. The 'plant hormone signal transduction' was the key pathway during peel coloring process, while the 'carbon metabolism' was the key pathway in core coloring process. Additionally, transient overexpression of AaMYB308like in N. tabacum confirmed the active role of photoresponse factor AaMYB308like in anthocyanin biosynthesis. This study will provide significant molecular resources for understanding the intricate mechanisms and pathways involved in anthocyanin biosynthesis. Furthermore, this study will help the researchers to improve the fruit appearance.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/biology10070648/s1, Figure S1: Distribution of unigenes and transcripts, Figure S2: Venn diagram and volcano plot of DEGs, Figure S3: KEGG enrichment of DEGs from WP110 vs. TP110 and WX110 vs. TX110, Figure S4: The expression patterns of 12 possible genes were involved in three KEGG pathways including anthocyanin biosynthesis, hormone signal transduction, and carbon metabolism, Table S1: The annotation information of unigenes in databases, Table S2: Information for six comparisons in peel, Table S3: Information for six comparisons in core, Table S4: The up/downregulated transcripts in WP110 vs TP110, Table S5: The up/down-regulated transcripts in WX110 vs TX110, Table S6: KEGG enrichment information of DEGs in WP110 vs TP110 and WX110 vs TX110, Table S7: The specific annotation information of a representative group of DEGs related to three main pathways including 'anthocyanin biosynthesis', 'plant hormone signal transduction' and 'carbon metabolism' involved in light-dependent fruit (peel and core) coloration, Table S8: MYB-typed DEGs in WP110 vs TP110 and WX110 vs TX110.