Exogenous Abscisic Acid Regulates Anthocyanin Biosynthesis and Gene Expression in Blueberry Leaves

: Blueberry ( Vaccinium corymbosum ) leaves have a positive influence on health because of their phenolic contents, including anthocyanins. Phytohormone abscisic acid (ABA) promotes anthocyanin accumulation, but the underlying mechanisms are unclear in blueberry leaves. In this study, we found that exogenous ABA promotes anthocyanin accumulation in blueberry leaves and we explored the global molecular events involved in these physiological changes by treating in vitro-grown blueberry seedlings with ABA and performing transcriptome deep sequencing (RNA-seq). We identified 6390 differentially expressed genes (DEGs), with 2893 DEGs at 6 h and 4789 at 12 h of ABA treatment compared to the control. Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways related to plant hormone signal transduction and phenylpropanoid and flavonoid biosynthesis were significantly enriched at both stages of the ABA treatment. Analysis of DEGs in plant hormone signal transduction pathways revealed that exogenous ABA affected the expression of genes from other plant hormone signaling pathways, especially brassinosteroid, auxin, and gibberellin signaling. To elucidate the mechanism driving anthocyanin biosynthesis in blueberry in response to ABA treatment, we screened anthocyanin biosynthesis structural genes (ASG) from the phenylpropanoid and flavonoid biosynthetic pathways, MYB transcription factor genes from R2R3-MYB subgroups 5, 6, and 7 and ABRE-binding factor (ABF) genes from the ABA signal transduction pathway. Pearson’s correlation coefficient ( r ) analysis indicated that the ABF s, MYB s, and structural genes form a network to regulate ABA-induced anthocyanin biosynthesis and MYBA1 is likely to play an important role in this regulatory network. These findings lay the foundation for improving anthocyanin biosynthesis in blueberry leaves.

Horticulturae 2024, 10,192 2 of 17 ABA regulates multiple aspects of plant physiology via its signal transduction pathway, but how ABA regulates anthocyanin biosynthesis via the signal transduction pathway in blueberry leaves is unclear.
Blueberry (Vaccinium corymbosum) is a popular fruit worldwide, in part due to the health benefits associated with its rich flavonoid compounds, particularly anthocyanins [27,28].The fruits of the Vaccinium species are the main commercial products, whereas the leaves have been used mainly for medicinal purposes [29][30][31][32].In the wild lowbush blueberry (V.angustifolium) and lingonberry (V.vitis-idaea), the total antioxidant capacities in terms of radical scavenging activity and reducing power were much higher in the leaves of both plants as compared to their fruits, which were in correlation with phenolic contents including total flavonoids, anthocyanins, and tannins [33].The total flavonoids from blueberry leaves had both hypolipidemic and antioxidant effects, which could help the cure and prevention of hyperlipidemia diseases [34].Thus, it is very important to study the mechanism of ABA induced anthocyanin biosynthesis in blueberry leaves.
In this study, to obtain differential expression genes (DEGs) from blueberry leaves in response to ABA treatment, we performed transcriptome deep sequencing (RNA-seq) analysis by treating in vitro-grown blueberry seedlings with ABA and annotated DEGs from KEGG pathways related to plant hormone signal transduction.To elucidate the mechanism driving anthocyanin biosynthesis in blueberry in response to ABA treatment, we assembled a regulatory network of ABFs-MYB-ASGs by gene co-expression analysis.Our findings increase understanding of ABA-induced anthocyanin biosynthesis and provide a basis for improving the value of blueberry leaves via ABA treatment.

Plant Materials and ABA Treatments
In vitro-grown seedlings were obtained from the blueberry (Vaccinium corymbosum) cultivar 'Northland'.The materials were cultured on modified woody plant medium (WPM) containing Murashige and Skoog vitamins under a 16-h-light/8-h-dark photoperiod of 2000 Lux from cool white fluorescent tubes at 25 • C. The seedlings were cultured on WPM containing 1.0 mg/L trans-Zeatin for 5 weeks, and then, seedlings with six to ten blades were transferred to the medium containing 100 µM ABA for 0, 6 or 12 h for RNA-seq and RT-qPCR analysis and measurement of the anthocyanin contents [35,36].

Transcriptome Sequencing
Total RNA was extracted from seedlings after 0, 6, or 12 h of ABA treatment using an RNAprep Pure Plant Kit (Tiangen, Beijing, China), and RNA concentration and purity were measured using a NanoDrop 2000 (Thermo Fisher Scientific, Wilmington, DE, USA).Sequencing libraries were generated using a Hieff NGS Ultima Dual-mode mRNA Library Prep Kit for Illumina (Yeasen Biotechnology (Shanghai) Co., Ltd., Shanghai, China).The libraries were sequenced on an Illumina NovaSeq platform to generate 150-bp pairedend reads.The raw reads were processed with a bioinformatic pipeline tool, BMKCloud (www.biocloud.net(accessed on 1 July 2023)).Clean data (clean reads) were obtained by removing reads containing adapters, reads containing ploy-N, and low-quality reads from raw data.The Q30, GC content, and sequence duplication level of the clean data were calculated.All downstream analyses were based on clean, high-quality data.The clean reads were mapped to the reference Vaccinium corymbosum cv.Draper V1.0 genome sequence (https://www.vaccinium.org/genomes(accessed on 1 August 2023)) using Hisat2 tools software version 2.0.4.
KEGG pathway enrichment analysis of the DEGs was performed using the KOBAS database and clusterProfiler software (version 4.4.4)[37].Heatmaps of the expression of DEGs in the plant hormone signal transduction pathways and the anthocyanin biosynthetic pathway under ABA treatment were assembled using log 2 (FPKM) values with Tbtools (v1.098761) software [38].The regulatory network of the anthocyanin biosynthetic pathway was constructed based on Pearson's correlation coefficients (r) between R2R3-MYB transcription factor genes and anthocyanin biosynthesis/ABF genes using IBM SPSS software (version 19.0).

Analysis of Abscisic Acid Responsiveness Element in the Promoter
The upstream sequences of screened MYB and anthocyanin biosynthesis genes were downloaded from the GDV.About 2000-bp upstream sequences were submitted to the online program PlantCARE (https://bioinformatics.psb.ugent.be/webtools/plantcare/html/ (accessed on 1 November 2023)) to predict abscisic acid responsiveness element (ABRE) or MYB binding sites.

Reverse-Transcription Quantitative PCR
Total RNA was extracted from each sample using an RNA Extraction Kit (Sangon Biotech, Shanghai, China), and first-strand cDNAs were synthesized using PrimeScriptTM RT Master Mix (TaKaRa, Kusatsu, Japan).RT-qPCR analysis was performed using an ABI 7900HT real-time PCR system.Nine ABA signal transduction pathway genes (ABF2d, ABF4b, PP2C-24a, PP2C-37a, PP2C-7a, PP2C-51b, PYL1, SnRK2A, and PYL4a) were selected to validate the RNA-seq data.Glyceraldehyde-3-phosphate dehydrogenase housekeeping (GAPDH; GenBank accession no.AY123769) was used as the reference gene, and the relative expression levels of each gene were calculated using the 2 −∆∆Ct method.Primer sequences for RT-qPCR are shown in Table S1.

Measurement of Anthocyanin Contents
Total anthocyanin contents were measured as described by Yang et al. [39].Each sample was ground to a powder, and 1 g was extracted in 6 mL of 1% (v/v) HCl in methanol and incubated at 4 • C for 16 h.After the mixtures were centrifuged (8000× g, 10 min, 4 • C), 4 mL of the supernatant was diluted with 4 mL ddH 2 O, and the absorbance was measured at 530 and 650 nm.Total anthocyanin content was calculated using the equation A 530 − 0.25 × A 650 [40].

Cluster Analysis and Sequence Alignment of MYB Transcription Factors
A phylogenetic tree was constructed through amino sequencing of 23 differentially expressed R2R3-MYB genes from blueberry and 124 R2R3-MYB proteins from Arabidopsis using the neighbor-joining method with MEGA X software (https://www.megasoftware.net(accessed on 30 October 2023)).The phylogenetic tree was divided into different subgroups based on MYB transcription factors from Arabidopsis [24,41].The MYBs from blueberry were named based on cluster and functional annotations from SwissProt databases.Sequence alignment and phylogenetic analysis between MYBA1, MYBA2 and VcMYBA (Gen-Bank accession no.MH105054) were performed using DNAMAN software, Version 8.192.

Statistical Analysis
The gene expression levels or anthocyanin contents were carried out using three independent biological replicates, and three technical replicates were performed for each biological replicate.One-way analysis of variance (ANOVA) was used to assess the differences in gene expression levels or anthocyanin contents during ABA treatment, and Tukey's test was used to identify significant differences at p-value ≤ 0.05 using SPSS 19.0 software.

ABA Treatment Promotes Anthocyanin Accumulation
To characterize the effect of exogenous ABA on anthocyanin accumulation in blueberry leaves, we treated in vitro-grown blueberry seedlings with exogenous ABA and without ABA for 0, 6, or 12 h (Figure 1a).The tip of the leaves turned red after 6 h of ABA treatment.We measured the anthocyanin contents of in vitro-grown blueberry seedlings under ABA treatment (Figure 1b).The anthocyanin contents significantly increased after 6 h and 12 h of ABA treatment compared to the 0 h control.At the same time, without ABA treatment, the content of anthocyanin did not change at 0 h, 6 h or 12 h.Thus, exogenous ABA promotes anthocyanin accumulation in vitro-grown blueberry seedlings.

Transcriptome Sequencing and Analysis of DEGs
To elucidate the molecular changes underlying the response of blueberry to exogenous ABA, we performed RNA-seq analysis using in vitro-grown blueberry seedlings under ABA treatment for 6 or 12 h, using untreated seedlings as the 0 h control.Table 1 summarizes the nine samples used for RNA-seq.More than 19.3 million clean reads were obtained, generating more than 5.78 Gb of clean bases for each sample.We mapped the clean reads to the blueberry reference genome at a rate of approximately 89.9% to 91.2% per sample.The percentage of Q30 bases was 90.45% or above, and the GC content ranged from 45.7% to 46.55%.
Genes with FDR < 0.01 and absolute log2 (fold-change) ≥ 1 between control samples (0 h) and samples under ABA treatment (6 and 12 h) were considered to be DEGs.We identified 2893 DEGs in the 0 h vs. 6 h comparison, consisting of 1765 upregulated and 1128 downregulated DEGs (Figure 2a).For the 0 h vs. 12 h comparison, the number of DEGs increased to 4789, with more downregulated DEGs (2984) than upregulated DEGs (1805) (Figure 2b).As shown in the Venn diagram in Figure 2c, across the two comparisons, we detected 6390 DEGs, with 1292 DEGs shared between 0 h vs. 6 h and 0 h vs. 12 h.

Transcriptome Sequencing and Analysis of DEGs
To elucidate the molecular changes underlying the response of blueberry to exogenous ABA, we performed RNA-seq analysis using in vitro-grown blueberry seedlings under ABA treatment for 6 or 12 h, using untreated seedlings as the 0 h control.Table 1 summarizes the nine samples used for RNA-seq.More than 19.3 million clean reads were obtained, generating more than 5.78 Gb of clean bases for each sample.We mapped the clean reads to the blueberry reference genome at a rate of approximately 89.9% to 91.2% per sample.The percentage of Q30 bases was 90.45% or above, and the GC content ranged from 45.7% to 46.55%.Genes with FDR < 0.01 and absolute log 2 (fold-change) ≥ 1 between control samples (0 h) and samples under ABA treatment (6 and 12 h) were considered to be DEGs.We identified 2893 DEGs in the 0 h vs. 6 h comparison, consisting of 1765 upregulated and 1128 downregulated DEGs (Figure 2a).For the 0 h vs. 12 h comparison, the number of DEGs increased to 4789, with more downregulated DEGs (2984) than upregulated DEGs (1805) (Figure 2b).As shown in the Venn diagram in Figure 2c, across the two comparisons, we detected 6390 DEGs, with 1292 DEGs shared between 0 h vs. 6 h and 0 h vs. 12 h.12 h; purple, DEGs for 0 h vs. 6 h and 0 h vs. 12h.

Functional Annotation and KEGG Pathway Analysis of DEGs under ABA Treatment
We annotated the functions of the DEGs under ABA treatment using the COG, GO, KEGG, KOG, NR, Pfam, Swiss-Prot, and eggNOG databases (Table S2).The total number of annotated DEGs increased from 2730 (94.37%) for the 0 h vs. 6 h comparison to 4623 (96.53%) for the 0 h vs. 12 h comparison; we annotated 1952 and 3453 of these DEGs using the KEGG database for each comparison, respectively.
To identify the biological functions of the DEGs in blueberry, we subjected all DEGs to KEGG pathway enrichment analysis and KEGG functional classification (Figure 3).Analysis of the top 20 enriched KEGG pathways showed that the plant hormone signal transduction pathway (ko04075), phenylpropanoid biosynthetic pathway (ko00940), and flavonoid biosynthetic pathway (ko00941) are significantly enriched at all stages (Figure 3a-c).At 6 and 12 h into ABA treatment, 291, 197, and 85 DEGs were enriched in the plant hormone signal transduction pathway, phenylpropanoid biosynthetic pathway, and

Functional Annotation and KEGG Pathway Analysis of DEGs under ABA Treatment
We annotated the functions of the DEGs under ABA treatment using the COG, GO, KEGG, KOG, NR, Pfam, Swiss-Prot, and eggNOG databases (Table S2).The total number of annotated DEGs increased from 2730 (94.37%) for the 0 h vs. 6 h comparison to 4623 (96.53%) for the 0 h vs. 12 h comparison; we annotated 1952 and 3453 of these DEGs using the KEGG database for each comparison, respectively.
To identify the biological functions of the DEGs in blueberry, we subjected all DEGs to KEGG pathway enrichment analysis and KEGG functional classification (Figure 3).Analysis of the top 20 enriched KEGG pathways showed that the plant hormone signal transduction pathway (ko04075), phenylpropanoid biosynthetic pathway (ko00940), and flavonoid biosynthetic pathway (ko00941) are significantly enriched at all stages (Figure 3a-c).At 6 and 12 h into ABA treatment, 291, 197, and 85 DEGs were enriched in the plant hormone signal transduction pathway, phenylpropanoid biosynthetic pathway, and flavonoid biosynthetic pathway, respectively (Figure 3c).KEGG functional classification showed that the DEGs are involved in cellular processes, environmental information processing, genetic information processing, metabolism, and organismal systems (Figure 3d-f).
The KEGG pathway with the greatest number of DEGs was the plant hormone signal transduction pathway for all stages, with more upregulated DEGs than downregulated DEGs at 6 h and more downregulated than upregulated DEGs at 12 h.Among metabolic pathways, those with the greatest number of DEGs were the phenylpropanoid biosynthetic pathway (197, 7.9%), the starch and sucrose metabolism pathway (158, 6.3%), and the flavonoid biosynthetic pathway (85, 3.4%) (Figure 3f).Thus, DEGs from the plant hormone signal transduction pathway, the phenylpropanoid biosynthetic pathway, and the flavonoid biosynthetic pathway play important roles in the response of blueberry to exogenous ABA.
Horticulturae 2024, 10, x FOR PEER REVIEW 7 of 17 flavonoid biosynthetic pathway, respectively (Figure 3c).KEGG functional classification showed that the DEGs are involved in cellular processes, environmental information processing, genetic information processing, metabolism, and organismal systems (Figure 3df).
The KEGG pathway with the greatest number of DEGs was the plant hormone signal transduction pathway for all stages, with more upregulated DEGs than downregulated DEGs at 6 h and more downregulated than upregulated DEGs at 12 h.Among metabolic pathways, those with the greatest number of DEGs were the phenylpropanoid biosynthetic pathway (197, 7.9%), the starch and sucrose metabolism pathway (158, 6.3%), and the flavonoid biosynthetic pathway (85, 3.4%) (Figure 3f).Thus, DEGs from the plant hormone signal transduction pathway, the phenylpropanoid biosynthetic pathway, and the flavonoid biosynthetic pathway play important roles in the response of blueberry to exogenous ABA.

Exogenous ABA Induces the Expression of Plant Hormone Signal Transduction Pathway Genes
To elucidate the functions of DEGs from the plant hormone signal transduction pathway, we identified 209 DEGs (Table S3).These DEGs were mapped to the ABA, GA, JA, auxin, SA, BR, ethylene, and CK biosynthetic pathways (Figure 4).Of these DEGs, 101, 35, 26, and 18 were mapped to the BR, ABA, auxin, and GA biosynthetic pathways, respectively.In the ABA signaling pathway, eight PYR/PYL family members were downregulated and one was upregulated.All PP2C gene family members were upregulated at 6 and 12 h of ABA treatment relative to the 0 h control.Of the five SnRK2 family members, one was downregulated and the remaining four were upregulated.All the ABFs were upregulated under ABA treatment.In the BR signal transduction pathway, 50 brassinosteroid insensitive 1 (BRI1) genes and 34 BRI1 associated protein kinase 1 (BAK1) genes were differentially expressed, including 9 that were upregulated and 75 that were downregulated.In the auxin signal transduction pathway, six AUX/IAA genes were upregulated and two were downregulated.In the GA signal transduction pathway, 17 genes encoding DELLA family members were identified, with 6 being upregulated and 11 being downregulated.These results indicate that ABA treatment not only induced the expression of ABA signal transduction pathway genes, but it also affected the expression of other plant hormone signal transduction pathway genes, especially the BR, auxin, and GA signal transduction pathways.
We selected nine genes from the ABA signal transduction pathway to validate the accuracy and reliability of the RNA-seq data by RT-qPCR (Figure S1a).Pearson's correlation coefficient (r) analysis showed that log2 (2 -ΔΔCt ) values from RT-qPCR and log2 (FC) from RNA-seq are significantly correlated for the 0 h vs. 6 h (0.943) and 0 h vs. 12 h (0.872) comparisons (Figure S1b).These results indicate that the RNA-seq data generated in this study are accurate and reliable.

Exogenous ABA Induces the Expression of Plant Hormone Signal Transduction Pathway Genes
To elucidate the functions of DEGs from the plant hormone signal transduction pathway, we identified 209 DEGs (Table S3).These DEGs were mapped to the ABA, GA, JA, auxin, SA, BR, ethylene, and CK biosynthetic pathways (Figure 4).Of these DEGs, 101, 35, 26, and 18 were mapped to the BR, ABA, auxin, and GA biosynthetic pathways, respectively.In the ABA signaling pathway, eight PYR/PYL family members were downregulated and one was upregulated.All PP2C gene family members were upregulated at 6 and 12 h of ABA treatment relative to the 0 h control.Of the five SnRK2 family members, one was downregulated and the remaining four were upregulated.All the ABFs were upregulated under ABA treatment.In the BR signal transduction pathway, 50 brassinosteroid insensitive 1 (BRI1) genes and 34 BRI1 associated protein kinase 1 (BAK1) genes were differentially expressed, including 9 that were upregulated and 75 that were downregulated.In the auxin signal transduction pathway, six AUX/IAA genes were upregulated and two were downregulated.In the GA signal transduction pathway, 17 genes encoding DELLA family members were identified, with 6 being upregulated and 11 being downregulated.These results indicate that ABA treatment not only induced the expression of ABA signal transduction pathway genes, but it also affected the expression of other plant hormone signal transduction pathway genes, especially the BR, auxin, and GA signal transduction pathways.
We selected nine genes from the ABA signal transduction pathway to validate the accuracy and reliability of the RNA-seq data by RT-qPCR (Figure S1a).Pearson's correlation coefficient (r) analysis showed that log 2 (2 −∆∆Ct ) values from RT-qPCR and log 2 (FC) from RNA-seq are significantly correlated for the 0 h vs. 6 h (0.943) and 0 h vs. 12 h (0.872) comparisons (Figure S1b).These results indicate that the RNA-seq data generated in this study are accurate and reliable.

Exogenous ABA Treatment Regulates the Expression of R2R3-MYB Transcription Factor Genes
To elucidate the functions of transcription factors in response to ABA treatment, we analyzed the differentially expressed R2R3-MYB transcription factor genes under ABA treatment according to Pfam, Swiss-Prot, and NR functional annotations.After removing the repeat MYB transcription factor genes based on sequence analysis, we identified 23 R2R3-MYB subfamily genes, of which 18 were upregulated and 5 were downregulated under ABA treatment (Table S4).To predict the potential functions of the proteins encoded by these R2R3-MYB genes, we reconstructed a phylogenetic tree based on their deduced protein sequences and 124 R2R3-MYB proteins from Arabidopsis (Figure 5).The blueberry R2R3-MYB transcription factors belong to 12 subgroups, which were previously defined for Arabidopsis R2R3-MYB proteins.MYBA1 and MYBA2 were clustered in subgroup 6; their encoding were upregulated under exogenous ABA treatment.Sequence alignment and phylogenetic analysis showed that MYBA1 shared 97.64% amino

Exogenous ABA Treatment Regulates the Expression of R2R3-MYB Transcription Factor Genes
To elucidate the functions of transcription factors in response to ABA treatment, we analyzed the differentially expressed R2R3-MYB transcription factor genes under ABA treatment according to Pfam, Swiss-Prot, and NR functional annotations.After removing the repeat MYB transcription factor genes based on sequence analysis, we identified 23 R2R3-MYB subfamily genes, of which 18 were upregulated and 5 were downregulated under ABA treatment (Table S4).To predict the potential functions of the proteins encoded by these R2R3-MYB genes, we reconstructed a phylogenetic tree based on their deduced protein sequences and 124 R2R3-MYB proteins from Arabidopsis (Figure 5).The blueberry R2R3-MYB transcription factors belong to 12 subgroups, which were previously defined for Arabidopsis R2R3-MYB proteins.MYBA1 and MYBA2 were clustered in subgroup 6; their encoding genes were upregulated under exogenous ABA treatment.Sequence alignment and phylogenetic analysis showed that MYBA1 shared 97.64% amino sequence identity with VcMYBA (GenBank accession no.MH105054) (Figure S2).Thus, MYBA1 is probably homologous to VcMYBA, which belongs to subgroup 6 and induces anthocyanin accumulation by heterologous expression in Nicotiana benthaminana [21].MYB5, MYB11, and MYB12 belong to subgroup 7, and their encoding genes were upregulated under ABA treatment.MYB23 and MYB123 clustered in subgroup 5, and their encoding genes were repressed in response to exogenous ABA treatment.In general, R2R3-MYB transcription factors of subgroups 5, 6, and 7 regulated flavonoid biosynthesis [41].Thus, these MYB transcription factors may be regulated ABA-induced anthocyanin accumulation.
Horticulturae 2024, 10, x FOR PEER REVIEW 10 of 17 sequence identity with VcMYBA (GenBank accession no.MH105054) (Figure S2).Thus, MYBA1 is probably homologous to VcMYBA, which belongs to subgroup 6 and induces anthocyanin accumulation by heterologous expression in Nicotiana benthaminana [21].MYB5, MYB11, and MYB12 belong to subgroup 7, and their encoding genes were upregulated under ABA treatment.MYB23 and MYB123 clustered in subgroup 5, and their encoding genes were repressed in response to exogenous ABA treatment.In general, R2R3-MYB transcription factors of subgroups 5, 6, and 7 regulated flavonoid biosynthesis [41].Thus, these MYB transcription factors may be regulated ABA-induced anthocyanin accumulation.

Exogenous ABA Induces the Expression of Anthocyanin Pathway Genes via MYB or ABF Transcription Factors
To clarify the regulatory network of ABA-induced anthocyanin biosynthesis, we screened the DEGs involved in anthocyanin biosynthesis from the phenylpropanoid and flavonoid biosynthetic pathway based on KEGG annotation.After removing repeat sequences, we retained seven DEGs (Table S5).Of these, C4H, CHI, CHI3, VcUFGT, and

Exogenous ABA Induces the Expression of Anthocyanin Pathway Genes via MYB or ABF Transcription Factors
To clarify the regulatory network of ABA-induced anthocyanin biosynthesis, we screened the DEGs involved in anthocyanin biosynthesis from the phenylpropanoid and flavonoid biosynthetic pathway based on KEGG annotation.After removing repeat sequences, we retained seven DEGs (Table S5).Of these, C4H, CHI, CHI3, VcUFGT, and UGT75C1 were upregulated and 4CL7 and LDOX were downregulated in response to ABA treatment (Figure 6 and Table S5).We collected the 2000 bp upstream sequences from the ATG start codons of the differentially expressed MYB and anthocyanin biosynthesis genes and looked for ABREs or MYB binding sites using the PlantCARE online tool (Figure S3).The promoter sequences of all the MYBs contain ABREs, including MYB23, which contains six ABREs.Among the anthocyanin biosynthesis genes, the C4H, 4CL7, CHI, and LDOX promoters contained one to four ABREs.However, the promoter sequences of anthocyanin biosynthesis genes contain three to nine MYB binding sites.In general, bZIP-type AREB/ABF transcription factors target ABREs and MYB transcription factors target MYB binding sites in the promoter regions of genes to regulate their expression.
type AREB/ABF transcription factors target ABREs and MYB transcription factors target MYB binding sites in the promoter regions of genes to regulate their expression.
To investigate the functions of MYB transcription factors from subgroups 5, 6, and 7 and ABF transcription factors in the anthocyanin biosynthetic pathway, we calculated Pearson's correlation coefficients (r) between ABFs and MYBs, between ABFs and anthocyanin biosynthetic genes, and between MYBs and anthocyanin biosynthetic genes according to FPKM values under ABA treatment (Table S6).The expression of MYB12, MYB23, MYB5, and MYBA1 was positively correlated with that of most anthocyanin metabolism genes.However, the expression of some MYBs was negatively correlated with that of anthocyanin pathway genes, such as MYB123 and CHI3, MYB23 and VcUFGT, and MYBA1 and 4CL7.The expression of ABFs was positively correlated with that of most anthocyanin metabolism genes.Correlation analysis between ABFs and MYBs showed that the expression of most ABFs is positively correlated with that of MYBs, and the expression of most ABFs and MYBs was positively correlated with that of anthocyanin metabolism genes.Based on the results of correlation analysis, we constructed a proposed regulatory network of ABA-induced anthocyanin biosynthetic pathway genes by MYB and ABF transcription factors (Figure 6).To investigate the functions of MYB transcription factors from subgroups 5, 6, and 7 and ABF transcription factors in the anthocyanin biosynthetic pathway, we calculated Pearson's correlation coefficients (r) between ABFs and MYBs, between ABFs and anthocyanin biosynthetic genes, and between MYBs and anthocyanin biosynthetic genes according to FPKM values under ABA treatment (Table S6).The expression of MYB12, MYB23, MYB5, and MYBA1 was positively correlated with that of most anthocyanin metabolism genes.However, the expression of some MYBs was negatively correlated with that of anthocyanin pathway genes, such as MYB123 and CHI3, MYB23 and VcUFGT, and MYBA1 and 4CL7.The expression of ABFs was positively correlated with that of most anthocyanin metabolism genes.Correlation analysis between ABFs and MYBs showed that the expression of most ABFs is positively correlated with that of MYBs, and the expression of most ABFs and MYBs was positively correlated with that of anthocyanin metabolism genes.Based on the results of correlation analysis, we constructed a proposed regulatory network of ABA-induced anthocyanin biosynthetic pathway genes by MYB and ABF transcription factors (Figure 6).

Exogenous ABA Promoter Anthocyanin Accumulation of Blueberry Leaves
Blueberries are the richest sources of polyphenolic antioxidants, especially anthocyanins, among all fruits and vegetables [42].Blueberry leaves, which have been used mainly for medicinal purposes, are considered essentially as waste or by-products.Various health benefits were reported for leaf extracts from the Vaccinium species [29][30][31][32].A large number of studies have shown that exogenous ABA stimulates anthocyanin accumulation of blueberry fruits [7,43,44].In this study, we also found that exogenous ABA also promotes anthocyanin biosynthesis in high blueberry leaves.Therefore, exogenous ABA treatment can improve the utilization value of blueberry leaves.

Exogenous ABA Regulated the Expression of ABA Signal Transduction Pathway Genes
The plant hormone ABA mediates various physiological responses in plants via the PYR/PYL/RCAR, PP2C, and SnRK2 family members; these three families constitute the core components of the ABA signal transduction pathway [5,9,10].ABA treatment regulates the expression of the genes encoding these family members.For example, in the primary roots of maize (Zea mays) seedlings, ABA treatment upregulated ZmPYL1-3 genes (encoding homologs of dimeric-type Arabidopsis ABA receptors) but downregulated ZmPYL4-10 (encoding homologs of monomeric-type Arabidopsis ABA receptors).However, the opposite trend was observed in the leaves of maize plants growing in medium containing ABA.At the same time, three PP2C and three SnRK2 family genes were upregulated in primary roots in response to ABA treatment [45].In this study, among the PYR/PYL/RCAR, PP2C, and SnRK2 family members, PYL1 was upregulated and PYL4a-g was downregulated in in vitro-grown blueberry seedlings in response to ABA treatment (PYL1 and PYL4 were named based on their homologs in Arabidopsis).In addition, the genes encoding all PP2C and most SnRK2 family members were upregulated in in vitro-grown blueberry seedlings under ABA treatment (Table S3 and Figure 4).These results indicate that the primary roots of maize seedlings and in vitro-grown blueberry seedlings employ similar mechanisms in response to exogenous ABA treatment.

Exogenous ABA Regulated the Expression of Other Plant Hormone Signal Transduction Pathway Genes
Exogenous ABA treatment not only induces the expression of ABA signal transduction pathway genes but also affects the expression of genes in other plant hormone signal transduction pathways, highlighting the complex crosstalk between ABA and other plant hormones in regulating plant defense responses and anthocyanin biosynthesis [1,46,47].Transcriptome analysis in mangrove (Kandelia obovata) showed that plant hormone signal transduction pathways were significantly enriched in DEGs related to response to ABA treatment during vivipary [48].Our data indicate that signal transduction pathways of other plant hormones are also significantly enriched during exogenous ABA treatment in in vitro-grown blueberry seedlings, particularly BR, followed by auxin and GA signaling (Table S3 and Figure 4).
BRs are steroidal plant hormones that are essential for plant growth and development [49,50].BRs have significant effects on plant secondary metabolism, particularly anthocyanin biosynthesis [51,52].However, several reports indicate that BR treatment inhibited red-flesh coloration and repressed high light-induced anthocyanin accumulation in red-fleshed apple (Malus domestica) seedlings [53,54].In this study, ABA treatment downregulated most BRI1 and BAK1 family genes.These family members form surface receptor kinase complexes that are responsible for perceiving BRs [55].Thus, exogenous ABA may inhibit endogenous BR biosynthesis to repress plant growth and promote flavonoid biosynthesis in blueberry.
Auxin, the first plant hormone discovered, promotes plant growth and regulates secondary metabolism [56,57].Auxin affects the levels of other phytohormones; in turn, other plant hormones affect the level of auxin.In sweet cherry (Prunus avium), 1-naphthaleneacetic acid (NAA) treatment altered ethylene production, which induced fruit ripening and enhanced anthocyanin production [58].However, in apple calli, NAA and 2,4-dichlorophenoxyacetic acid (2,4-D) treatment inhibited anthocyanin accumulation [59,60].Aux/IAAs repress auxin responses and play central roles in the auxin signal transduction pathway [61,62].Our transcriptome analysis showed that most Aux/IAA genes were upregulated under exogenous ABA treatment in in vitro-grown blueberry seedlings.Thus, exogenous ABA also increased endogenous auxin contents in blueberry.
GA promotes plant growth and cell division and regulates anthocyanin biosynthesis [2,63].DELLA proteins are repressors of GA signaling that modulate all aspects of GA-induced growth and development [4,64].ABA attenuates GA biosynthesis and increases DELLA gene expression, leading to a reduction in bioactive GA levels [2].DELLAs play a major role in mediating the crosstalk among plant hormonal signals [65].Our transcriptome analysis showed that 17 DELLA genes were regulated by ABA treatment in in vitro-grown blueberry seedlings, indicating that DELLAs modulate the crosstalk between the ABA and GA signaling pathways.

The regulatory Network of ABA-Induced Blueberry Anthocyanin Accumulation
ABA plays important roles in promoting anthocyanin production and fruit ripening [37,66].In general, ABA induces gene expression by binding to conserved cis-acting ABREs in the promoters of their target genes.ABF transcription factors regulate the ABREmediated transcription of downstream target genes [12,14,67].R2R3-MYB transcription factors from subgroups 5, 6, and 7 regulate flavonoid biosynthesis via the transcriptional regulation of flavonoid structural genes [21,23,68].In Arabidopsis, MYB49 interacts with ABSCISIC ACID-INSENSITIVE5 (ABI5) to regulate ABA-repressed cadmium accumulation [69].Here, we showed that exogenous ABA promoted anthocyanin accumulation in in vitro-grown blueberry seedlings (Figure 1).The promoter sequences of all R2R3-MYB genes from subgroups 5, 6, and 7 and the ASGs (C4H, 4CL7, CHI, and LDOX) contain one or more ABREs, and all the ASGs contain more than three MYB binding sites (Figure S3).The expression of ABF genes was upregulated under ABA treatment and also positively or negatively correlated with that of MYBs and ASGs.At the same time, MYBs also positively or negatively correlated with ASGs in expression levels under ABA treatment.These findings indicate that ABFs, MYBs, and ASGs form a regulatory network to co-regulate ABA-induced anthocyanin biosynthesis.
Most studies showed that R2R3-MYB genes from subgroup 6 control anthocyanin biosynthesis [19,21].In this study, the expression of MYBA1 and MYBA2 genes (subgroup 6) were upregulated under exogenous ABA treatment (Figure 5).In which, the MYBA2 gene expression was upregulated, was positively correlated with the expression of ASGs under UV-B radiation and played an important role in the UV-B-induced anthocyanin biosynthetic pathway [70].MYBA1 has been reported as VcMYB437and the expression of VcMYB437 was higher in the skin of ripe blueberry fruit than in the pulp, implying the potential roles in anthocyanin biosynthesis [71].VcMYBA (probably homologous to MYBA1) promoted anthocyanin accumulation by heterologous expression in Nicotiana benthaminana [21].Our study found that MYBA1 expression was significantly correlated with that of 4CL7 and VcUFGT, At the same time, 4CL7 and VcUFGT contain multiple MYB binding sites in their promoter sequences (Figure 6 and Figure S3).However, the expression of VcMYBA2 was not significantly correlated with that of ASGs.Thus, MYBA1 is likely to play a core role in the network of ABA-induced anthocyanin biosynthesis; however, its authentic roles still need to be documented in future studies.

Figure 1 .
Figure 1.ABA induces anthocyanin accumulation in in vitro-grown blueberry seedlings.Representative phenotypes (a) and anthocyanin contents (b) of in vitro-grown blueberry seedlings under 100 µM ABA treatment and without ABA treatment for 0 h, 6 h and 12 h.Values are means ± standard error (SE) of three independent biological replicates; different letters represent significant differences at p < 0.05 among samples by one-way ANOVA and Tukey's test.Bar, 1 cm.

Figure 1 .
Figure 1.ABA induces anthocyanin accumulation in in vitro-grown blueberry seedlings.Representative phenotypes (a) and anthocyanin contents (b) of in vitro-grown blueberry seedlings under 100 µM ABA treatment and without ABA treatment for 0 h, 6 h and 12 h.The top and bottom plants showed 100 µM ABA treatment and without ABA treatment, respectively.Values are means ± standard error (SE) of three independent biological replicates; different letters represent significant differences at p < 0.05 among samples by one-way ANOVA and Tukey's test.Bar, 1 cm.

Figure 2 .
Figure 2. Differentially expressed genes in response to ABA treatment in blueberry, as revealed by transcriptome deep sequencing (RNA-seq).Volcano plots showing all the genes at 6 h (a) or 12 h (b) compared to 0 h under ABA treatment.Each point represents a single gene.The red and green dots indicate upregulated and downregulated differentially expressed genes (DEGs), respectively, and the black dots indicate genes with unchanged expression.(c) Venn diagram showing the number of DEGs across pairwise comparisons.Pink, DEGs for 0 h vs. 6 h; light blue, DEGs for 0 h vs. 12 h; purple, DEGs for 0 h vs. 6 h and 0 h vs. 12h.

Figure 2 .
Figure 2. Differentially expressed genes in response to ABA treatment in blueberry, as revealed by transcriptome deep sequencing (RNA-seq).Volcano plots showing all the genes at 6 h (a) or 12 h (b) compared to 0 h under ABA treatment.Each point represents a single gene.The red and green dots indicate upregulated and downregulated differentially expressed genes (DEGs), respectively, and the black dots indicate genes with unchanged expression.(c) Venn diagram showing the number of DEGs across pairwise comparisons.Pink, DEGs for 0 h vs. 6 h; light blue, DEGs for 0 h vs. 12 h; purple, DEGs for 0 h vs. 6 h and 0 h vs. 12h.

Figure 3 .
Figure 3. KEGG pathway enrichment analysis and functional classification of differentially expressed genes in blueberry under ABA treatment.Analysis of the top 20 significantly enriched DEGs after 6 h (a), 12 h (b), or 6 and 12 h (c) of ABA treatment compared to 0 h.The number of DEGs is indicated by the size of the circle.The solid upright and inverted triangles represent upregulated and downregulated DEGs, respectively.KEGG functional classification of DEGs after 6 h (d), 12 h (e), or 6 and 12 h (f) of ABA treatment compared to 0 h.

Figure 3 .
Figure 3. KEGG pathway enrichment analysis and functional classification of differentially expressed genes in blueberry under ABA treatment.Analysis of the top 20 significantly enriched DEGs after 6 h (a), 12 h (b), or 6 and 12 h (c) of ABA treatment compared to 0 h.The number of DEGs is indicated by the size of the circle.The solid upright and inverted triangles represent upregulated and downregulated DEGs, respectively.KEGG functional classification of DEGs after 6 h (d), 12 h (e), or 6 and 12 h (f) of ABA treatment compared to 0 h.

Figure 4 .
Figure 4. Differentially expressed genes of plant hormone signal transduction KEGG pathways in blueberry under ABA treatment.Green indicates low expression, and red indicates high expression, based on log2 (FPKM) values.

Figure 4 .
Figure 4. Differentially expressed genes of plant hormone signal transduction KEGG pathways in blueberry under ABA treatment.Green indicates low expression, and red indicates high expression, based on log 2 (FPKM) values.

Figure 6 .
Figure 6.Proposed model of the ABF-regulated anthocyanin biosynthetic pathway under ABA treatment.Green squares indicate low expression and red squares indicate high expression, based on log 2 (FPKM).Solid lines indicate the known biosynthetic pathway; dotted lines indicate the putative signaling pathway based on co-expression analysis.Arrows indicate positive correlation; short lines indicate negative correlation.Red and green text indicates upregulated and downregulated DEGs, respectively.

Table 1 .
Summary of sequencing data.

Table 1 .
Summary of sequencing data.