White loose blueberry (Vaccinium corymbosum) calli were obtained from the cultivar ‘Northland’, which was cultured on a modified woody plant medium with Murashige and Skoog vitamins containing 3.0 mg/L 2,4-dichlorophenoxyacetic acid (2,4-D), 30 g/L sucrose, and 7 g/L agar, pH 5.4 ± 0.2. All calli were cultured under a 16-h light/8-h dark photoperiod at 25°C and subcultured every 21 days. UV-B was applied by means of narrow band lamps (TL20/01; 311-nm Philips, Netherlands) positioned above the calli at the height of about 10 cm. The calli were harvested right before (0 h) and after treatment consisting of 1, 3, 6, 12, and 24 h of UV-B radiation, frozen in liquid nitrogen and stored at –80°C. All samples were collected as three independent biological replicates and were labeled 0h_1, 0h_2, 0h_3, 1h_1, 1h_2, 0h_3, 3h_1, 3h_2, 3h_3, 6h_1, 6h_2, 6h_3, 12h_1, 12h_2, 12h_3, 24h_1, 24h_2, and 24h_3, respectively, for RNA-seq analysis. For reverse transcription quantitative PCR (RT-qPCR) and measurements of flavonol, anthocyanin, and proanthocyanidin contents (using the 0 h and 24 h samples), each biological replicate was assessed as three technical replicates.

Total RNAs were extracted with TRIzol reagent (Invitrogen, USA) and RNA concentration and purity was measured using a NanoDrop 2000 (Thermo Fisher Scientific, Wilmington, DE). RNA integrity was assessed using an RNA Nano 6000 Assay Kit for the Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, USA). Sequencing libraries were generated using a NEBNext UltraTM RNA Library Prep Kit for Illumina (NEB, USA). The libraries were sequenced on an Illumina HiSeq 2500 platform, and 150 bp paired-end reads were generated. Raw reads in fastq format were processed through in-house perl scripts to obtain clean data. Q20, Q30, GC-content, and the sequence duplication level of the clean reads were calculated. The clean reads were mapped to the reference Vaccinium corymbosum cv. Draper V1.0 genome sequence (https://www.vaccinium.org/genomes) using Hisat2 software.

Gene expression levels were quantified as fragments per kilobase of transcript per million fragments mapped (FPKM) values. Differentially expressed genes (DEGs) resulting from the comparison of the 0 h sample to the 1, 3, 6, 12, and 24 h samples were identified using the criteria of absolute log₂(fold change) ≥1 and a false discovery rate (FDR) < 0.01 by DEGSeq2. Gene function was annotated based on the following databases: Clusters of Orthologous Groups (COG) (Tatusov et al., 2000), Gene Ontology (GO) (Ashburner et al., 2000), Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa et al., 2004), Eukaryotic Orthologous Groups (KOG) (Koonin et al., 2004), NCBI non-redundant protein sequences (NR) (Deng et al., 2006), Protein family (Pfam) (Finn et al., 2014), Swiss-Prot (A manually annotated and reviewed protein sequence database) (Apweiler et al., 2004) and evolutionary genealogy of genes: Non-supervised Orthologous Groups (eggNOG) (Huerta-Cepas et al., 2017). KEGG pathway enrichment analysis of DEGs was implemented with KOBAS software (Mao et al., 2005). The phenylpropanoid and flavonoid biosynthetic pathways were mapped against the phenylpropanoid, flavonoid, flavonol, and anthocyanin KEGG pathways using the DEGs identified above. Heatmap representations of gene expression levels were drawn using log₁₀(FPKM) with the Tbtools (v1.098761) software (Chen et al., 2020).

WGCNA was performed with the WGCNA package in R (Langfelder and Horvath, 2008). First, the WGCNA algorithm assumes that the gene network follows a scale-free distribution, defines the gene co-expression correlation matrix and the adjacency function formed by the gene network, and then calculates the correlation coefficients of different nodes, based on which WGCNA builds a hierarchical clustering tree. The different branches of the clustering tree represent different gene modules. The co-expression degree of genes in significant individual modules is high, while the co-expression degree of genes belonging to different modules is low. Genes with different expression levels were assigned to various modules via a dynamic tree cut. There were at least 30 genes per co-expression module. Correlations among various modules were calculated using 0.25 as the similarity threshold. For genes in each module, KEGG pathway enrichment analysis was conducted to reveal the biological functions of each module. The genes within the WGCNA kMEblue module were selected to screen for enrichment of the plant hormone signal transduction pathway, transcription factor genes, and flavonoid biosynthetic pathway genes, and to draw the corresponding expression heatmap and regulatory network of the flavonoid pathway. Pearson’s correlation coefficients were compared using SPSS 19.0 software.

Total RNA was extracted from each sample using an RNA Extraction Kit (Sangon Biotech, Shanghai, China). RT-qPCR was performed on an ABI 7900HT real-time PCR system. Ten genes of interest (4CL2, VcCHI3, DFR, VcF3’5’H, VcF3H-2, VcLAR, MYB44, MYB114, VcMYBPA1, and ARF18) involved in flavonoid biosynthesis and belonging to the WGCNA kMEblue module were selected for analysis, using GAPDH (AY123769) as the reference transcript. Primer sequences are shown in Table S1 . The relative expression levels of each gene were calculated using the 2^–ΔΔCt method.

The contents of flavonoid compounds including flavonols, anthocyanins, and proanthocyanidins for calli at 0 h and 24 h UV-B treatment were determined as described by Yang et al. (2022a). Each sample was ground to a powder and 0.5 g was extracted in 5 mL 80% (v/v) methanol at 4°C for 2 h to isolate the flavonols. The mixtures were centrifuged (8,000 g, 10 min, 4°C), 1 mL of the supernatant was removed and mixed with 1 mL methanol, and then 0.1 mL of 10% (w/v) aluminum chloride, 0.1 mL 1 M KOAc, and 2.8 mL water were added and the mixtures was incubated for 30 min at 25°C. Rutin was used as a master standard and the absorbance at 415 nm was measured. To measure total anthocyanin contents, samples (0.5 g) were extracted in 3 mL of 1% HCl in methanol and incubated at 4°C for 16 h. After centrifugation (8,000 g, 10 min, 4°C), a 2-mL supernatant was diluted with 2 mL water and the absorbance was measured at 530 nm and 650 nm. Total anthocyanin content was calculated using a previously published formula (Rabino and Mancinelli, 1986). Proanthocyanidin was detected using the DMACA method. Briefly, samples (0.5 g) were extracted in 5 mL of 70% (v/v) acetone containing 0.1% (w/v) ascorbic acid at 4°C for 30 min. After centrifugation (8,000 g, 10 min, 4°C), a 3-mL supernatant aliquot was extracted with 3 mL of ether at −20°C for 1 h. Then, 2 mL of the lower phase of the extracted liquid was removed and mixed with 1 mL of methanol and 0.5 mL 2% (w/v) DMACA solution. The mixture was incubated for 20 min at 25°C, and the absorbance was then measured at 643 nm. Catechin was used as the master standard (Wang et al., 2017).

To reveal the regulatory network underlying the blueberry response to UV-B radiation, we performed an RNA-seq analysis using total RNA extracted from blueberry calli exposed to UV-B radiation for 1, 3, 6, 12, or 24 h. Table 1 summarizes the details of all RNA-seq samples; we obtained 5.85–8.66 Gb of clean bases for each sample, with a Q30 ranging from 93.70 to 94.97% and a GC content ranging from 46.35 to 46.90%. We successfully mapped approximately 90.88 to 91.90% of all clean reads per sample to the blueberry reference genome.

To investigate changes in gene expression under UV-B radiation, we identified differentially expressed genes (DEGs) between the control samples (0 h) and each time point of the UV-B treatment, which returned a total of 16,899 DEGs ( Table S2 ). We detected 2,706 DEGs in 0h_vs_1h, of which 2,428 were upregulated and only 278 were downregulated. The number of DEGs increased gradually with longer UV-B exposure of up to 6 h, with downregulated DEGs remaining scarcer than upregulated DEGs. The number of DEGs further increased to 10,146 in 0h_vs_12h and to 10,718 in 0h_vs_24h, respectively, but now with more downregulated DEGs than upregulated DEGs ( Figure 1A ). A Venn diagram comparing all the lists of DEGs shows that 806 DEGs are shared by all pairwise comparisons ( Figure 1B ). These results indicate that the number of DEGs increases with the duration of UV-B radiation treatment.

We annotated the function of DEGs in response to UV-B radiation with the COG, GO, KEGG, KOG, NR, Pfam, Swiss-Prot, and eggNOG databases ( Table S3 ). During UV-B treatment, the number of annotated DEGs increased, starting at 2,596 annotated DEGs out of 2,706 after 1 h, and rising to 9,953 annotated DEGs out of 10,718 after 24 h. The number of annotated DEGs also differed by database, with the highest number obtained with the NR database and the lowest with the COG database.

We then undertook a functional classification of DEGs that were annotated in the COG, eggNOG, and KOG databases ( Table S4 ). We observed that DEGs are most enriched in terms related to signal transduction mechanisms, secondary metabolite biosynthesis, transport and catabolism, as well as transcription, especially in the eggNOG and KOG functional annotations. Furthermore, we performed a GO classification of shared annotated DEGs for the three categories biological process, cellular component, and molecular function ( Figure S1 ; Table S5 ). In the biological process category, cellular processes and metabolic processes were the most enriched; in the cellular components category, cell and cell part were the most enriched; and in the molecular functions category, binding and catalytic activity were the most enriched.

To further elucidate the biological functions of the DEGs, we carried out a KEGG pathway enrichment analysis. We identified 136 KEGG pathways whose encoded proteins are affected during UV-B treatment in blueberry calli ( Table S6 ). Using a Q-value < 0.01 as cutoff, we determined that the KEGG pathway circadian rhythm-plant pathway was significantly affected at all time points of UV-B treatment ( Figure 2 ). Similarly, the plant hormone signal transduction pathway was significantly enriched at almost all stages, suggesting that plant hormones are involved in the response to UV-B stress. We also established that many DEGs are involved in the biosynthesis of phenylpropanoid-derived compounds. Among them, genes participating in the isoflavonoid biosynthetic pathway were induced after 1 and 3 h UV-B treatment ( Figures 2A, B ). Likewise, genes from the phenylpropanoid biosynthetic pathway began to be enriched after 3 h of UV-B treatment ( Figures 2B–E ), followed by genes from the anthocyanin and flavonoid biosynthetic pathways, starting at 6 h of UV-B treatment ( Figures 2C–E ). We conclude that genes involved in the biosynthesis of phenylpropanoid-derived compounds exhibit a precise spatial pattern: the expression of isoflavonoid biosynthetic genes is induced first, followed by genes involved in the phenylpropanoid pathway, and finally genes involved in the anthocyanin and flavonoid biosynthetic pathways. In agreement with this result, the top 20 enriched pathways identified when comparing the control samples (0 h) to any sample exposed to UV-B included anthocyanin biosynthesis, phenylpropanoid biosynthesis, and flavonoid biosynthesis ( Figure 2F ).

To elucidate the effects of UV-B radiation on the biosynthesis of flavonoid compounds, we measured the flavonoid contents of blueberry calli treated by UV-B radiation for 24 h ( Figure 3A ). The contents of flavonols, proanthocyanidins, and anthocyanins increased 3.1-, 7.0-, and 2.1-fold, respectively, after 24 h of UV-B irradiation relative to the control samples ( Figures 3B–D ). Thus, UV-B radiation promotes the accumulation of flavonoid compounds.

We screened for genes involved in the phenylpropanoid, flavonoid, flavonol, and anthocyanin KEGG pathways among the DEGs. We identified 73 genes encoding 32 types of enzymes, based on Enzyme Commission (EC) annotated data ( Table 2 ). We then plotted the expression levels of these genes as a heatmap. Again, the genes selected here clustered as a function of the duration of UV-B exposure, with the 0, 1, and 3 h treatments forming one group and the 6, 12, and 24 h treatments forming another ( Figure S2 ). Notably, the later induction of these DEGs suggested that the accumulation of phenylpropanoid metabolites mainly occurs after a longer UV-B exposure.

We mapped the 73 DEGs identified above onto the KEGG phenylpropanoid and flavonoid biosynthetic pathways ( Figure 4 ). Among them, only three genes were induced after 1 h of UV-B radiation, and many genes were induced after 3 h. Most of the genes, including phenylalanine ammonia lyase (PAL1), PAL3, 4-coumarate CoA ligase (4CL2), chalcone synthase synthase (CHS), chalcone isomerase (CHI3), VcFLS, and VcUFGT, were induced under UV-B radiation, and their expression reached the highest levels at 24 h ( Table 2 ). Most genes involved in the phenylpropanoid biosynthetic pathway were rapidly upregulated, especially PAL1 and 4CL2 from 6 h of UV-B radiation treatment onwards. By contrast, 4CL6 was rapidly downregulated from 3 h of UV-B radiation onwards. Most genes involved in the flavonoid biosynthetic pathway were upregulated; CHS expression was induced from 3 h of UV-B radiation treatment onwards, reaching a log₂(FC) value of 6.33 at 24 h. The expression of F3’5’H, F3’H, and FLS, which are involved in the flavonol biosynthetic pathway, was upregulated under UV-B radiation. Among genes involved in the anthocyanin biosynthetic pathway, VcUFGT expression was induced from 1 h of UV-B exposure onwards, and then rapidly upregulated from 3 h onwards, reaching the highest level at 24 h (7.86 for log₂(FC) value). VcLAR, which is involved in the proanthocyanin biosynthesis pathway, was also upregulated from 6 h onwards, reaching a log₂(FC) value of 2.13 at 24 h ( Table 2 ). These results indicate that the expression of most genes involved in phenylpropanoid metabolite biosynthesis is induced under UV-B radiation and that PAl1, 4CL2, CHS, VcFLS, VcUFGT, and VcLAR play important roles in UV-B-induced flavonoid accumulation.

To reveal the regulatory network underlying flavonoid metabolism under UV radiation, we conducted a weighted gene co-expression network analysis (WGCNA) using the DEGs defined above. WGCNA clustered the DEGs into three modules, namely, kMEblue, kMEbrown, and kMEturquoise, which contained 2,174, 783, and 2,171 genes, respectively ( Figures 5A, B ). KEGG enrichment analysis showed that the kMEblue module comprises many of the genes involved in the phenylpropanoid, flavonoid, and anthocyanin metabolism pathways ( Figure S3 ). We thus characterized the connection between these genes and those associated with plant hormone signal transduction, or encoding transcription factors by looking for their annotations in the KEGG, Swiss-Prot, and NR databases ( Table S7 ). We identified 37 genes fulfilling the above criteria, and their expression levels fell into two groups, as evidenced by a heatmap ( Figure 5C ). Genes from group I were downregulated, whereas genes from group II were upregulated under UV-B radiation; the latter group included all flavonoid metabolism genes. We also identified ten MYB transcription factor genes in the kMEblue module; these genes belonged to subgroups 4 (VcMYB4a), 11 (MYB102), 8 (MYB20), 5 (MYBA2.1 and VcMYBPA), 2 (VcMYB14), 7 (MYBF), 6 (MYBA2 and VcMYB144), and 22 (MYB44) ( Figure S4 ).

To investigate the functions of these MYB transcription factor genes in the flavonoid pathway, we calculated the Pearson’s correlation coefficients (r) between their expression levels and those of genes involved in flavonoid biosynthesis during UV-B radiation ( Table 3 ). We observed that the expression of MYBA2, MYBF, MYB102, MYBPA1, MYBPA2.1, and MYB114 is positively correlated with that of most genes involved in the flavonoid pathway. Conversely, the expression levels of MYB20, VcMYB14, MYB44, and VcMYB4a were negatively correlated with those of most flavonoid metabolism genes. We also observed a positive correlation between the expression levels of bHLH74, VcMYB20, and VcMYB44 and between bHLH25 and VcMYB4a, and a negative correlation between the expression levels of bHLH74 and MYB114, and between the expression levels of bHLH25 and those of MYB114, MYBA2, and MYBF. The expression of bHLH74 and bHLH25 was also negatively correlated with the expression levels of genes involved in flavonoid metabolism.

To elucidate the roles of plant hormones in regulating the flavonoid biosynthesis pathway, we calculated the correlation between MYB expression levels and the expression of plant-hormone-related genes ( Table 3 ). In the auxin biosynthetic pathway, the expression levels of indole-3-acetic acid inducible 14 (IAA14) and auxin-response factor 18 (ARF18) were negatively correlated with those of MYB114 and positively correlated with those of MYB20 and MYB44. In the gibberellic acid biosynthetic pathway, chitin-inducible gibberellin-responsive protein 1 (CIGR1) expression was positively correlated with that of MYB20 and MYB44 and negatively correlated with that of MYB114. Scarecrow-like transcription factor (PAT1) expression was positively correlated with that of MYB20, MYB44, and VcMYB14 and negatively correlated with that of VcMYBPA1. The expression level of Ga insensitive dwarf 2 (GID2) was positively correlated with those of CIGR1 and PAT1. However, among genes involved in the brassinosteroid biosynthetic pathway, only brassinosteroid-signaling kinase KINASE 7 (BSK7) expression levels appeared to be positively correlated with MYB20 and MYB44 expression and negatively correlated with MYB114 expression ( Table 3 ). Figure 6 summarizes the regulatory network of UV-B-induced flavonoid biosynthetic pathway through hormone signal transduction and transcriptional regulation pathways based on DEGs identified from the WGCNA kMEblue module in this study.

To validate the accuracy and reliability of the RNA-seq data, we selected ten genes of interest involved in flavonoid biosynthesis from the WGCNA kMEblue module for RT-qPCR analysis ( Figure S5 ). A linear regression analysis showed a positive correlation between the RNA-seq and RT-qPCR results, with correlation coefficients (R ²) of 0.8926 (0 h vs 1 h), 0.8338 (0 h vs 3 h), 0.8849 (0 h vs 6 h), 0.8856 (0 h vs 12 h), and 0.8336 (0 h vs 24 h). This observation confirmed that the RNA-seq data in this study are accurate and reliable.