Skip to main content
Download PDF

Genome-wide analysis of MYB transcription factors of Vaccinium corymbosum and their positive responses to drought stress

BMC Genomics volume 22, Article number: 565 (2021) Cite this article

Abstract

Background

Blueberry (Vaccinium corymbosum L.) is an important species with a high content of flavonoids in fruits. As a perennial shrub, blueberry is characterized by shallow-rooted property and susceptible to drought stress. MYB transcription factor was reported to be widely involved in plant response to abiotic stresses, however, the role of MYB family in blueberry responding to drought stress remains elusive.

Results

In this study, we conducted a comprehensive analysis of VcMYBs in blueberry based on the genome data under drought stress, including phylogenetic relationship, identification of differentially expressed genes (DEGs), expression profiling, conserved motifs, expression correlation and protein-protein interaction prediction, etc. The results showed that 229 non-redundant MYB sequences were identified in the blueberry genome, and divided into 23 subgroups. A total of 102 MYB DEGs with a significant response to drought stress were identified, of which 72 in leaves and 69 in roots, and 8 differential expression genes with a > 20-fold change in the level of expression. 17 DEGs had a higher expression correlation with other MYB members. The interaction partners of the key VcMYB proteins were predicted by STRING analysis and in combination with physiological and morphological observation. 10 key VcMYB genes such as VcMYB8, VcMYB102 and VcMYB228 were predicted to be probably involved in reactive oxygen species (ROS) pathway, and 7 key VcMYB genes (VcMYB41, VcMYB88 and VcMYB100, etc..) probably participated in leaf regulation under drought treatment.

Conclusions

Our studies provide a new understanding of the regulation mechanism of VcMYB family in blueberry response to drought stress, and lay fundamental support for future studies on blueberry grown in regions with limited water supply for this crop.

Background

Blueberry is an important perennial shrub within the genus Vaccinium of the Ericaceae and its fruit is rich in anthocyanin, which has significant value to human health [1, 2]. In recent years, drought has become a major threat to crop production with the climate change [3]. Blueberries are shallow-rooted plants with slender and underdeveloped roots [4], which is therefore more vulnerable to drought stress. Across the globe, drought stress led to reduced blueberry yield by about 25–30 % [5], which has evidently become a crucial factor in limiting the blueberry supply and production chain. Nevertheless, it is still unclear about the drought-tolerance mechanisms for blueberry seedlings.

Plants are universally confronted with many extreme environmental events during its lifetime. In order to survive, plants have evolved a set of elaborate and complicated self-regulation networks. Transcription factors (TFs) usually function as key regulators in plants responding to abiotic and biotic stresses [6,7,8]. Among them, MYB TF family was reported to be widely involved in a range of abiotic stresses [9].

MYB TF family exists in all eukaryotes and the MYB domain is generally composed of 1 ~ 4 imperfect amino acid sequence repeats (R) with about 52 amino acids in plants. MYB TFs are classified into four subfamilies (1R-MYB, R2R3-MYB, 3R-MYB, and a typical-MYB) based on the number and position of these repeats in their DNA-binding domains [10, 11]. R2R3-MYB is the main type of TF identified in the plant up to now, with more than 130 members in Arabidopsis and 90 in Oryza sativa [12]. To date, studies on MYB TF family responding to abiotic stress have been reported in many species. A total of 196, 155 and 265 MYB genes were identified in Arabidopsis, Oryza sativa, and M. truncatula, respectively [11,12,13]. In Arabidopsis, 51 % of AtMYB genes are up-regulated and 41 % are down-regulated under drought condition [14]; in Oryza sativa, 65 % of the expressed MYB genes were differentially regulated by drought stress in seedlings [15]. Multiple complex biological processes are involved in MYB TFs-regulated drought tolerance in plants such as ROS signaling pathway, stomatal regulation, cell strtuctue and component regulation and phytohormones-mediated signaling pathway. For example, in Arabidopsis, AtMYB44 responds to drought stress by participating in stomatal regulation and ROS accumulation [16,17,18]. AtMYB15 and AtMYB2 act as positive regulators under drought stress by activating the dehydration response genes such as AtRD22; overexpression of AtMYB52 can improve drought tolerance of transgenic Arabidopsis by affecting the cell wall structure [19, 20]. AtMYB96 regulated cuticular wax biosynthesis and contributed to drought resistance in Arabidopsis [21, 22]. In woody palnts, MdS1MYB1 particapated in positive regulation of drought stress by activating the auxin signaling pathway in apples (Malus × domestica) [23]. In addition, MYB TFs can simultaneously respond to a variety of adverse circumstances. In soybean (Glycine max), GmMYB118 maintains cell homeostasis by regulating osmotic and oxidative substances, thus improving the tolerance of soybean to drought and salt stress [24]. Overexpression of PtsrMYB in tobacco (Nicotiana tabacum L.) confers enhanced salt, dehydration and cold tolerance with lower levels of ROS in transgnic tabacco compared to wild type plants [25].

In this study, genome-wide analysis of the MYB family in blueberry under drought stress was performed. Two hundred and twenty-nine non-redundant MYB sequences were identified and divided into 23 subgroups. A total of 102 MYB DEGs responding to drought stress were identified. Furthermore, 23 important potential VcMYBs were determined based on the expression correlation analysis and 20-fold DEGs identification. Among them, 10 key genes were probably involved in the ROS pathway, and 7 key genes were likely to be involved in leaf regulation responding to droughts stress. The key VcMYBs and the possible interation proteins identified in this study provide the candidate genes for genetic improvement and lay the basis for the future research to explore the molecular mechanisms of blueberry grown in regions with limited water supply for this crop.

Results

Identification and classification of VcMYBs from blueberry

299 cDNA sequences were identified with full-length ORFs as putative MYB genes from the transcriptome database on the free online platform of Majorbio Cloud Platform (https://cloud.majorbio.com/), and all potentially redundant VcMYBs sequences were discarded by ClustalX software. We performed a BLAST search against the blueberry genome database (https://www.vaccinium. org/) to verify the validity of the 299 cDNA sequences, then designated them as VcMYB1 to VcMYB229 for further analysis (Additional file 1). Based on the number of repeat units at the MYB-domain, the screened VcMYBs were classified into four groups, namely, “1R-MYB,” “R2R3-MYB,” “3R-MYB,” and “a typical-MYB”. We found that the number of the four types of VcMYB genes was 19 (8.30 %), 191 (83.41 %), 10 (4.37 %), and 9 (3.93 %), respectively. The subfamily R2R3-MYB had the largest number of MYB proteins in blueberry (Table 1).

Table 1 The MYB-domain analysis of VcMYB genes based on GRAVY and molecular weight
Full size table

The physicochemical characteristics analysis of these VcMYB proteins showed that the 229 predicted VcMYB proteins varied from 71 (VcMYB40) to 1686 (VcMYB137) amino acids in length, with an average length of 328 aa, suggesting the functional diversity for VvMYB family in blueberry. The molecular weight of these VcMYB proteins ranged from 7.96 kD for VcMYB40 to 191.08 kD for VcMYB137. The mean value of GRAVY and pI was − 0.69 and 7.12, respectively (Table 1, Additional file 1 and Additional file 2).

Phylogenetic analysis of VcMYB genes

In order to predict the potential function of VcMYB, we performed phylogenetic reconstruction between 229 VcMYB genes and 133 AtMYB genes using the NJ method (Fig. 1 and Additional file 1). An unrooted composite phylogenetic tree was created, in which VcMYBs were classified into 23 major groups (C1-C23) with supported bootstrap values. We found that subgroup C11 had the largest number of MYB proteins with 38, whereas C21 has the lowest number with only 4 members. In Arabidopsis, 126 AtMYB genes were divided into 25 subfamilies [26, 27], and these defined clades were compared and labeled with the composite evolutionary tree in our study (Fig. 1). In the evolutionary tree, 14 subgroups (C1, C2, C3, C8, C10, etc.) in our study were matched with 23 Arabidopsis subgroups (S1, S2, S3, S4, S5, etc.) reported in previous study [27]. However, S8 and S17 Arabidopsis subgroups were not retrieved from the composite evolutionary tree and no AtMYBs genes in Arabidopsis were matched with the subgroups C15 and C23 in blueberry (Fig. 1), suggesting that some divergence occurred in MYB family among different plant species during the long evolutionary process.

Fig. 1
figure 1

Phylogenetic analysis of MYB proteins between blueberry and Arabidopsis. Complete sequence alignments of 229 MYB amino acid sequences in blueberry and 133 MYB amino acid sequences in Arabidopsis were performed to construct the phylogenetic tree. The black filled circle denotes VcMYB genes; the hollow circle denotes AtMYB genes. The red font indicates differentially expressed genes

Full size image

Identification and cluster analysis of differently expressed VcMYB genes under drought stress

The expression level of each gene was calculated according to the FPKM (Fragments per kb per million reads method) on the free online platform of Majorbio Cloud Platform and RSEM (http://deweylab.biostat.wisc.edu/rsem/) was used to quantify gene abundances. The differential expression analyses of VcMYBs were performed with the R package DESeq2 (http://bioconductor.org/ packages/stats/bioc/DESeq2/). As shown in Fig. 2, 69 and 72 differentially expressed genes (DEGs) were screened out (at a fold ratio with 2 and P-adjusted value with 0.05) in leaf and root, respectively (Fig. 2 A and B). A total of 102 DEGs were identified in leaves and roots, among which 39 DEGs were found to be co-expressed both in leaves and roots, whereas 33 and 30 genes specially expressed in leaves and roots, respectively (Fig. 2 C).

Fig. 2
figure 2

Identification and cluster analysis of differently expressed VcMYB genes in blueberry under drought stress. Venn diagram of DEGs in response to drought in leaves (A) and in roots (B). (C) represents the DEGs co-expressed in leaf and root. (D) and (E) represent the hierarchical cluster and heatmap of DEGs in leaf or root, respectively. Red indicates up-regulated genes, and blue indicates down-regulated genes. CK_L represents control in leaf; CK_R represents control in root; MD_L represents moderate drought treatment in leaf; MD_R represents moderate drought treatment in root; SD_L represents severe drought treatment in leaf; SD_R represents severe drought treatment in root; MYB_L represents the DEGs in the leaf of blueberry; MYB_R represents the DEGs in root of blueberry. The diamond indicates 20-fold VcMYB DEGs. L1 represents the genes with high expression of MD_L. L2 represents the genes with high expression of CK_L. L3 represents the genes with high expression of MD_L and SD_L. L4 represents the genes with high expression of SD_L. R1 represents the genes with high expression of MD_R. R2 represents the genes with high expression of CK_R. R3 represents the genes with high expression of MD_R and SD_R. R4 represents the genes with high expression of SD_R

Full size image

In order to facilitate the analysis of these DEGs, a heatmap was constructed based on log10 FPKM (Fig. 2D and E). Under drought stress, the expression of DEGs in leaves and roots can be clustered into 2 groups, namely the up-regulated group (L1, L3, L4, R1, R3 and R4) and the down-regulated group (L2 and R2). The number of up-regulated DEGs was 43 in leaves and 48 in roots, and down-regulated DEGs in leaves and roots was 29 and 21, respectively. Meanwhile, eight 20-fold differential genes were identified, and clustered into C13 (VcMYB2, VcMYB8, VcMYB14, VcMYB29, VcMYB48, VcMYB102, VcMYB108 and VcMYB227) and C15 (VcMYB29) subfamilies (Fig. 1).

In order to verify the accuracy of RNA-seq data and the reliability of the data analysis, we randomly selected 16 VcMYB TFs for qRT-PCR analysis in leaves and roots. We found that the qRT-PCR results were basically consistent with the RNA-seq data, except that the expression level of VcMYB2 in leaves showed slight difference. Our results indicated that the RNA-Seq data in this study were accurate and reliable (Additional file 3).

Phylogenetic and conserved motifs analysis of VcMYB DEGs

In order to further analyze the characteristics of these 102 DEGs in 23 subgroups, the conserved motifs of VcMYB proteins were analyzed using the MEME online program. 10 motifs, named motif 1–10 (Additional file 4) were determined for these MYB DEGs proteins (Fig. 3), among which motif 1, 6 and 7 were identified as the core conserved domains and constituted the SANT domain of MYB (Additional file 4). Among these conserved motifs, motif 3 and 7 correspond to the most conserved genes with 88, and motif 9 corresponds to the least conserved genes with only 6 based on the prediction. Motif 6 and 8 contained 50 amino acids, but motif 7 contained only 8 amino acids (Additional file 4). It should be noted that no differential genes were identified in C12, C19 and C21 subfamilies. VcMYB proteins in the same cluster subgroup were observed to have the similar motif composition (Fig. 3). For example, subgroup C20 contained motif 1, 3 and 7, indicating the functional similarities in the same subgroup [28, 29].

Fig. 3
figure 3

Phylogenetic tree and conserved motifs analysis of VcMYB DEGs. The phylogenetic tree was constructed using MEGA6.0 software, and was decorated by TBtools software. MEME (v4.11.1) online program was employed to predict the motif

Full size image

In order to elucidate the sequence feature of these VcR2R3-MYB proteins in blueberry, multiple sequence alignment was perfomed by using Clustal X. 71 highly conserved VcR2R3-MYB amino acid regions were identified among 102 differently expressed VcMYB proteins (such as VcMYB2, VcMYB8, VcMYB14 etc.) (Additional file 5). Subsequently, ggseqlogo was used to generate sequence logos. The R2 and R3 MYB repeats of the VcR2R3-MYBs contained many conserved amino acids, such as the typical tryptophan (Trp), which is crucial to the sequence-specific binding of DNA[28]. Five conserved Trp residues were identified in the R2 and R3 repeats (Fig. 4 A and 4B). As with its counterparts in the Chinese Pear (Pyrus bretschneideri Rehd.) [28], the first conserved Trp residue in the R3 repeat was generally replaced by other amino acid (Fig. 4B), meanwhile, the second and third conserved Trp residue were identified in the R3 repeat in bluebrerry. Some other highly conserved amino acids were also discovered, such as Gly (G), Glu (E), Asp (D), Cys (C), Arg (R), Leu (L), Iie (I), Thr (T), Asn (N) and Lys (K).

Fig. 4
figure 4

The sequence logos of the R2 (A) and R3 (B) VcMYB repeats. These logos were based on multiple full-length alignments of all blueberry R2R3-MYB domains. The bit score represents the information content for each position in the sequence. Asterisks represent the conserved residues that are identical among all R2R3-MYB domains, and triangles denote the typically conserved residues (Trp) in the R2R3-MYB domains. Gly (G), Glu (E), Asp (D), Cys (C), Arg (R), Leu (L), Iie (I), Thr (T), Asn (N) and Lys (K)

Full size image

GO annotation analysis of VcMYB DEGs

The potential functions of DEGs for VcMYB in blueberry were predicted by GO annotation analysis. In this study, the 102 VcMYB DEGs were assigned to 16 GO categories (Additional file 6 and Additional file 7). The results showed that 102 DEGs participated in the biological process (38), cellular component (7) and molecular function (21). The “metabolic process” gene (11) was the dominant category in the biological process category accounting for 16.92 %. VcMYB48, VcMYB88, VcMYB108, VcMYB229, VcMYB129 and VcMYB228 were categorized into responding to stimulus. In the cellular component category, VcMYB genes were mainly located in the “cell,” “organelle,” and “cell part.” Regarding the molecular function category, “binding” (14) was the most dominant group accounting for 21.54 % (Additional file 6).

GO functional enrichment analysis was further carried out by Goatools based on the Majorbio Cloud Platform. A total of 6 GO terms were considered to be significantly enriched among these DEGs, in which 4 and 2 GO terms belonged to “biological process” (BP) and “molecular function” (MF), respectively (Fig. 5 A). By analyzing the genes involved in the 6 GO terms, we found that VcMYB108 participated in the most of GO terms (5), and VcMYB228, VcMYB48, VcMYB88, VcMYB229 and VcMYB129 participated in 4 GO terms, respectively (Fig. 5B).

Fig. 5
figure 5

Gene ontology enrichment analysis of DEGs. (A) Bubble diagram of DEGs based on GO functional-enrichment analysis. MF represents molecular function, and BP represents biological process. (B) String diagram of DEGs on the basis of important GO terms

Full size image

Expression correlation analysis of VcMYB DEGs in blueberry

The possible co-expression relationship between genes can be reflected by the correlation analysis of expression levels of genes. Therefore, we constructed two expression correlation networks of 72 DEGs in leaf and 69 DEGs in root in blueberry based on the spearman correlation algorithm. The results showed that 12 VcMYB genes were correlated in the expression correlation network of leaves and roots (Fig. 6 A and 6B), among which 7 genes were co-expressed including VcMYB88, VcMYB100, VcMYB105, VcMYB108, VcMYB129, VcMYB146 and VcMYB228 (Fig. 6 C). Notably, VcMYB88 and VcMYB228 had the largest node whether in root or leaf, indicating that the two genes were significantly correlated with other genes and may play a crucial role in blueberry response to drought stress.

Fig. 6
figure 6

Expression correlation analysis of DEGs in blueberry under drought stress. DEGs’ expression correlation network under drought stress in leaf (A) and in root (B). Lines represent expression correlation between genes, and each node represents an MYB gene. Large nodes represent that this gene has an expression correlation with a large number of other genes. The blue font indicates that this gene coexists in the expression correlation network of the leaf and the root. The black font indicates that this gene is specifically present in the expression correlation network of the leaf or the root. (C) Venn diagram of DEGs that were expressed in leaf or root

Full size image

10 VcMYB DEGs probably participated in drought response through ROS-mediated signaling pathway

According to the above screening, eight 20-fold differential genes including VcMYB2, VcMYB8, VcMYB14, VcMYB48, VcMYB29, VcMYB102, VcMYB108 and VcMYB207 were identified (Fig. 2). The expression correlation analysis of MYB TFs showed that 17 genes (Fig. 6) were significantly correlated in the expression correlation network of blueberry leaves and roots. Thus, based on the orthologs of Arabidopsis, we further predicted the potential interaction proteins of these 10 key VcMYB proteins by STRING analysis to explore the possible functioin of VcMYB protein in blueberry response to drought stress (Fig. 7 A). The results showed that multiple candidate key proteins interact with each other, i.e., MYB44 (homolog of VcMYB228) could interact with EIN2, MPK3, VIP1, WRKY33 and STZ (homolog of VcEIN2, VcMPK3, VcVIP1, VcWRKY33 and VcSTZ, respectively), MYB106 (homolog of VcMYB100 and VcMYB193) could interact with NOA1 (homolog of VcNOA1), MYB121 (homolog of VcMYB2, VcMYB8, VcMYB14, VcMYB48, VcMYB102, VcMYB108 and VcMYB227) could interact with AT5G47390 (homolog of VcKUA1) (Fig. 7 A and Additional file 8). Previous studies reported that overexpression of AtMYB44, AtEIN2, AtMPK3, AtVIP1, AtWRKY33, AtSTZ, AtMYB106, AtNOA1, AtMYB121 and AtKUA1 could enhance drought tolerance by regulating the accumulation of ROS and cell homeostasis [16, 18, 30,31,32,33]. Therefore, we speculate that the 10 key VcMYB DEGs were likely to be involved in drought response through ROS-mediated signaling pathway in blueberry.

Fig. 7
figure 7

Effect of drought stress on SOD and H2O2 of blueberry. (A) Protein-protein interaction analysis of 10 VcMYB DEGs by STRING. (B) The expression profile of VcMYB8, VcMYB102, VcMYB193 and VcMYB288 using qRT-PCR under drought stress. Each value is represented as the mean value ± standard error of three independent determinations, and different letters indicate significant differences at p < 0.05 by Duncan’s multiple range test. CK: control group; MD: moderate drought stress; SD: severe drought stress. (C) The expression heatmap of the interacting proteins in leaf and root. (D) In situ visualization of H2O2 accumulation by DAB staining of the blueberry leaf under drought stress. (E) SOD activity and H2O2 content in leaf or root. Each value is indicated as the mean value ± standard error of three independent determinations, and different letters represent significant differences at p < 0.05 by Duncan’s multiple range tests

Full size image

The qRT-PCR results showed that the relative expression level of VcMYB8, VcMYB102 and VcMYB228 increased significantly and VcMYB193 decreased significantly with the duration of drought stress in leaves and roots (Fig. 7B). Furthermore, the relative expression levels of the other 6 potentially key VcMYBs (VcMYB2, VcMYB14, VcMYB48, VcMYB100, VcMYB108 and VcMYB227) were also significantly induced under drought stress in leaves and roots (Additional file 3), implying the functional divergence and significant role of these VcMYBs genes in blueberry under drought stress. Under drought stress, among the DEGs proteins that could interact with key VcMYBs in blueberry, VcMPK3, VcWRKY33 and VcSTZ showed trends upward, VcNOA1, VcEIN2 and VcVIP1 showed the downward trend (Fig. 7 C).

Furthermore, we examined the change of superoxide dismutase (SOD) activity and the content of H2O2 in the leaf and root of blueberry under drought stress. The results indicated that with the duration of drought stress, the content of H2O2 in leaves and roots were significantly increased, and more and more reddish-brown plaques were observed in leaves (Fig. 7Dand E), indicating that drought stress increased the amount of oxygen released by H2O2. Nevertheless, the activity of SOD in leaves and roots increased significantly with the aggravation of drought treatment (Fig. 7E).

7 VcMYB DEGs probably participated in drought response through the leaf regulation pathway

In this study, based on the orthologs analysis of Arabidopsis, we predicted the interactions of 7 key VcMYB proteins (VcMYB41, VcMYB88, VcMYB100, VcMYB152, VcMYB171, VcMYB193 and VcMYB202) of blueberry by STRING analysis (Fig. 8 A). The results showed that many candidate proteins interact with each other, for example, MYB41 (homolog of VcMYB171 and VcMYB202) could interact with WAKL4 (homolog of VcWAKL4), MYB106 (homolog of VcMYB100 and VcMYB193) could interact with SHN1 and GL3 (homolog of VcSHN1 and VcGL3, respectively), and MYB63 (homolog of VcMYB41) could interact with NAC073, IRX3, IRX1, NAC010, VND7, NST1, KNAT7, NAC012 and NAC066 (homolog of VcNAC73, VcIRX3, VcIRX1, VcNAC10, VcNAC30, VcNST1, VcKNAT7, VcNAC12 and VcNAC66, respectively) (Fig. 8 A and Additional file 8). Previous reports showed that overexpression of the interacted proteins, such as AtMYB106, AtSHN1, AtGL3, AtMYB41, AtWALK4, AtMYB63 enhanced the drought tolerance of Arabidopsis by affecting leaf morphology[34], cell wall biosynthesis [19, 34, 35], wax biosynthesis [36, 37] and cuticle biosynthesis [38]. The leaf morphology and structure changes occurred due to the damaged biosynthesis of leaf cell wall, wax, and cuticle [39]. Thus, we speculated that the 7 key VcMYB DEGs might participate in drought response through the regulation of leaf morphology and structure in blueberry.

Fig. 8
figure 8

Morpho-anatomical trait and structural parameter changes of blueberry under drought stress. (A) Protein-protein interaction analysis of 7 VcMYB DEGs by STRING. (B) The expression profile of VcMYB41, VcMYB88, VcMYB100 and VcMYB171 under drought stress. Each value is represented as the mean value ± standard error of three independent determinations, and different letters indicate significant differences at p < 0.05 by Duncan’s multiple range test. (C) The expression heatmap of interacting proteins in leaf. (D) Morpho-anatomical trait changes of blueberry under drought stress. (E) Total thickness, upper epidermis thickness, lower epidermis thickness and tissue structure tense ratio of the leaf. Each value is indicated as the mean value ± standard error of three independent determinations. Different letters represent significant differences at p < 0.05 by Duncan’s multiple range test

Full size image

The qRT-PCR results showed that the relative expression level of VcMYB41, VcMYB88, VcMYB100 and VcMYB171 in leaves changed significantly with the duration of drought stress (Fig. 8B). Meanwhile, the relative expression levels of the other 3 VcMYBs (VcMYB152, VcMYB193 and VcMYB202) of the regulatory genes potentially involved in leaf morphology and structure were also significantly induced under drought stress in leaves (Fig. 7B, Additional file 3). The expression heatmap of these interaction proteins showed that these genes could respond to drought stress (Fig. 8 C), implying their possible role of these genes in response to drought stress. Moreover, we observed that the blueberry leaf cells shortened, shrunk and damaged under drought treatments (Fig. 8D). The total thickness of leaves, the thickness of the upper epidermis and lower epidermis of leaves decreased with drought stress. The tissue structure tense ratio of the leaf (CTR) showed an obvious increase with the aggravation of drought degrees. Under SD treatment, the CTR increased by 35.31 and 8 % compared to CK and MD treatment, respectively (Fig. 8E).

Discussion

In this study, 229 MYB TFs were identified in blueberry genome and on the free online platform of Majorbio Cloud Platform. Our previous study showed that blueberry genome was closely related to Vitis vinifera L. [40]. In order to reveal the divergence and differences of MYB gene family between different plant species, we compared the numbers of MYBs from Vitis vinifera L. (134 MYBs) [41], Pyrus bretschneideri Rehd. (129 MYBs) [28], Arabidopsis (197 MYBs) [15], Brassica rapa ssp. pekinensis (104 MYBs) [42], and O. sativa (155 MYBs) [15] with the number of VcMYBs in blueberry. It was found that the number of MYB genes in blueberry was obviously higher than that in other plant species, indicating that the large number of VcMYB genes in blueberry, which is tetraploid, was probably realted to the occurrence of multiple whole-genome duplication (WGD) events during the species evolution. It was also reported that the divergence of subfamily or gene number in same gene family of different species is related to their genome size or evolutionary differences [43]. Interestingly, we found that C15 and C23 subgroups belonged to two newly expanded subfamilies in the phylogenetic tree (Fig. 1), which may have important significance for blueberry.

Generally, there is a quite conserved DNA domain (MYB domain) at the N-terminus of MYB protein. The domain usually includes 1–4 incomplete repeats (R), and each R domain contains 3 α-helix structures [44, 45]. Six conserved tryptophan residues were evenly distributed in the R2 and R3 domains of MYB proteins [15, 46]. In this study, we found five conserved Trp residues of the VcR2R3-MYBs conservative domain in blueberry and the first Trp residue of R3 repeats was variable. Actually, a similar phenomenon was also reported in other species, the first Trp residue of R3 repeats was also variable in Alfalfa [11]. Substitution at the first Trp residue may be responsible for recognizing novel target genes or may lead to loss of DNA-binding activity against target genes [28]. Compared to the relatively conserved N terminal, the C-terminal domain of MYB protein has high variability [26, 47]. In the present study, by using the MEME program, we found that the conserved motifs located at the C-terminal of VcMYB proteins were quite different, such as VcMYB4, VcMYB113 and VcMYB226. Therefore, variable C-terminal domains and relatively conserved N-terminal domains formed the MYB gene family with functional diversity.

The expression correlation is the embodiment of the correlation coefficient between genes obtained by the spearman correlation algorithm, and is also a useful analysis method for predicting important putative genes [48]. Based on these methods, we identified 17 highly expressive genes in blueberry under drought stress. In particular, VcMYB88 and VcMYB228, which had the largest node, have more correlation with other MYB genes, probably playing an important role in blueberry response to drought stress. In ramie (Boehmeria nivea), Bn4CL3 was predicted to be involved in the synthesis of lignin through the expression correlation, and further research verified that Bn4CL3 negatively regulates lignin synthesis in the experiment of recombinant protein [49].

Although the genome sequencing of blueberry has been completed, the functional annotation of many genes was relatively imperfect [50, 51]. Therefore, it is a useful supplementary solution to perform systematic evolution analysis in order to understand and predict the potential functions of the key genes [52]. In the differential gene analysis, we identified 8 VcMYB genes with 20-fold differences, and these genes were mainly classified into C13 subfamily. Furthermore, in the expression correlation analysis, 17 VcMYB genes with higher correlation were obtained, which were mainly classified into the C11 and C18 subfamilies. Notably, the three subgroups C11, C13 and C18 contained many Arabidopsis MYB genes, such as AtMYB12, AtMYB41 and AtMYB108, and these genes have been reported to participate in plant responses to drought stress [20, 53, 54]. Therefore, we speculated that the VcMYB genes belonged to C18, C11and C13 subgroups, such asVcMYB2, VcMYB48 and VcMYB108, etc., probably play crucial roles in blueberry responding to drought stress. By using qRT-PCR analysis, we confirmed that these genes could be induced under drought stress.

A body of evidence demonstrated that MYB TFs participate in drought response by regulating the ROS pathway. For example, in Arabidopsis, overexpression of TaMYB33 gene improves tolerance to drought stress through the ROS pathway [55]; the GbMYB5 overexpression in tobacco increased the tolerance to drought stress by elevating the antioxidant enzyme activities and inducing the expression of stress-responsive genes (SOD, CAT and GST) [56]. According to this study, a total of 23 key VcMYB genes were identified (Figs. 2 and 6). Based on the orthologs analysis of Arabidopsis, the interaction partners of the key VcMYB proteins were predicted by STRING analysis. 7 interaction proteins of 10 key VcMYB proteins were identified (Fig. 7) and its orthologs in Arabidopsis such as AtEIN2, AtMPK3, AtKNAT7, AtSHN1 and AtSTZ, etc. [31, 33, 57] were reported to function through the ROS-mediated signaling and leaf regulation pathway under drought stress. Thus, we speculated that as the interaction partners of these proteins, VcMYB228, VcMYB100, VcMYB193, VcMYB2, VcMYB8, VcMYB14, VcMYB48, VcMYB102, VcMYB108 and VcMYB227, probably play crucial roles in ROS accumulation and elimination under drought stress (Figs. 7 and 9). Most abiotic stresses can cause oxidative stress and changes of ROS in plants [58, 59]. In the present study, we also observed that the reddish-brown patches of blueberry leaves by diaminobenzidine (DAB) staining increased under drought stress, indicating that drought stress increased the amount of oxygen released by the H2O2 in cells of blueberry leaves, which in turn further increase the degree of oxidation in cells (Fig. 7 D and E). With the aggravation of drought treatment, the activity of SOD in leaves and roots increased significantly (Fig. 7E). In fact, similar results have been reported in other species. For instance, in Manihotesculenta, the SOD activity of leaves increased significantly under drought stress [60]. In Brassica campestris, SOD activity and the content of H2O2 increased obviously when plants were subjected to drought stress [61]. Nevertheless, further study are needed for the 10 key VcMYB DEGs we screened to reveal the accurate signaling pathways of its involvement.

Fig. 9
figure 9

A prediction model of the mechanism of blueberry VcMYB genes in response to drought.

Full size image

The leaf is indispensable for regulating plant water transport and the changes of leaf morphology and structure is crucial for plants responding to drought stress [62]. CTR and leaf thickness are key indicators for the leaf morphology and structure [63, 64]. The leaf morphology and structure changes usually occurred due to the damaged biosynthesis of leaf cell walls, wax, and cuticle [39]. Here, we found that with the degree of drought increased, the leaf thickness of blueberries seriously decreased, while the CTR increased (Fig. 8D and E), which may be due to the disruption of xylem water supply and lead to the retardation of cell elongation and expansion [65]. In Quercus ilex, the leaves tended to get smaller under drought stress, while rehydration promotes larger leaves [66]. Previous studies showed that MYB TFs participates in drought response by leaf morphology and structure pathway. In Arabidopsis, the AtMYB96 gene overexpression improves tolerance to drought stress by changing leaf morphology [67]. In this study, based on the orthologs analysis in Arabidopsis, the interactions of key VcMYB proteins were predicted by STRING analysis. 14 interaction proteins of 7 key VcMYB proteins were identified and its orthologs in Arabidopsis such as AtKNAT7, AtSHN1 and AtNST1, etc. [68,69,70] were reported to function through the leaf regulation pathway under drought stress. Therefore, we speculated that as the important interaction partners of these proteins, VcMYB41, VcMYB88, VcMYB100, VcMYB152, VcMYB171, VcMYB193 and VcMYB202 are probably involved in the blueberry response to drought stress by regulating the leaf morphology and structure (Figs. 8 and 9).

In order to verify the accuracy of RNA-seq data, we randomly selected 16 VcMYB DEGs of blueberry under drought stress, and analyzed the RNA-seq and qRT-PCR data of these genes. The results showed that the expression trend of the selected gene were basically similar for RNA-seq and qRT-PCR data, suggesting the transcriptome sequencing was accurate and validity. Nevertheless, we also found that the folds of the expression changes were different. We speculated that two different sequencing methods probably lead to a slight difference. Zhu et al. [71] found that there were some differences between the two sequencing methods under different operating conditions. Interestingly, we found that some genes in blueberries have different expression trends in leaves and roots under drought treatment. For example, the expression levels of VcMYB20, VcMYB157 and VcMYB226 in leaves were significantly down-regulated under drought stress, while up-regulated in roots. Meanwhile, some tissue-specific genes, such as VcMYB48, were not expressed in leaves, but roots, indicating that VcMYB48 was a root-specific gene. Our study provide new promising candidates and insights for blueberry in response to drought stress.

Conclusions

229 non-redundant VcMYB sequences were identified in blueberry and can be divided into 23 subfamilies. 102 DEGs that significantly responded to drought treatment were identified, including 72 in leaves and 69 in roots, and 8 differential expression genes with a > 20-fold change in the level of expression. Among 102 DEGs, 17 DEGs had higher expression correlation with other MYB members. In addition, 15 key VcMYB genes including VcMYB8, VcMYB14, VcMYB88, VcMYB102, VcMYB108 and VcMYB228 etc.., which were probably involved in ROS and leaf morphology and structure pathways under drought treatment, were identified by STRING analysis combined with physiological and morphological observation. These findings provide new candidates and a new understanding of the regulatory mechanism of VcMYB genes family in blueberry response to drought stress.

Methods

Plant materials and experimental treatments

In this study, tissue-culture seedlings of blueberry variety ‘Bluecrop’ with a height of about 18 cm were used as the test materials. These blueberry seedlings were purchased from Dalian Senmao Modern Agriculture Co., Ltd., Dalian, China. They are very common fruit trees in China and are not endangered species. No specific permits were required for the sample collection.

In March 2018, the 3-month-old uniform tissue culture seedlings were transplanted into small plastic pots with the size of 15 × 15 cm at the top and 10 × 10 cm at the bottom. The nutrient soil included natural peat-soil and vermiculite and the volume ratio of peat-soil to vermiculite was 1:1. The available nutrients in the soil contained nitrogen with 140 mg·kg-1, phosphorus with 100 mg·kg-1, potassium with 180 mg·kg-1, soil organic matte with 45 % and pH value was 5.5. The seedlings were grown in the greenhouse with 16-h light/8-h dark photoperiod and watered with tap water modulated to pH 5.2–5.5 by dilute sulphuric acid. The greenhouse temperature was 25 ± 1℃ and the relative humidity was 55–65 %. All pots were irrigated with an equal amount of water every 3 days until the drought treatment was applied by withholding irrigation. The experiment was designed as a complete randomized block with three replicates per block, and nine pots were employed in each replicate. Three uniform plants were grown in each pot and three treatments were desighed as follows: CK, control group with soil water content 75–80 %; MD, moderate drought group with soil water content 55–60 %; SD, severe drought group with soil water content 30–35 %. Under the condition of greenhouse, three treatments of plants were applied by withholding irrigation 15 days after transplanting. The soil water content (SWC) reached the predetermined levels without water for 10, 20 and 40 days, respectively[7]. The gravimetric method was employed to control the SWC [72]. Every pot was weighed at 15:00 Beijing Standard Time (BST) each day. For RNA transcriptome sequencing, samples were taken when the extent of drought treatment reached the three predetermined levels, respectively. For physiological index detection, another 10-day treatment were conducted for each drought treatment after the three predetermined levels reached and then sampled. The pots were destroyed between 15:00 and 17:00 BST.

For sampling, the treated blueberry seedlings were cut off at soil level and leaves (L) and roots (R) were washed with distilled water and then wiped dry with paper towel. All samples were immediately frozen in liquid nitrogen and stored at -80℃ for further analysis. The experimental materials used for RNA-seq and qRT-PCR were single plant sample with three biological replicates. The materials used for physiological and biochemical analysis were mixed samples, and six plants were chosen as one biological replicate, and three biological replicates were included for each assay.

Identification of the MYB gene family in blueberry

Samples of the three treatments were sent to Shanghai Majorbio Bio-pharm Biotechnology Co., Ltd. (Shanghai, China) for library construction and sequenced by means of Illumina HiSeq xten system (www.majorbio.com) [7]. Total RNA was extracted from the sample using TRIzol® Reagent (Plant RNA Purification Reagent for plant tissue) according the manufacturer’s instructions (Invitrogen) and genomic DNA was removed using DNase I (TaKara). RNA-seq transcriptome library was prepared following the TruSeqTM RNA sample preparation Kit from Illumina (San Diego, CA) using 5 µg of total RNA. After quantified by TBS380, the paired-end RNA-seq sequencing library was sequenced with the Illumina HiSeq xten (PE150). The raw reads were trimmed and quality controlled by SeqPrep (https://github.com/jstjohn/SeqPrep ) and Sickle (https://github.com/najoshi/sickle ) with default parameters. Then clean reads were separately aligned to a reference genome using Hisat2 with default parameters (Http://tophat.cbcb.umd.edu/,version2.0.0) software. The original data were stored in the free online platform of Majorbio Cloud Platform (https://cloud.majorbio.com/) and the National Genomics Data Center, China National Center for Bioinformation / Beijing Institute of Genomics, Chinese Academy of Sciences (https://bigd.big.ac.cn/gsub/) (CRA003258).

The sequences of VcMYB protein sequences were obtained and analyzed from the free online platform of Majorbio Cloud Platform (https://cloud.majorbio.com/) and the blueberry genome (https://www.vaccinium.org/), and the sequences of AtMYB protein sequences in Arabidopsis were downloaded from the Arabidopsis genome TAIR database (https://www.arabidopsis.org/). Subsequently, all candidate VcMYBs were verified using SMART3 and Pfam Program (http://pfam.xfam.org/) to confirm that they contained the core domains [28]. Moreover, the remaining sequences were analyzed by ClustalX software, and all potentially redundant VcMYBs sequences were discarded. Non-redundant sequences were renamed as predicted VcMYB genes. The theoretical isoelectric point (PI), molecular protein weight (kDa) and an average of hydropathicity (GRAVY) index values of each VcMYB were obtained using the ProtParam tool (https://web.expasy.org/protparam/) [73].

Conserved motif and phylogenetic analysis of VcMYB genes

Multiple sequence alignments of the identified MYB amino acid sequences in Vaccinium corymbosum and Arabidopsis were performed using Clustal X (http://www.clustal.org/clustal2/) with default parameters [74]. The phylogenetic dendrogram was constructed by the neighbor-joining (NJ) method using MEGA6.0 software (https://www.megasoftware.net/) [75] with 1000 bootstrap value, and the phylogenetic tree was decorated by EvolView (https://evolgenius.info//evolview-v2/#login). In addition, the conserved motifs among VcMYB members were analyzed using MEME (v4.11.1) online program (http://meme-suite.org/index.html) [76], and the conserved motifs were confirmed according to the following parameters: the number of motifs was 10, and other parameters used the default settings. The graphics were decorated by TBtools software [77], and the phylogenetic analysis of VcMYB was the same as above. In order to analyze the sequences features of the R2 and R3 domains of VcR2R3-MYBs, the amino acid sequences of R2 and R3 repeats in Vaccinium corymbosum were extracted. Multiple sequence alignments of these identified R2R3-MYB proteins were aligned using Clustal X with default parameters. Multiple alignment files for R2 and R3 repeats were submitted to ggseqlogo [78] using the default settings to acquire sequence logos.

Differential expression analysis, functional enrichment and data analysis

The data was analyzed and manipulated on the free online platform of the Majorbio Cloud Platform. To identify DEGs between two different samples, the expression level of each gene was calculated according to the FPKM. RSEM (http://deweylab.biostat.wisc.edu/rsem/) [79] was used to quantify gene abundances. R statistical package software DESeq2 (http://bioconductor.org/ packages/stats/bioc/DESeq2/) was utilized for differential expression analysis [80]. The DEGs were defined with the p-adjust < 0.05 and fold change ≥ 2. Veen diagrams of DEGs were generated using the Majorbio Cloud Platform. Heatmaps were generated according to the log10 fold-change values at MD_R/SD_R/MD_L/SD_L compared with CK_R and CK_L[81]. GO functional-annotation analysis was carried out by BLAST2GO [82] and GO(Gene Ontology, http://www.geneontology.org) database. GO functional-enrichment analysis was performed to identify which DEGs were significantly enriched in GO terms at Bonferroni-corrected p-value ≤ 0.05 compared with the whole- genome background. GO functional enrichment analysis was carried out by Goatools (https://github. com/tanghaibao/Goatools) [83]. GO functional-annotation, functional-enrichment and string diagram were performed using the free online platform of the Majorbio Cloud Platform. Cluster analysis of the DEGs based on the log10 fold-change values, was generated using Expander 7 [84] software with the K-means algorithm [85] on the Majorbio Cloud Platform. Finally, expression correlation analysis of VcMYBs was analyzed based on the spearman correlation algorithm with the significant differences at q < 0.05 and the correlation coefficient for 0.5. The expression correlation network was constructed by the free online platform of the Majorbio Cloud Platform[7]. 15 VcMYB protein sequences were used as the targets and the STRING website (https://string-db.org/cgi/input. pl) was used to predict protein-protein interactions. Their orthologs in Arabidopsis thaliana were appointed as references [86]. All analysis settings were set at their default values unless previously mentioned.

Measurement of morphological and physiological indexes

The content of H2O2 and SOD activity were measured using commercial kits provided by Jiancheng Bioengineering Institute (Nanjing, China) according to the manufacturer’s instructions. According to Zhang et al, the histochemical staining of H2O2 was performed by DAB (3,3’- diaminobenzidine) method [87].

3 ~ 5 leaves below the top of the branch of seedlings were sampled and the tissue in the middle of the leaves with a width of 0.5 cm was fixed by FAA (ethanol: formaldehyde: glacial acetic acid = 90: 5:5). Paraffin section and safranine-fast green staining were performed according to the method by Miksche [88], and the slice thickness was 8 μm. The results were analyzed and photographed under a DM 3000 optical microscope [89], and leaf anatomic structures were measured and calculated using ImageJ and Nano Measurer 1.2 software. Total leaf thickness, upper epidermis thickness, palisade tissue thickness and lower epidermis thickness were observed. The CTR were calculated as follows: CTR = alisade tissue thickness/blade thickness ×100 %. Each treatment chose 3 leaves, and all experiments were carried out for three biological replicates.

Validation of RNA-seq sequencing data by RT-qPCR

Total RNA was extracted from leaves and roots of three group seedlings using TRIzol® Reagent (OmniPlant RNA Kit(DNase I))according to the manufacturer’s instructions (Invitrogen). The first-strand cDNA was synthesized with FastQuant cDNA First-Strand Synthesis Kit (Tiangen Biotechnology, Beijing). After adjusting the cDNA template concentration, qRT-PCR assays were detected by Fluorescence Quantification PCR Kit (Tiangen Biotechnology, Beijing) on StepOnePlusTM Real-Time PCR System (ABI, USA). The primer sequences of 20 DEGs were shown in Additional file 9 and VcUBC28 was used as the reference gene for analyzing expression data (Vashisth et al., 2011; Liang et al., 2019) [7, 90]. The RNA-seq data were displayed by log10 (FPKM + 1). All the qRT- PCR experiments were performed for three biological replicates.

Statistical analysis

The IBM-SPSS version 23 was used for statistical analysis. The statistical difference analysis was performed by Dunnett’s test with p < 0.05 designated as significant difference.

Availability of data and materials

All data and materials used in this study are publicly available. The datasets supporting the conclusions of this article are available in the NCBI repository (https://0-www-ncbi-nlm-nih-gov.brum.beds.ac.uk/bioproject/PRJNA737006/) with Bio project number PRJNA737006. The other datasets supporting the conclusions of this article are included within the article (and its additional files).

Abbreviations

aa:

Amino acid

MYB:

V-myb avian myeloblastosis viral oncogene homolog

TFs:

Transcription factors

DEGs:

Differentially expressed genes

GO:

Gene Ontology

ROS:

Reactive oxygen species

RNA-Seq:

RNA sequencing

SOD:

Superoxide dismutase

CTR:

Tissue structure tense ratio of the leaf

DAB:

3,3’-diaminobenzidine

qRT-PCR:

Real-time Quantitative polymerase chain reaction

References

  1. Plunkett BJ, Espley RV, Dare AP, Warren BA, Grierson ER, Cordiner S, et al. MYBA from blueberry (Vaccinium section Cyanococcus) is a subgroup 6 type R2R3MYB transcription factor that activates anthocyanin production. Frontiers in Plant Science. 2018;9:1300.

    Article  PubMed  PubMed Central  Google Scholar 

  2. Yang J, Li B, Shi W, Gong Z, Chen L, Hou Z. Transcriptional activation of anthocyanin biosynthesis in developing fruit of blueberries (Vaccinium corymbosum L.) by preharvest and postharvest UV irradiation. Journal of Agricultural Food Chemistry. 2018;66(42):10931–42.

    Article  CAS  PubMed  Google Scholar 

  3. FAO. The state of agricultural commodity markets 2018. Agricultural Trade. Climate Change: Rome; 2018.

    Google Scholar 

  4. Chen X, Qiu L, Guo H, Wang Y, Yuan H, Yan D, et al. Spermidine induces physiological and biochemical changes in southern highbush blueberry under drought stress. Brazilian Journal of Botany. 2017;40(4):841–51.

    Article  Google Scholar 

  5. Drobek M, Frc M, Cybulska J. Plant biostimulants: importance of the quality and yield of horticultural crops and the improvement of plant tolerance to abiotic stress—a review. Agronomy. 2019; 6(9).

  6. Golldack D, Li C, Mohan H, Probst N. Tolerance to drought and salt stress in plants: unraveling the signaling networks. Frontiers in Plant Science. 2014; 5.

  7. Liang K, Wang A, Sun Y, Yu M, Zhang L. Identification and expression of NAC transcription factors of Vaccinium corymbosum L. in response to drought stress. Forests. 2019;10(12):1088.

    Article  Google Scholar 

  8. Wang M, Huang Q, Lin P, Zeng Q, Li Y, Liu Q, et al. The overexpression of a transcription factor gene VbWRKY32 enhances the cold tolerance in Verbena bonariensis. Frontiers in Plant Science. 2019;10:1746.

    Article  PubMed  Google Scholar 

  9. Zhao Y, Cheng X, Liu X, Wu H, Bi H, Xu H. The wheat MYB transcription factor TaMYB is involved in drought stress responses in Arabidopsis. Frontiers in Plant Science. 2018;9:1426.

    Article  PubMed  PubMed Central  Google Scholar 

  10. Zhang K, Logacheva M, Meng Y, Hu J, Wan D, Li L, et al. Jasmonate-responsive MYB factors spatially repress rutin biosynthesis in Fagopyrum tataricum. Journal of Experimental Botany. 2018;69(8):1955–66.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  11. Zhou Q, Jia C, Ma W, Cui Y, Jin X, Luo D, et al. MYB transcription factors in alfalfa (Medicago sativa): genome-wide identification and expression analysis under abiotic stresses. PeerJ. 2019;7:e7714.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  12. Jiang CZ, Gu X, Peterson T. Identification of conserved gene structures and carboxy-terminal motifs in the MYB gene family of Arabidopsis and Oryza sativa L. ssp indica. Genome Biology. 2004; 57(7).

  13. Chen YH, Yang XY, He K, Liu MH, Li JG, Gao ZF, et al. The MYB transcription factor superfamily of Arabidopsis: expression analysis and phylogenetic comparison with the rice MYB family. Plant Molecular Biology. 2006;60(1):107–24.

    Article  CAS  Google Scholar 

  14. Zimmermann P, Hirschhoffmann M, Hennig L, Gruissem W. GENEVESTIGATOR. Arabidopsis microarray database and analysis toolbox. Plant Physiology. 2004;136(1):2621–32.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  15. Katiyar A, Smita S, Lenka SK, Rajwanshi R, Chinnusamy V, Bansal KC. Genome-wide classification and expression analysis of MYB transcription factor families in rice and Arabidopsis. BMC Genomics. 2012;13(1):544–4.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  16. Jung C, Seo JS, Han SW, Koo YJ, Kim CH, Song SI, et al. Overexpression of AtMYB44 enhances stomatal closure to confer abiotic stress tolerance in transgenic Arabidopsis. Plant Physiology. 2008;146(2):623–35.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  17. Seo JS, Sohn HB, Noh K, Jung C, An JH, Donovan CM, et al. Expression of the Arabidopsis AtMYB44 gene confers drought/salt-stress tolerance in transgenic soybean. Molecular Breeding. 2012;29(3):601–8.

    Article  CAS  Google Scholar 

  18. Persak H, Pitzschke A. Dominant repression by Arabidopsis transcription factor MYB44 causes oxidative damage and hypersensitivity to abiotic stress. International Journal of Molecular Sciences. 2014;15(2):2517–37.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  19. Park MY, Kang J, Kim SY. Overexpression of AtMYB52 confers ABA hypersensitivity and drought tolerance. Molecules Cells. 2011;31(5):447–54.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  20. Baldoni E, Genga A, Cominelli E. Plant MYB transcription factors: their role in drought response mechanisms. International Journal of Molecular Sciences. 2015;16(7):15811–51.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  21. Seo P, Lee S, Suh M, Park M, Go Y, Park C. The MYB96 transcription factor regulates cuticular wax biosynthesis under drought conditions in Arabidopsis. The Plant Cell. 2011;23(3):1138–52.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  22. Ambawat S, Sharma P, Yadav N, Yadav R. MYB transcription factor genes as regulators for plant responses: an overview. Physiology Molecular Biology of Plants. 2013;19(3):307–21.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  23. Wang RK, Cao ZH, Hao YJ. Overexpression of a R2R3 MYB gene MdSIMYB1 increases tolerance to multiple stresses in transgenic tobacco and apples. Physiologia Plantarum. 2014;150(1):76–87.

    Article  CAS  PubMed  Google Scholar 

  24. Du Y, Zhao M, Wang C, Gao Y, Wang Y, Liu Y, et al. Identification and characterization of GmMYB118 responses to drought and salt stress. BMC Plant Biology. 2018;18(1):1–18.

    Article  CAS  Google Scholar 

  25. Sun PP, Zhu XF, Huang XS, Liu JH. Overexpression of a stress-responsive MYB transcription factor of Poncirus trifoliata confers enhanced dehydration tolerance and increases polyamine biosynthesis. Plant Physiology Biochemistry. 2014;78:71–9.

    Article  CAS  PubMed  Google Scholar 

  26. Stracke R, Werber M, Weisshaar B. The R2R3-MYB gene family in Arabidopsis thaliana. Current Opinion in Plant Biology. 2001;4(5):447–56.

    Article  CAS  PubMed  Google Scholar 

  27. Dubos C, Stracke R, Grotewold E, Weisshaar B, Martin C, Lepiniec L. MYB transcription factors in Arabidopsis. Trends in Plant Science. 2010;15(10):573–81.

    Article  CAS  PubMed  Google Scholar 

  28. Cao Y, Han Y, Li D, Lin Y, Cai Y. MYB Transcription factors in Chinese pear (Pyrus bretschneideri Rehd.): genome-wide identification, classification, and expression profiling during fruit development. Frontiers in Plant Science. 2016;7:577.

    Article  PubMed  PubMed Central  Google Scholar 

  29. Fan F, Wang Q, Wen X, Ding G. Transcriptome-wide identification and expression profiling of Pinus massoniana MYB transcription factors responding to phosphorus deficiency. Journal of Forestry Research. 2019;31(3):909–19.

    Article  CAS  Google Scholar 

  30. Xingwang L, Bin L, Shudan X, Yanlinq C, Wenzhu Q, Chen J, et al. Cucumber (Cucumis sativus L.) nitric oxide synthase associated gene1 (CsNOA1) plays a role in chilling stress. Frontiers in Plant Science. 2016;7:1652.

    Google Scholar 

  31. Kim SH, Kim HS, Bahk S, An J, Yoo Y, Kim JY, et al. Phosphorylation of the transcriptional repressor MYB15 by mitogen-activated protein kinase 6 is required for freezing tolerance in Arabidopsis. Nucleic Acids Research. 2017;45(11):6613–27.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  32. Li KQ, Xing CH, Yao ZH, Huang XS. PbrMYB21, a novel MYB protein of Pyrus betulaefolia, functions in drought tolerance and modulates polyamine levels by regulating arginine decarboxylase gene. Plant Biotechnology Journal. 2017;15(9):1186–203.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  33. Nascimento FX, Glick JRM. BRJFiPe. Ethylene and 1-aminocyclopropane-1-carboxylate (ACC) in plant–bacterial interactions. Frontiers in Plant Science. 2018;9:114.

  34. Kumar A, Yogendra KN, Karre S, Kushalappa AC, Dion Y, Choo TM. WAX INDUCER1 (HvWIN1) transcription factor regulates free fatty acid biosynthetic genes to reinforce cuticle to resist Fusarium head blight in barley spikelets. Journal of Experimental Botany. 2016;67(14):4127–39.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  35. Zhang H, Ying YQ, Wang J, Zhao XH, Zeng W, Beahan C, et al. Transcriptome analysis provides insights into xylogenesis formation in moso bamboo (Phyllostachys edulis) shoot. Scientific Reports. 2018; 8.

  36. MacMillan CP, Birke H, Chuah A, Brill E, Tsuji Y, Ralph J, et al. Tissue and cell-specific transcriptomes in cotton reveal the subtleties of gene regulation underlying the diversity of plant secondary cell walls. Bmc Genomics. 2017; 18.

  37. Sajeevan RS, Nataraja KN, Shivashankara KS, Pallavi N, Gurumurthy DS, Shivanna M. Expression of Arabidopsis SHN1 in indian mulberry (Morus indica L.) increases leaf surface wax content and reduces post-harvest Water Loss. Frontiers in Plant Science. 2017;8:418.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  38. Oshima Y, Mitsuda N. The MIXTA-like transcription factor MYB16 is a major regulator of cuticle formation in vegetative organs. Plant Signaling Behavior. 2013;8(11):e26826.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  39. Tsukaya H, Beemster GT. Genetics, cell cycle and cell expansion in organogenesis in plants. Journal of Plant Research. 2006;119(1):1–4.

    Article  PubMed  Google Scholar 

  40. Li L, Zhang H, Liu Z, Cui X, Zhang T, Li Y, et al. Comparative transcriptome sequencing and de novo analysis of Vaccinium corymbosum during fruit and color development. BMC Plant Biology. 2016;16(1):223.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  41. Wong DCJ, Schlechter R, Vannozzi A, Holl J, Hmmam I, Bogs J, et al. A systems-oriented analysis of the grapevine R2R3-MYB transcription factor family uncovers new insights into the regulation of stilbene accumulation. DNA Research. 2016;23(5):451–66.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  42. Saha G, Park J, Ahmed NU, Kayum A, Kang K, Nou I. Characterization and expression profiling of MYB transcription factors against stresses and during male organ development in Chinese cabbage (Brassica rapa ssp. pekinensis). Plant Physiology Biochemistry. 2016;104:200–15.

    Article  CAS  PubMed  Google Scholar 

  43. Liu J, Chen N, Chen F, Cai B, Dal Santo S, Tornielli GB, et al. Genome-wide analysis and expression profile of the bZIP transcription factor gene family in grapevine (Vitis vinifera). BMC Genomics. 2014;15(1):281.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  44. Kranz H, Denekamp M, Greco R, Jin H, Leyva A, Meissner R, et al. Towards functional characterisation of the members of the R2R3-MYB gene family from Arabidopsis thaliana. The Plant Journal: for Cell Molecular Biology. 1998;16(2):263–76.

    Article  CAS  Google Scholar 

  45. Li Q, Zhang CJ, Li J, Wang L, Ren ZH. Genome-wide identification and characterization of R2R3MYB family in cucumis sativus. Plos One. 2012; 7(10).

  46. Wang F, Li X. Erratum to: Genome-wide characterization and expression analysis of MYB transcription factors in Lotus japonicas and Medicago truncatula. Genes Genomics. 2017;39(8):831–42.

    Article  CAS  Google Scholar 

  47. Du H, Yang S, Liang Z, Feng B, Liu L, Huang Y, et al. Genome-wide analysis of the MYB transcription factor superfamily in soybean. BMC Plant Biology. 2012;12:106.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  48. Bageritz J, Willnow P, Valentini E, Leible S, Boutros M, Teleman AA. Gene expression atlas of a developing tissue by single cell expression correlation analysis. Nature Methods. 2019;16(8):750–6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  49. Tang Y, Liu F, Xing H, Mao K, Chen JJG. Correlation analysis of lignin accumulation and expression of key genes involved in lignin biosynthesis of Ramie (Boehmeria nivea). 2019; 10(5):389.

  50. Vikas G, Estrada AD, Ivory B, Rob R, Ketan P, Meyer MD, et al. RNA-seq analysis and annotation of a draft blueberry genome assembly identifies candidate genes involved in fruit ripening, Biosynthesis of Bioactive Compounds, and Stage-Specific Alternative Splicing. 2015; 4(1):5.

  51. Colle M, Leisner CP, Wai CM, Ou S, Bird KA, Wang J, et al. Haplotype-phased genome and evolution of phytonutrient pathways of tetraploid blueberry. GigaScience. 2019;8(3):12.

    Article  CAS  Google Scholar 

  52. Irvine D, Tuerk C, Gold L. Systematic evolution of ligands by exponential enrichment with integrated optimization by non-linear analysis. Selexion. 1991;222(3):739–61.

    CAS  Google Scholar 

  53. Lippold F, Sanchez DH, Musialak M, Schlereth A, Scheible WR, Hincha DK, et al. AtMyb41 regulates transcriptional and metabolic responses to osmotic stress in Arabidopsis. Plant Physiology. 2009;149(4):1761–72.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  54. Schenke D, Bottcher C, Scheel D. Crosstalk between abiotic ultraviolet-B stress and biotic (flg22) stress signalling in Arabidopsis prevents flavonol accumulation in favor of pathogen defence compound production. Plant Cell Environment. 2011;34(11):1849–64.

    Article  CAS  Google Scholar 

  55. Qin YX, Wang MC, Tian YC, He WX, Han L, Xia GM. Over-expression of TaMYB33 encoding a novel wheat MYB transcription factor increases salt and drought tolerance in Arabidopsis. Molecular Biology Reports. 2012;39(6):7183–92.

    Article  CAS  PubMed  Google Scholar 

  56. Chen TZ, Li WJ, Hu XH, Guo JR, Liu AM, Zhang BL. A cotton MYB transcription factor, GbMYB5, is positively involved in plant adaptive response to drought stress. Plant Cell Physiology. 2015;56(5):917–29.

    Article  CAS  PubMed  Google Scholar 

  57. Park HC, Han HJ, Lee SM, Yun DJ, Chung WS. ASYMMETRIC LEAVES1 is phosphorylated by MPK3/6 in Arabidopsis thaliana. Journal of Plant Biology. 2013;56(4):208–15.

    Article  CAS  Google Scholar 

  58. Zhu J. Abiotic stress signaling and responses in plants. Cell. 2016;167(2):313–24.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  59. Choudhury FK, Rivero RM, Blumwald E, Mittler R. Reactive oxygen species, abiotic stress and stress combination. Plant Journal. 2017;90(5):856–67.

    Article  CAS  Google Scholar 

  60. Shan Z, Luo X, Wei M, Huang T, Khan A, Zhu Y. Physiological and proteomic analysis on long-term drought resistance of cassava (Manihot esculenta Crantz). Scientific Reports. 2018;8(1):17982.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  61. Rezayian M, Niknam V, Ebrahimzadeh H. Penconazole and calcium ameliorate drought stress in canola by upregulating the antioxidative enzymes. Functional Plant Biology: FPB. 2020;47(9):825–39.

    Article  CAS  PubMed  Google Scholar 

  62. Sack L, Cowan PD, Jaikumar N, Holbrook NM. The ‘hydrology’ of leaves: co-ordination of structure and function in temperate woody species. Plant Cell Environment. 2003;26(8):1343–56.

    Article  Google Scholar 

  63. Lin WU, Yan H, Xiao-Lan N, Li-Jun W. Studies on leaf tissue parameters and their relations to drought resistance of sea buckthorn. Journal of Jilin Agricultural University. 2003;25(4):390–3.

    Google Scholar 

  64. Ogburn RM, Edwards EJ. Quantifying succulence: a rapid, physiologically meaningful metric of plant water storage. Plant Cell Environment. 2012;35(9):1533–42.

    Article  Google Scholar 

  65. KM J, LS H, WS Y. Cork oak (Quercus suber L.) forest decline in Tunisia: a linkage between physiological adaptation and stress. Scientific Research Essays. 2011;6(6):1143–6.

    Google Scholar 

  66. PeA-Rojas K, Aranda X, Joffre R, Fleck I. Leaf morphology, photochemistry and water status changes in resprouting Quercus ilex during drought. Plant Function Evolutionary Biology. 2005;32(2):117–30.

    Article  Google Scholar 

  67. Seo PJ, Xiang FN, Qiao M, Park JY, Lee YN, Kim SG, et al. The MYB96 transcription factor mediates abscisic acid signaling during drought stress response in Arabidopsis. Plant Physiology. 2009;151(1):275–89.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  68. Zhang J, Xie M, Tuskan G, Muchero W, Chen J. Recent advances in the transcriptional regulation of secondary cell wall biosynthesis in the woody plants. Frontiers in Plant Science. 2018;9:1535.

    Article  PubMed  PubMed Central  Google Scholar 

  69. Laubscher M, Brown K, Tonfack L, Myburg A, Mizrachi E, Hussey S. Temporal analysis of Arabidopsis genes activated by Eucalyptus grandis NAC transcription factors associated with xylem fibre and vessel development. Sci Rep-Uk. 2018;8(1):10983.

    Article  CAS  Google Scholar 

  70. Oshima Y, Mitsuda N. The MIXTA-like transcription factor MYB16 is a major regulator of cuticle formation in vegetative organs. Plant Signaling Behavior. 2013;8(11):26826.

    Article  CAS  Google Scholar 

  71. Zhu C, Shi F, Chen Y, Wang M, Zhao Y, Geng G. Castanea mollissima transcriptome analysis of Chinese chestnut (Blume) in response to yasumatsu infestation. International Journal of Molecular Sciences. 2019; 20(4).

  72. Luo YZ, Liu H, Yan G, Li G, Turner NC. Roots of lucerne seedlings are more resilient to a water deficit than leaves or stems. Agronomy. 2019;9(3):123.

    Article  CAS  Google Scholar 

  73. Elisabeth G, Alexandre G, Christine H, Ivan I, Appel RD, Amos B. ExPASy: the proteomics server for in-depth protein knowledge and analysis. Nuclc Acids Research. 2003;31(13):3784–8.

    Article  Google Scholar 

  74. Larkin MA, Blackshields G, Brown NP, Chenna R, Mcgettigan PA, Mcwilliam H, et al. Clustal W and Clustal X version 2.0. Bioinformatics. 2007;23(21):2947–8.

    Article  CAS  PubMed  Google Scholar 

  75. Kumar S, Stecher G, Tamura K. MEGA7: molecular evolutionary genetics analysis version 7.0 for bigger datasets. Molecular Biology Evolution. 2016;33(7):1870–4.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  76. Bailey T, Boden M, Buske F, Frith M, Grant C, Clementi L, et al. MEME SUITE: tools for motif discovery and searching. Nucleic Acids Research. 2009;37:W202–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  77. Chen C, Chen H, Zhang Y, Thomas H, Frank M, He Y, et al. TBtools: an integrative toolkit developed for interactive analyses of big biological data. Molecular Plant. 2020;13(8):1194–202.

    Article  CAS  PubMed  Google Scholar 

  78. Wagih O. Ggseqlogo: a versatile R package for drawing sequence logos. Bioinformatics. 2017;33(22):3645–7.

    Article  CAS  PubMed  Google Scholar 

  79. Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-seq data with or without a reference genome. BMC Bioinformatics. 2011;12(1):323.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  80. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology. 2014;15(12):550.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  81. Li W, Li X, Chao J, Zhang Z, Wang W, Guo Y. NAC family transcription factors in tobacco and their potential role in regulating leaf senescence. Frontiers in Plant Science. 2018; 9:1900.

  82. Conesa A, Gotz S, Garcia-Gomez JM, Terol J, Talon M, Robles MJB. Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research. Bioinformatics. 2005;21(18):3674–6.

    Article  CAS  PubMed  Google Scholar 

  83. Klopfenstein D, Zhang L, Pedersen BS, Ramírez F, Vesztrocy AW, Naldi A, et al. GOATOOLS: a python library for gene ontology analyses. Scientific Reports. 2018;8(1):10872.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  84. Ulitsky I, Maron-Katz A, Shavit S, Sagir D, Linhart C, Elkon R, et al. Expander: from expression microarrays to networks and functions. Nature Protocols. 2010;5(2):303.

    Article  CAS  PubMed  Google Scholar 

  85. Shamir R, Maron-Katz A, Tanay A, Linhart C, Steinfeld I, Sharan R, et al. EXPANDER–an integrative program suite for microarray data analysis. BMC Bioinformatics. 2005;6(1):232.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  86. Meng X, Li Y, Yuan Y, Zhang Y, Liu M. The regulatory pathways of distinct flowering characteristics in Chinese jujube. Horticulture Research. 2020;7(1):123.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  87. Zhang H, Cui X, Guo Y, Luo C, Zhang L. Picea wilsonii transcription factor NAC2 enhanced plant tolerance to abiotic stress and participated in RFCP1-regulated flowering time. Plant Molecular Biology. 2018;98(6):471–93.

    Article  CAS  PubMed  Google Scholar 

  88. De Micco V, Aronne G. Combined histochemistry and autofluorescence for identifying lignin distribution in cell walls. Biotechnic Histochemistry. 2007;82(4):209–16.

    Article  PubMed  CAS  Google Scholar 

  89. Griffiths H, Weller G, Toy LF, Dennis RJ. You’re so vein: bundle sheath physiology, phylogeny and evolution in C3 and C4 plants. Plant Cell Environment. 2013;36(2):249–61.

    Article  CAS  Google Scholar 

  90. Vashisth T, Johnson LK, Malladi A. An efficient RNA isolation procedure and identification of reference genes for normalization of gene expression in blueberry. Plant Cell Reports. 2011;30(12):2167–76.

    Article  CAS  PubMed  Google Scholar 

Download references

Acknowledgements

We thank Majorbio Company (www.majorbio.com) for the RNA-seq and Dr. Xiehai Song in Beijing Forestry University for the graphical design of Fig. 1 and 4.

Funding

This study was funded by the grant from the National Natural Science Foundation of China ( grant number 31872942), the Natural Science Foundation of Beijing ( grant number 6172026) and the project from Agricultural Ministry of China (grant number 2016ZX08009-003). These funding bodies had no role in the study’s design, sample collection, analysis or interpretation of data, and writing the manuscript.

Author information

Authors and Affiliations

  1. Research & Development Center of Blueberry, Key Laboratory of Forest Silviculture and Conservation of the Ministry of Education, College of Forestry, Beijing Forestry University, 35 QinghuaEast Road, 100083, Beijing, China

    Aibin Wang, Kehao Liang, Shiwen Yang, Yibo Cao, Lei Wang, Ming Zhang, Jing Zhou & Lingyun Zhang

Authors
  1. Aibin Wang
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Kehao Liang
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Shiwen Yang
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Yibo Cao
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Lei Wang
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Ming Zhang
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. Jing Zhou
    View author publications

    You can also search for this author in PubMed Google Scholar

  8. Lingyun Zhang
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

LZ and AW designed the experiment. AW, LW, MZ, and SY carried out the experiment. AW, YC, LW, JZ, and KL performed the statistical analyses. AW and LZ wrote the manuscript. LZ and YC revised the manuscript. All authors commented on the manuscript. All authors contributed to the article and approved the submitted version.

Corresponding author

Correspondence to Lingyun Zhang.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare no conflict of interest.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1: Table S1.

The MYB transcription factors in blueberry and Arabidopsis used in this study.

Additional file 2: Table S2.

The detailed information of VcMYB proteins.

Additional file 3: Fig. S1.

Validation of RNA-seq sequencing data by RT-qPCR. Expression of 16 VcMYB DEGs in leaf (A) or root (B) under drought stress. The left y-axis (black bars) is the relative expression level of qRT-PCR. The right y-axis (gray bars) is log10 (FPKM + 1) of RNA-sEq. Each value is indicated as the mean value ± standard error of three independent determinations, and different letters represent significant differences at p < 0.05 by Duncan’s multiple range tests.

Additional file 4: Fig. S2.

Sequence logos of MYB transcription factor domains in blueberry.

Additional file 5: Table S3.

Seventy-one highly conserved VcR2R3-MYB genes from blueberry.

Additional file 6: Table S4.

Functional annotation (gene ontology) of VcMYB proteins.

Additional file 7: Fig. S3.

Functional annotation of DEGs on the basis of gene ontology.

Additional file 8: Table S5.

Summary of the homologous genes in blueberry and Arabidopsis used in this study.

Additional file 9: Table S6.

Primer sequences were used in the experiment.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Wang, A., Liang, K., Yang, S. et al. Genome-wide analysis of MYB transcription factors of Vaccinium corymbosum and their positive responses to drought stress. BMC Genomics 22, 565 (2021). https://0-doi-org.brum.beds.ac.uk/10.1186/s12864-021-07850-5

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://0-doi-org.brum.beds.ac.uk/10.1186/s12864-021-07850-5

Keywords

BMC Genomics

ISSN: 1471-2164

Contact us