Next Article in Journal
Unraveling the Role of AtSRT2 in Energy Metabolism, Stress Responses, and Gene Expression during Osmotic Stress in Arabidopsis thaliana
Previous Article in Journal
MADS-Box Family Genes in Lagerstroemia indica and Their Involvement in Flower Development
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Analyzing the Diversity of MYB Family Response Strategies to Drought Stress in Different Flax Varieties Based on Transcriptome Data

1
School of Life Science, Inner Mongolia University, Hohhot 010020, China
2
Inner Mongolia Academy of Agricultural & Animal Husbandry Sciences, Hohhot 010031, China
3
Key Laboratory of Black Soil Protection and Utilization, Ministry of Agriculture and Rural Areas, Inner Mongolia Key Laboratory of Degradation Farmland Ecological Remediation and Pollution Control, Inner Mongolia Conservation Tillage Engineering Technology Research Center, Hohhot 010031, China
4
Agricultural College, Inner Mongolia Agricultural University, Hohhot 010019, China
*
Authors to whom correspondence should be addressed.
Plants 2024, 13(5), 710; https://doi.org/10.3390/plants13050710
Submission received: 28 November 2023 / Revised: 20 February 2024 / Accepted: 28 February 2024 / Published: 2 March 2024
(This article belongs to the Section Plant Genetics, Genomics and Biotechnology)

Abstract

:
The MYB transcription factor family has numerous members, and is involved in biological activities, such as ABA signaling, which plays an important role in a plant’s resistance to abiotic stresses such as drought. However, the diversity of MYB members that respond to drought stress and their regulatory mechanisms in different flax varieties were unclear. In this study, we obtained 855.69 Gb of clean data from 120 flax root samples from 20 flax (Linum usitatissimum L.) varieties, assembled 92,861 transcripts, and identified 434 MYB family members in each variety. The expression profiles of the MYB transcription factor family from 20 flax varieties under drought stress were analyzed. The results indicated that there are four strategies by which the MYB family responds to drought stress in these 20 flax varieties, each of which has its own specific processes, such as development, reproduction, and localization processes. The four strategies also include common biological processes, such as stimulus responses, metabolic processes, and biological regulation. The WGCNA method was subsequently employed to identify key members of the MYB family involved in response strategies to drought stress. The results demonstrated that a 1R-MYB subfamily gene co-expression network is significantly related to the gibberellin response and cytokinin-activated signaling pathway processes in the ‘Strategy 4’ for MYB family response to drought, identifying core genes such as Lus.scaffold70.240. Our results showed a diversity of MYB family responses to drought stress within flax varieties, and these results contribute to deciphering the mechanisms of the MYB family regulation of drought resistance. This will promote the more accurate breeding development of flax to adapt to agricultural production under drought conditions.

1. Introduction

Drought is one of the most significant abiotic stresses that affect plant growth and development, causing a series of reactions, such as the wilting of plant leaves, premature flowering, growth stagnation, or even death [1,2]. The MYB family of transcription factors (V-myb avian myeloblastosis viral oncogene homolog) is one of the largest families of transcription factors in plants [3,4,5,6]. According to the location and number of repetitive sequences in the MYB domain, MYB transcription factors can be divided into four subfamilies as follows: 1R-MYB (MYB-Related), 2R-MYB (R2R3-MYB), 3R-MYB, and 4R-MYB [7,8,9,10]. When subjected to drought stress, MYB transcription factors can play a role in conducting ABA signals, regulating the transcription of functional genes and other functions, and enhancing the drought resistance of plants [11,12]. The results of multiple studies support MYB as a potential transcription factor for plant breeding and improvement [13,14].
MYB transcription factors are involved in plant drought resistance in a variety of ways. Firstly, MYBs mainly depend on ABA signals to participate in the regulation of stomatal open–close movement [15]. For example, MYB44 can prevent the inhibition of RCA 1/PYL9 on ABI1 phosphatase activity by binding to ABA receptors such as PYL8, PYL9, and so on, thus participating in ABA signaling [16]. AtMYB77 can interact with MYB44 to form a heterodimer to work together under drought stress, and the myb44–myb77 double mutants have shown strong drought resistance [17]. The MYB transcription factor is also involved in the drought-resistant process associated with leaf permeability, by regulating the synthesis of flavonoids and the stratum corneum. The R2R3-MYB subfamily members of Arabidopsis thaliana, MIXTA-like genes MYB16 and MYB106, can regulate the formation of the stratum corneum in trichomes [18]. The transcription factor MYB12 in Arabidopsis thaliana is a key regulator of flavonoid biosynthesis, and the overexpression of MYB12 and MYB75 in transgenic plants significantly increases flavonoid accumulation, thus improving drought tolerance [19]. In addition, the drought resistance of the MYB family has also been widely described in other plants. For example, the expression of the AtMYB68 gene in Brassica napus enhanced its drought resistance [20]. Under drought stress, the expression of the TdMYB4A063 gene in durum wheat in the root system was significantly up-regulated after 2 h [21].
However, even for homologous genes, members of the MYB family have different expression patterns in the life activities of different species, and often even in different varieties of the same species. For example, the expression levels of transcription factors such as MYB15 in response to cold stress were significantly lower in tea varieties with high cold tolerance than those of sensitive varieties [22]. This indicates that the differences in varieties also diversify the MYB gene expression profile.
Flax (Linum usitatissimum L.) can be categorized based on its development into fiber flax, oil flax, and oil–fiber flax. Oil flax is one of the most important oil crops in the world; it is rich in nutrients that are beneficial to human health, such as omega-3 fatty acids [23,24,25]. Flax has good drought tolerance, making it a relevant material for studying plant drought resistance [26]. In this study, by analyzing the expression profiles of MYB family members in 20 flax varieties under drought stress, we identified combination patterns between the biological processes that MYB families participate in under drought stress. And then we summarized the response strategies of MYB members of different flax varieties to drought stress, and dug out the preference of MYB family members in these strategies. These results have significant reference value for shedding light on the cooperative drought resistance mechanism among members of the MYB family, and improving the efficiency of molecular breeding for drought resistance in crops such as flax.

2. Results

2.1. Flax MYB Family Member Distribution

By screening the genes containing the MYB domain in the amino acid sequence of the flax genome, 434 MYB genes were obtained, among which 235 genes were members of the 1R-MYB subfamily, 194 genes were members of the 2R-MYB subfamily, 3 genes were members of the 3R-MYB subfamily, and 1 gene belonged to the 4R-MYB subfamily. In addition, the Lus.scaffold127.115 gene contained five MYB domains. The gene IDs and domains of the MYB family members are given in Supplementary Table S1.

2.2. Quality Control of the Transcriptome Sequencing of 20 Flax Varieties

This sequencing involved three biological replicates of 20 flax varieties in two treatment groups (drought group and control group), with a total of 120 samples, from which 855.69 Gb of clean data were obtained, and 92,861 transcripts were assembled. The average error rates of the sequenced bases corresponding to data from each sample ranged from 0.0251% to 0.0268%, that of Q20 ranged from 97.33% to 97.98%, that of Q30 ranged from 92.61% to 94.11%, and that of GC content ranged from 45.99% to 49.65% (Supplementary Table S2).

2.3. Expression Levels and Difference Analysis of the MYB Family

DEGs in the MYB family of 20 flax varieties were screened and counted (Table 1, Figure 1). The results showed that all varieties had members of the MYB family that were differentially expressed under drought stress. Variety B (Zhangya No.1) had the most DEGs in the MYB family (193). In addition, the number of DEGs in the MYB families of the M (Ningya No.19, 110), H (Dingxi No.17, 151), E (Yiya No.4, 102), and P (Baya No.15, 99) varieties was also high (≥99), suggesting that the members of the MYB family played a larger role in the process of drought resistance within these varieties. The number of DEGs in the MYB family of variety C (Ningya No.15) was the least (15). In addition, the numbers of DEGs in the MYB family of R (Longya No.10, 22), T (Neiya No.9, 24), I (BGOLDXREDWING44X3, 26), L (Ningya No.17, 27), and Q (Gaolanbai, 28) were also small (<30).
Among the 20 flax varieties selected in this study, only A (Longza No.1) and T (Neiya No.9) were water-sensitive varieties. A’s (Longza No.1) MYB family had 64 DEGs, which was considered a medium number among the 20 flax varieties. The MYB family of T (Neiya No.9) had only 24 DEGs, which is considered small; thus, it was speculated that the expression level of the MYB family might affect the drought resistance of T (Neiya No.9).

2.4. The Diversity of Expression Patterns

The difference in the expression patterns of MYB family members in different varieties may have reflected the diversity of their regulatory pathways. According to the expression difference (log2FoldChange) of multiple MYB family members in the drought group and the control group in each flax variety, the MYB family distance of different flax varieties was calculated, hierarchical clustering was performed on the flax varieties, and a heat map was drawn (Figure 2). According to the clustering results, the MYB family expression patterns of the 20 flax varieties were divided into eight groups, and all MYB genes in the flax genome were divided into four groups.
In the expression pattern clustering results, R (Longya No.10) and L (Ningya No.17) were clustered into Group 1. A (Longza No.1), G (Zhangya No.2), and F (Dingya No.15) were clustered into Group 2, and T (Neiya No.9) was not clustered with other varieties, thus becoming Group 3 alone. D (Longya No.8) and C (Ningya No.15) were clustered into Group 4, and N (Baxuan No.3), E (Yiya No.4), and P (Baya No.15) were clustered into Group 5. S (Yiya No.5), J(R43), B (Zhangya No.1), O (Sha Che Zao Shu Zhong Hong), and H (Dingxi No.17) were clustered into the largest group in this variety clustering, Group 6. A large number of flax varieties in this group were down-regulated in Group 4 of the MYB gene clustering. Varieties M (Ningya No.19) and Q (Gaolanbai) were in the seventh group, while I (BGOLDXREDWING44X3) and K (Lixian) were in the eighth group (Figure 2).

2.5. Analysis of Potential Biological Significance of Expression Patterns

According to the expression value clustering of MYB family members, the 20 flax varieties were divided into eight groups (Figure 2). Next, we selected one flax variety from each group as the representative, and conducted the functional annotation of its DEGs in the MYB family to explore the underlying biological significance of different expression patterns. In addition to the water-sensitive variety A (Longza No.1, 64 DEGs) in the second group and the water-sensitive variety T (Neiya No.9, 24 DEGs) in the third group, the varieties with the most DEGs were selected in other groups in this study, including L (Ningya No.17, 27 DEGs), D (Longya No.8, 43 DEGs), E (Yiya No.4, 102 DEGs), B (Zhangya No.1, 193 DEGs), M (Ningya No.19, 110 DEGs), and K (Lixian, 52 DEGs). The KEGG annotations (Table 2) and GO analysis (Table 3 and Supplementary Table S3) were performed on DEGs of the MYB family in these eight varieties, and the response strategies of the MYB family to drought stress in different flax varieties were summarized by analyzing the combination patterns of biological processes annotated by MYB family DEGs in different flax varieties.
The KEGG annotation of the MYB family DEGs was performed on the eight flax varieties, and the results showed that Lus.scaffold440.7, a member of the flax 2R-MYB subfamily, was annotated with “Pentose and glucuronide conversions” under the “Metabolism category” in five flax varieties (A, D, E, L, T). Lus.scaffold111.105, a member of the 1R-MYB subfamily, was annotated with the “Circadian rhythm-plant” pathway in four varieties (B, K, L, M). “Plant Hormone Signal Transfer” was a unique KEGG pathway of variety B (Zhangya No.1), with which three 1R-MYB subfamily members (Lus.scaffold70.63, Lus.scaffold45.334 and Lus.scaffold34.70) were annotated (Table 3, Supplementary Table S1).
The GO molecular functions annotated for the MYB families of the eight varieties included transcription regulator activity and binding. The greater the number of DEGs present in the varieties, the more obvious the gene amplification of the binding function. In addition, the catalytic activity in molecular functions was annotated by DEGs within the MYB family of seven flax varieties (A, T, L, D, E, B, and M), indicating the presence of MYB transcription factors with catalytic activity in multiple flax varieties (Table 2 and Supplementary Table S3).
The GO annotation of biological processes showed that the biological processes involved in all eight varieties included the response to stimulation, metabolic processes, and biological regulation. These biological processes were a common element of the strategies of the MYB family in responding to drought stress across all eight flax varieties. The MYB family response to drought stress in the eight groups was divided into four strategies, based on specific biological processes where the annotated DEGs were not present in all varieties. The groups represented by the D and L varieties were categorized as ‘Strategy 1’, and did not contain specific biological processes. The group represented by the T variety was ‘Strategy 2’, and the unique biological process was the developmental process. The group represented by A, E, and K was ‘Strategy 3’, and the specific processes were the developmental process, cellular component organization or biogenesis, the reproductive process, and growth. The group represented by varieties B and M was categorized as ‘Strategy 4’, and the unique biological processes were the developmental process, cellular component organization or biogenesis, the reproductive process, growth, and the localization process. In each strategy, common biological processes and unique biological processes were combined into a combination module able to respond to drought stress.

2.6. MYB Family Co-Expression Network and Core Drought Resistance Members

Among the four strategies of the MYB family of the 20 flax varieties in response to drought stress, ‘Strategy 4’ had the most biological processes annotated by DEGs, suggesting that this strategy is more likely to require coordination and interaction among MYB family members. Therefore, according to the expression of MYB family members and the results of the functional annotations, we selected the five following flax varieties belonging to Strategy 4 of the MYB family response to drought stress: S (Yiya No.5), J (R43), B (Zhangya No.1), O (Sha Che Zao Shu Zhong Hong), and H (Dingxi No.17). We took them as an example to mine the co-expression module of MYB family members in Strategy 4 in response to drought stress, and mined the preferred MYB gene of this strategy. There were three biological replicates for each variety in the drought and control groups, resulting in a total of 30 samples. The gene co-expression modules and preferred genes of the flax MYB family related to drought resistance were screened and identified based on the WGCNA (weighted gene co-expression network analysis) algorithm.
After the background correction and standardization of the expression data of the MYB family members, abnormal and small-mutation genes were filtered. Subsequently, 255 genes exhibiting a related action intensity, most in line with the scale-free distribution condition, were selected for analysis from the 434 MYB family members. The soft threshold of β = 8 was selected for module identification and gene co-expression network construction, because the scale-free network reached a value above 0.9 for the first time in these conditions, indicating good connectivity of the co-expression network (Figure 3).
The blockwiseModules function was used for module identification and module analysis (Figure 4). In this study, a total of two co-expression modules, ‘Turquoise’ and ‘Blue’, were identified in Strategy 4. There were 152 MYB genes in the ‘Turquoise’ module, and 46 MYB genes in the ‘Blue’ module. The results showed that most genes in the ‘Turquoise’ module had a positive correlation with the control samples (b, h, j, o, and s) and a negative correlation with the drought samples (B, H, J, O, and S), while most genes in the ‘Blue’ module had a positive correlation with the drought-treated samples (B, H, J, O, and S) and a negative correlation with the control samples (b, h, j, o, and s). Therefore, it was considered that most genes in the ‘Turquoise’ module were down-regulated in expression under drought stress, while most genes in the ‘Blue’ module were up-regulated in expression under drought stress.
The MM-GS scatter plots of the ‘Turquoise’ and ‘Blue’ modules were constructed (Figure 5). The results showed that, in the ‘Turquoise’ module, there was a significant correlation between GS and MM (cor = 0.469, p < 0.00894), while in the ‘Blue’ module, the correlation between GS and MM was low (cor = −0.25, p < 0.183). Therefore, we judged that there were more genes in the ‘Turquoise’ module that were highly correlated with the samples. And there were also more MYB family members within its modules that were highly correlated with drought resistance. Therefore, we further annotated and enriched the MYB family members in the ‘Turquoise’ module.
All MYB family members in the ‘Turquoise’ module were annotated and enriched for GO analysis, and a GO enrichment bubble map was constructed (Figure 6). The results showed that some genes were significantly enriched in the gibberellin response and the signaling pathway processes of cytokinin activation. Therefore, the ‘Turquoise’ module could conduct signal transduction by participating in the process of the response to plant hormones, thereby helping flax varieties resist drought stress.
The gene nodes in the top 30 for node connectivity (degree) in the ‘Turquoise’ module were screened to construct a gene co-expression network (Figure 7). Located in the inner circle were the five genes with the highest node connectivity in the module. The five MYB family members in the inner ring (Lus.scaffold4.116, Lus.scaffold381.22, Lus.scaffold8.186, Lus.scaffold144.21, Lus.scaffold70.240) belonged to the 1R-MYB subfamily, where the gene Lus.scaffold70.240 was a homologous gene of AtMYB88, a member of the Arabidopsis thaliana MYB family (Supplementary Table S5). Therefore, these five genes play an important role in the ‘Turquoise’ module, and belong to the core genes of the module. Based on this, these five genes are the preferred genes of the flax MYB family in response to drought stress ‘Strategy 4’.
In order to clarify the expression condition of the ‘Strategy 4’ preference gene, which serves as the core gene of the ‘Turquoise‘ module, and to verify the accuracy of the transcriptome data, we selected variety S (Yiya No.5) (Figure 2, Table 4) to represent the flax variety following ‘Strategy 4’. Subsequently, q-PCR experiments were conducted using the Lus.scaffold70.240 gene, which is a homologous gene of Arabidopsis thaliana AtMYB88, as an example (Figure 8, Supplementary Table S5). The results showed that the expression trend of Lus.scaffold70.240 gene was consistent in q-PCR and transcriptome sequencing (Figure 8). The expression of this gene was down-regulated under drought stress, which was also consistent with our conclusion that most genes in this module were down-regulated under drought stress by drawing the heat map (Figure 4) of the ‘Turquoise’ module.

3. Discussion

The inability of plants to move leads to the need to mobilize gene expression in order to resist the growth retardation and oxidative damage caused by a drought environment, in which the transcription-factor-mediated gene expression regulatory network plays an important role [27,28]. The MYB transcription factor family is one of the largest transcription factor families in plants, and plays an important role in regulating plant organ development, secondary metabolism, the hormone signaling network, the cell cycle, and other biological processes during drought stress [29,30]. In this study, we obtained the transcriptome data of the MYB family in 20 flax varieties through field water treatment experiments. By analyzing the expression profiles of the MYB family under drought stress, we summarized the four strategies of the MYB family in these 20 flax varieties in response to drought stress and the preferred MYB family members in Strategy 4. This study enabled us to understand the combination models of the flax MYB transcription factor family that regulate biological processes under drought stress, and will help us to decode the synergistic drought resistance mechanism of MYB family members in the future. This will allow us to improve the efficiency of molecular breeding for the drought resistance of flax and other crops.

3.1. The Distribution of DEGs in the MYB Family of 20 Flax Varieties Was Significantly Different

Because different varieties have different paternal and maternal parents, and have been subjected to selection pressures from different places of origin, there is a high degree of genetic diversity between the genomes of different varieties of the same species [31]. This genomic diversity indicates that the expression patterns of different varieties of the same species will also be diverse under different growth stages or abiotic stress conditions [32]. A large number of transcriptome analysis results support this claim [33,34,35,36]. For example, the expression profiles of two potato varieties with contrasting drought tolerance showed significant differences during both the early and late stages of drought stress, as well as after rehydration [34]. Among various sorghum varieties, drought stress induced distinct changes in the expression patterns of many homologous genes related to cutin monomer and cutin wax biosynthesis [35].
In our previous studies, we confirmed that there is a wide genetic diversity in the genomes of different varieties of flax, which is influenced by many factors, such as parental and geographical origin [37]. In this study, we confirmed that the distribution of DEGs in the MYB family of different flax varieties is also significantly different under drought stress (Table 1, Figure 1). For example, there was a large gap between the 20 varieties in the number of DEGs (15~193). Five varieties (Zhangya No.1, Ningya No.19, Dingxi No.17, Yiya No.4, and Baya No.15) had more DEGs in the MYB family (≥99), while another five varieties (Ningya No.15, Longya No.10, Neiya No.9, Ningya No.17, Gaolanbai, and BGOLDXREDWING44X3) had fewer DEGs in the MYB family (<30) (Table 1, Figure 1).

3.2. The Response Strategies of the MYB Family to Drought Stress Were Diverse among Different Flax Varieties

The first step in exploring MYB’s response strategies to drought stress is to clarify the expression of MYB genes under drought stress. Similar expression patterns in the MYB family will correspond to similar strategies in response to drought stress [38]. In order to further clarify the similarity of the expression patterns of the MYB family in 20 flax varieties, we clustered the 20 flax varieties based on changes in the expression levels (log2FoldChange) of MYB family members under drought stress (Figure 2). The results showed that the 20 flax varieties could be divided into eight groups, and the expression patterns of the MYB family members in each group were similar (Figure 2).
Since the expression patterns of the flax MYB family in each group were similar, we screened the expression patterns of the flax MYB family in each group under drought stress, functionally annotating the DEGs. The KEGG annotation results showed that Lus.scaffold440.7, a member of the 2R-MYB subfamily of flax, was annotated with pentose and glucuronic acid conversion pathways under the metabolism category in five flax varieties (A, D, E, L, and T) (Table 2). The impact of drought stress on the expression of genes related to sugar metabolism has also been reflected in a large number of previous studies. For example, following drought stress, the genes related to starch synthesis in tea trees were down-regulated, while those related to starch decomposition were up-regulated, with 74 DEG annotations in the pentose and glucuronic acid conversion pathways [39]. The MYB28, MYB29, and MYB76 transcription factors have been identified as important regulatory factors of Aliphatic Glucosinolate (AG) biosynthesis [40,41,42]. The transcription factors MYB34, MYB51, and MYB122 are important regulators of Indolic Glucosinolate (IG) biosynthesis [43,44].
According to the biological processes in the GO functional annotation of the DEGs in the MYB family, eight varieties could be classified into four strategies in response to drought (Table 3). The common biological processes in these four strategies included the response to stimulation, metabolic processes, and biological regulation. These biological processes have been proven by a large number of studies to be the basic regulatory processes initiated by transcription factors under drought stress [45,46,47]. The unique biological process of ‘Strategy 2’ was the developmental process. The functional prediction results of the MYB family of Cajanus cajan also showed that CcR2R3-MYBs was involved in the developmental process and responded to abiotic stress, and the overexpression of CcMYB107 significantly improved the drought resistance of plants [48]. The specific processes of ‘Strategy 3’ were the reproductive process and growth. Nine R2R3-MYB members that were constitutively expressed in male reproductive tissues were also identified in Chinese cabbage, based on organ-specific expression data [49]. The biological processes specific to ‘Strategy 4’ were localization and others. TvCyP1 cyclophilin has been proven to regulate the nuclear translocation of Myb1 and Myb3. It can directly bind to Myb3, and it can mediate a series of the conformational switches of Myb3 from the membrane chamber to the nucleus [50]. The pathways by which the MYB family regulates drought resistance have been extensively studied, and our work identified the pattern of the combinations of these biological processes involved in the regulation of the MYB family in the 20 flax varieties [44,51]. In the future, we will explore which model of the MYB family is the most efficient in regulating drought resistance, which will have important reference value for the breeding work of flax and other crops.
In this study, Lus.scaffold440.7 (2R-MYB subfamily), a homolog of the AtMYB93 gene, was simultaneously annotated in two GO terms (metabolic process and catalytic activity function) in five flax varieties (A, D, E, L, and T). Therefore, it was speculated that this gene product had catalytic activity, and could play a catalytic role in the metabolic process of multiple flax varieties (Supplementary Table S1). Several studies have shown that AtMYB93 is a negative regulator of lateral root development in Arabidopsis thaliana; it is strongly expressed in early lateral root endodermal cells, and it is induced by auxin in the meristem of the primary root [52,53]. In the future, we plan to conduct more detailed bioinformatics and molecular biological analyses of this gene in order to determine how effective it is in improving the drought resistance of different flax varieties.

3.3. WGCNA Mines the Preferred MYB Family Members of Strategy 4 Such as Lus.scaffold70.240

In this study, we found a 1R-MYB subfamily gene co-expression module (‘Turquoise’ module) that was highly related to drought resistance in Strategy 4, and most of the genes in the module were down-regulated under drought stress (Figure 4). The GO annotation and enrichment results indicated that some genes in this module were significantly enriched in the gibberellin response process and the cytokinin activation signaling pathway (Figure 6). Research has shown that the NtMYB330 transcription factor in tobacco can regulate the expression of genes related to abscisic acid (ABA) and gibberellic acid (GA) signaling, affecting both ABA-GA crosstalk and seed germination [54]. Studies have shown that the transcription of cytokinin target genes is regulated by type-b response regulatory proteins (RRBs), a transcription factor of the MYB family [55,56]. In addition, a transgenic experiment on the VyMYB24 gene of a highly drought-resistant grape variety showed that the gene may inhibit plant development by regulating gibberellin (GA) metabolism [57].
WGCNA screened the five following genes, which had the highest connectivity in this module: Lus.scaffold4.116, Lus.scaffold381.22, Lus.scaffold8.186, Lus.scaffold144.21, and Lus.scaffold70.240 (Figure 7). In the qRT-PCR experiment, the expression of the Lus.scaffold70.240 gene was down-regulated under drought stress, which was consistent with the trend of transcriptome sequencing (Figure 8). The Lus.scaffold70.240 gene is a homologue of the AtMYB88 gene of Arabidopsis thaliana, and the AtMYB88 gene can regulate the late division of stomatal cell lineage [58]. qRT-PCR showed that PsFLP, the homologous gene of AtMYB88 in Pisum sativum, was involved in the response of Pisum sativum to drought stress [59]. Lus.scaffold4.116 and Lus.scaffold144.21 are the homologous genes of the Arabidopsis thaliana MYB family TBP-like type genes—AT1G07540 and AT3G12560 (Supplementary Table S5)—which play an important role in regulating the structure and function of eukaryotic telomeres. Severe drought stress can cause DNA damage in plants [60]. RT-PCR experiments showed that some genes encoding TBP (Telomeric DNA-binding proteins), such as AtTRP2, were down-regulated when the Arabidopsis thaliana DNA was damaged [61]. This result is consistent with the conclusion of this study that most genes in the ‘Turquoise’ module are down-regulated under drought stress, which can confirm the conclusion of this study (Figure 4). Due to the highest connectivity, we speculated that these five genes are the preferred genes in the MYB family of flax for coping with drought stress under ‘Strategy 4’. These five genes will be the focus of our future research. We will use a large number of molecular biology experiments, such as q-PCR, sub-cellular localization, immunoprecipitation, gene knockout, etc., to explore the specific biological processes in which these genes participate in regulation under drought stress, as well as their impact on crop drought resistance. This will have important reference value for the molecular breeding of crops such as flax.

4. Materials and Methods

4.1. Identification of MYB Family Members

HMMER-3.0 software [62] was used to identify the MYB domain (PF00249) in the flax amino acid sequence. The criteria for determining the characteristics of the conservative domain were as follows: (1) The coefficient of similarity with the confirmed conservative sequence of the MYB domain, or conditional E-value value, was less than 1 × 10−5; and (2) the conservative sequence comprised at least one conservative tryptophan (W) [63]. ID simplification was performed on the obtained genes, and the genes obtained after the deletion of repeat sequences according to the GFF files were counted as members of the MYB family. According to the number of conservative domains of MYB, it was divided into four subfamilies as follows: 1R-MYB (MYB-Related), 2R-MYB (R2R3-MYB), 3R-MYB, and 4R-MYB [63,64]. The flax amino acid sequence, GFF, and other data were obtained from the Figshare database (https://figshare.com/articles/dataset/Annotation_files_for_Longya-10_genome/13614311, accessed from 15 October 2022).

4.2. Plant Materials

Based on previous research, such as the identification and cluster analysis of drought resistance in early field experiments, combined with the recommendations of drought-resistant varieties from the national oil crop industry technical system, 18 drought-resistant flax varieties and two water-sensitive varieties as control varieties were screened in this study to conduct the drought stress test (Table 4, Supplementary Table S4) [37].

4.3. MYB Gene Expression Induction Method

In this experiment, a split-plot design was adopted, and the main treatments were the two following kinds of water treatments: drought treatment and control treatment. The secondary treatment included 20 flax varieties, of which 18 were drought resistant and two were water sensitive. Each variety of each treatment had three repeats.
The 20 flax varieties grew under natural conditions in the field until the water-sensitive varieties grew to 20 cm, at which time, a dry shed was set up for the water treatment, and all the flax varieties had entered the full bloom stage. The water treatment was set as two gradients, which were divided into (1) control treatment (the soil water content was maintained at 70~75% of the maximum field water holding capacity) and (2) drought treatment (the soil water content was maintained at 25~30% of the maximum water holding capacity).
When the water treatment reached a significant level, where 80% of the water-sensitive varieties wilted or even lodged under drought treatment, and about 50% of the varieties in the drought-resistant group exhibited leaf wilting (at this time, the water treatment had lasted for 30 days, Figure 9), the roots of the flax were sampled, and three replicates were randomly selected from each group, totaling 120 samples (20 varieties × 2 treatments × 3 replicates). The samples were immediately frozen in a −80 °C refrigerator for cryopreservation until transcriptome sequencing at Majorbio company (Beijing, China). At the same time, we also took three samples of variety S (Yiya No.5) in the control group and treatment group (1 varieties × 2 treatments × 3 replicates) for the qRT-PCR experiment to verify the transcriptome data.

4.4. Determination and Analysis of the MYB Gene Expression Level

4.4.1. RNA Extraction

The total RNA in the samples was obtained via the Invitrogen Method, and genomic DNA was removed using DNase I (TaKara, Kyoto, Japan) [65]. The concentration and purity of the extracted RNA were detected using NanoDrop 2000 (NanoDrop Technologies, Wilmington, NC, USA). RNA integrity was detected using agarose gel electrophoresis, and the RIN value was determined by Agilent 2100 (Agilent Technologies, Santa Clara, CA, USA) to ensure the transcriptome sequencing was performed using qualified samples (OD 260/280 = 1.8–2.2, OD260/230 ≥ 2.0, RIN ≥ 6.5, 28S: 18S ≥ 1.0, >1 μg).

4.4.2. cDNA Library Construction and Computer Sequencing

The RNA library was constructed using the TruseQTM RNA Sample Preparation Kit (IILUMINA, San Diego, CA, USA). MRNAs were first isolated and enriched from the total RNA by A-T base pairing with Oligo(dT)-bearing magnetic beads with the 3′ terminal ployA of mRNA. The mRNA was randomly fragmented by adding fragmentation buffer, and small fragments of about 300 bp were isolated via magnetic bead screening. Next, using a Superscript Double-Stranded cDNA Synthesis Kit (Invitrogen, Carlsbad, CA, USA), six-base random hexamers (IILUMINA, San Diego, CA, USA) were added, and mRNA was used as the template to invert and synthesize one-strand cDNA, followed by two-strand synthesis in order to form a stable double-strand structure. The double-stranded cDNA had a sticky end. End Repair Mix was added to make it a flat end, followed by an “A” base at the 3′ end for a connection to the Y-junction. After the enrichment of the cDNA by PCR, DNA Clean Beads screened the 200–300 bp band. After quantification via TBS380 (PicoGreen) (Molecular Probes, Eugene, OR, USA), the library was high-throughput-sequenced using the Illumina NovaSeq 6000 sequencing platform with a sequencing read length of PE150.

4.4.3. Transcriptome Data Quality Control

Fastp-0.19.5 software [66] was used to remove the adapter sequence from the reads, and remove the reads without the inserted fragments due to linker self-ligation, etc. We pruned low-quality sequences with a 5′ end mass value of less than 20 and a 3′ end mass value of less than 3, and removed reads containing ambiguous bases. Sequences with a length of less than 30 bp after the removal of the adapter and mass pruning were discarded.

4.4.4. Analysis of MYB Family Gene Expression

TPM transformation was performed on the Read Counts of each sample gene to obtain standardized gene expression levels. Subsequently, the expression levels of MYB family genes in all varieties and treatments were extracted to form an MYB family genes expression matrix, and expression difference analysis was performed using the DESeq2- 1.24.0 package in R language to calculate the differential expression multiple log2FoldChange of each variety in the drought group and the control group. The DEGs (differentially expressed genes) of the MYB family in each variety under drought stress were identified using p-value < 0.05, |log2FoldChange| > 1.

4.4.5. Expression Pattern Cluster Analysis

The expression pattern clustering analysis was conducted based on the log2FoldChange matrix of the MYB genes of all 20 varieties using RSEM-1.3.3 software [67]. The clustering method of the genes and varieties was hierarchical clustering, and the clustering algorithm was Spearman. The expression patterns of the 20 flax varieties were grouped according to the variety clustering results. The MYB family DEGs of one variety for each group were selected for functional annotation, including GO annotation and KEGG annotation, to identify the functional diversity of the MYB family DEGs in different varieties.

4.4.6. Analysis of the Gene Co-Expression Network

A WGCNA (weighted gene co-expression network analysis) algorithm was used to construct a gene co-expression network, and the code was written with reference to the WGCNA-1.63 package in R language for analysis [68]. The specific analysis process was as follows:
  • Data preprocessing. Background correction and standardization were performed on the MYB gene expression data, and the genes with abnormal and small variations were filtered. Then, sample clustering and scale-free curve analysis were performed to determine the soft threshold β of the recognition module.
  • Module identification. Module identification was performed using the blockwiseModules function, with the parameters set to minModuleSize = 30, minKMEtoStay = 0.3, and mergeCutHeight = 0.25.
  • Module analysis. A target module that was highly correlated with the drought resistance of flax was found, and the genes in this module were analyzed for functional enrichment by using the Goatools-0.6.5 software. The analysis method was Fisher’s exact test, p-adjust < 0.05.
  • Mining core genes of the target module. The degree of the connectivity of each node in the target module was calculated using Cytoscape-3.9.0 software, and the top 30 nodes were screened to construct a gene co-expression network [69].

4.4.7. qRT-PCR Experiment for Validation

Primer-5 software was used to design the primer sequence, and the output production was 150–300 bp. The Actin gene of flax was used as the internal reference, and the primers and target gene sequences were shown in Supplementary Table S6. The reaction system was SYBR Premix Ex Taq (Tli RNaseH Plus) (2×) 10 μL, PCR Forward Primer 0.5 μL, PCR Reverse Primer 0.5 μL, ROX Reference Dye or Dye (50×) 0.4 μL, cDNA solution 1 μL, ddH2O 7.6 μL. The PCR procedure was 95 °C 30 s; 95 °C 5 s; 60 °C 31 s for a total of 40 cycles. Three biological replicates and three technical replicates were set up, and the relative gene expression was calculated using 2−∆∆CT.

5. Conclusions

This study investigated the expression profile diversity of the MYB transcription factor family in 20 flax varieties subjected to drought stress. The results indicated that there are four strategies by which the MYB family responds to drought stress in 20 flax varieties, each of which has its own specific processes, such as development, reproduction, and localization processes. The four strategies also involve common biological processes, such as stimulus responses, metabolic processes, and biological regulation. A 1R-MYB subfamily gene co-expression network that was significantly related to the gibberellin response and the cytokinin-activated signaling pathway processes was identified in the ‘Strategy 4’ for MYB response to drought, and core genes such as Lus.scaffold70.240 were obtained by WGCNA. In conclusion, the response strategies of the MYB families of different flax varieties to drought stress are diverse.

Supplementary Materials

The following supporting information can be downloaded at https://www.mdpi.com/article/10.3390/plants13050710/s1, Supplementary Table S1: Gene ID and domain of flax MYB family members; Supplementary Table S2: Statistical table of sequencing quality of flax varieties; Supplementary Table S3 Annotated MYB gene statistics for each GO term; Supplementary Table S4 Basic information of flax varieties; Supplementary Table S5: Comparison results of ‘Strategy 4’ preference genes and Arabidopsis thaliana MYB family homologous genes using the Blastp algorithm; Supplementary Table S6: Primer sequences used in the q-PCR experiment.

Author Contributions

Conceptualization, Z.L.; data curation, X.Z and L.Y.; formal analysis, F.Z., J.M. and Y.L.; project administration, S.S. and Y.C.; software, F.Z. and S.B.; visualization, F.Z. and L.C.; writing—original draft, F.Z.; writing—review and editing, F.Z., Y.L., J.M. and X.Z. All authors have read and agreed to the published version of the manuscript.

Funding

This study was funded by the National Natural Science Foundation of China (32260457, U23A20195), Scientific and Technological Projects of Grassland Talents in Inner Mongolia Autonomous Region (Lu., Zhao.), the Leading Talent Project of “Science and Technology Leading Talent Team Project of Inner Mongolia Autonomous Region” (2022LJRC0010).

Data Availability Statement

The transcriptome (RNA-seq) data used in this study were deposited in the Sequence Read Achieve (SRA) of the NCBI database (BioProject ID: PRJNA1046211).

Conflicts of Interest

The authors declare no conflicts of interest.

References

  1. Sayed, M.A.; Schumann, H.; Pillen, K.; Naz, A.A.; Léon, J. AB-QTL analysis reveals new alleles associated to proline accumulation and leaf wilting under drought stress conditions in barley (Hordeum vulgare L.). BMC Genet. 2012, 13, 61. [Google Scholar] [CrossRef] [PubMed]
  2. Singhal, P.; Jan, A.T.; Azam, M.; Haq, Q.M.R. Plant abiotic stress: A prospective strategy of exploiting promoters as alternative to overcome the escalating burden. Front. Life Sci. 2015, 9, 52–63. [Google Scholar] [CrossRef]
  3. Wang, X.; Niu, Y.; Zheng, Y. Multiple Functions of MYB Transcription Factors in Abiotic Stress Responses. Int. J. Mol. Sci. 2021, 22, 6125. [Google Scholar] [CrossRef] [PubMed]
  4. Riechmann, J.; Heard, J.; Martin, G.; Reuber, L.; Jiang, C.; Keddie, J.; Adam, L.; Pineda, O.; Ratcliffe, O.; Samaha, R.; et al. Arabidopsis transcription factors: Genome-wide comparative analysis among eukaryotes. Science 2000, 290, 2105–2110. [Google Scholar] [CrossRef] [PubMed]
  5. Wang, D.; Lipsick, J. Mutational analysis of the transcriptional activation domains of v-Myb. Oncogene 2002, 21, 1611–1615. [Google Scholar] [CrossRef]
  6. Persson, M.; Andrén, Y.; Mark, J.; Horlings, H.; Persson, F.; Stenman, G. Recurrent fusion of MYB and NFIB transcription factor genes in carcinomas of the breast and head and neck. Proc. Natl. Acad. Sci. USA 2009, 106, 18740–18744. [Google Scholar] [CrossRef]
  7. Nakagoshi, H.; Nagase, T.; Ueno, Y.; Ishii, S. Transcriptional trans-repression by the c-myb proto-oncogene product. Nucleic Acids Res. 1989, 17, 7315–7324. [Google Scholar] [CrossRef]
  8. Bussolari, R.; Candini, O.; Colomer, D.; Corradini, F.; Guerzoni, C.; Mariani, S.A.; Cattelani, S.; Silvestri, C.; Pecorari, L.; Iacobucci, I.; et al. Coding sequence and intron–exon junctions of the c-myb gene are intact in the chronic phase and blast crisis stages of chronic myeloid leukemia patients. Leuk. Res. 2007, 31, 163–167. [Google Scholar] [CrossRef]
  9. Dubos, C.; Stracke, R.; Grotewold, E.; Weisshaar, B.; Martin, C.; Lepiniec, L. MYB transcription factors in Arabidopsis. Trends Plant Sci. 2010, 15, 573–581. [Google Scholar] [CrossRef]
  10. Antje, F.; Katja, M.; Braun, E.L.; Erich, G. Evolutionary and comparative analysis of MYB and bHLH plant transcription factors. Plant J. 2011, 66, 94–116. [Google Scholar]
  11. Agarwal, M.; Hao, Y.; Kapoor, A.; Dong, C.H.; Fujii, H.; Zheng, X.; Zhu, J.K. A R2R3 type MYB transcription factor is involved in the cold regulation of CBF genes and in acquired freezing tolerance. J. Biol. Chem. 2006, 281, 37636–37645. [Google Scholar] [CrossRef]
  12. Nakabayashi, R.; Yonekura-Sakakibara, K.; Urano, K.; Suzuki, M.; Yamada, Y.; Nishizawa, T.; Matsuda, F.; Kojima, M.; Sakakibara, H.; Shinozaki, K.; et al. Enhancement of oxidative and drought tolerance in Arabidopsis by overaccumulation of antioxidant flavonoids. Plant J. Cell Mol. Biol. 2014, 77, 367–379. [Google Scholar] [CrossRef] [PubMed]
  13. Rodrigues, J.A.; Espley, R.V.; Allan, A.C. Genomic analysis uncovers functional variation in the C-terminus of anthocyanin-activating MYB transcription factors. Hortic. Res. 2021, 8, 77. [Google Scholar] [CrossRef]
  14. Liu, X.; Zhang, S.; Sun, M.; Guo, Y.; Zhao, S.; Zhou, X.; Bai, X.; Dai, K.; Li, H.; Yuan, X.; et al. SiMYBS3, Encoding a Setaria italica Heterosis-Related MYB Transcription Factor, Confers Drought Tolerance in Arabidopsis. Int. J. Mol. Sci. 2023, 24, 5418. [Google Scholar] [CrossRef]
  15. Xie, Z.; Li, D.; Wang, L.; Sack, F.D.; Grotewold, E. Role of the stomatal development regulators FLP/MYB88 in abiotic stress responses. Plant J. 2010, 64, 731–739. [Google Scholar] [CrossRef] [PubMed]
  16. Li, D.; Li, Y.; Zhang, L.; Wang, X.; Zhao, Z.; Tao, Z.; Wang, J.; Wang, J.; Lin, M.; Li, X.; et al. Arabidopsis ABA receptor RCAR1/PYL9 interacts with an R2R3-type MYB transcription factor, AtMYB44. Int. J. Mol. Sci. 2014, 15, 8473–8490. [Google Scholar] [CrossRef] [PubMed]
  17. Jung, C.; Seo, J.S.; Han, S.W.; Koo, Y.J.; Kim, C.H.; Song, S.I.; Nahm, B.H.; Choi, Y.D.; Cheong, J.J. Overexpression of AtMYB44 enhances stomatal closure to confer abiotic stress tolerance in transgenic Arabidopsis. Plant Physiol. 2008, 146, 623–635. [Google Scholar] [CrossRef]
  18. Oshima, Y.; Mitsuda, N. The MIXTA-like transcription factor MYB16 is a major regulator of cuticle formation in vegetative organs. Plant Signal Behav. 2013, 8, e26826. [Google Scholar] [CrossRef]
  19. Mehrtens, F.; Kranz, H.; Bednarek, P.; Weisshaar, B. The Arabidopsis transcription factor MYB12 is a flavonol-specific regulator of phenylpropanoid biosynthesis. Plant Physiol. 2005, 138, 1083–1096. [Google Scholar] [CrossRef]
  20. Deng, M.; Wang, Y.; Kuzma, M.; Chalifoux, M.; Tremblay, L.; Yang, S.; Ying, J.; Sample, A.; Wang, H.M.; Griffiths, R.; et al. Activation tagging identifies Arabidopsis transcription factor AtMYB68 for heat and drought tolerance at yield determining reproductive stages. Plant J. 2020, 104, 1535–1550. [Google Scholar] [CrossRef]
  21. Blanco, E.; Curci, P.L.; Manconi, A.; Sarli, A.; Zuluaga, D.L.; Sonnante, G. R2R3-MYBs in Durum Wheat: Genome-Wide Identification, Poaceae-Specific Clusters, Expression, and Regulatory Dynamics Under Abiotic Stresses. Front. Plant Sci. 2022, 13, 896945. [Google Scholar] [CrossRef]
  22. Wang, L.; Yao, L.; Hao, X.; Li, N.; Wang, Y.; Ding, C.; Lei, L.; Qian, W.; Zeng, J.; Yang, Y.; et al. Transcriptional and physiological analyses reveal the association of ROS metabolism with cold tolerance in tea plant. Environ. Exp. Bot. 2019, 160, 45–58. [Google Scholar] [CrossRef]
  23. Allaby, R.G.; Peterson, G.W.; Merriwether, D.A.; Fu, Y.B. Evidence of the domestication history of flax (Linum usitatissimum L.) from genetic diversity of the sad2 locus. Theor. Appl. Genet. 2005, 112, 58–65. [Google Scholar] [CrossRef]
  24. Sadak, M.S.; Bakry, B.A. Alleviation of drought stress by melatonin foliar treatment on two flax varieties under sandy soil. Physiol. Mol. Biol. Plants 2020, 26, 907–919. [Google Scholar] [CrossRef]
  25. Wang, N.; Lin, Y.; Qi, F.; Xiaoyang, C.; Peng, Z.; Yu, Y.; Liu, Y.; Zhang, J.; Qi, X.; Deyholos, M.; et al. Comprehensive Analysis of Differentially Expressed Genes and Epigenetic Modification-Related Expression Variation Induced by Saline Stress at Seedling Stage in Fiber and Oil Flax, Linum usitatissimum L. Plants 2022, 11, 2053. [Google Scholar] [CrossRef] [PubMed]
  26. Wang, W.; Wang, L.; Wang, L.; Tan, M.; Ogutu, C.O.; Yin, Z.; Zhou, J.; Wang, J.; Wang, L.; Yan, X. Transcriptome analysis and molecular mechanism of linseed (Linum usitatissimum L.) drought tolerance under repeated drought using single-molecule long-read sequencing. BMC Genom. 2021, 22, e109. [Google Scholar] [CrossRef]
  27. Zhang, H.; Zhu, J.; Gong, Z.; Zhu, J.-K. Abiotic stress responses in plants. Nat. Rev. Genet. 2021, 23, 104–119. [Google Scholar] [CrossRef] [PubMed]
  28. Chen, Q.; Hu, T.; Li, X.; Song, C.-P.; Zhu, J.-K.; Chen, L.; Zhao, Y. Phosphorylation of SWEET sucrose transporters regulates plant root:shoot ratio under drought. Nat. Plants 2021, 8, 68–77. [Google Scholar] [CrossRef]
  29. Du, H.; Zhang, L.; Liu, L.; Tang, X.F.; Yang, W.J.; Wu, Y.M.; Huang, Y.B.; Tang, Y.X. Biochemical and molecular characterization of plant MYB transcription factor family. Biochemistry 2009, 74, 1–11. [Google Scholar] [CrossRef]
  30. Cominelli, E.; Galbiati, M.; Vavasseur, A.; Conti, L.; Sala, T.; Vuylsteke, M.; Leonhardt, N.; Dellaporta, S.L.; Tonelli, C. A guard-cell-specific MYB transcription factor regulates stomatal movements and plant drought tolerance. Curr. Biol. 2005, 15, 1196–1200. [Google Scholar] [CrossRef]
  31. Nergui, K.; Jin, S.; Zhao, L.; Liu, X.; Xu, T.; Wei, J.; Chen, X.; Yang, Y.; Li, H.; Liu, Y.; et al. Comparative analysis of physiological, agronomic and transcriptional responses to drought stress in wheat local varieties from Mongolia and Northern China. Plant Physiol. Biochem. 2022, 170, 23–35. [Google Scholar] [CrossRef]
  32. Li, Y.; Shi, L.C.; Pei, N.C.; Cushman, S.A.; Si, Y.T. Transcriptomic responses to drought stress among natural populations provide insights into local adaptation of weeping forsythia. BMC Plant Biol. 2021, 21, e273. [Google Scholar] [CrossRef]
  33. Rahimi, Y.; Ingvarsson, P.K.; Bihamta, M.R.; Alipour, H.; Taleei, A.; Khoshnoodi Jabar Abadi, S. Characterization of Dynamic Regulatory Gene and Protein Networks in Wheat Roots Upon Perceiving Water Deficit Through Comparative Transcriptomics Survey. Front. Plant Sci. 2021, 12, e710867. [Google Scholar] [CrossRef]
  34. Ponce, O.P.; Torres, Y.; Prashar, A.; Buell, R.; Lozano, R.; Orjeda, G.; Compton, L. Transcriptome profiling shows a rapid variety-specific response in two Andigenum potato varieties under drought stress. Front. Plant Sci. 2022, 13, 1003907–1003927. [Google Scholar] [CrossRef]
  35. Zhang, X.; Ni, Y.; Xu, D.; Busta, L.; Xiao, Y.; Jetter, R.; Guo, Y. Integrative analysis of the cuticular lipidome and transcriptome of Sorghum bicolor reveals cultivar differences in drought tolerance. Plant Physiol. Biochem. 2021, 163, 285–295. [Google Scholar] [CrossRef]
  36. Niu, Y.; Li, J.; Sun, F.; Song, T.; Han, B.; Liu, Z.; Su, P. Comparative transcriptome analysis reveals the key genes and pathways involved in drought stress response of two wheat (Triticum aestivum L.) varieties. Genomics 2023, 115, 110688. [Google Scholar] [CrossRef] [PubMed]
  37. Zhang, F.; Wu, J.; Bateer, S.; Yi, L.; Zhao, X.; Lu, Z. Genetic Diversity Analysis of Drought Tolerant Flax Varieties Based on Genome-wide SNP. In Proceedings of the 2022 International Conference on Biomedical and Intelligent Systems, Chengdu, China, 26 August 2022; p. 124580Q. [Google Scholar]
  38. Lammers, A.B.N.C.; Garcia, H.G.; Eisen, M.B. Unified bursting strategies in ectopic and endogenous even-skipped expression patterns. bioRxiv 2023. [Google Scholar] [CrossRef]
  39. Liu, S.C.; Jin, J.Q.; Ma, J.Q.; Yao, M.Z.; Ma, C.L.; Li, C.F.; Ding, Z.T.; Chen, L. Transcriptomic Analysis of Tea Plant Responding to Drought Stress and Recovery. PLoS ONE 2016, 11, e0147306. [Google Scholar] [CrossRef] [PubMed]
  40. Gigolashvili, T.; Yatusevich, R.; Berger, B.; Müller, C.; Flügge, U.-I. The R2R3-MYB transcription factor HAG1/MYB28 is a regulator of methionine-derived glucosinolate biosynthesis inArabidopsis thaliana. Plant J. 2007, 51, 247–261. [Google Scholar] [CrossRef] [PubMed]
  41. Malitsky, S.; Blum, E.; Less, H.; Venger, I.; Elbaz, M.; Morin, S.; Eshed, Y.; Aharoni, A. The Transcript and Metabolite Networks Affected by the Two Clades of Arabidopsis Glucosinolate Biosynthesis Regulators. Plant Physiol. 2008, 148, 2021–2049. [Google Scholar] [CrossRef] [PubMed]
  42. Sønderby, I.E.; Burow, M.; Rowe, H.C.; Kliebenstein, D.J.; Halkier, B.A. A Complex Interplay of Three R2R3 MYB Transcription Factors Determines the Profile of Aliphatic Glucosinolates in Arabidopsis. Plant Physiol. 2010, 153, 348–363. [Google Scholar] [CrossRef] [PubMed]
  43. Frerigmann, H.; Gigolashvili, T. MYB34, MYB51, and MYB122 distinctly regulate indolic glucosinolate biosynthesis in Arabidopsis thaliana. Mol. Plant 2014, 7, 814–828. [Google Scholar] [CrossRef]
  44. Frerigmann, H.; Piślewska-Bednarek, M.; Sánchez-Vallet, A.; Molina, A.; Glawischnig, E.; Gigolashvili, T.; Bednarek, P. Regulation of Pathogen-Triggered Tryptophan Metabolism in Arabidopsis thaliana by MYB Transcription Factors and Indole Glucosinolate Conversion Products. Mol. Plant 2016, 9, 682–695. [Google Scholar] [CrossRef]
  45. Chen, H.; Lai, L.; Li, L.; Liu, L.; Jakada, B.H.; Huang, Y.; He, Q.; Chai, M.; Niu, X.; Qin, Y. AcoMYB4, an Ananas comosus L. MYB Transcription Factor, Functions in Osmotic Stress through Negative Regulation of ABA Signaling. Int. J. Mol. Sci. 2020, 21, 5727. [Google Scholar] [CrossRef]
  46. Ke, Y.J.; Zheng, Q.D.; Yao, Y.H.; Ou, Y.; Chen, J.Y.; Wang, M.J.; Lai, H.P.; Yan, L.; Liu, Z.J.; Ai, Y. Genome-Wide Identification of the MYB Gene Family in Cymbidiumensifolium and Its Expression Analysis in Different Flower Colors. Int. J. Mol. Sci. 2021, 22, 13245. [Google Scholar] [CrossRef]
  47. Oh, J.E.; Kwon, Y.; Kim, J.H.; Noh, H.; Hong, S.-W.; Lee, H. A dual role for MYB60 in stomatal regulation and root growth of Arabidopsis thaliana under drought stress. Plant Mol. Biol. 2011, 77, 91–103. [Google Scholar] [CrossRef] [PubMed]
  48. Li, H.; Yang, J.; Ma, R.; An, X.; Pan, F.; Zhang, S.; Fu, Y. Genome-wide identification and expression analysis of MYB gene family in Cajanus cajan and CcMYB107 improves plant drought tolerance. Physiol. Plant. 2023, 175, e13954. [Google Scholar] [CrossRef] [PubMed]
  49. Saha, G.; Park, J.-I.; Ahmed, N.U.; Kayum, M.A.; Kang, K.-K.; Nou, I.-S. Characterization and expression profiling of MYB transcription factors against stresses and during male organ development in Chinese cabbage (Brassica rapa ssp. pekinensis). Plant Physiol. Biochem. 2016, 104, 200–215. [Google Scholar] [CrossRef]
  50. Chu, C.H.; Huang, Y.H.; Liu, H.W.; Hsu, H.M.; Tai, J.H. Membrane localization of a Myb3 transcription factor regulated by aTvCyP1 cyclophilin in the parasitic protozoanTrichomonas vaginalis. FEBS J. 2018, 285, 929–946. [Google Scholar] [CrossRef]
  51. Hussain, Q.; Asim, M.; Zhang, R.; Khan, R.; Farooq, S.; Wu, J. Transcription Factors Interact with ABA through Gene Expression and Signaling Pathways to Mitigate Drought and Salinity Stress. Biomolecules 2021, 11, 1159. [Google Scholar] [CrossRef]
  52. Gibbs, D.; Coates, J. AtMYB93 is an endodermis-specific transcriptional regulator of lateral root development in arabidopsis. Plant Signal Behav. 2014, 9, e970406. [Google Scholar] [CrossRef]
  53. Gibbs, D.; Voß, U.; Harding, S.; Fannon, J.; Moody, L.; Yamada, E.; Swarup, K.; Nibau, C.; Bassel, G.; Choudhary, A.; et al. AtMYB93 is a novel negative regulator of lateral root development in Arabidopsis. New Phytol. 2014, 203, 1194–1207. [Google Scholar] [CrossRef]
  54. Zhao, L.; Song, Z.; Wang, B.; Gao, Y.; Shi, J.; Sui, X.; Chen, X.; Zhang, Y.; Li, Y. R2R3-MYB Transcription Factor NtMYB330 Regulates Proanthocyanidin Biosynthesis and Seed Germination in Tobacco (Nicotiana tabacum L.). Front. Plant Sci. 2022, 12, 819247. [Google Scholar] [CrossRef]
  55. Leuendorf, J.E.; Schmülling, T. Meeting at the DNA: Specifying Cytokinin Responses through Transcription Factor Complex Formation. Plants 2021, 10, 1458. [Google Scholar] [CrossRef] [PubMed]
  56. Aki, S.S.; Morimoto, T.; Ohnishi, T.; Oda, A.; Kato, H.; Ishizaki, K.; Nishihama, R.; Kohchi, T.; Umeda, M. R2R3-MYB transcription factor GEMMA CUP-ASSOCIATED MYB1 mediates the cytokinin signal to achieve proper organ development in Marchantia polymorpha. Sci. Rep. 2022, 12, 21123. [Google Scholar] [CrossRef] [PubMed]
  57. Zhu, Z.; Quan, R.; Chen, G.; Yu, G.; Li, X.; Han, Z.; Xu, W.; Li, G.; Shi, J.; Li, B. An R2R3-MYB transcription factor VyMYB24, isolated from wild grape Vitis yanshanesis J. X. Chen., regulates the plant development and confers the tolerance to drought. Front. Plant Sci. 2022, 13, 966641. [Google Scholar] [CrossRef] [PubMed]
  58. Lai, L.B.; Nadeau, J.A.; Lucas, J.; Lee, E.-K.; Nakagawa, T.; Zhao, L.; Geisler, M.; Sack, F.D. The Arabidopsis R2R3 MYB Proteins FOUR LIPS and MYB88 Restrict Divisions Late in the Stomatal Cell Lineage. Plant Cell 2005, 17, 2754–2767. [Google Scholar] [CrossRef] [PubMed]
  59. Ning, C.; Yang, Y.; Chen, Q.; Zhao, W.; Zhou, X.; He, L.; Li, L.; Zong, D.; Chen, J. An R2R3 MYB transcription factor PsFLP regulates the symmetric division of guard mother cells during stomatal development in Pisum sativum. Physiol. Plant. 2023, 175, e13943. [Google Scholar] [CrossRef]
  60. Tuteja, N.; Singh, M.; Misra, M.; Bhalla, P.; Tuteja, R. Molecular mechanisms of DNA damage and repair: Progress in plants. Crit. Rev. Biochem. Mol. Biol. 2001, 36, 337–397. [Google Scholar] [CrossRef]
  61. Hwang, M.G.; Kim, K.; Lee, W.-K.; Cho, M.H. AtTBP2 and AtTRP2 in Arabidopsis encode proteins that bind plant telomeric DNA and induce DNA bending in vitro. Mol. Genet. Genom. 2005, 273, 66–75. [Google Scholar] [CrossRef] [PubMed]
  62. Mistry, J.; Finn, R.D.; Eddy, S.R.; Bateman, A.; Punta, M. Challenges in homology search: HMMER3 and convergent evolution of coiled-coil regions. Nucleic Acids Res. 2013, 41, e121. [Google Scholar] [CrossRef]
  63. Yanhui, C.; Xiaoyuan, Y.; Kun, H.; Meihua, L.; Jigang, L.; Zhaofeng, G.; Zhiqiang, L.; Yunfei, Z.; Xiaoxiao, W.; Xiaoming, Q.; et al. The MYB transcription factor superfamily of Arabidopsis: Expression analysis and phylogenetic comparison with the rice MYB family. Plant Mol. Biol. 2006, 60, 107–124. [Google Scholar] [CrossRef]
  64. Du, H.; Wang, Y.B.; Xie, Y.; Liang, Z.; Jiang, S.J.; Zhang, S.S.; Huang, Y.B.; Tang, Y.X. Genome-Wide Identification and Evolutionary and Expression Analyses of MYB-Related Genes in Land Plants. DNA Res. 2013, 20, 437–448. [Google Scholar] [CrossRef] [PubMed]
  65. Rio, D.; Ares, M.; Hannon, G.; Nilsen, T. Purification of RNA using TRIzol (TRI reagent). Cold Spring Harb. Protoc. 2010, 2010, prot5439. [Google Scholar] [CrossRef]
  66. Chen, S.; Zhou, Y.; Chen, Y.; Gu, J. fastp: An ultra-fast all-in-one FASTQ preprocessor. Bioinformatics 2018, 34, i884–i890. [Google Scholar] [CrossRef] [PubMed]
  67. Li, B.; Dewey, C. RSEM: Accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinform. 2011, 12, e323. [Google Scholar] [CrossRef] [PubMed]
  68. Yang, T.; Yang, S.; Chen, Z.; Tan, Y.; Bol, R.; Duan, H.; He, J. Global transcriptomic analysis reveals candidate genes associated with different phosphorus acquisition strategies among soybean varieties. Front. Plant Sci. 2022, 13, 1080014. [Google Scholar] [CrossRef] [PubMed]
  69. Shannon, P.; Markiel, A.; Ozier, O.; Baliga, N.S.; Wang, J.T.; Ramage, D.; Amin, N.; Schwikowski, B.; Ideker, T. Cytoscape: A Software Environment for Integrated Models of Biomolecular Interaction Networks. Genome Res. 2003, 13, 2498–2504. [Google Scholar] [CrossRef]
Figure 1. Distribution diagram of DEGs of the MYB family in 20 flax varieties.
Figure 1. Distribution diagram of DEGs of the MYB family in 20 flax varieties.
Plants 13 00710 g001
Figure 2. Cluster analysis of the expression levels in the drought group of 20 flax varieties. Each row in the figure represents a variety, and each column represents an MYB gene. The color in the heat map represents the value of the log2FoldChange. On the left is the tree diagram of gene clustering, and on the top is the tree diagram and variety code of sample clustering. Varieties marked in red were involved in the subsequent functional analysis of this study, and the color of the tree represents their clustering grouping.
Figure 2. Cluster analysis of the expression levels in the drought group of 20 flax varieties. Each row in the figure represents a variety, and each column represents an MYB gene. The color in the heat map represents the value of the log2FoldChange. On the left is the tree diagram of gene clustering, and on the top is the tree diagram and variety code of sample clustering. Varieties marked in red were involved in the subsequent functional analysis of this study, and the color of the tree represents their clustering grouping.
Plants 13 00710 g002
Figure 3. The scaling free curve of the gene co-expression module. (A): no scale tolerance curve, (B): scale-free mean connectivity curve.
Figure 3. The scaling free curve of the gene co-expression module. (A): no scale tolerance curve, (B): scale-free mean connectivity curve.
Plants 13 00710 g003
Figure 4. Flax MYB family WGCNA module classification tree. The top half of the figure (A) is a gene-level cluster tree, with each branch representing an MYB gene. The middle part of the figure (B) is the co-expression module, to which this MYB gene belongs. Each color represents a module, and gray modules represent genes that were not divided into specific modules. The lower part of the figure (C) is a heat map of the correlation between this MYB gene and the sample. Each row represents a set of treatments for a flax variety. The darker the color, the greater the correlation between this MYB gene and the sample, with red indicating a positive correlation and green indicating a negative correlation.
Figure 4. Flax MYB family WGCNA module classification tree. The top half of the figure (A) is a gene-level cluster tree, with each branch representing an MYB gene. The middle part of the figure (B) is the co-expression module, to which this MYB gene belongs. Each color represents a module, and gray modules represent genes that were not divided into specific modules. The lower part of the figure (C) is a heat map of the correlation between this MYB gene and the sample. Each row represents a set of treatments for a flax variety. The darker the color, the greater the correlation between this MYB gene and the sample, with red indicating a positive correlation and green indicating a negative correlation.
Plants 13 00710 g004
Figure 5. Group 6 flax variety MM-GS scatter plot. (A): Turquoise module MM-GS scatter diagram; (B): Blue module MM-GS scatter diagram. Module membership (MM): also known as kME (eigengene-based connectivity) correlation between the expression pattern of a gene and module characteristic genes; gene significance (GS): a measure of the importance of a gene.
Figure 5. Group 6 flax variety MM-GS scatter plot. (A): Turquoise module MM-GS scatter diagram; (B): Blue module MM-GS scatter diagram. Module membership (MM): also known as kME (eigengene-based connectivity) correlation between the expression pattern of a gene and module characteristic genes; gene significance (GS): a measure of the importance of a gene.
Plants 13 00710 g005
Figure 6. GO enrichment bubble diagram of the ‘Turquoise’ module. The vertical axis represents the GO Term, and the horizontal axis represents the Rich factor, which is the ratio of enriched genes to annotated genes in that term. The size and color of each dot indicate the number of genes and p-adjust ranges, respectively. Only results with a p-adjust < 0.5 are displayed for the top 20 enrichments.
Figure 6. GO enrichment bubble diagram of the ‘Turquoise’ module. The vertical axis represents the GO Term, and the horizontal axis represents the Rich factor, which is the ratio of enriched genes to annotated genes in that term. The size and color of each dot indicate the number of genes and p-adjust ranges, respectively. Only results with a p-adjust < 0.5 are displayed for the top 20 enrichments.
Plants 13 00710 g006
Figure 7. Turquoise module gene co-expression network. Each node in the figure represents an MYB gene and is labeled with the gene ID. The size of a node’s circle represents its degree of connectivity. The larger the circle, the higher the degree of connectivity and the more radiation edges. The inner circle represents the five nodes with the highest connectivity.
Figure 7. Turquoise module gene co-expression network. Each node in the figure represents an MYB gene and is labeled with the gene ID. The size of a node’s circle represents its degree of connectivity. The larger the circle, the higher the degree of connectivity and the more radiation edges. The inner circle represents the five nodes with the highest connectivity.
Plants 13 00710 g007
Figure 8. Expression levels of Lus.scaffold70.240 of the S (Yiya No.5) variety in the control group and the drought treatment group in q-PCR experiments and transcriptome sequencing. The left Y-axis represents the expression of this gene in the control group and the drought group in the q-PCR experiment, and the expression value is calculated by 2−∆∆CT method; the right Y-axis represents the expression of this gene in the transcriptome sequencing in the control group and the drought group, and the expression level is expressed by the TPM value.
Figure 8. Expression levels of Lus.scaffold70.240 of the S (Yiya No.5) variety in the control group and the drought treatment group in q-PCR experiments and transcriptome sequencing. The left Y-axis represents the expression of this gene in the control group and the drought group in the q-PCR experiment, and the expression value is calculated by 2−∆∆CT method; the right Y-axis represents the expression of this gene in the transcriptome sequencing in the control group and the drought group, and the expression level is expressed by the TPM value.
Plants 13 00710 g008
Figure 9. The drought shed and plants used in the experiment. (A): flax moisture treatment test in the dry shed; (B): flax plants during transcriptome sampling. The plant on the left is the Zhangya No.1 variety under drought treatment, and the plant on the right is a sample of Zhangya No.1 under control treatment.
Figure 9. The drought shed and plants used in the experiment. (A): flax moisture treatment test in the dry shed; (B): flax plants during transcriptome sampling. The plant on the left is the Zhangya No.1 variety under drought treatment, and the plant on the right is a sample of Zhangya No.1 under control treatment.
Plants 13 00710 g009
Table 1. Quantity distribution of the DEGs of 20 flax varieties.
Table 1. Quantity distribution of the DEGs of 20 flax varieties.
Variety CodeDEGsSig-UpSig-DownVariety CodeDEGsSig-UpSig-Down
A_vs_a643430K_vs_k522230
B_vs_b19387106L_vs_l27216
C_vs_c1569M_vs_m1104961
D_vs_d431627N_vs_n452421
E_vs_e1025448O_vs_o844836
F_vs_f39309P_vs_p994059
G_vs_g665115Q_vs_q281513
H_vs_h1519259R_vs_r22139
I_vs_i26188S_vs_s622240
J_vs_j883058T_vs_t241311
Note: in variety code, capital letters indicate the drought group, and lowercase letters indicate the control group. Sig-up: significantly differentially expressed genes of up-regulation; Sig-down: significantly differentially expressed genes of down-regulation.
Table 2. MYB family KEGG annotation classification statistical table.
Table 2. MYB family KEGG annotation classification statistical table.
First CategorySecond CategoryPath IDDescribeGene IDVarieties
MetabolismCarbohydrate metabolismmap00040Pentose and glucuronate interconversionsLus.scaffold440.7A, D, E, L, T
Organismal systemsEnvironmental adaptationmap04712Circadian rhythm-plantsLus.scaffold111.105B, K, L, M
Lus.scaffold292.3B, K, M
Lus.scaffold97.69B, K
Lus.scaffold80.178K
Environmental information processingSignal transductionmap04075Plant hormone signal transductionLus.scaffold70.63B
Lus.scaffold45.334B
Lus.scaffold34.70B
Note: the “Varieties” column shows flax varieties with differential expression of this gene. The green background indicates annotated KEGG information related to ‘Metabolism’. The orange background indicates annotated KEGG information related to ‘Organismal systems’. The blue background indicates annotated KEGG information related to ‘Environmental information processing’.
Table 3. GO annotation classification of the MYB family DEGs of different varieties.
Table 3. GO annotation classification of the MYB family DEGs of different varieties.
GO TypeGO TermDLTAEKBM
Biological_processBiological regulation7371724144828
Cellular process 1461752615
Response to stimulus122375148
Metabolic process12435163
Developmental process 258586
Cellular component organization or biogenesis 14172
Reproductive process 12131
Reproduction 2
Growth 22132
Localization 11
Cellular_componentCell part4327246210052191109
Organelle422623619852189108
Membrane 1 21
Membrane part 1 74
Organelle part 1 6
Protein-containing complex 1 5
Extracellular region 1 11
Molecular_functionBinding26111436553311361
Transcription Regulator Activity9251014133421
Catalytic activity12233 93
Note: The numbers in the table represent the number of MYB DEGs annotated to the GO term. The yellow background represents the flax varieties belonging to ‘Strategy 1’, the green background represents the flax variety belonging to ‘Strategy 2’, the orange background represents the flax varieties belonging to ‘Strategy 3’, and the blue background represents the flax varieties belonging to ‘Strategy 4’.
Table 4. Names and codes of the flax varieties.
Table 4. Names and codes of the flax varieties.
Variety CodeVariety NameVariety CodeVariety NameVariety CodeVariety NameVariety CodeVariety Name
A *Longza No.1FDingya No.15KLixianPBaya No.15
BZhangya No.1GZhangya No.2LNingya No.17QGao Lan Bai
CNingya No.15HDingxi No.17MNingya No.19RLongya No.10
DLongya No.8IBGOLDXREDWING44X3NDaxuan No.3SYiya No.5
EYiya No.4JR43OSha Chen Zao Shu Zhong HongT *Neiya No.9
Note: * indicates that the variety is water-sensitive.
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Zhang, F.; Liu, Y.; Ma, J.; Su, S.; Chen, L.; Cheng, Y.; Buter, S.; Zhao, X.; Yi, L.; Lu, Z. Analyzing the Diversity of MYB Family Response Strategies to Drought Stress in Different Flax Varieties Based on Transcriptome Data. Plants 2024, 13, 710. https://doi.org/10.3390/plants13050710

AMA Style

Zhang F, Liu Y, Ma J, Su S, Chen L, Cheng Y, Buter S, Zhao X, Yi L, Lu Z. Analyzing the Diversity of MYB Family Response Strategies to Drought Stress in Different Flax Varieties Based on Transcriptome Data. Plants. 2024; 13(5):710. https://doi.org/10.3390/plants13050710

Chicago/Turabian Style

Zhang, Fan, Ying Liu, Jie Ma, Shaofeng Su, Liyu Chen, Yuchen Cheng, Siqin Buter, Xiaoqing Zhao, Liuxi Yi, and Zhanyuan Lu. 2024. "Analyzing the Diversity of MYB Family Response Strategies to Drought Stress in Different Flax Varieties Based on Transcriptome Data" Plants 13, no. 5: 710. https://doi.org/10.3390/plants13050710

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop