**A MYB-Related Transcription Factor from** *Lilium lancifolium* **L. (LlMYB3) Is Involved in Anthocyanin Biosynthesis Pathway and Enhances Multiple Abiotic Stress Tolerance in** *Arabidopsis thaliana*

#### **Yubing Yong, Yue Zhang and Yingmin Lyu \***

Beijing Key Laboratory of Ornamental Germplasm Innovation and Molecular Breeding, National Engineering Research Center for Floriculture, College of Landscape Architecture, Beijing Forestry University, Beijing 100083, China

**\*** Correspondence: luyingmin@bjfu.edu.cn

Received: 13 May 2019; Accepted: 27 June 2019; Published: 29 June 2019

**Abstract:** Most commercial cultivars of lily are sensitive to abiotic stresses. However, tiger lily (*Lilium lancifolium* L.), one of the most widely distributed wild lilies in Asia, has strong abiotic stresses resistance. Thus, it is indispensable to identify stress-responsive candidate genes in tiger lily for the stress resistance improvement of plants. In this study, a MYB related homolog (LlMYB3) from tiger lily was functionally characterized as a positive regulator in plant stress tolerance. LlMYB3 is a nuclear protein with transcriptional activation activity at C-terminus. The expression of *LlMYB3* gene was induced by multiple stress treatments. Several stress-related cis-acting regulatory elements (MYBRS, MYCRS, LTRE and DRE/CRT) were located within the promoter of *LlMYB3*; however, the promoter activity was not induced sufficiently by various stresses treatments. Overexpressing *LlMYB3* in *Arabidopsis thaliana* L. transgenic plants showed ABA hypersensitivity and enhanced tolerance to cold, drought, and salt stresses. Furthermore, we found *LlMYB3* highly co-expressed with *LlCHS2* gene under cold treatment; yeast one-hybrid (Y1H) assays demonstrated LlMYB3 was able to bind to the promoter of *LlCHS2*. These findings suggest that the stress-responsive LlMYB3 may be involved in anthocyanin biosynthesis pathway to regulate stress tolerance of tiger lily.

**Keywords:** MYB; CHS; anthocyanin biosynthesis; abiotic stresses; lily

#### **1. Introduction**

Environmental stresses, such as drought, high salinity, and extreme temperatures, adversely affect the growth, development, and productivity of plants. To adapt to these environmental stressors, plants have evolved complex signaling cascades to regulate the expression of stress-related genes that can improve stress tolerance either directly or indirectly [1]. Many proteins and genes in the complex signaling networks are regulated by multiple transcription factor (TF) families. As one of the largest TF groups in plants, the MYB family has been proven to be essential for responding to abiotic stresses [2,3].

MYB proteins are characterized by their highly conserved DNA-binding domains [2,4]. According to the number of imperfect repeats of the SANT (for SWI3, ADA2, N-CoR, and TFIIIB) domain (50–53 amino acids) inMYB DNA-binding domains, plantMYB proteins can be mainly sub divided into three subfamilies: MYB-related (one single SANT domain), the R2R3-MYB (two SANT domains) and R1R2R3-MYB (three SANT domains) [5]. Accordingly, research on *MYB* genes has mainly focused on the R2R3-MYB gene family, which has been shown to play important roles in many plant-specific processes including the response to abiotic stress in past decades [6]. In *Arabidopsis thaliana* L., AtMYB14 and AtMYB15 enhance freezing tolerance by regulating *CBF* and its downstream target genes [7,8]; AtMYB20 and AtMYB44 confer salt and drought resistance respectively by downregulating the expression of *PP2Cs*[9,10]; AtMYB2, AtMYB15 and AtMYB96 function in the ABA-mediated drought and salt stress response [11–13]. In rice (*Oryza sativa* L.), OsMYB4 was reported to be a positive regulator in transgenic Arabidopsis, tomato (*Solanum lycopersicum* L. cv. Tondino), and apple (*Malus pumila* Mill.) to cold and drought tolerance [14–16]; Ectopic expression of *OsMYB2* facilitated salt, cold, dehydration tolerance in rice [17]. In wheat (*Triticum aestivum* L.), overexpression of *TaMYBsm1*, *TaMYB33*, *TaMYB2A* and *TaMYB30-B* have been shown to improve the drought tolerance in Arabidopsis [18–21]; TaMYB73 can improve salinity stress tolerance in Arabidopsis [22]; TaMYB19 has been found to participate in responses to abiotic stress in transgenic Arabidopsis [23]. However, compared to 2R-MYB genes, there are few reports of functional studies of other MYB subfamilies in abiotic stress response in plants [5].

MYB TFs have also been identified to be the major determinant regulators in anthocyanin biosynthesis [24]. Interestingly, a cross-talk exists between abiotic stresses responses and anthocyanin biosynthesis. For instance, low temperature induced and high temperature suppressed anthocyanin biosynthesis in Arabidopsis, which involved the altered regulation of *AtMYB3*, *AtMYB6* and *AtMYBL2* [25,26]. Overexpression of *AtPAP1* or *AtMYB12*, two flavonol synthesis regulators, enhances oxidative and drought tolerance in Arabidopsis [27]. *MaMYB10* in apple, *PcMYB10* in pear (*Pyrus communis* L.) and *BrMYB2-2* in *Brassica rapa* are responsible for the temperature affected anthocyanin accumulation [28–30]. Thus, it is pertinent to propose that some MYB proteins may be involved in the correlation between anthocyanin accumulation and abiotic stress tolerance.

As a wild stress-resistant plant, tiger lily (*Lilium lancifolium* L.) has been shown to have capacity for resisting multiple abiotic stresses [31,32], which could be an ideal model to study stress tolerance mechanisms and signaling regulation of a stress-resistant plant. Based on our RNA-seq data published previously [32], we isolated a MYB-related type gene, *LlMYB3*, from *L. lancifolium*, to further study the role of LlMYB3 in plant stress response. Our work showed that *LlMYB3* was induced by cold, drought, salt and exogenous ABA treatments. Ectopic expression of *LlMYB3* could improve cold, drought and salt tolerance by upregulating the transcription of several stress-related genes in Arabidopsis. Moreover, LlMYB3 TF might be involved in anthocyanin biosynthesis, which can bind to the promoter of *LlCHS2*. These results provide valuable insights into the role of the LlMYB3 in regulating plant stress response.

#### **2. Results**

#### *2.1. Isolation of LlMYB3 and Sequence Analysis*

*LlMYB3* gene comprises 1205 bp nucleotides with 462 bp open reading frame, 341 bp 5 UTR, 402 bp 3 UTR. It encodes a putative protein of 153 amino acids with a calculated molecular mass of 17.71kD and a pI of 9.54. Amino acid analysis revealed that the LlMYB3 is a MYB-related type protein with one single conserved SANT domain at the N-terminus between 39 and 86 amino acids (Figure 1a). A phylogenetic tree based on the amino acid sequences of some well-studied MYB proteins was constructed [33], which revealed that LlMYB was clustered closely to OsMYB4, OsMYB2 and AtMYB41 (Figure 1b). Furthermore, we BLAST the DNA sequence of *LlMYB3* to the whole-genome sequence of Arabidopsis thaliana using The Arabidopsis Information Resource (TAIR) to locate the chromosomal location. The BLAST result showed that the highest score (bits) significant alignment of *LlMYB3* was *AT5G62320.1* which was located in No.5 chromosome (Figure S1).

**Figure 1.** Characterization of tiger lily LlMYB3 protein. (**a**) Alignment of LlMYB3 with *Antirrhinum majus* AmMYB308, *Arabidopsis* AtMYB3, AtMYB4 and AtMYB5. The conserved MYB domain is marked and identical amino acids are shaded in black. (**b**) Phylogenetic tree analysis of LlMYB3 with other known stress-responsive MYB proteins. The LlMYB3 is marked with red triangle. (**c**) GFP and GFP-LlMYB3 fusion proteins were transiently expressed in tobacco leaves under control of the CaMV35S promoter and observed under a laser scanning confocal microscope. Scale bars for 35S-GFP, 50 μm; for 35S-LlMYB3-GFP, 25 μm. (**d**) Full-length protein (LlMYB3-A), N-terminal fragment (LlMYB3-N) and C-terminal fragment (LlMYB3-C) were fused with GAL4 DNA binding domain and expressed in yeast strain Y2HGold. The pGBDKT7 vector was used as a negative control. Transformed yeasts were dripped on the SD/-Trp, SD/-Trp-His-x-gal and SD/-Trp-His-Ade, after being cultured for 3 days in the growth chamber.

#### *2.2. Subcellular Localization and Transactivation Assay of LlMYB3*

The GFP-LlMYB3 fusion construct and the GFP control in pBI121-GFP vector driven by CaMV35S promoter were transiently expressed in tobacco epidermal cells and visualized under a laser scanning confocal microscope to determine the subcellular localization of LlMYB3. Results showed that the fluorescence signals from GFP alone were widely distributed throughout the cells, whereas, the GFP-LlMYB3 fusion protein fluorescence signal was mainly detected in the nucleus (Figure 1c). Thus, these results demonstrated that LlMYB3 is a nuclear protein.

To investigate the transcriptional activity of LlMYB3 protein, the entire coding region, N-terminal and C-terminal domain coding sequence were inserted into the pGBDKT7 vector, which contains the GAL4 DNA-binding domain. The transactivation results showed that all transformed yeast cells grew well on SD/-Trp medium (Figure 1d). The yeast strain containing the full-length LlMYB3 (LlMYB3-A) and the C-terminus of LlMYB3 (LlMYB3-C) could grow well on the selection medium SD/-Trp/-His/-Ade, while the cells with the N-terminus of LlMYB3 (LlMYB3-N) and pGBDKT7 empty vector could not grow normally (Figure 1d). Furthermore, the yeast cells that grew well on the SD/-Trp/-His-x-α-gal medium appeared blue in the presence of α-galactosidase, indicating the activation of the reporter gene *Mel1* (Figure 1d). These results indicated that LlMYB3 is a transcriptional activator, and its transactivation domain locates in the C-terminal region.

#### *2.3. Expression Patterns of LlMYB3 under Multiple Stresses and ABA*

To explore the possible involvement of *LlMYB3* in abiotic stress response, we analyzed the expression patterns of *LlMYB3* in tiger lily plants after treatment with abiotic stresses and ABA. The qRT-PCR analyses revealed that the *LlMYB3* gene has relatively high expression levels in bulb and flower petal, while its expression was low in the leaf and stem (Figure 2a). The expression of *LlMYB3* was significantly and rapidly induced within 1 h after cold, salt and drought treatment, leading to sevenfold to ninefold increase, but only in the salt treatment, the transcript level of *LlMYB3* increased again during 6–24 h (Figure 2c–e). Treatment of tiger lily plants with ABA induced the expression of *LlMYB3*, showing fivefold to sixfold increases at 24 h after treatment (Figure 2b). These data indicate that *LlMYB3* is a stress-responsive MYB-related gene in tiger lily, and its expression is sensitive to cold, drought and salt signaling molecules.

**Figure 2.** Expression patterns of *LlMYB3* in tiger lily seedlings under different stress treatments. Expression patterns of *LlMYB3* in leaves, bulbs, roots, stems, and flower petal (**a**), and after ABA (**b**), cold (**c**), drought (**d**) and NaCl (**e**) treatments in leaves by qRT-PCR analysis. Transcript levels were normalized to *LlTIP1*. Values are means ± SD of three replicates. Three independent experiments were performed. Asterisks indicate a significant difference (\*\* *p* < 0.01) compared with the corresponding controls.

#### *2.4. Promoter Analysis of LlMYB3 under Multiple Stresses and ABA*

To clarify the mechanism underlying the stress-inducible expressions of *LlMYB3*, the 1374 bp upstream of ATG start codon *LlMYB3* promoter sequence was cloned and used to drive the *GUS* expression in Arabidopsis. Sequences of various putative stress-related cis-acting regulatory elements were identified, including MYB, MYC recognition sites, LTRE, DRE/CRT, TGA and ERE elements (Table 1).



By histochemical GUS staining, we only observed prominent GUS staining in leaf veins of over 2-week old transgenic *LlMYB3* promoter plants (Figure 3a). There was no obvious difference in GUS staining between stress-treated and non-treated transgenic plants. Thus, the stress inducible activity of the *LlMYB3* promoter was revealed by measuring *GUS* gene expression levels in the transgenic Arabidopsis through qRT-PCR analyses. The result showed that *GUS* gene transcript level could be induced by cold, ABA, salt and drought treatments with a maximal level at 2 and 12 h, respectively (Figure 3b). By GUS enzyme activity assay, however, only an extremely weak fluorescence signal was detected. These results indicated that the activity of *LlMYB3* promoter can be induced by multiple stress treatments, while it is not strong enough to mediate the GUS enzyme activity.

**Figure 3.** GUS activity of the *LlMYB3* promoter in transgenic Arabidopsis plants. (**a**) Beta-glucuronidase (GUS) expression in untreated transgenic Arabidopsis plants. (**b**) The *GUS* transcript levels in the leaves of the transgenic Arabidopsis under cold (4 ◦C), drought, salt and ABA treatments. The untreated transformants served as controls. There were 12 transgenic lines acquired. Values are means ±SD of three replicates. Three independent experiments were performed. Asterisks indicate a significant difference (\*\* *p* < 0.01) compared with the corresponding controls.

#### *2.5. Overexpression of LlMYB3 in Arabidopsis Improves Tolerance to Cold and Drought Stresses*

To explore the function of LlMYB3 in providing tolerance to abiotic stress in plants, transgenic Arabidopsis plants overexpressing *LlMYB3* driven by the CaMV35S promoter were generated. Two independent homozygous lines LlMYB3-5 (L5) and LlMYB3-8 (L8) with relatively high expression levels (Figure S2) were selected for the analysis.

To study the effect of *LlMYB3* overexpression on cold stress, *LlMYB3* transgenic lines and wild-type (WT) plants were grown in equal amounts of potting soil for 4 weeks under normal conditions, and cold stress was applied by being exposed to −4, −6, −8 ◦C for 12 h. The results showed that all plants grew well under −4 ◦C treatment as the same as under normal temperature 22 ◦C (Figure 4a). When the temperature decreased to −6 ◦C, most of WT plants were dead with a survival rate at approximately 20%, but over half of transgenic plants survived (Figure 4a,b). Furthermore, all WT plants were dead, whereas the survival rate for transgenic plants was observed at 30–35% under −8 ◦C treatment (Figure 4a,b). In a further experiment, 4-week-old plants were treated at 4 ◦C for 3 h, and the relative electrolyte leakage and soluble sugars were measured after treatment. As a result, the electrolyte leakage was lower in transgenic plants relative to WT plants (Figure 4c); and transgenic plants produced remarkably higher levels of soluble sugars under a chilling condition compared to WT plants (Figure 4d).

**Figure 4.** Overexpression of *LlMYB3* in Arabidopsis improves the freezing and drought tolerance. Performance of wild-type (WT) and *LlMYB3* transgenic lines after freezing (**a**) and drought (**e**) treatments. (**b**) Survival rate of plants in (**a**) under freezing temperatures (blue region) and in (**e**) after drought treatment for 30 days (red region). Each data point is the mean of four experiments, and each experiment comprises 30 plants. Relative electrolyte leakage (**c**) and soluble sugar content (**d**) in WT and *LlMYB3* transgenic lines after 4 ◦C and 16.1% PEG6000 treatments for 3 h. (**f**) Water loss rate of leaves from WT and transgenic Arabidopsis. The data represent the means from three experiments. The bars show the SD. Significant differences between the transgenic and WT lines are indicated as \* 0.01 < *p* < 0.05 and \*\* *p* < 0.01.

Similarly, to study drought stress tolerance, after withholding water for 30 days, WT plants showed visible symptoms of drought-induced damage, such as drying, wilting, and even death while some transgenic plants remained green with expanded leaves (Figure 4e). Further analyses showed that after re-watering, few WT plants survived, whereas about 78–82% of transgenic plants continued to grow (Figure 4e,f). Additionally, after being treated with 16.1% PEG6000 (−0.5 MPa) for 3 h, transgenic plants showed lower electrolyte leakage and higher levels of soluble sugars compared to WT plants (Figure 4c,d). The water-loss rates were also slightly lower in transgenic plants (L5) than in WT plants after 3 h treatment (Figure 4f).

#### *2.6. Overexpression of LlMYB3 in Arabidopsis Increases Seed Sensitivity to ABA and Tolerance to NaCl*

The salt tolerance and ABA sensitivity of *LlMYB3* transgenic plants was assessed. NaCl significantly inhibited Arabidopsis germination when the seeds were cultivated on MS medium supplemented with 50 mM NaCl (Figure 5a). Only about 30% of the WT seeds germinated in MS medium containing 50 mM NaCl while about 60–65% of transgenic plants seeds were able to germinate (Figure 5a,b). In contrast, except that WT seeds germinated slower in MS medium containing 2 μM ABA, no obvious difference was observed between the germination of WT seeds cultivated on MS medium supplemented with 0 and 2 μM ABA. However, the germination ratio of transgenic plants seeds was remarkably lower than that of WT seeds in MS medium containing 2 μM ABA (Figure 5a,b). Therefore, we suggested that *LlMYB3* transgenic plants are more tolerant to salt stresses and more hypersensitive to ABA than WT plants.

**Figure 5.** Hypersensitivity and enhancing tolerance of *LlMYB3* transgenic lines to ABA and NaCl. Germination of WT seeds of Col-0 and 35S:LlMYB3 on MS supplemented with 50 mM NaCl and 2 μM ABA after 9 days of incubation at 22 ◦C (**a**). Germination rate of seeds counted after 5, 7 and 9 days after sowing. The data represent the means from three experiments. The bars show the SD. Significant differences between the transgenic and WT lines are indicated as \*\* *p* < 0.01.

#### *2.7. Altered Expression of Stress-Responsive Genes in LlMYB3 Transgenic Plants*

The *LlMYB3* transgenic plants exhibited an improved tolerance to freezing, drought and salt stresses. We then measured the expression levels of genes involving stress response in the transgenic plants under normal conditions. Except for *AtCOR47*, transcripts of *AtRD29A*, *AtRD29B*, *AtRD20*, *AtABI5*, *AtGolS1*, *AtLEA14* and *AtAPX2* genes (NCBI accession numbers are shown in Table S1) accumulated in *LlMYB3* transgenic plants compared to WT plants (Figure 6). The enhanced expression of these genes in transgenic plants might contribute to the stronger stress tolerance, which also implied that LlMYB3 TF may confer stress tolerance through regulating various stress-responsive genes.

**Figure 6.** Expression levels of the stress-associated genes under normal condition in WT and *LlMYB3* transgenic plants. Gene-specific primers were used to detect the relative transcript levels of the stress-related genes. Values are means ± SD of three replicates. Three independent experiments were performed. Asterisks indicate a significant difference (\*\* *p* < 0.01) compared with the corresponding controls.

#### *2.8. LlMYB3 Can Bind to the Promoter of LlCHS2*

In our previous study, through the analysis of gene co-expression networks involved in cold resistance of tiger lily, we found that the *LlMYB3* was highly co-expressed with genes involved in anthocyanin biosynthesis pathway, including phenylalanine ammonia-lyase (*PAL*), cinnamic acid 4-hydroxylase (*C4H*), 4-hydroxycinnamoyl-CoA ligase (*4CL*), chalcone synthases (*CHS*) and flavonol synthase (*FLS*) [32]. In this study, the results of qRT-PCR and Pearson's correlation coefficient (*r*) confirmed that the expression pattern of *LlMYB3* was significantly similar to *LlCHS2*'s (chalcone synthase2) (*r* > 0.8) under continuous cold treatment (Figure 7a and Table S4). The *LlCHS2* gene information is shown in Figure S3. Thus, we performed the Y1H assay to explore whether there is an interaction between LlMYB3 protein and *LlCHS2* promoter. The 932 bp upstream of ATG start codon *LlCHS2* promoter sequence was cloned, and the fragment (−820 to −553) of the *LlCHS2* promoter containing four MYB binding sites was isolated (Figure S3 and Table S3). The minimal inhibitory concentration of Aureobasidin A (AbA) for bait yeast strains was found to be 200 ng·mL−<sup>1</sup> (Figure S4). Yeast cells transformed with pGADT7-LlMYB3 and pAbAi-proLlCHS2 grew well on SD/Leu plates with 200 and 250 ng·mL−<sup>1</sup> AbA (Figure 7b). This showed that LlMYB3 could bind to the promoter of *LlCHS2*; suggesting LlMYB3 is involved in the regulatory pathway of *LlCHS2*.

**Figure 7.** Yeast one-hybrid analysis of LlMYB3 binding to *LlCHS2* promoter. (**a**) Correlation analysis of expression patterns of *LlMYB3* and some structural genes involved in anthocyanin biosynthesis pathway under continuous cold treatment. (**b**) Yeast one-hybrid analysis and schematic representation ofMYB binding sites in the *LlCHS2* promoter. Yeast strain Y1HGold was co-transformed with bait (pAbAi-proLlCHS2) and a prey (pGADT7 or pGADT7-LlMYB3) construct. Interaction between bait and prey was determined by cell growth on SD medium lacking Leu in the presence of 200 and 250 ng·mL−<sup>1</sup> AbA.

#### **3. Discussion**

Considerable studies indicate that plant MYB family members play critical roles in response to abiotic stresses. However, the evidence of improved stress resistance for most *MYB* genes is mainly from model species as Arabidopsis and rice. In this study, a cold stress-responsive gene *LlMYB3* was cloned and characterized from tiger lily (*Lilium lancifolium* L.), to investigate the role of this *MYB* gene response to various abiotic stresses. Sequence analysis shows that LlMYB3 is a MYB-related type protein with one single conserved SANT domain and displays high identity with reported stress responsive MYB member OsMYB4, OsMYB2 and AtMYB41. As predicted that LlMYB3 is a TF, LlMYB3 protein is localized in the nucleus and both the C-terminal and full-length *LlMYB3* have high transactivation ability in yeast. Our previous transcriptome data analysis identified a unigene contig 10,499 coding for LlMYB3, which showed significant changes in expression in tiger lily under cold stress [31,32]. Here, we further confirmed that *LlMYB3* was also up-regulated by drought, salt and ABA treatments.

Several stress-related cis-elements are present in the promoter of *LlMYB3*. The *LlMYB3* promoter activity is shown to be induced by multiple stress treatments; however, it is not strong enough to mediate the GUS activity significantly. It indicates that the *LlMYB3* promoter is not an effective stress-inducible promoter, and the expression of *LlMYB3* responding to abiotic stresses might mainly be regulated by the upstream regulatory factors. Furthermore, the promoter of *LlMYB3* shows vascular vein specific expression in transgenic Arabidopsis leaves. Usually, an expression pattern of a gene by promoter analysis could reflect its function [34]. The gene expression in leaf veins is usually regulated by the alteration of environments and physiological metabolic signals of this tissue during leaf development and growth [34]. Thus, it is assumed that *LlMYB3* gene might function in vascular tissues in leaves during stress response.

Compared to R2R3-MYB genes, few studies of the MYB-related genes in abiotic stress response have been reported in plants [23]. For instance, AtMYBC1, a 1R-MYB protein, was reported to be a repressor of freezing tolerance in a CBF independent pathway in Arabidopsis [35]. In rice, MYBS3 was shown to be essential for cold stress tolerance [36]; and overexpression of *OsMYB48-1* enhanced drought and salinity tolerance in rice [23]. In potato, the single MYB domain TF StMYB1R-1 has been shown to involve in drought tolerance by activation of drought-related genes [37]. Overexpression of a single-repeat MYB TF *AmMYB1* from grey mangrove conferred salt tolerance in transgenic tobacco [38]. In the present study, we generated transgenic Arabidopsis plants overexpressing *LlMYB3* gene under the control of the constitutive CaMV35S promoter. Both morphological and physiological evidence strongly demonstrated that transgenic lines had more pronounced tolerances to cold, drought and salt than WT. At gene transcription level, qRT-PCR analysis showed that the expression level of 7 of picked 8 stress-responsive genes, including *AtRD29A*, *AtRD29B*, *AtRD20*, *AtGSTF6*, *AtGolS1*, *AtLEA14* and *AtAPX2* genes were higher in the transgenic plants compared with those of WT. It suggests that these genes might be transcriptionally regulated directly or indirectly by LlMYB3, and overexpressed *LlMYB3* gene may enhance stress tolerance by regulating downstream stress-responsive genes in Arabidopsis.

Moreover, we showed overexpression of *LlMYB3* resulted in enhanced ABA sensitivity, which was also observed on many MYB TFs from Arabidopsis, such as AtMYB2, AtMYB15, AtMYB96 and AtDIV2 [10–12,39]. Given that the expression level and promoter activity of *LlMYB3* can be induced by ABA treatment, *LlMYB3* might be involved in ABA signaling pathway to response to stresses. However, promoter analysis showed that there are no known ABA responsive related cis-acting elements located in the promoter of *LlMYB3*. We thus assume that there are two possible reasons. The first is that novel ABA responsive related cis-elements might exist in the promoter of *LlMYB3*; the second is that the transport of ABA signaling molecules to *LlMYB3* is through a complex signaling network rather than by directly recognizing ABA responsive related cis-elements in *LlMYB3* promoter.

On the other hand, MYB proteins play essential roles by regulating the expression of a large number of anthocyanin biosynthesis genes. For example, in Arabidopsis, the anthocyanin regulators MYB75/PAP1, MYB90/PAP2, MYB113, and MYB114 control the expression of the late biosynthetic genes *DFR* and *LDOX*/*ANS* [2,40,41]; the flavonol regulators MYB12/PFG1, MYB11/PFG2, and MYB111/PFG3 regulate expression of the four early biosynthetic genes *CHS*, *CHI*, *F3H*, *F3'H* and *FLS* [42,43]. More importantly, some anthocyanin biosynthetic genes are even the direct targets of MYB proteins in response to abiotic stresses [3]. MYB1 in carrot can bind to the box-L-like sequences of phenylalanine ammonia-lyase 1 (*PAL1*) promoter specifically and activates *PAL1* under UV-B irradiation [44]. MYB134 in poplar, which is essential for wound and UV-B tolerance, regulates stress-responsive proanthocyanidin biosynthesis by binding to the promoter of proanthocyanidin biosynthetic genes, such as *ANR2* [45]. In rice, OsC1-MYB protein is shown to bind to the MYB responsive elements in the promoters of stress-induced flavonoid pathway genes *OsDFR* and *OsANS* [46]; overexpression of *OsMYB4* in transgenic Arabidopsis increases chilling and freezing tolerance by transactivating *PAL2* and other cold inducible genes [15]. Here, we found LlMYB3 can bind to the MYB binding sites in the promoter of cold-responsive *LlCHS2* gene, suggesting LlMYB3 protein may also function in the correlation between anthocyanin accumulation and cold stress tolerance.

In conclusion, LlMYB3 is a nucleus-localized transcriptional activator which is regulated by cold, drought and salt stresses and sensitive to ABA. Overexpressing *LlMYB3* in Arabidopsis showed ABA hypersensitivity and enhancing tolerance of transgenic plants to freezing, dehydration and salt conditions by up-regulate many stress-responsive genes. Furthermore, cold-responsive *LlCHS2* is a direct target of LlMYB3 TF in response to abiotic stresses. Therefore, our findings provide a novel MYB-related gene, which plays a positive role in plant stress resistance and might be involved in anthocyanin biosynthesis pathway in response to cold stress. Our future efforts will be focused on investigating the role of upstream regulatory factors in regulating expression and modulating the function of LlMYB3 under various stress conditions.

#### **4. Materials and Methods**

#### *4.1. Plant Materials*

The tiger lily seedlings preparation method is described in our previous study [32]. The bulbs of tiger lily were cleaned, disinfected, and then stored at 4 ◦C; in March, the bulbs were box-cultivated in a greenhouse (116.3◦ E, 40.0◦ N) under controlled conditions. The model plant *Arabidopsis thaliana* L. Columbia-0 (Col-0) was selected for the transgenic study of LlMYB3. Arabidopsis plants were grown in 8 cm × 8 cm plastic pots containing a 1:1 mixture of sterile peat soil and vermiculite under controlled conditions (22/16 ◦C, 16 h light/8 h dark, 65% relative humidity, and 1000lx light intensity). Seeds of *Nicotiana benthamiana* L. were planted and cultured under the same conditions.

#### *4.2. Cloning and Sequence Analysis of LlMYB3*

The complete sequence cDNA of *LlMYB3* gene was obtained from the transcriptome data of cold-treated tiger lily leave in our laboratory. Primer pairs (Table S1) were designed to amplify the coding sequence (CDS). The PCR products were cloned into pEASYT1-Blunt vector (TransGen Biotech, Beijing, China). After confirmation by sequencing, plasmid pEASYT1-LlMYB3 was used as a template for all experiments. The homolog genes of *LlMYB3* were searched through BLAST (http://www.ncbi. nlm.nih.gov/BLAST/) database [47]. Multiple sequence alignments were performed using DNAMAN (version 7). Phylogenetic tree analysis was performed using neighbor-joining method in MEGA5 software with 1000 replications [48]. The NCBI accession numbers of genes used in multiple sequence alignments and phylogenetic tree analysis are shown in Table S2. The theoretical molecular weight and isoelectric point were calculated using ExPASy (http://expasy.org/tools/protparam.html) [49].

Genomic DNA was extracted from tiger lily leaves using the DNeasy Plant Mini Kit (Qiagen, Valencia, CA, USA). The promoter of *LlMYB3* gene was isolated using a Genome Walker Kit (Clontech, Mountain View, CA, USA) with nest PCR according to the manufacturer's instructions. Conserved cis-element motifs of *LlMYB3* promoter were predicted using PLACE (http://www.dna.affrc.go.jp/PLACE/signalscan.html) databases [50].

#### *4.3. Abiotic Stresses Treatment and Quantitative Real-Time PCR Analysis*

For expression analysis of *LlMYB3* in response to abiotic stress and ABA treatment, 8-week-old tiger lily seedlings were treated with 4 ◦C, 16.1% PEG6000 (−0.5 MPa), 100 mM NaCl and 100 μM exogenous ABA for 0, 1, 3, 6, 12 and 24 h, respectively. Leaf samples were collected and immediately frozen with liquid nitrogen and stored at −80 ◦C for RNA isolation.

Total RNA was isolated using an RNAisomate RNA Easyspin isolation system (Aidlab Biotech, Beijing, China). First-strand cDNA synthesis was performed using Prime Script II 1st strand cDNA Synthesis Kit (Takara, Shiga, Japan). The qRT-PCR was performed using a Bio-Rad/CFX Connect™ Real-Time PCR Detection System (Bio-Rad, San Diego, CA, USA) with SYBR® qPCR mix (Takara, Shiga, Japan). Relative mRNA content was calculated using the 2−ΔΔ*C*<sup>t</sup> method [51] against the internal reference gene encoding tiger lily tonoplast intrinsic protein 1 (LlTIP1) [31] and Arabidopsis *Atactin* gene (NCBI accession No. NM\_112764). The primers used in this study were designed with Primer Premier 5 and are listed in Table S1. All reactions were performed in three biological replicates. Student's t-test was performed for all statistical analysis in this study. The heat-map was generated using MeV4.9 and clustered by hierarchical clustering (HCL) with default parameters [52]. Pearson's correlation coefficient (*r*) to define similarity of expression levels between *LlMYB3* and structural genes involved in anthocyanin biosynthesis pathway.

#### *4.4. Subcellular Localization and Transactivation Assay*

To determine its subcellular localization, the whole *LlMYB3* coding region without stop codon was amplified and cloned into pBI121-GFP at *Xho*I and *Sal*I by using ClonExpressII One Step Cloning Kits (Vazyme, Piscataway, NJ, USA) to produce LlMYB3-GFP fusion construct driven by CaMV35S promoter. The recombinant constructs and empty GFP vector were transformed into *Agrobacterium tumefaciens* GV3101 and infiltrated separately into tobacco (*N. benthamiana*) epidermal cells. After agro-infiltration for 32–48 h, GFP fluorescence signals were excited at 488nm and detected under Leica TCS SP8 Confocal Laser Scanning Platform (Leica SP8, Leica, Wetzlar, Germany) using a 500–530 nm emission filter.

The transactivation experiment was carried out according to the manual of Yeast Protocols Handbook (Clontech). The full-length coding region and truncated fragments N-terminus (1–231 bp) and C-terminus (232–462) of *LlMYB3* generated by PCR amplification were fused in frame to the GAL4 DNA binding domain in the vector of pGBKT7 (Invitrogen, Carlsbad, CA, USA). These constructs and negative control pGBKT7 were transformed into yeast strain Y2HGold by using Quick Easy Yeast Transformation Mix (Clontech). The transformed yeast strains were screened on the selective medium plates SD/-Trp, SD/-Trp/-His-Ade and SD/-Trp-His-x-α-gal plates. The transactivation activity was detected according to their growth status and α-galactosidase activity.

#### *4.5. Yeast One-Hybrid (Y1H) Assays*

Y1H assay was carried out using the Matchmaker™ Gold Yeast One-Hybrid System (Clontech). The *LlCHS2* promoter was amplified by genome walking nested PCR method described previously for *LIMYB3* promoter, and the fragment (−820 to −553) of *LlCHS2* promoter containing four MYB binding sites was isolated and cloned into the pAbAi (bait) vector (shown in Figure 7b and Figure S3). Full-length *LlMYB3* was inserted into pGADT7 (prey) vector yielding plasmid pGADT7-LlMYB3. The bait plasmids were linearized and transformed into the yeast strain Y1HGold. Positive yeast cells were then transformed with pGADT7-LlMYB3 plasmid. The DNA-protein interaction was determined based on the growth ability of the co-transformants on SD/-Leu medium with Aureobasidin A (AbA) according to the manual.

#### *4.6. Generation of Transgenic Arabidopsis*

The *LlMYB3* open read frame (ORF) was cloned into pBI121 vector under control of a CaMV35S promoter; the *LlMYB3* promoter region was inserted into CaMV35S-GUS vector by replacing the CaMV35S promoter. The recombinant vectors and empty GUS vector were transformed into 5-week-old Arabidopsis ecotype Col-0 plants by using *Agrobacterium tumefaciens* GV3101 and the floral-dip method [53]. Transformed seeds were selected on MS medium containing 50 mg/L kanamycin. T3-generation homozygous lines were selected for gene functional analysis.

#### *4.7. Histochemical Staining and Fluorometric GUS Assay*

Histochemical staining and fluorometric GUS assay analysis for GUS activity was carried out as described before [54]. To understand the effects of different stresses on GUS expression mediated by the *LlMYB3* promoter, transgenic *LlMYB3* Arabidopsis plants were treated with 4 ◦C, 16.1% PEG6000 (−0.5 MPa), 100 mM NaCl and 100 μM exogenous ABA for different durations before sampling. The leaves of stress-treated transgenic *LlMYB3* Arabidopsis were incubated in GUS reaction buffer (3 mg/mL X-gluc, 40 mM sodium phosphate pH 7, 10 mM EDTA, 0.1% Triton X-100, 0.5 mM ferricyanatum kalium, 0.5 mM ferrocyanatum kalium and 20% methanol). After overnight incubation at 37 ◦C, the stained samples were bleached with 70% (*v*/*v*) ethanol to remove chlorophyll. Photos of those stained samples were obtained by a Leica TL3000 Ergo microscope under white light. Leaves of stress-treated transgenic Arabidopsis were also used to exam *GUS* gene expression level by qRT-PCR, and determine GUS enzyme activity and measuring the fluorescence of 4-methylumbelliferone produced by GUS cleavage of 4-methylumbelliferyl-β-d-glucuronide (4-MUG). Protein amount was determined using a Protein Assay kit (Bio-Rad, Hercules, CA, USA) using bovine serum albumin as a standard.

#### *4.8. Evaluation of Transgenic Plants Abiotic Stress Tolerance and ABA Sensitivity*

The seeds of *LlMYB3* T3-generation homozygous lines and the wild type (WT) were sown on vermiculite soil in pots for freezing and drought treatment. There were 3-week-old seedlings at 4 ◦C for 3 h, then at −4, −6 or −8 ◦C, respectively, for 12 h. After that, the plants were kept at 4 ◦C for 3 h before transferring to a normal condition at 22 ◦C. For the drought treatment, the water intake of 3-week-old potted Arabidopsis plants in water-saturated substrate was withheld for 30 days, followed by rehydrating the seedlings for 7 days. The survival rates of transgenic and WT seedlings were statistical analyzed.

For determining the salt tolerance and ABA sensitivity in transgenic plants, Arabidopsis seeds were cultivated on MS medium supplemented with 0 and 2 μM ABA or 50 mM NaCl, respectively, under continuous light at 22 ◦C in a growth chamber. The germination rate was scored on the 9th day after planting on the plates.

#### *4.9. Measurements of Relative Electrolyte Leakage, Soluble Sugar, and Water Loss Rate*

The relative electrolyte leakage, soluble sugar content and water loss rate were evaluated following the method described previously [55,56]. The relative electrolyte leakage was evaluated by determining the relative conductivity of fresh leaves (100 mg) in solution using a conductivity detector. The anthrone-sulfuric acid colorimetry was used for determining the soluble sugar. The water loss rate was calculated related to the initial fresh weight of the leaf samples; the samples were placed on the lab bench (20−22 ◦C, humidity 45−60%) and weighed at designated time points. All the measurements were performed with ten plants in triplicate.

**Supplementary Materials:** Supplementary materials can be found at http://www.mdpi.com/1422-0067/20/13/3195/s1. Table S1 Primers used in this study. Table S2 NCBI accession numbers of genes used in multiple sequence alignments and phylogenetic tree analysis. Table S3 MYB binding sites identified in the promoter region of *LlCHS2*. Table S4 Expression pattern correlation between *LlMYB3* and anthocyanin biosynthesis structural genes under continuous cold stress. Asterisks indicate a significant difference (0.01 < \* *p* < 0.05). Figure S1 Chromosomal location

of *AT5G62320.1*. The highest score (bits) significant alignment of *LlMYB3* in the whole-genome sequence of *Arabidopsis thaliana* was *AT5G62320.1* which located in No.5 chromosome. Figure S2 Overexpression of *LlMYB3* confirmed by qRT-PCR. 12 independent homozygous were selected for the analysis. Wild type Arabidopsis is served as negative control. Values are means ± SD of three replicates. Three independent experiments were performed. Asterisks indicate a significant difference (\*\* *p* < 0.01) compared with the corresponding controls. The lines 5 and 8, which showed relative high expression levels of *LlMYB3* transcripts, were selected for further study. Figure S3 Sequence information about chalcone synthase 2 gene *LlCHS2* from tiger lily. The fragment (−820 to −553) of the *LlCHS2* promoter containing four MYB binding sites used in Y1H assay is underlined. The MYB binding sites described in Table S3 are highlighted in yellow. The partial ORF sequence of *LlCHS2* is shaded in grey. Figure S4 Minimal inhibitory concentration of Aureobasidin A (AbA) selected for bait yeast strains. 200 ng·mL−<sup>1</sup> of Aureobasidin A (AbA) was shown to be the minimal inhibitory concentration for bait (pAbAi-proLlCHS2) yeast strains.

**Author Contributions:** Y.Y. designed of the study and performed the statistical analysis and molecular experiments. Y.Z. conceived of the study, and helped to draft the manuscript. All authors read and approved the final manuscript.

**Funding:** This research was funded by China National Key Research & Development Project, grant number 2018YFD1000402 and China National Natural Science Foundation (grant no. 31672190, 31872138, 31071815 and No. 31272204).

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **Abbreviations**


#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **QTL Analysis of Resistance to High-Intensity UV-B Irradiation in Soybean (***Glycine max* **[L.] Merr.)**

**Min Young Yoon 1,**†**, Moon Young Kim 1,2, Jungmin Ha 1,2, Taeyoung Lee 1, Kyung Do Kim <sup>3</sup> and Suk-Ha Lee 1,2,\***


Received: 4 June 2019; Accepted: 3 July 2019; Published: 4 July 2019

**Abstract:** High-intensity ultraviolet-B (UV-B) irradiation is a complex abiotic stressor resulting in excessive light exposure, heat, and dehydration, thereby affecting crop yields. In the present study, we identified quantitative trait loci (QTLs) for resistance to high-intensity UV-B irradiation in soybean (*Glycine max* [L.]). We used a genotyping-by-sequencing approach using an F6 recombinant inbred line (RIL) population derived from a cross between Cheongja 3 (UV-B sensitive) and Buseok (UV-B resistant). We evaluated the degree of leaf damage by high-intensity UV-B radiation in the RIL population and identified four QTLs, *UVBR12-1*, *6-1*, *10-1*, and *14-1*, for UV-B stress resistance, together explaining 20% of the observed phenotypic variation. The genomic regions containing *UVBR12-1* and *UVBR6-1* and their syntenic blocks included other known biotic and abiotic stress-related QTLs. The QTL with the highest logarithm of odds (LOD) score of 3.76 was *UVBR12-1* on Chromosome 12, containing two genes encoding spectrin beta chain, brain (SPTBN, Glyma.12g088600) and bZIP transcription factor21/TGACG motif-binding 9 (bZIP TF21/TGA9, Glyma.12g088700). Their amino acid sequences did not differ between the mapping parents, but both genes were significantly upregulated by UV-B stress in Buseok but not in Cheongja 3. Among five genes in *UVBR6-1* on Chromosome 6, Glyma.06g319700 (encoding a leucine-rich repeat family protein) had two nonsynonymous single nucleotide polymorphisms differentiating the parental lines. Our findings offer powerful genetic resources for efficient and precise breeding programs aimed at developing resistant soybean cultivars to multiple stresses. Furthermore, functional validation of the candidate genes will improve our understanding of UV-B stress defense mechanisms.

**Keywords:** soybean; UV-B stress resistance; quantitative trait loci; spectrin beta chain, brain; bZIP transcription factor21/TGA9; leucine-rich repeat family protein; stress defense signaling

#### **1. Introduction**

Increased solar ultraviolet-B radiation (UV-B, 280–315 nm) since the late 1980s is considered a serious environmental issue owing to the lengthy expected time of the recovery of the destroyed stratospheric ozone layer to pre-1980 levels [1,2]. High-intensity UV-B radiation beyond the level of positive stimulus to sessile plants exerts multiple stresses, such as strong light, high temperatures, and dehydration, causing physiological and morphological damages including reduced photosynthetic capacity, leaf discoloration, and reduced biomass and seed yields [3,4]. According to the United Nations Environment Programme (UNEP) annual report, increased UV-B radiation in terrestrial areas reduces plant productivity by about 6% [5]. To prevent such yield losses in crop plants, genetic studies of resistance to high-intensity UV-B radiation as a complex abiotic stress and the identification of genetic elements involved in the UV-B defense response are needed for major crops across the world.

Soybean (*Glycine max* [L.] Merr.) is one of the most important crops for food, feed, energy production, and industrial resources worldwide. Under supplemental UV-B radiation, soybean cultivars show differences in physiological damage and morphological changes; UV-B-sensitive cultivars show significant yield reductions [6–8]. Quantitative trait loci (QTLs) associated with resistance to supplementary UV-B treatment have been identified on chromosomes (Chrs) 3, 6, 7, and 19 using a recombinant inbred line (RIL) population derived from a cross between Keunol (UV-B sensitive) and Iksan10 (UV-B resistant) [9]. Using the same population, the QTL *qUVBT1* on Chr 7 was identified with the 180K AXIOM SoyaSNP array and a gene encoding a UV excision repair protein was identified as a candidate gene involved in UV-B tolerance [10]. In our previous study, we identified the most resistant and most sensitive genotypes to elevated UV-B, Buseok and Cheongja 3, among 140 soybean genotypes, including 94 *G. max* and 46 *G. soja* accessions [11]. Transcriptome profiling of these two genotypes differing in UV-B resistance has revealed differentially expressed genes involved in stress defense signaling, immune responses, and reactive oxygen species metabolism [12].

In the present study, we identified QTLs associated with resistance to high-intensity UV-B irradiation using an F6 RIL population of Cheongja 3 (UV-B sensitive) × Buseok (UV-B resistant). Furthermore, we investigated nucleotide variation in genes located in the QTLs in the mapping parents and their expression levels in response to UV-B treatment.

#### **2. Results**

#### *2.1. Phenotypic Evaluation of UV-B Stress Resistance in a RIL Population of Cheongja 3* × *Buseok*

UV-B-resistant Buseok and UV-B-sensitive Cheongja 3 showed leaf damage of 26.8% and 62.4%, respectively, in response to high-intensity UV-B irradiation (Figure 1). The difference in UV-B resistance between the two parents was consistent with the results of previous studies [11,12]. The RIL population of Cheongja 3 × Buseok showed high phenotypic variation in leaf damage, ranging from 10% to 100%, with a mean value of 50.3% (Figure 1). The normal distribution of the degree of UV-B leaf damage in the RIL population (Shapiro-Wilk: W = 0.993, P = 0.665), with transgressive segregation, indicates that UV-B resistance is quantitatively regulated by multiple genes.

**Figure 1.** Distribution of the degree of leaf damage under high-intensity ultraviolet-B (UV-B) irradiation in the recombinant inbred line (RIL) population of Cheongja 3 × Buseok. X-axis indicates the degree of

leaf damage from 10% to 100%. For each degree of damage, RGB (red-green-blue) and two-color-transformed images of a representative injured leaf are shown. Blue and red in the two-color-transformed image indicate intact and damaged parts of the leaf exposed to UV-B irradiation, respectively. \*\* indicates a significant difference between Cheongja 3 and Buseok at *p* < 0.01 based on a Student's *t*-test.

#### *2.2. Genotyping-by-Sequencing Analysis and Genetic Linkage Map Construction*

Out of about 562 million genotyping-by-sequencing (GBS) reads generated from two libraries of the two parents and 174 RILs, 441,142,257 (78.4%) high-quality clean reads were obtained (Table S1). The number of reads per line ranged from 1,147,211 to 10,657,712 with an average of 3,197,053 reads. From these reads, we identified 271,254 unfiltered single nucleotide polymorphisms (SNPs) differentiating Cheongja 3 and Buseok. We used 4604 high-quality single nucleotide polymorphisms after filtering to construct a genetic linkage map of the RIL population. The final linkage map contained 3136 SNP markers evenly distributed on 20 chromosomes according to their physical locations on the soybean reference genome (Figure S1 and Table S2). The map spanned 3607.3 cM with an average interval of 1.15 cM between adjacent markers (Table S2). On average, each chromosome contained 156 markers spanning an average length of 180.4 cM. The chromosomes ranged from 65.8 cM (Chr 16) to 269.2 cM (Chr 15). The shortest chromosome, Chr 16, was most saturated, containing 95 SNP markers with an average marker density of 0.7 cM. Chr 4 had the largest intervals of 2.3 cM between adjacent markers. The longest chromosome, Chr 15, harbored 319 markers, covering 269.2 cM with an average marker interval of only 0.8 cM.

#### *2.3. Identification of QTLs for UV-B Resistance*

Based on the constructed genetic map, we identified four QTLs underlying UV-B resistance in soybean on Chr 12, 14, 10, and 6 (in order of logarithm of odds [LOD] score), together explaining 20.1% of phenotypic variation (Table 1). The QTL *UVBR12-1* was located at Chr12:7,261,406 ... 7,273,570 and had the highest LOD score of 3.8, explaining 9.5% of the phenotypic variation. The QTLs *UVBR14-1* and *UVBR10-1* on Chr 14 and 16 had LOD scores of 2.2 and 1.1 and explained 5.3% and 2.7% of phenotypic variation, respectively. On Chr 6, a minor QTL for UV-B resistance between markers Chr06:50,843,417 and Chr06:50,873,593 had the lowest *R*<sup>2</sup> value (2.6%) among all QTLs detected (Table 1).

We searched for genes located within the four QTLs *UVBR6-1*, *10-1*, *12-1*, and *14-1* according to the physical locations of SNP markers associated with the QTLs (Table 1 and Figure 2). Two genes encoding spectrin beta chain, brain (*SPTBN*, Glyma.12g088600) and bZIP transcription factor21 (*bZIP TF21*/*TGA9*, Glyma.12g088700) were located within *UVBR12-1* of the 121.6 kb region flanked by the marker positions Chr12:7273570 and Chr12:72761406 (Table 1 and Figure 2). *UVBR6-1* harbored five protein-coding genes anchored in the 30.2 kb region between Chr06:50,843,417 and Chr06:50,873,593 (Table 1 and Figure 2), including two genes (Glyma.06g319600 and Glyma.06g319700) encoding leucine-rich repeat (LRR) family proteins, one gene (Glyma.06g319800) for alfin-like 1, and two genes (Glyma.06g319900 and Glyma.06g320000) encoding family with sequence similarity 136, member A (FAM136A)-like protein (Table 1 and Figure 2). The 1.33 Mb genomic region of *UVBR10-1* (Chr10:41,185,273 ... Chr10:42,517,624) included 140 genes, and 60 genes were located in the 647.6 kb genomic region corresponding to *UVBR14-1* (Chr14:47,368,499 ... Chr14:46,720,930) (Table 1 and Table S3).



Genetic position of a QTL peak in the linkage map constructed in the present study. b Maximum-likelihood logarithm of odds (LOD) score for the individual QTL. c Allelic effect.Percent of phenotypic variance explained by the QTL. e Number of protein-coding genes within marker intervals on the basis of G. max gene models ver. 1.1. f Known stress-relatedQTLswithin2MbsurroundingtheQTLsidentifiedinthisstudy.

d *Int. J. Mol. Sci.* **2019** , *20*, 3287

**Figure 2.** LOD peaks and chromosomal locations of two QTLs *UVBR12-1* and *UVBR6-1* for resistance to high-intensity UV-B irradiation on chromosome (Chr) 12 and Chr 6. Red curves present LOD score distribution of detected QTLs on Chr 12 and Chr 6. These loci contain two and five protein-coding genes, respectively.

#### *2.4. SNPs and Gene Expression Di*ff*erences in the QTLs in the Mapping Parents*

To identify candidate genes in the QTLs associated with UV-B resistance, we investigated SNPs in two and five genes within the two UV-B resistance QTLs *UVBR12-1* and *UVBR6-1*, respectively, by a sequence analysis of the parental genotypes Cheongja 3 and Buseok (Table). In *UVBR12-1*, the two genes *SPTBN* (Glyma.12G088600) and *bZIP TF21*/*TGA9* (Glyma.12G088700) had nine and 25 SNPs between Cheongja 3 and Buseok, respectively (Table S4). In *SPTBN* (Glyma.12G088600), we detected one and five SNPs in the up- and downstream regions, respectively. We found three additional SNPs in the genic regions of *SPTBN*, one of which was a synonymous SNP in an exon (Figure 3). Among 25 total SNPs in *bZIP TF21*/*TGA9* (Glyma.12G088700), we found 15 SNPs in the genic region, consisting of three in 3 UTR, 10 in introns and two in exons (Figure 3 and Table S5). Each of the 2 kb up- and downstream regions of *bZIP TF21*/*TGA9* possessed five SNPs. Two SNPs in the coding sequence of *bZIP TF21*/*TGA9* were synonymous (Figure 3).

Among five genes within *UVBR6-1*, three genes had nucleotide differences between Cheongja 3 and Buseok (Tables S4 and S5). In Glyma.06g319700, encoding a LRR family protein, we detected two nonsynonymous missense mutations in exons and two SNPs in the 2 kb upstream region (Figure 3). In particular, we discovered amino acid changes from Phe to Leu and Arg to Gly at positions 247 and 254 of Glyma.06g319700, respectively. The two genes Glyma.06g319900 and Glyma.06g320000 (encoding FAM136A-like proteins) had eight and five SNPs in non-coding regions, respectively (Table S5).

Among five genes within *UVBR6-1*, three genes had nucleotide differences between Cheongja 3 and Buseok (Tables S4 and S5). In Glyma.06g319700, encoding a LRR family protein, we detected two nonsynonymous missense mutations in exons and two SNPs in the 2 kb upstream region (Figure 3). In particular, we discovered amino acid changes from Phe to Leu and Arg to Gly at positions 247 and 254 of Glyma.06g319700, respectively. The two genes Glyma.06g319900 and Glyma.06g320000 (encoding FAM136A-like proteins) had eight and five SNPs in non-coding regions, respectively (Table S5).

We further investigated transcriptional differences of *SPTBN* (Glyma.12G088600) and *bZIP TF21*/*TGA9* (Glyma.12G08870) on *UVBR12-1* between Cheongja 3 and Buseok by qRT-PCR (Figure 4) owing to the lack of protein sequence differences between the mapping parents (Figure 3). Both of the genes were upregulated in Buseok subjected to 6 h of UV-B treatment compared with levels in the control. Relative to UV-B-sensitive Cheongja 3, Buseok showed significantly higher expression levels of the two genes in response to 6 h of UV-B exposure.

#### *2.5. Comparisons of UV-B Stress Resistance QTLs with Known Stress-Related QTLs Based on Synteny*

The QTL with the highest LOD, *UVBR12-1*, on Chr 12 was located near *SCN39-4*, a QTL for the reaction to the soybean cyst nematode (SCN; *Heterodera glycines*) linked to Sctt009 (Table 1). The genomic region (Chr12:6,657,383 ... 11,713,748) harboring *UVBR12-1* and *SCN39-4* had a duplicated block on the same chromosome (Chr12:34,009,158 ... 36,919,294), carrying QTLs for drought and flood tolerance (*Drought tolerance 6-4* and *Flood tolerance 7-1*, respectively) linked to Sat 175 (Figure 5). Three homeologous blocks that show syntenic relationships with the two duplicated regions on Chr 12 were also detected on Chr 6, 11, and 13. The homeologous region (Chr11:24,411,218 ... 26,359,492) on Chr 11 with a median Ks value of 0.135 probably resulted from the recent whole genome duplication (WGD) event 13 million years ago [13] and only had collinearity with the beginning of the genomic region (including only *UVBR12-1*) on Chr 12. This syntenic block contained six QTLs for SCN resistance, *SCN 17-1*, *18-2*, *20-1*, *23-1*, *24-1*, and *32-2*, associated with Satt583. On Chr 13, the marker Satt554, associated with the QTL *Asian soybean rust 2-3* for resistance to Asian soybean rust, resided in another duplicated region (Chr13:39,150,251 ... 41,460,183) resulting from the recent WGD (median Ks, 0.133). In *UVBR12-1*, *bZIP TF21*/*TGA9* was retained in all three duplicated regions but *SPTBN* was only retained in the homeologous block on Chr 11 and was lost in the other blocks (Figure 5). On Chr 6, the remaining duplicated block (Chr06:47,897,433 ... 50,898,694) was mainly conserved with the middle part of the genomic region of *UVBR12-1*. This block with the UV-B resistance QTL *UVBR6-1* included another UV-B-related QTL (*Plant damage UV-B induced 1-2*) as well as biotic stress-related QTLs, including *SCN 17-3* and *20-2*, and *SDS 7-5* for sudden death syndrome resistance (reaction to *Fusarium solani* f. sp. *glycines*), which were associated with Satt371.

The UV-B resistance QTL *UVR14-1* on Chr 14 was mapped near three QTLs for iron efficiency, *Fe-e*ff*ect 3-2*, *9-2*, and *10-3* (Table 1). On Chr 10, *UVR10-1* was co-localized with several QTLs related to biotic as well as abiotic stresses, including *Flood tolerance 4-7*, *Drought tolerance 6-3*, *Sclero* (*Reaction to Sclerotinia sclerotiorum infection*) *2-23*, *4-10*, and *3-18*, and *Phytoph* (*Reaction to Phytophthora sojae infection*) *5-3*.

**Figure 3.** Structures of three genes encoding SPTBN (Glyma.12G088600), bZIP TF21/TGA9 (Glyma.12G088700), and LRR family protein (Glyma.06G319700), and single nucleotide polymorphisms (SNPs) in the mapping parents Cheongja and Buseok. Exons and untranslated regions (UTRs) are indicated by red- and blue-filled boxes, respectively. SNP positions, nucleotide replacements, and amino acid substitutions between Cheongja 3 and Buseok are presented. One and two SNPs in exons of *SPTBN* (Glyma.12G088600) and *bZIP TF21*/*TGA9* (Glyma.12G088700), respectively, were synonymous. The SNPs T/G and C/G in exons of Glyma.06G319700 (LRR family protein) caused amino acid replacements of Phe to Leu and Arg to Gly at the 247th and 254th residues, respectively.

**Figure 4.** Expression levels of *SPTBN* (Glyma.12G088600) (**a**) and *bZIP TF21*/*TGA9* (Glyma.12G088700) (**b**) in Cheongja 3 (green) and Buseok (orange). Y-axis represents the relative transcript abundance determined by qRT-PCR. Control, 0.5, and 6 h on the X-axis refer to 0, 0.5, and 6 h UV-B irradiation, respectively. Error bars represent the standard error for three independent replicates. \*, significant at *p* < 0.05.

#### **3. Discussion**

Crop plants under abiotic stresses, such as drought, salinity, and extreme temperatures, induce common cellular signaling mechanisms associated with osmotic stress, which disrupts homeostasis and alters the ion balance in the cell [14–16]. Similarly, high-intensity UV-B radiation above ambient levels induces cellular responses to nonspecific (genotoxic) damage, similar to various stress defense mechanisms but distinct from photomorphogenesis to low-dose UV-B, indicating that it is a complex environmental stressor [17–22]. We previously found that, compared with gene expression in UV-B-sensitive Cheongja 3, UV-B-resistant Buseok exhibits the upregulation of genes involved in phosphatidic acid-dependent signaling pathways as well as several downstream pathways, such as ABA signaling, mitogen-activated protein kinase cascades, and reactive oxygen species overproduction [12]. Therefore, high-density UV-B irradiation is a useful stress treatment for molecular biological and genetic studies of genes involved in extensive resistance to multiple abiotic stresses caused by climate change.

To identify the genetic elements responsible for UV-B resistance by QTL mapping, we developed a RIL population from a cross between UV-B-sensitive Cheongja 3 and UV-B-resistant Buseok. High-intensity UV-B irradiation turns most leaves yellow with red spots in UV-B-sensitive Cheongja 3, resulting in a dramatic loss in aerial dry weight, whereas Buseok retains more healthy leaves [11,12]. Our RIL population showed a wide distribution in the degree of leaf damage caused by UV-B stress (Figure 1). We identified four QTLs associated with UV-B resistance on Chr 12, 14, 10, and 6 (in order of LOD score) in this population, and all LOD scores were less than four (Table 1). Using a different UV-B-resistant soybean genotype, Iksan 10, previous QTL studies based on SSR genotyping and a SoyaSNP assay identified two major loci for UV-B resistance on Chr 19 and 7, respectively [9,10]. These previously reported QTLs do not overlap with our UV-B resistance QTLs. We identified a novel QTL, *UVBR12-1*, with a LOD score of 3.76 on Chr 12, explaining about 10% of phenotypic variation (Table 1). Thus, the genetic determinants of UV-B resistance in the two soybean genotypes Buseok and Iksan 10 are presumed to be different. However, a minor QTL, *UVBR6-1*, on Chr 6 was mapped in the vicinity of a known minor QTL *UV-B induced plant damage 1-2* associated with Satt371 [9] (Table 1 and Figure 5), which was only detected by multiple regression in the previous study [9].

The newly identified UV-B resistance QTLs were co-localized with other known QTLs for resistance to biotic as well as abiotic stresses (Table 1 and Figure 5). These biotic stresses include infections with *H. glycines*, *F. solani* f. sp. *glycines*, *S. sclerotiorum*, and *P. sojae*, and the abiotic stresses include drought, flood, and iron deficiency. Plants subjected to different combinations of multiple stresses show extensive overlap and crosstalk between stress-response signaling pathways, together with specific responses to individual stresses, indicating a high degree of complexity in plant molecular responses to external stresses [23]. Our results indicate that some genetic elements mediating resistance to UV-B stress may be shared with mechanisms underlying responses to other stresses, consistent with previous findings [12,24].

We identified candidate genes controlling resistance to high-intensity UV-B irradiation among two genes on *UVBR12-1* and six genes on *UVBR6-1* (Figures 3 and 4). The gene Glyma.12G088600 on *UVBR12-1* is a homolog of *Arabidopsis SPTBN* (AT5G22450), likely encoding a β-spectrin (βI to βV), which are actin-binding proteins in mammals [25]. Though plants, including *Arabidopsis*, lack spectrins or spectrin-like proteins [26,27], sequences with spectrin repeats and N-terminal calponin homology domains for actin binding are present in the *Arabidopsis* genome [28]. Additionally, spectrins or spectrin-like proteins are localized in plant cellular organelles, such as the Golgi apparatus [29], endoplasmic reticulum [30], and nucleus [31], and in the border plasma membrane of elongating plant cells [32,33]. In mammalian cells, actins bound to β-spectrins and other actin-binding proteins, such as Protein 4.1, Adducin, and Dematin, are connected to the junctional complex at the intracellular side of the plasma membrane [34–36]; these are implicated in signal targeting as well as the maintenance of cell shape and structure [25,36]. β-Spectrins interact with membrane phosphoinositides (PtdIns) via Pleckstrin homology (PH) domains present in a number of proteins involved in cellular signaling [37–39]. Direct evidence for a relationship between the spectrin-based membrane cytoskeleton and plant stress

signaling is very limited; however, levels of PtdIns, such as PtdIns4P and PtdIns(4,5)P2, and other phospholipid molecules in plant cells are increased by osmotic stress from salinity and dehydration [40]. Several genes encoding phosphatidylinositol phosphate kinases and phosphatidylinositol-specific phospholipase C involved in phospholipid signaling are upregulated in Buseok after exposure to high-intensity UV-B [12]. In this study, despite a lack of amino acid differences between the mapping parents (Figure 4 and Table S6), elevated *SPTBN* expression in response to UV-B stress was observed in Buseok but not in Cheongja 3 (Figure 4). Further investigations are necessary to characterize the cellular function of SPTBN in soybean using knockout and overexpression mutants.

HY5 (elongated hypocotyl 5) and HYH (HY5 homolog) are bZIP TFs induced by UV-B. They are key components of low-level UV-B photomorphogenic signaling mediated by β-propeller protein and E3 ubiquitin ligase, encoded by *UVR8* (*UV Resistance Locus 8*) and *COP1* (*Constitutively Photomorphogenic 1*), respectively [19,41,42]. The bZIP TF gene Glyma.12g088700 on *UVBR12-1* is an ortholog of *Arabidopsis* bZIP21/TGA9 (AT1G08320), which belongs to Clade IV of the TGACG motif-binding (TGA) protein family and is essential for anther development [43,44]. In a phylogenetic tree of bZIP genes from *G. max* and *A. thaliana*, the bZIP21/TGA9 orthologs Glyma.12g088700 and AT1G08320 belong to a different cluster from that including HY5 (AT5g11260) (Figure S2), corresponding to previous classifications in *A. thaliana* [43]. The expression of Glyma.12g088700 on *UVR12-1* was significantly increased by exposure to high-intensity UV-B in Buseok (Figure 4), but direct evidence for the role of this type of bZIP gene in resistance to abiotic stresses is lacking. bZIP21/TGA9 and TGA10 are only known to be involved in plant pathogen-associated molecular pattern-triggered immunity stimulated by the immunogenic peptide of bacterial flagellin (flg22) [45,46].

In *UVBR6-1*, Glyma.06g319700 showed amino acid differences between the mapping parents; the gene is orthologous to the *Arabidopsis* gene AT1G33590, which encodes a LRR family protein. The LRR motif plays a central role in recognizing different pathogen-associated molecules in the innate host defense of plants and animals [47]. Recent studies have shown that several genes encoding LRR-containing proteins are involved in the abiotic stress response [48], including *LP2* (*LEAF PANICLE 2*), *TaPRK2697* (*Triticum aestivum PROTEIN OF RECEPTOR KINASES 26697*), *AtPXL1* (*Arabidopsis thaliana PHLOEM INTERCALATED WITH XYLEM-LIKE 1*), and *LRR-RLK-VIII* (*LEUCINE-RICH REPEAT RECEPTOR-LIKE KINASE VIII*), which are upregulated by drought, salt, cold, and toxic metals, respectively, in diverse plant species [49–52].

Three candidate genes responsible for UV-B resistance in soybean were identified from two QTLs, *UVBR12-1* and *UVBR6-1*. These genes encode SPTBN, bZIP TF21/TGA9, and a LRR family protein and are likely involved in stress defense signaling; they should be experimentally verified using overexpression or knockout mutants. Further functional studies will improve our understanding of UV-B stress defense mechanisms. Furthermore, our results provide powerful genetic tools for more efficient and precise breeding programs aimed at the development of highly adaptable soybean cultivars under various abiotic stresses caused by global climate change.

#### **4. Materials and Methods**

#### *4.1. Plant Materials and DNA Extraction*

Two soybean genotypes were used as parents to develop a RIL population for genetic map construction and the QTL analysis. Buseok, a paternal line, is a UV-B-resistant genotype, and Cheongja 3, a maternal line, is a UV-B-sensitive genotype [11,12]. Artificial crossing was performed in summer 2012, and 176 F6 RILs were generated from F2 seeds using the single seed descent method from winter 2012 to spring 2016. Healthy young leaves from two parental genotypes and RILs were collected, and high-quality genomic DNA was extracted using the GeneAll® Exgene Plant SV Kit (GeneAll Biotechnology, Seoul, Korea). DNA quality was assessed by the 260/280 nm ratio using a Nanodrop 3000 spectrometer (Thermo Scientific, Wilmington, DE, USA). DNA was quantified using the Invitrogen Quant-iT PicoGreen® dsDNA Assay Kit (Life Technologies, Burlington, ON, Canada) and adjusted to 20 ng/μL.

#### *4.2. UV-B Treatment and Phenotypic Evaluation*

Soybean seeds of the mapping parents and their RILs were planted (three seeds per pot) in 3 L pots containing a 1:1 mixture of desalinated sand and commercial potting soil (Baroker, Seoul Bio Co., Seoul, Korea) in a greenhouse at the experimental farm of Seoul National University, Suwon, Korea in August 2016. Supplemental UV-B irradiation was used to treat soybean plants at the V3 growth stage two to three weeks after emergence according to a previous study [11]. Supplemental UV-B radiation was applied using G40T10E UV-B lamps (Sankyo Denki, Nagano, Japan) with a mean UV-B intensity of 5.68 <sup>±</sup> 0.4 Wm−<sup>2</sup> at the plant level under the lamps. The plants were exposed to high-intensity UV-B stress for 1 h every day at 11:00 am from 17 August to 20 August. The intensity of 1 h of UV-B irradiation was equivalent to two times 11.5 kJ/m2, the daily soybean UV-B biological effective dose [12,53]. Three sets of first and second trifoliate leaves above unifoliate leaves from three plants after UV-B treatment were collected with three replications per line. To investigate leaf color changes caused by high-intensity UV-B exposure, the collected leaf samples were scanned with an EPSON Perfection V33 scanner (Epson Inc., Long Beach, CA, USA) and scanned images were analyzed using the WinDIAS 3 Leaf Image Analysis System (DELTA-T DEVICES LTD., Cambridge, UK). The degree of damage was scored on a scale of 1–10 (where 1 = 0%–10%, 2 = 11%–20%, and up to 10 = 91%–100%; damaged area/healthy area × 100 (%)) and these scores was used in QTL mapping for UV-B resistance. Phenotypic distribution was tested for deviation from normality with the Shapiro–Wilk test.

#### *4.3. Genotyping-by-Sequencing*

A GBS technique was used to detect SNPs to genotype the RIL population. A GBS library was constructed following the protocol described by [54]. DNAs from the mapping parents and 176 RILs were digested individually with ApeKI, which recognizes a degenerate 5 bp sequence (GCWGC, where W is either A or T). Barcoded adapters were ligated with digested DNA fragments (Table S6), and ligated DNA was amplified with appropriate primers. Then, two separate libraries were constructed by pooling amplified DNA samples of 88 RILs for each library. Single-end sequencing of the GBS libraries was performed on two lanes of an Illumina HiSeq2000 instrument (Illumina Inc., San Diego, CA, USA).

#### *4.4. Sequence Analysis, Genetic Map Construction, and QTL Analysis*

Raw GBS reads were mapped against the *G. max* reference genome (Wm82.a2) downloaded from Phytozome (https://phytozome.jgi.doe.gov/pz/portal.html) using Bowtie v2.1 after sequence quality control, such as the removal of barcode, adapter, and ApeKI overhang sequences as well as reads with Phred scores of <15 for at least 80% of bases [55,56]. SNPs were called using in-house python scripts with the following filtering criteria: read depth ≥3, Q value ≥30, and missing error rate ≤10%. All SNPs were deposited in The European Variation Archive (https://www.ebi.ac.uk/ eva,PRJEB32685). Using identified SNP markers, a genetic map was constructed with JoinMap 4.1 (https://www.kyazma.nl/index.php/mc.JoinMap). SNP markers were grouped using the Kosambi mapping function, and segregation distortion of individual markers was calculated using the χ<sup>2</sup> test in JoinMap 4.1. The SNP genotyping data and degrees of leaf damage by UV-B treatment for 176 RILs of Cheongja 3 × Buseok were evaluated using QTL IciMapping v.4.1.0.0, and QTLs for resistance to high-intensity UV-B were identified by inclusive composite interval mapping (https: //www.integratedbreeding.net/386/breeding-services/more-software-tools/icimapping). To determine the statistically significant threshold for the LOD score, 1000 permutation test was performed at the 5% significance level.

#### *4.5. SNP Survey and qRT-PCR Analysis of Genes in the QTLs for UV-B Stress Resistance*

To investigate SNPs in genes located in the QTLs for UV-B resistance in Cheongja 3 and Buseok, whole genome sequences of these genotypes reported in a previous study were used [11]; the sequences are available from the website of the Crop Genomics Laboratory at Seoul National University (http://plantgenomics.snu.ac.kr/). To detect high-confidence SNPs, the following cut-off values were used: minimum read depth of 5, maximum read depth of 20, and Phred-scaled probability score above 30. A maximum read depth of 20 was applied to remove false positive SNP calls that may result from duplicated or repetitive sequences. The SNP positions were determined based on the Glyma 2.0 gene models (Wm82.a2.v1). To compare gene expression levels within the QTL *UVBR12-1* for UV-B resistance between Cheongja 3 and Buseok, gene-specific primers were designed for qRT-PCR using Primer3 (http://bioinfo.ut.ee/primer3-0.4.0/) (Table S7). Total RNA from each sample for control, 0.5 h, and 6 h UV-B treatments in Cheongja 3 and Buseok was used to synthesize cDNA using a Bio-Rad iScript™ cDNA Synthesis Kit (Hercules, CA, USA). UV-B treatments for qRT-PCR analysis followed as described in the previous RNA-seq study [12]: 0.5 h exposure (equivalent to 11.5 KJ/m2 soybean UV-BBE) induced photomorphogenic (non-damaging) cellular response to low-dose UV-B radiation, while 6 h exposure (12-times higher than daily UV-BBE) resulted in nonspecific (genotoxic) damage induced by abiotic stress. qRT-PCR was performed using the synthesized cDNA as a template and a Bio-Rad iQ™ SYBR Green Supermix Kit on a LightCycler® 480 (Roche Diagnostics, Laval, QC, Canada). The qRT-PCR mixture in a total volume of 20 μL contained 100 ng of cDNA, each primer at 300 μM, 8 μL of sterile water, and 10 μL of Bio-Rad iQ™ SYBR Green Supermix. The amplification conditions were as follows: 5 min of denaturation at 95 ◦C followed by 40 cycles of 95 ◦C for 10 s and 60 ◦C for 1 min. The samples were analyzed in triplicate, and actin was used as a reference gene for the normalization of target gene expression in soybean. Data were analyzed based on the stable expression level of the reference gene according to the method of Livak and Schmittgen [57].

#### *4.6. Co-Localization of Stress-Related QTLs with UV-B Resistance QTLs*

Abiotic and biotic stress-related QTLs for soybean were obtained from the SoyBase website (http://soybase.org/). The chromosomal positions of QTLs on soybean chromosomes were determined using marker information from version 4.0 of the soybean map from SoyBase [58]. Syntenic blocks in the soybean reference genome sequence were analyzed using MCScanX with default parameters [59], and syntenic blocks overlapping with genomic regions of approximately 2 to 5 Mbp encompassing SNP markers associated with UV-B resistance were identified.

**Supplementary Materials:** Supplementary materials can be found at http://www.mdpi.com/1422-0067/20/13/3287/s1.

**Author Contributions:** Conceptualization, M.Y.Y. and M.Y.K.; methodology, M.Y.Y. and M.Y.K.; software, M.Y.Y. and T.L.; validation, M.Y.Y.; formal analysis, M.Y.Y.; investigation, M.Y.Y.; resources, M.Y.Y.; data curation, J.H.; writing—original draft preparation, Y.M.Y. and M.Y.K.; writing—review and editing, M.Y.K., J.H., K.D.K and S.-H.L.; visualization, M.Y.Y., T.L. and M.Y.K.; supervision, S.-H.L.; project administration, S.-H.L.; funding acquisition, S.-H.L.

**Funding:** This work was supported by the Next Generation BioGreen 21 Program (Code No. PJ01322401), Rural Development Administration, Republic of Korea

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Abbreviations**



#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

International Journal of *Molecular Sciences*

#### *Article*
