To search for functional SNPs and genes for ML in Block 1, we performed a sequence alignment analysis of the 7 nonsynonymous SNPs in the LD block. In the full population there were four main haplotypes (A–D) involving more than three accessions. Three (B, C, and D) and all four haplotypes were detected in japonica and indica, respectively (**Figure 5B**). Accessions carrying haplotype C had shorter ML than those carrying B or D, with significant differences in indica. The peak non-synonymous SNP (Chr7\_13611491) was the only locus with a consistent allele (A) in long mesocotyl haplotypes B and D, whereas allele (T) was present in short mesocotyl haplotype C (**Figure 5B**). We suggest that the peak non-synonymous SNP might be the functional locus for ML in Block 1.

To check the inference we divided the population into groups carrying alleles (A) and (T) of the peak locus, and reexamined the relationship between ML and each of other 6 nonsynonymous SNPs above threshold in Block 1. In 222 and 151 accessions carrying allele (A) in indica and japonica, respectively,

Independent-sample *T*-tests for ML differences between alleles of each of the 6 non-synonymous SNPs in Block 1 using accessions carrying allele (A) of the peak signal (Chr7\_13611491). Data in parentheses show number of accessions. Manhattan plots for ML in Block 1 using (D) 373 accessions and (E) 222 *indica* accessions carrying allele (A) of the peak signal. Red points show 6 non-synonymous SNPs in Block 1. (F) Expression analysis of *OsML2* in the mesocotyls from 6 accessions with long- and short-mesocotyls. Names and ML of accessions are plotted on the *X*-axis. Red and blue bars show long- and short-mesocotyl accessions, respectively. Data represent means ± s.d (*n* = 3). Asterisks represent significant differences between accessions with long- and short-mesocotyls (\*\**p* < 0.01, Student's *t* test). (G) Haplotype analysis (left) of 12 non-synonymous SNPs in *OsML2* and comparison of ML (right) among haplotypes of *OsML2* in *indica* and *japonica*, respectively. Different letters above the violins indicate significant differences (*p* < 0.05) detected by Duncan's test.

there was no significant difference in ML between alleles of each non-synonymous SNP (**Figure 5C**). In 104 and 2 accessions carrying allele (T) in indica and japonica, we detected no allelic diversity in the other 6 non-synonymous SNPs (**Figure 5B**). We performed another complementary GWAS using accessions carrying allele (A) of the peak SNP in the full population and indica, respectively. There was no signal for ML in Block 1, and the 6 non-synonymous SNPs showed p-values <−log(p) = 1 in the full population and indica (**Figures 5D,E**). These results showed that the other 6 non-synonymous SNPs were not functional loci associated with ML, despite the higher −log(p) values in the original Non-Syn GWAS of the full population and indica (**Tables 1** and Table S8). We suggest that the signals in the other 6 non-synonymous SNPs were produced in accessions carrying allele (T) of the peak SNP (Chr7\_13611491).

We checked the expression level of the candidate LOC\_Os07g24010 including non-synonymous SNP Chr7\_13611491. There was about a 15.9-fold higher expression in short-mesocotyl accessions compared to long-mesocotyl accessions (**Figure 5F**). We renamed LOC\_Os07g24010 as OsML2; it also encoded a hypothetical protein. There were 12 non-synonymous SNPs in OsML2. Among them, six (13611166, 13611223, 13611304, 13611463, 13611491, and 13611535) showed high signals with −log(p) > 4 in the Non-Syn GWAS using the full population (Table S14). Based on the 12 nonsynonymous SNPs in OsML2, four and two haplotypes were identified in indica and japonica, respectively. There were clear differences in ML among accessions carrying haplotypes A, B, C, or D (**Figure 5G**). The results showed that OsML2 was an important expressed gene in regulating mesocotyl elongation, and several non-synonymous SNPs were associated with ML.

#### Cross-Validation of Two Genes for Mesocotyl Length

For efficient utilization of OsML1 and OsML2 in breeding for deep-sowing tolerance, it is important to uncover the genotypic effect of each haplotype and the combined haplotypes of two genes. Based on the above sequence alignment analysis of each gene, we focused on OsML1-A, OsML2-A, OsML2-B, and OsML2-D as long ML haplotypes, whereas OsML1-B, OsML1- C, and OsML2-C were considered to be short ML haplotypes. Among 12 possible allelic combinations of the two genes, eleven could be detected in 621 rice accessions. There were more genetic diversity in indica, including 11 allelic combinations, but only 3 allelic combinations were identified in japonica (**Figure 6A**). Considering the complexity of population structure and genetic background, we performed statistical analysis in indica by one-way ANOVA. The MLs were the shortest in accessions without long ML genotypes, OsML1-C + OsML2- C (0.62 cm) and OsML1-B + OsML2-C (0.87 cm), followed by accessions with only one long ML genotypes, OsML1-C + OsML2-A (1.12 cm), OsML1-C + OsML2-B (1.15 cm), OsML1-B + OsML2-A (1.35 cm), OsML1-B + OsML2-B (1.66 cm), OsML1- B + OsML2-D (1.91 cm), OsML1-C + OsML2-D (2.27 cm), while accessions with both long ML genotypes had longest ML, OsML1- A + OsML2-B (3.07 cm), OsML1-A + OsML2-A (3.69 cm), and OsML1-A + OsML2-D (3.76 cm) (**Figure 6A**). Haplotypes OsML1-A and OsML2-D formed the best combination with longest ML and should be the best genotype for use in molecular marker assisted breeding for long ML.

For validation, we identified accessions with long ML haplotypes for each gene from 3,024 cultivated rice accessions (Li J. Y. et al., 2014; Alexandrov et al., 2015). Three hundred and sixty seven and 982 accessions were detected with OsML1-A and OsML2-D, respectively, and 110 had both alleles. We selected 7 lines with OsML2-D, 2 lines with OsML1-A and OsML2-A and 1 line with OsML1-A and OsML2-D, and grew them in a plastic box with 10 cm of soil cover. The mean MLs were 4.57, 4.70, and 3.74 cm, respectively (**Figure 6B**). Unfortunately no seed of other accessions was available, but the results supported the hypothesis that OsML1 and OsML2 are functional genes and that the long ML haplotypes acting in an additive manner can be used in molecular marker assisted breeding. Additionally, we demonstrated a powerful strategy for efficient cloning of complex trait QTL that combines Non-Syn GWAS, known linkage map information, allelic frequencies differences between phenotypic pools, expression analyses and haplotype analyses.

#### Evidence of Positive Selection on OsML1 in Japonica

To explore whether the two genes were targeted by natural or human selection during rice domestication, we performed signature identification of selection using 621 cultivated and 446 wild rice accessions (Huang et al., 2012a). We extracted the sequencing data including the gene and a 2 kb promoter region, respectively, within OsML1 and OsML2. Selective signal scan were performed within the both genes using the ratio of the genetic diversity in wild rice to that in indica and japonica (πW/πI and πW/πJ), respectively (**Table 2**). High selective signal was detected within OsML1 in japonica (πW/πJ = 4.59). The Tajima's D was lowest at the value of −2.19, suggesting strong positive selection across OsML1 in japonica.

# DISCUSSION

#### Mesocotyl Elongation Is a Key Trait for Deep-Sowing Tolerance

Rice accessions for dry direct-seeding need have the ability to withstand deep planting for access to moisture or to tolerate

TABLE 2 | Nucleotide diversity and Tajima's *D* Test.

FIGURE 6 | Functional validation of haplotype combinations of *OsML1* and *OsML2*. (A) ML of different combinations of haplotypes in *OsML1* and *OsML2*. Different letters indicate significant differences (*p* < 0.05) detected by one-way ANOVA. (B) Cross-validation by phenotyping of 10 accessions representing different genotypes of *OsML1* and *OsML2*.

variable planting depth associated with machine sowing. Many studies have been carried out to identify deep-sowing tolerant cereal genotypes and to examine the relationship between mesocotyl and/or coleoptile length and deep-sowing tolerance (Turner et al., 1982; Kirby, 1993; Luo et al., 2007; Chung, 2010; Wu et al., 2015; Lu et al., 2016). Here a diverse set of 621 cultivars was tested for deep-sowing tolerance under field conditions and in subsequent depth-of-sowing experiments in plastic boxes. Our results, like previous studies, confirmed two points: (1) ML is induced by deeper soil covering (Wu et al., 2015; Lu et al., 2016), (2) and variation in ML fits a negative binomial distribution and there are few elite cultivated rice genotypes with long mesocotyls (Luo et al., 2007). Moreover, our results showed that ER was not affected greatly until the sowing depth reached 3 cm. Therefore, we suggest that new accessions for dry direct-seeding should have ML of more than 3 cm.

Comparison of ML between the two main rice subgroups showed that average ML in indica was significantly longer than in japonica. Given that indica rice was originally developed from crosses between japonica and local wild rice (Huang et al., 2012a) and that weedy and wild rice have longer ML (Cai and Morishima, 2002; Chung, 2010; Eizenga et al., 2016), we suggest that short mesocotyls are probably the result of an evolutionary change due to the long period of traditional rice cultivation by transplantation. Weedy and wild rice could be elite genetic resources for solving the problem of poor seedling establishment in direct-seeded rice (Chung, 2010). However, in the present work we also identified several germplasm with long mesocotyls and good seedling emergence under deep sown conditions and those elite germplasm can be used as parents in breeding rice accessions for dry direct-seeding.

# Non-syn GWAS Was Conducive in Identifying Functional Loci and Alleles

GWAS has become a common method in searching for candidate genes underlying target traits in rice (Huang X. et al., 2010; Huang et al., 2012b; Si et al., 2016; Yano et al., 2016). The Mixed Linear Model (MLM) with PC and kinship is efficient in decreasing false positives due to population structure and relationships among individuals (Yu et al., 2006). However, for complex traits affected by population structure, MLM weakens the real association and reduces signal strength of known genes. CMLM improves the statistical power by using group kinship (Zhang et al., 2010a). Constructing a trait-specific kinship derived from weighted SNP analysis is another strategy to improve MLM, as achieved by SUPER and FarmCPU (Wang et al., 2014; Liu et al., 2016). We derived two sets of SNPs (groups I and II) from 15 million polymorphisms according to annotated gene locations and possible influence on gene function, and then we obtained 572,511 non-synonymous SNPs and 1,801,421 SNPs located in the 5′ sequences of genes (≤1 Kb upstream of the first ATG) by this approach. Finally, we performed and compared GWAS using these groups to search for candidate genes affecting ML. this provides an alternative solution for constructing a trait-specific kinship derived from polymorphisms that can be screened out according to biological function of target traits. Additionally, this greatly reduced the computational burden by performing GWAS using these SNPs that accounted for amino acid sequence and expression.

There were two differences between Non-Syn GWAS and GWAS using a random sample of genetic markers (group III), one was the loci, the other was PC and kinship matrix derived from these loci. In our studies, we assessed their impact on the GWAS results. By comparison of GWAS between ML and SNPs in group I using different PC and kinship, Non-Syn GWAS detected the same SNPs for ML as GWAS with PC and kinship derived from group III. The results suggested that CMLM using PC and kinship derived from group I had same efficiency as CMLM using PC and kinship from group III. Further comparison of Non-Syn GWAS and GWAS using group III indicated that there were less loci in QTL regions (local LD block of peak SNP) identified by Non-Syn GWAS. Thus, we were able to quickly and conveniently identify functional loci by the integrated analysis without confounding by other loci with little likelihood of biological function.

#### OsML1 and OsML2 Are Natural Variants for Mesocotyl Length, and Could Be Useful in Breeding for Deep-Sowing Tolerance

Four of 13 QTLs in Non-Syn GWAS of the full population overlapped with loci identified in previous linkage mapping studies (Cai and Morishima, 2002; Ouyang et al., 2005; Huang C. et al., 2010; Lee H. S. et al., 2012). Among them, qFML3-2 located in chromosome 3 was previously detected four times in bi-parental mapping studies. However, none QTL in Non-Syn GWAS overlapped with association signal in previous GWAS (Wu et al., 2015; Lu et al., 2016). The most likely reason for this finding was that different depths of soil cover were applied in these GWAS (Wu et al., 2015; Lu et al., 2016).

Interestingly, both genes showed higher expression in shortmesocotyl accessions compared to long-mesocotyl accessions. Furthermore, long ML haplotypes of both genes (OsML1-A, OsML2-A, OsML2-B, and OsML2-D) could be loss-of-function or partial loss-of-function variations. Nucleotide diversity analysis and Tajima's D Test showed that OsML1 was strongly directly selected in japonica accessions. These results supported that short ML was a target of selection during the long term process of conversion from wild rice to japonica under the traditional transplanting system.

Sequence alignment of two genes also showed that there was less genetic diversity of OsML1 and OsML2 in japonica, and most of japonica rice had haplotypes OsML1-B and OsML2-B associated with short mesocotyls. However, there were significant differences among OsML1 and OsML2 haplotypes in indica. Superior haplotypes with long mesocotyls (OsML1-A and OsML2-D) can be used for molecular marker assistant selection, and those two genes can also be further cloned and analyzed, then used for gene transformation and breeding of rice accessions for deep-sowing tolerance.

# AUTHOR CONTRIBUTIONS

YZ, WZ, CJ, JL, and ZL designed the research, performed most of experiments and analyzed the data. XiaW, HX, ET, ZY, YC, and XinW performed part of the experiments. JX, YP, MR, and HZ conceived and supervised the project. YZ and ZL conceived the experiment and wrote the manuscript.

#### ACKNOWLEDGMENTS

This work was supported by Program of Science and Technology cooperation (2015DFG31900), Key Program of Hainan Department of Science and Technology (ZDYF2016217),

#### REFERENCES

and the China National Key Technologies Research and Developments Program (2015BAD02B01, 2013BAD01B02-15 and 2016YFD0100101).

#### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2018. 00332/full#supplementary-material

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Zhao, Zhao, Jiang, Wang, Xiong, Todorovska, Yin, Chen, Wang, Xie, Pan, Rashid, Zhang, Li and Li. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Prioritization of Candidate Genes in QTL Regions for Physiological and Biochemical Traits Underlying Drought Response in Barley (Hordeum vulgare L.)

Kornelia Gudys 1,2, Justyna Guzy-Wrobelska1 \*, Agnieszka Janiak 1 , Michał A. Dziurka3 , Agnieszka Ostrowska3 , Katarzyna Hura4 , Barbara Jurczyk 4 , Katarzyna Zmuda ˙ 4 , Daria Grzybkowska1 , Joanna Sróbka ´ 1 , Wojciech Urban1 , Jolanta Biesaga-Koscielniak 3 , Maria Filek 3 , Janusz Koscielniak 4 , Krzysztof Mikołajczak 5 , Piotr Ogrodowicz 5 , Karolina Krystkowiak 5,6, Anetta Kuczynska ´ 5 , Paweł Krajewski 7 and Iwona Szarejko1

#### Edited by:

Hanwei Mei, Shanghai Agrobiological Gene Center, China

#### Reviewed by:

Enrico Francia, Università degli Studi di Modena e Reggio Emilia, Italy Shabir Hussain Wani, Sher-e-Kashmir University of Agricultural Sciences and Technology, India

#### \*Correspondence:

Justyna Guzy-Wrobelska<EMAIL_ADDRESS>
#### Specialty section:

This article was submitted to Plant Breeding, a section of the journal Frontiers in Plant Science

Received: 31 October 2017 Accepted: 18 May 2018 Published: 12 June 2018

#### Citation:

Gudys K, Guzy-Wrobelska J, Janiak A, Dziurka MA, Ostrowska A, Hura K, Jurczyk B, Zmuda K, Grzybkowska D, ˙ Sróbka J, Urban W, ´ Biesaga-Koscielniak J, Filek M, Koscielniak J, Mikołajczak K, Ogrodowicz P, Krystkowiak K, Kuczynska A, Krajewski P and ´ Szarejko I (2018) Prioritization of Candidate Genes in QTL Regions for Physiological and Biochemical Traits Underlying Drought Response in Barley (Hordeum vulgare L.). Front. Plant Sci. 9:769. doi: 10.3389/fpls.2018.00769 1 Department of Genetics, Faculty of Biology and Environmental Protection, University of Silesia, Katowice, Poland, 2 Department of Botany and Nature Protection, Faculty of Biology and Environmental Protection, University of Silesia, Katowice, Poland, 3 Department of Developmental Biology, Institute of Plant Physiology, Polish Academy of Sciences, Krakow, Poland, 4 Department of Plant Physiology, Faculty of Agriculture and Economics, University of Agriculture, Krakow, Poland, 5 Department of Biotechnology, Institute of Plant Genetics, Polish Academy of Sciences, Poznan, Poland, 6 Department of Plant Functional Metabolomics, Institute of Bioorganic Chemistry, Polish Academy of Sciences, Poznan, Poland, 7 Department of Biometry and Bioinformatics, Institute of Plant Genetics, Polish Academy of Sciences, Poznan, Poland

Drought is one of the most adverse abiotic factors limiting growth and productivity of crops. Among them is barley, ranked fourth cereal worldwide in terms of harvested acreage and production. Plants have evolved various mechanisms to cope with water deficit at different biological levels, but there is an enormous challenge to decipher genes responsible for particular complex phenotypic traits, in order to develop drought tolerant crops. This work presents a comprehensive approach for elucidation of molecular mechanisms of drought tolerance in barley at the seedling stage of development. The study includes mapping of QTLs for physiological and biochemical traits associated with drought tolerance on a high-density function map, projection of QTL confidence intervals on barley physical map, and the retrievement of positional candidate genes (CGs), followed by their prioritization based on Gene Ontology (GO) enrichment analysis. A total of 64 QTLs for 25 physiological and biochemical traits that describe plant water status, photosynthetic efficiency, osmoprotectant and hormone content, as well as antioxidant activity, were positioned on a consensus map, constructed using RIL populations developed from the crosses between European and Syrian genotypes. The map contained a total of 875 SNP, SSR and CGs, spanning 941.86 cM with resolution of 1.1 cM. For the first time, QTLs for ethylene, glucose, sucrose, maltose, raffinose, α-tocopherol, γ-tocotrienol content, and catalase activity, have been mapped in barley. Based on overlapping confidence intervals of QTLs, 11 hotspots were identified that enclosed more than 60% of mapped QTLs. Genetic and physical map integration allowed the identification of 1,101 positional CGs within the confidence intervals of drought response-specific QTLs. Prioritization resulted in the designation of 143 CGs, among them were genes encoding antioxidants, carboxylic acid biosynthesis enzymes, heat shock proteins, small auxin up-regulated RNAs, nitric oxide synthase, ATP sulfurylases, and proteins involved in regulation of flowering time. This global approach may be proposed for identification of new CGs that underlies QTLs responsible for complex traits.

Keywords: barley, drought tolerance, function map, QTL, QTL hotspot, CG prioritization, GO enrichment

#### INTRODUCTION

Drought is one of the most devastating abiotic stresses which limits strongly crops growth and productivity. It varies in occurrence, duration and severity, and also from location to location, and in the same location from year to year. Increases in the frequency, severity, and affected areas of droughts are projected due to global climate changes (Baum et al., 2007; Cattivelli et al., 2011; Fang and Xiong, 2015). Plant water deficit occurs when the rate of transpiration exceeds water uptake, and is a component of several different stresses including drought (Bray, 1997). In order to survive and reduce stress damage, plants respond to drought stress with a variety of defense reactions, including stomata closure, limitation of transpiration, repression of photosynthesis and cell growth, activation of respiration, accumulation of osmoprotectants, or antioxidants, although their growth, development, and yield are usually negatively affected (Shinozaki and Yamaguchi-Shinozaki, 2007; Farooq et al., 2009; Cattivelli et al., 2011).

To adapt to water deficit conditions, plants have evolved complex strategies of drought resistance, which integrate multiple adaptations, from the cellular to the whole plant level, and often are specific for particular genotype × environment relationships. Thanks to extensive research, four major drought resistance mechanisms have been described and their indicators dissected: drought escape manifested by a short life cycle, photoperiod sensitivity and plasticity of development; drought avoidance through reduced water loss or increased water uptake; drought tolerance through osmotic adjustment and antioxidant capacity; and drought recovery related to capability to resume the growth after a complete loss of turgor pressure and leaf dehydration (Bray, 2007; Farooq et al., 2009; Chen et al., 2010; Fang and Xiong, 2015). Despite the efforts, much slower progress has been made in understanding genetic basis of drought resistance and in developing more tolerant genotypes, owing to the genetic complexity of this trait and its quantitative inheritance with hundreds of genes of small, often epistatic and/or pleiotropic effects and low heritability (Tuberosa and Salvi, 2007; Cattivelli et al., 2011; Fan et al., 2015).

Barley (Hordeum vulgare L.) is the fourth most widely cultivated cereal in terms of harvested acreage and production (FAO, 2016). It is an excellent model plant to decipher genetic background of drought resistance as it has considerable genetic adaptability to a wide range of environments and a high level of drought tolerance showed by local landraces and wild barley (H. vulgare spp. spontaneum). Being a valuable resource of alleles for adaptive traits, drought tolerant genotypes have been increasingly exploited in quantitative trait loci (QTL) studies to uncover the genetic control of multiple adaptations to drought stress (Baum et al., 2007; Cattivelli et al., 2011). Over last two decades numerous QTLs controlling agronomic performance and yield components under drought stress have been identified for barley (Teulat et al., 2001b; Baum et al., 2003; Talamé et al., 2004; von Korff et al., 2008; Cuesta-Marcos et al., 2009; Kalladan et al., 2013; Mansour et al., 2014; Tondelli et al., 2014). While these studies were oriented to decipher genomic regions important for barley breeding in order to maintain crop yield and grain quality under drought stress, they have seldom been used to elucidate the mechanisms of drought resistance at the genetic and molecular level (Baum et al., 2007; Fang and Xiong, 2015).

To fulfill this purpose, various indicator traits of physiological/biochemical processes occurring under drought have been investigated in barley using diverse mapping populations and drought stress conditions. Most studies concentrated on the parameters related to plant water status, resulting primarily in the identification of QTLs for relative water content (RWC), and then, for osmotic adjustment (OA), osmotic potential (OP), water content (WC), or carbon isotope discrimination (Teulat et al., 1997, 1998, 2001a, 2003; Diab et al., 2004; Chen et al., 2010; Szira et al., 2011; Wójcik-Jagła et al., 2013; Honsdorf et al., 2014; Fan et al., 2015; Mora et al., 2016). In contrast, only individual studies have detected QTLs for water-soluble carbohydrates (WSC, Teulat et al., 2001a; Diab et al., 2004) and proline content (Sayed et al., 2012; Fan et al., 2015). Additionally, some studies on the photosynthetic efficiency under drought stress conditions have identified QTLs related to chlorophyll fluorescence parameters, chlorophyll content or PSII (photosystem II) photochemical activity (Guo et al., 2008; von Korff et al., 2008; Wójcik-Jagła et al., 2013; Honsdorf et al., 2014; Mora et al., 2016) and plasma membrane integrity (Wójcik-Jagła et al., 2013). Recently, QTLs for multiple metabolites accumulation under drought stress (Piasecka et al., 2017), including QTLs for the content of fat-soluble antioxidants: α- and γ-tocopherols (Templer et al., 2017), have been detected.

These papers present a valuable source of data on the chromosomal regions potentially involved in drought response in barley, however, the comprehensive integration of the majority of these results and the anchoring of the identified QTLs on barley genome to decipher genes and networks for QTLs underlying particular physiological traits could be challenging. The difficulties arise not only from inherent limitations associated with particular mapping populations and applied methods but foremost from a lack of the information on the proper genomic localization of these QTLs which preclude their precise positioning on physical maps (Mir et al., 2012; Salvi and Tuberosa, 2015). Today, more comprehensive and effective approach to decipher the genetic basis of QTLs (Monclus et al., 2012; Bargsten et al., 2014; Kumar et al., 2017) may be applied also in barley. The availability of SNP (single nucleotide polymorphism) markers based on within-gene polymorphisms (Close et al., 2009) together with high-throughput genotyping methods enable production of function maps with high density genome coverage for QTL mapping (Potokina et al., 2008; Mammadov et al., 2012). Such maps provide straight-forward connection of the identified QTL intervals to the reference barley genome sequence (Mayer et al., 2012), which could be followed by the analysis of large candidate gene assemblies and their biological interpretation using public databases (e.g., Gene Ontology; GO) and bioinformatics high-throughput enrichment tools.

This study present a first comprehensive approach for elucidation of genetic basis of physiological mechanisms of drought response/tolerance in barley based on the identification of the positional candidate gene (CG) assembly within the QTL confidence intervals, followed by the exploration of their putative functions related to drought tolerance. The study aimed particularly (1) to identify QTLs for a wide range of physiological/biochemical traits representative for plant water status, photosynthetic efficiency, osmoprotectant and hormone content, as well as activity or accumulation of antioxidants under drought stress on a created high-density function map; (2) to project QTL confidence intervals on physical barley genome map; and (3) to retrieve a set of potentially causative CGs, underlying the analyzed traits using a Gene Ontology (GO) enrichment approach.

#### MATERIALS AND METHODS

#### Plant Material

Barley (H. vulgare L.) population of 100 recombinant inbred lines (RIL), named MCam, produced by the single-seed descent (SSD) method from the cross between spring genotypes "Maresi" and Cam/B1/CI08887//CI05761 (here after referred as CamB), was used in this study to phenotype a set of physiological and biochemical traits under drought stress and control conditions, followed by their QTL mapping on a high-density function map. "Maresi" is a German advanced cultivar (semidwarf, two-rowed, malting type), of high and stable yielding under the European environmental conditions (pedigree and seed source: http://genbank.vurv.cz/genetic/resources/asp2/default\_a. htm), whereas CamB is a Syrian breeding line (two-rowed, early heading), kindly obtained from ICARDA (International Center for Agricultural Research in the Dry Areas), adapted to dry environments (detailed characteristics of parental genotypes: Górny, 2001; Filek et al., 2015, 2016; Chmielewska et al., 2016).

#### Plant Growth Conditions and Drought Stress Treatment

The experiment was carried out in a phytotron (growth chamber). Seeds of each RIL were sterilized and sown separately in a box (22 L; 10 RIL per day), filled with a mixture of sandy loam and sand (7/2, v/v). In this substrate, a pF range of 2.2–3.0 indicated easily available water whereas pF > 4.2 was the permanent wilting point, as calculated based on the water retention curve (Filek et al., 2015). The pF-value is defined as a logarithm of the pressure p (expressed in centimeters of water head) necessary for removal of water from soil capillaries (Mikołajczak et al., 2017). Soil moisture was regulated two times a day. After germination (4 days, at 25◦C), the boxes containing 30 plants were kept at 5◦C (day/night) during the next 10 days, then the following parameters were settled in phytotron growth chamber: photoperiod 16/8 h (day/night), the temperature of 20/17◦C (day/night), irradiance of 520 µmol(photon)m−2 s −1 . After 10 days, the number of plants in a box was reduced to 25. The substrate humidity was determined by monitoring boxes weight, and it was kept at 11% water content (VWC), i.e., pF = 2.8. After 13 days (at the moment of the emergence of the 4th leaf) the temperature was set to 25/16◦C and soil drought (3.65% VWC, i.e., pF = 4.0) was applied to plants for 10 days. Plants grown in boxes with 11% VWC were used as the control (optimal condition growth). The boxes were randomly distributed in the growth chamber, and their location was changed twice a day. Measurements of the physiological/biochemical traits were performed at the seedling stage, on the third leaf for the plants from both (control and drought stress) conditions after 10 days of drought. Depending on the trait, the measurements were performed in 4–25 biological replications.

#### Phenotyping

Altogether, 40 different physiological/biochemical traits were measured at the seedling stage for RILs and the parents (described below and in **Table 1**). To facilitate the analysis, the traits were grouped into four categories according to their roles in drought stress response: plant water status (5 traits), photosynthetic efficiency (19 traits), osmoprotectant and hormone content (7 traits), activity and accumulation of antioxidants (9 traits).

#### Plant Water Status

Relative water content (RWC), water content (WC), water loss (WL) rate, water use efficiency (WUE), and electrolyte leakage (EL) were determined in this category. RWC was calculated according to Barrs (1968), where RWC = (FW – DW)/(TW – DW) × 100%, and where FW, DW and TW, respectively, are fresh, dry and turgid weight. To measure TW, leaf samples were placed in water (7◦C) in darkness for 24 h, for the complete rehydration. WC was calculated as (DW/FW) × 100%. WL rate was calculated on the basis of the depletion of leaf FW during 24 h under the appropriate growing conditions. WUE was calculated as Pn/E, where Pn and E are net photosynthesis and transpiration rates, respectively, estimated by the measurement of gas exchange (next section). EL test was used to determine the plasma membrane integrity. For each genotype, samples were prepared as described by Płazek et al. (2014).

#### Photosynthetic Efficiency

Photochemical efficiency was estimated by means of chlorophyll a fluorescence registration. Measurements were taken using a fast chlorophyll fluorescence induction kinetics fluorometer Handy PEA and modulated fluorescence system FMS2 (Hansatech,

#### TABLE 1 | Physiological and biochemical traits used in QTL analysis.

DW, dry weight; FW, fresh weight; U, arbitrary unit of enzyme activity per mg of total soluble protein.

Kings Lynn, UK). The induction of a chlorophyll fluorescence signal was measured after 30 min of leaf dark adaptation in clips (Hansatech). Before measurements, the LED-light source of the fluorometer was calibrated using an SQS light meter (Hansatech, Kings Lynn, UK). The conditions used to determine polyphasic chlorophyll a fluorescence transients were: excitation irradiance of 3,000 µmol m−2 s −1 , a pulse duration of 1 s, and fixed gain of 0.7 (Handy PEA ver. 1.3 software, Hansatech). Fourteen parameters were calculated (Kalaji et al., 2011): ABS/RC and ABS/CS—light absorption flux, respectively, per PSII reaction center (RC) and per excited leaf cross-section at t = 0 (CS); ET0/RC and ET0/CS (electron transport flux per RC and CS); TR0/RC and TR0/CS (trapped energy flux per RC and CS); DI0/RC and DI0/CS (dissipated energy flux per RC and CS); RC/CS (the maximum number of active RC per CS); φpo (TR0/ABS—maximal quantum yield of primary photochemistry); ψo (ET0/TR0—exciton transfer efficiency to the electron transport chain); φeo (ET0/ABS—electron transport yield); PIabs (the performance index per absorption); (1-B)av (the average fraction of open RC during the time needed to complete the closure of all RCs). FMS2 measurements were made after light adaptation of the leaf (usually 2–5 min at 500 µmol (quanta) m−2 s −1 ) when the fluorescence signal (Fs) became constant. Three parameters: the photochemical quantum yield of PSII (8PSII), the photochemical quenching (qP) and the efficiency of excitation energy capture by open PSII RC (Fv ′ /Fm ′ ) were calculated according to Genty et al. (1989). The gas exchange rate (Pn; E) was measured using an infrared gas analyzer (Ciras-1, PP Systems, Hitchin, UK) and Parkinson leaf chamber (PLC6). The controlled measuring conditions were: CO2 concentration of 400 µmol (CO2) mol−1 (air), 30% relative humidity, irradiance of 800 µmol (quanta) m−2 s −1 and the leaf temperature of 25◦C. The measurements of chlorophyll fluorescence and gas exchange were carried out from 10 a.m. to 2 p.m. (beginning of day 6 a.m., day length 16 h), and until stable measurement values were obtained. The order of plants measurements (photosynthesis and fluorescence) during each day of was random.

#### Osmoprotectant and Hormone Content

Seven compounds were quantified in this category. Soluble carbohydrates (glucose, fructose, sucrose, raffinose, and maltose) were identified by their retention times and quantified by integrating peak areas against the internal standard according to the procedure by Janeczko et al. (2010). Measurements were made using a high performance liquid chromatography consisting of the following modules: a gradient pump (Agilent 1200, Santa Clara, CA, USA), autosampler (Agilent 1200), the thermostat STH 585 (Dionex, Sunnyvale, CA, USA), ESA detector Coulochem II Analitical Cell 5040 with a gold working electrode and a palladium reference electrode, an analog/digital converter (Agilent), program control and data collector software ChemStation Rev.B.04.01 (Agilent).

The content of free proline was determined spectrophotometrically (at λ = 520 nm) according to the methods by Bates et al. (1973) and Marin et al. (2009), with the use of UV/VIS spectrophotometer (UV-1800, Rayleigh, Beijing, China) and it was calculated from the standard curve.

For ethylene production, each sample (0.1 g FW) was homogenized in 70% ethanol. Following the centrifugation, the supernatant was evaporated (at 40◦C) and the pellet was dissolved in H2O (1 ml), and 1 pmol HgCl2 was added. Then, the sample was transferred into a new vial and incubated on ice (for 5 min) with a mixture of 5% NaOCl and saturated NaOH (2/1, v/v). The gas sample was taken from the vial, and ethylene was quantified by gas chromatography (Hewlett Packard 5890 Series II, Palo Alto, CA, USA) with Porapak R column (80/100 mesh, Agilent, Santa Clara, CA, USA) and detector (flame ionization) as described by Grzesiak et al. (2013).

#### Activity and Accumulation of Antioxidants

Antioxidant enzymes (3), fat-soluble antioxidants (5) and βcarotene were quantified in this category. To measure the activity of antioxidant enzymes: superoxide dismutase (SOD), catalase (CAT), and peroxidase (POD), the plant material was homogenized at 4◦C in 0.05 mM potassium phosphate buffer (KP) with 0.1 mM EDTA. The supernatant was divided into three subsamples after centrifugation. The activity of SOD [EC <IP_ADDRESS>.] was determined spectrophotometrically by the cytochrome reduction method (McCord and Fridovich, 1969) with the modifications by Szechynska-Hebda et al. ´ (2012). Changes in absorbance were followed with a Biochrom Ultrospec II spectrophotometer (LKB, Sweden) at λ = 550 nm. The inhibition of oxidized cytochrome c absorbance was monitored. The activity of CAT [EC <IP_ADDRESS>] was measured spectrophotometrically according to Aebi (1984) with the modifications by Wojtania et al. (2016). The decrease in H2O2 absorbance was measured at λ = 240 nm. The activity of POD [EC <IP_ADDRESS>] was measured spectrophotometrically using the method by Lück (1962) with the modifications by Wojtania et al. (2016). The increase in oxidized p-phenylenediamine absorbance was monitored at λ = 485 nm.

Tocochromanols (α- γ- δ-tocopherols and α- γ-tocotrienols) and β-carotene were measured based on the modified method by Surówka et al. (2016). Briefly, lyophilized samples were extracted in 0.1% solution of butylated hydroxytoluene (BHT) in ethanol/acetone/methanol/2-propanol (8/3/3/1 v/v) at 70◦C in shaking water bath for 15 min, then 80% KOH was added, and the extraction was continued for 30 min. Next, samples were diluted with H2O (1/1 v/v) and cleaned on Chromabond XTR cartridge (Macherey-Nagel, Germany). Compounds of interest were eluted by n-hexane, vacuum evaporated (Rota Vapor, Switzerland) and reconstituted in 1% BHT in methanol/dichloromethane (3/1 v/v) prior HPLC separation. The Agilent 1260 (Santa Clara, CA, USA) UHPLC binary system with diode array (DAD) and fluorescence (FLD) detectors was used. Separation was achieved on Ascentis Express C-18 (Supelco Analytical, Sigma Aldrich, USA) analytical columns at 0.8 ml/min, 60◦C and linear gradient of A) 0.5% formic acid (FA) in acetonitrile (ACN)/H2O (6/4 v/v) and B) 0.5% FA in 2-propanol/ACN (9/1 v/v), from 40 to 100% of B in 15 min. Tocochromanols were detected by FLD, whereas β-carotene by DAD. Identity and quantity of compounds were confirmed by the comparison with data obtained for the pure standards under identical conditions.

#### Statistical Analysis

All physiological traits data were analyzed with Statistica 12.0 software (Stat. Soft Inc., USA). The normality of trait distributions was verified using the Shapiro–Wilk test. Linear correlation coefficients (Pearson's) were calculated between all the analyzed traits separately for drought and control conditions. The F-test was used to assess the homogeneity of variance and the Student's t-test to compare the statistical significance of differences of the analyzed traits between control and drought conditions.

In the QTL analysis absolute values of all physiological/biochemical traits were used which described the drought response/tolerance level of RILs, and drought stress indices (DSI), which are relative values, were also calculated for all measured traits to compare among RILs (Bouslama and Schapaugh, 1984; Wójcik-Jagła et al., 2013) according to the formula: DSI (%) = (d/c) × 100%, where d and c are the absolute values obtained under drought and control conditions, respectively. DSI data were included into QTL mapping as a set of additional phenotypic traits.

# Genotyping

Mapping Functional CGs on a Consensus Barley Map A consensus barley linkage map of SNP markers from BOPA1 (barley oligonucleotide pool assay 1; Close et al., 2009) and simple sequence repeat (SSR) markers, published earlier (Mikołajczak et al., 2016), was used in this study for the map enrichment with functional candidate genes (CG) and the construction of a new high-density, function barley map for the QTL analysis. Briefly, the consensus map was constructed for the MCam population and two other bi-parental RILs: LCam (derived from "Lubuski" × CamB) and GH (derived from "Georgie" × "Harmal"), produced from the crosses between European (drought susceptible) and Syrian (drought tolerant) spring barley cultivars. The map consisted of 819 markers, spanned 953.8 cM (an average resolution of 1.2 cM) and comprised of 13 linkage groups attributed to barley chromosomes 1H-7H. The map was uniformly covered with the markers, and was in good agreement with other integrated barley maps (Varshney et al., 2007; Close et al., 2009). The number of linkage groups exceeding the number of barley chromosomes was probably the effect of a diverse genetic background of parental genotypes, but it did not affect the marker order within the linkage groups, nor the CGs mapping and the overall quality of the map.

For the enrichment of the map, a group of 41 CGs potentially involved in drought response was selected based on the transcriptome analysis of "Maresi" and CamB genotypes under drought stress (Janiak et al., 2018). In addition, 23 genes encoding barley orthologs of drought tolerance-related genes, described in model species, were chosen (Supplementary Material S1). Publicly available genomic sequence information in GenBank (http://www.ncbi.nlm.nih.gov/genbank/) was utilized to identify barley genomic sequences of these CGs, so called functional CGs (Pflieger et al., 2001). The information on the structure and the chromosomal localizations of the CGs was derived from Ensembl Plants (http://plants.ensembl.org/index.html; ver. 082214v1).

For polymorphisms identification in the functional CG sequences, total genomic DNAs of two RIL populations, MCam and GH, and their parental forms, were extracted using a method by Doyle and Doyle (1987) with minor modifications. Primer pairs (Supplementary Material S2) for the PCR amplifications of CG fragments (800–1,200 bp) from genomic DNA of parental genotypes were designed with the software Primer 3 (http:// bioinfo.ut.ee/primer3-0.4.0/). PCR products were sequenced from both ends at Genomed company (Warsaw, Poland; www.genomed.pl). The sequences assembly was done with the CodonCode Aligner software (CodonCode Corporation, Centerville, MA, USA).

Sequence polymorphisms (SNPs or insertions/deletions) identified in the CG sequences between the parental genotypes of RILs ("Maresi" vs. CamB or "Georgie" vs. "Harmal") were subsequently genotyped within the adequate RIL population using one of the following methods: CAPS (cleaved amplified polymorphisms) or dCAPS (derived CAPS) for SNPs, and PCR for indels. In case of CAPS, the same primer pairs were used for the PCR amplification, then PCR products were cleaved with the appropriate restriction enzymes and visualized by 2.0% agarose gel electrophoresis. In the dCAPS, mismatch primers were designed with the dCAPS Finder 2.0. software (http:// helix.wustl.edu/dcaps/dcaps.html) for nested-PCR amplification step, which followed the standard PCR amplification, then PCR products were cleaved with the appropriate restriction enzymes and visualized by 4.0% agarose gel electrophoresis. The indels were genotyped using the PCR amplifications with a newly designed primer pairs, followed by 4.0% agarose gel electrophoresis of the PCR products. Primer pairs designed for the methods of genotyping polymorphisms are given in Supplementary Material S2.

The functional CG linkage analysis was performed on the individual maps of the adequate mapping populations (MCam or GH), and was followed by construction of a new consensus, function map with the use all three maps and the software JoinMap 3.0 (Van Ooijen and Voorrips, 2001), according the procedure described previously (Mikołajczak et al., 2016). Marker order within each linkage group was check for the accordance with the other consensus barley maps (Varshney et al., 2007; Close et al., 2009). Final chromosome maps were drawn with the MapChart software (Voorrips, 2002).

#### QTL Analysis

QTL analysis was performed with the MapQTL 5.0 software (Van Ooijen, 2004) for 40 physiological/biochemical traits measured under both water regimes (optimal water supply and drought) and for their DSI. QTLs were first mapped by interval mapping (IM), followed by multiple-QTL model (MQM), using the marker closest to the peak at each putative QTL as a cofactor. After performing a genome-wide permutation test with 1,000 iterations (Churchill and Doerge, 1994), a LOD (logarithm of the odds) thresholds from 2.8 to 3.5 (depending on the trait) were used to establish the presence of significant QTLs (p < 0.05). Confidence intervals for the QTLs were estimated based on two-LOD support interval, by taking two positions around the peak of the LOD profile, which had LOD-values by 2.0 lower than the maximum. The percentage of phenotypic variation (R 2 ) explained by each QTL was calculated and a QTL was considered as major when it explained >10% of the phenotypic variation. The additive genetic effects were also calculated for the QTLs, and positive value indicated that the "Maresi" allele increased the trait value, whereas negative value indicated the decrease in the trait value caused by the "Maresi" allele. The QTLs which showed overlapping confidence intervals were clustered into hotspots.

#### Identification of Positional Candidate Genes Within QTL Confidence Intervals

In order to identify the positional candidate genes (CG; i.e., closely linked genes localized within QTL regions; Pflieger et al., 2001), genetic and physical map integration was performed based on the positions of markers defining the boundaries of QTL confidence intervals (or the closest to them) in the genome. Nucleotide sequences of markers were mapped to the barley reference genomic sequence deposited in Ensembl Plants database (version 082214v1; http://plants.ensembl.org/index. html) using BLAST tool in order to project QTLs on the physical map (maximum E-value = 1E-100, minimum 95% identity of the sequence). Gene models found within the physical QTL intervals were retrieved using BioMart tool and grouped into the same four gene categories as the physiological/biochemical traits and corresponding QTL regions. This classification was based on the assumption that positional CGs identified within QTL regions related to a given trait category should be involved in the biological processes relevant to this trait category under the analyzed water regime.

#### GO Enrichment Analysis

With the aim to identify GO terms (Biological Processes; BPs) associated with positional CGs and to determine the overrepresentation of a given GO term in an analyzed set of genes compared to the genome-wide background frequency, the GO enrichment analysis was performed using the PLAZA Monocots database version 3.0 (http://bioinformatics.psb.ugent. be/plaza/versions/plaza\_v3\_monocots; Proost et al., 2015). The significance of over-representation was determined using the hyper geometric distribution followed by the Bonferroni correction for multiple testing (corrected p ≤ 0.05). Furthermore, PLAZA was used to barley/Arabidopsis cross-species analysis in order to identify putative barley orthologs in the Arabidopsis thaliana genome. It was motivated by the extensive annotation features available for this model species.

# RESULTS

#### Phenotypic Variation and Correlations Among Traits

Descriptive statistics for physiological/biochemical traits evaluated in the MCam RIL population and parental genotypes is shown in Supplementary Material S3. The mean values of most of the measured parameters under both conditions: optimal water supply (C) and drought stress (D) as well as their stress indices (DSI), significantly differentiated both parental genotypes and varied among the RILs. Almost all traits were distributed normally within RIL population, and lines which exceeded the parental range of variation for different traits were observed. Significant transgressive segregation suggests a broad genetic diversity of parental genotypes and the polygenic inheritance of the investigated traits. Drought stress conditions caused a significant decrease of mean values for the majority of physiological parameters describing plant water status, specific energy fluxes, quantum efficiency ratios, photochemical activity of PSII and efficiency of gas exchange. On the other hand, an increase in proline and ethylene content, as well as in the activity of antioxidant enzymes and fat-soluble antioxidants, was observed. Analyzed genotypes varied in terms of carbohydrates content, both under control and drought conditions. Drought stress increased the accumulation of glucose, fructose and sucrose, whereas maltose and raffinose contents were diminished.

Under control conditions, strong significant correlations (|r| ≥ 0.5, p < 0.05) were exclusively observed among photosynthesis-related traits and parameters describing plant water status (**Figure 1**). The strongest positive correlations were revealed between: net photosynthesis (Pn) and transpiration (E); the photochemical quantum yield of PSII (8PSII) and the photochemical quenching (qP); the performance index per absorption (PIabs) and quantum efficiency ratios (φpo, ψo, φeo). Furthermore, three groups of parameters were found with significant internal correlations: specific energy fluxes (ABS, TR0, ET0, DI0) per excited leaf cross-section (CS), specific energy fluxes per reaction centre (RC), and quantum efficiency ratios (φpo, ψo, φeo). The strongest negative correlations were observed between: transpiration (E) and water use efficiency (WUE); and between ABS/RC, TR0RC, DI0/RC and the efficiency of excitation energy capture by open PSII reaction centres (Fv'/Fm'), PIabs, φpo, and φeo.

Interestingly, the increasing number and magnitude of significant correlations, both among physiological parameters and between physiological and biochemical traits, was revealed under drought stress conditions (**Figure 1**). Photosynthetic efficiency parameters and plant water status traits were positively correlated in the majority of cases. The strongest relationships, not identified under control conditions, were observed between: the RWC and 8PSII, Pn, E, Fv ′ /Fm ′ ; as well as 8PSII and the number of active RCs per CS (RC/CS), specific energy fluxes (ABS, TR0, ET0,) per CS, PIabs, and (1-B)av. Only three parameters, i.e., electrolyte leakage (EL), ABS/RC, and DI0/RC were strongly and negatively correlated with all of the other physiological traits. This indicates drought-induced increase of the permeability of cell membranes, inactivation of a part of RC pool and increase of the antenna size. Under drought stress, in contrast to optimal water supply, we observed numerous and mainly negative correlations between the physiological and biochemical parameters (in particular, for SOD, POX, γ- and α-tocopherols activities). Taken together, these results revealed the changes in cellular metabolism during water deficit which activated the processes enhancing barley adaptation capacity to unfavorable conditions and prevented from a negative impact of drought stress on photosynthesis efficiency.

#### Quantitative Trait Loci Identification

QTL analysis was performed using our previously published high-density SNP and SSR-based consensus genetic map (Mikołajczak et al., 2016) which was enriched with 64 functional CGs including differentially expressed genes (DEGs) derived from the transcriptome analysis of parental genotypes under drought stress and barley orthologs of drought tolerance-related genes described in model species (Supplementary Material S1). The newly constructed barley function map consists of 875 loci and spans 941.86 cM with an average resolution of 1.1 cM. Based on comparison with physical map we estimated that it covers about 95% of the barley genome. This provides an excellent framework for QTL identification. In order to extensively evaluate the response of the analyzed genotypes to drought, both, the direct values of the 40 physiological and biochemical parameters measured under control conditions and drought stress, as well as the relative values (DSI) for each trait, were used in QTL mapping.

A total of 64 QTLs for 25 drought tolerance-related traits that describe plant water status, photosynthetic efficiency, osmoprotectant, and hormone content, as well as activity and accumulation of antioxidants were identified among all of the chromosomes, except 4H (**Figure 2**). The number of detected QTLs varied from 1 to 4, depending on the trait and the experiment variant (C/D/DSI). The highest number of QTLs were positioned on chromosomes 2H and 3H (18 and 17, respectively), and 5H (12). Eight regions for QTLs were mapped on chromosomes 6H and 7H, whereas chromosome 1H contained a single QTL. The maximum LOD scores estimated in the QTL confidence intervals ranged from 3.0 to 20.76, while the phenotypic variation explained by an individual QTL varied from 6.7 to 87.5%. Chromosome 3H contained regions that accounted for most of the observed phenotypic variation in the investigated parameters (14 QTLs with at least 30% of R 2 ). Favorable alleles which had positive effects on the variation in analyzed traits came from both parental genotypes, "Maresi" and CamB.

Under control conditions, the QTL analysis revealed 32 chromosomal regions for 17 traits (**Table 2**). Most (22) QTLs were found to determine photosynthetic efficiency. The others were related to the plant water status traits (4), raffinose content (2), maltose content (2), sucrose content (1), and catalase activity (1). About a half of the identified QTLs under control conditions explained a large proportion of the observed phenotypic variation (R 2 ≥ 30%). Chromosomal regions with the highest R 2 values (exceeding 60%) were mapped for RWC, TR0/CS, DI0/CS and maltose content. Furthermore, two QTLs located on chromosome 7H, which explained more than 85% of the observed phenotypic variation in maltose content, were characterized by the highest LOD scores (>20.0) among all the QTLs revealed in the analyzed RIL population. Three functional CGs, co-segregating with the maximum LOD scores

stress (red), control conditions (green), stress indices (blue), and QTL hotspots (patterned).

TABLE 2 | Summary of QTLs detected for 17 physiological and biochemical traits in the MCam population under control conditions.

Functional candidate genes that co-segregated with the maximum LOD scores within confidence intervals of mapped QTLs are underlined; LOD ≥ 10 and R2 ≥ 30 are in bold. R 2 , percentage of phenotypic variation explained by individual QTL.

within confidence intervals of the QTLs mapped under control conditions, were found. X71362.1 gene, encoding a member of the dehydrin family, was associated with the QTL for ET0/RC on chromosome 5H. On the same chromosome, HvABI5 gene, involved in plant response to abscisic acid, was located in the QTL confidence interval for TR0/RC. Finally, the QTL for maltose content mapped on chromosome 7H co-segregated with barley ortholog of ADC1 gene which encodes arginine decarboxylase.

The QTL analysis performed for the same physiological/biochemical parameters of RILs grown under drought stress conditions showed the presence of 19 QTL regions for 11 traits (**Table 3**). All detected QTLs were declared as the major ones, as they explained more than 10% of the phenotypic variation of the particular traits. The highest number of QTLs was identified for photosynthetic efficiency parameters (7), as well as osmoprotectant content (6). The remaining QTLs were found to determine plant water status (4), ethylene content (1), and α-tocopherol activity (1). The highest R 2 values and LOD scores described the chromosomal regions involved in the proline and maltose accumulation. Drought-induced proline level was found to be controlled by three major QTLs mapped on chromosomes 2H and 3H, which explained between 73.6 and

Functional candidate genes that co-segregated with the maximum LOD scores within confidence intervals of mapped QTLs are underlined; LOD ≥ 10 and R2 ≥ 30 are in bold. R2 , percentage of phenotypic variation explained by individual QTL.

77.5% of the observed phenotypic variation, but an additional locus with a smaller effect was identified on chromosome 2H. In addition, two QTLs located on chromosome 3H, involved in maltose content regulation, were characterized by LOD scores above 18.0 and R 2 values that exceeded 85%. Two differently expressed functional CGs were found to be associated with the chromosomal regions for fast kinetics of chlorophyll a fluorescence parameters, located on chromosome 2H. Gene AK356764.1 that encodes transketolase co-localized with the QTL for TR0/RC, whereas gene AK353596.1 encoding a member of the peroxidase family showed coincidence with the LOD peak region for ET0/RC.

All 13 QTLs mapped for stress indices (DSI) of 11 parameters were classified as the major ones (**Table 4**). Most of them were detected for photosynthetic efficiency and osmoprotectant content (6 and 4, respectively), similarly to the drought stress conditions. The other QTLs were identified for leaf water content (1) and tocochromanols activity (2). The highest proportion of observed phenotypic variation, ranging from 60.8 to 82.1%, was explained by two QTLs positioned on chromosome 3H: for sucrose and for γ-tocotrienol stress indices. Two functional CGs from the transcriptome analysis of "Maresi" and CamB co-segregated with maximum LOD scores within the QTL confidence intervals, identified on chromosome 5H. Dehydrin gene (X71362.1) co-localized with the QTLs for leaf water content stress index, while AK356764.1 gene that encodes cysteine protease was associated with the QTL for TR0/RC stress index.

#### Co-location of QTLs for Different Traits

The positioned QTLs were not distributed evenly in the barley genome, and they clearly tend to be clustered in the particular chromosome regions. Based on their overlapping confidence intervals, 11 hotspots were identified that contained together more than 60% of mapped QTLs for different traits (**Table 5**, **Figure 2**). The highest number of hotspots (4), that combined the QTLs from all investigated trait categories analyzed in all experiment variants (C/D/DSI), was located on chromosome 2H, followed by three hotspots which contained 2–6 overlapping QTL confidence intervals detected on chromosome 3H. Two hotspots were localized on chromosome 5H, each of them for two physiological parameters. Then, a single genomic region on chromosome 6H enclosed 7 QTLs for photosynthetic efficiency under control water conditions. The last hotspot, that included 2 overlapping QTLs for DSI of fast kinetics of chlorophyll a fluorescence parameters, was identified on chromosome 7H.

#### QTL Projection on Barley Genome

One of the main goals of this study was to identify the positional candidate genes (CGs) within QTL confidence intervals, followed TABLE 4 | Summary of QTLs for stress indices (DSI) for 11 physiological and biochemical traits in the MCam population.

Functional candidate genes that co-segregated with the maximum LOD scores within confidence intervals of mapped QTLs are underlined; LOD ≥ 10 and R2 ≥ 30 are in bold. R2 , percentage of phenotypic variation explained by individual QTL.

by the exploration of their putative functions related to drought tolerance. We used the availability of the barley reference genome sequence to search for positional CGs in the chromosomal regions corresponding to the confidence intervals of the QTLs determined on the consensus map. The integration of the genetic and physical map led to the identification of overall 3,198 positional candidate genes that could be responsible for the variation in the analyzed traits (Supplementary Material S4; **Figure 3A**). Among them, there were the CGs (1,101) underlying the drought response-specific QTLs, i.e., those QTLs which were detected under drought stress conditions and for the stress indices. The number of genes included in a particular QTL region for all analyzed traits varied between 1 and 550, and was correlated with the confidence interval size, whereas the total ratio between physical and genetic distances was 2.11 Mbp/cM.

Comparative analysis of the number of the positional CGs classified into four trait categories, as well as three analyzed experiment variants (C/D/DSI) was shown on **Figures 3B–E**. When considering plant water status and antioxidant activity, 699 and 172 genes, respectively, were specific to the analyzed experiment variants. Altogether, 2,004 positional CGs were identified for photosynthetic efficiency, and approximately 30% of CGs involved in this process under drought stress (192 out of 665) were also detected in the control conditions. Even larger overlap was observed among 933 genes corresponding to the osmoprotectant and hormone content. Approximately 30% of positional CGs identified under drought stress (88 out of 289) were also detected under control conditions, while as much as 66% (192) were common with stress indices of the analyzed traits. Only 9 genes were identified exclusively under drought stress within this trait category.

# QTL Candidate Genes Prioritization

Due to the limited resolution of QTL mapping which resulted in a high number of identified positional candidate genes, we applied Gene Ontology (GO) enrichment approach to: (1) predict the functional ontologies associated with positional CGs, (2) to determine the over-representation of certain GO terms in the analyzed gene sets compared to the genome-wide background frequency, (3) to reduce the number of CGs to potentially casual genes responsible for the variation in the considered traits.

As for the gene sets detected for the analyzed trait categories under optimal water supply, the functional annotation showed their participation in Biological Processes (BP) clearly relevant for the proper shaping of plant growth and development at the morphological, anatomical and physiological levels. They included e.g., post-embryonic organ development, carbohydrate mediated signaling, pentose metabolic process, protein phosphorylation, stomatal movement, defense response, as well as biosynthesis and transport of various compounds (Supplementary Material S5). In order to gain an insight into molecular mechanisms of drought stress response in barley, the special attention was paid to 1,101 positional CGs underlying drought response-specific QTLs. The GO enrichment analysis resulted in narrowing down this number to 143 prime candidates involved in the significantly overrepresented BPs which were closely related to the exposition of the analyzed genotypes to drought stress (**Figure 4**). Some of them are described as unknown in existing annotations and assessment of their potential engagement in drought response in barley requires further exploration of their biological functions.

D, drought stress; C, control conditions; DSI, stress index.

#### CGs for Plant Water Status

The GO enrichment in the set of positional CGs identified for plant water status under drought stress, revealed 17% of them to be involved in the oxidation-reduction process (Supplementary Material S6). A detailed analysis showed that two barley genes, MLOC\_5505 and MLOC\_69302, encoding enzymes belonging to the dehydrogenase and peroxidase families, respectively, were identified in qHS2.1 hotspot which included the overlapping QTLs for leaf water content and water loss under drought, as well as glucose content stress index and catalase activity under control water conditions. The other genes of this GO category are: MLOC\_6299 (an ortholog of Arabidopsis AT5G65110, involved in biosynthesis of long chain fatty acid), MLOC\_4579 (flavonol synthase), and MLOC\_71498 (an ortholog of Arabidopsis gene AT4G25650, encoding translocation channel at the inner envelope membrane of chloroplasts). All these genes were detected in the QTL confidence interval for water loss, on chromosome 7H. The remaining highly enriched

group of genes were annotated as intermembrane transport (MLOC\_52403 encoding vesicle-associated membrane protein), homoserine metabolic process and "de novo" L-methionine biosynthetic process (MLOC\_71910, cystathionine gamma), cyanidin 3-O-glucoside metabolic process and beta-glucoside metabolic process (MLOC\_36976, UDP-glucosyl transferase), and regulation of brassinosteroid mediated signaling pathway (MLOC\_72613, L-type lectin-domain containing receptor kinase III). Among the CGs for water content stress index, the most enriched Biological Process was carboxylic acid biosynthesis. A detailed analysis of the genes classified into this BP category revealed two barley genes: MLOC\_18300 (a member of NCEDrelated gene family) encoding the key enzyme involved in ABA biosynthesis, and MLOC\_53947 for putative plastid pyruvate dehydrogenase. The others significantly over-represented GO terms were: proline biosynthesis process, ornithine catabolic process, arginine catabolic process to glutamate (MLOC\_35821, ornithine aminotransferase), regulation of ATPase activity, and response to heat (MLOC\_37449, protein DnaJ).

#### CGs for Photosynthetic Efficiency

The functional annotation analysis of the positional CGs identified for photosynthetic efficiency under drought stress showed that the majority of genes belong to the GO categories of protein folding, response to heat, response to high light intensity, and response to hydrogen peroxide (Supplementary Material S7). They include genes located mostly on chromosomes 2H and 7H that encoded heat shock proteins (MLOC\_568, MLOC\_31567, MLOC\_41281) and chaperonins (MLOC\_51927 and MLOC\_54083). Some of them are involved in folding of RuBisCO, the major enzyme at the first step of carbon fixation in the Calvin cycle. A particularly interesting GO category comprised three genes mapped on chromosomes 2H and 5H (MLOC\_15256, MLOC\_19075, MLOC\_38928), encoding ATP sulfurylases involved in the sulfate assimilation pathway. The other highly enriched Biological Process among CGs for photosynthetic efficiency in response to drought was detection of abiotic stimulus. Two genes, one encoding phytochrome A (MLOC\_824) and the other involved in the chromatin remodeling in the response to environmental cues (MLOC\_60235) were found within this category. The most enriched GO terms of the CGs for photosynthetic efficiency DSI included three groups. The first group consisted of the over-represented BPs related to DNA dealkylation and DNA demethylation (MLOC\_53029 and MLOC\_11707, encoding DNA repair proteins). Another highly represented GO category was lipid transport. Here, MLOC\_51456 encoding a nonspecific lipid-transfer protein was found. The BPs belonging

FIGURE 4 | Biological processes significantly over-represented in the gene sets for the analyzed trait categories: under drought stress (D) and for stress indices (DSI).

to the last group were related to short-day photoperiodism that included MLOC\_6879 (an ortholog of Arabidopsis gene AT3G24440, VERNALIZATION INSENSITIVE 3-LIKE 1) and MLOC\_19228 (an ortholog of Arabidopsis gene AT5G46210, CULLIN4). Red or far-red light signaling pathway was another overrepresented BP from this group, including MLOC\_20045, encoding nucleoside diphosphate kinase and MLOC\_53845, an ortholog of Arabidopsis gene PFT1, encoding phytochrome and flowering time regulatory protein).

#### CGs for Osmoprotectant and Hormone Content

The majority of drought-specific CGs for the osmoprotectant and hormone content category were common for drought stress conditions and stress indices (DSI). According to the GO enrichment, they are involved in tripeptide and modified amino acid transport, glutathione transport (MLOC\_174, a member of Crt-like transporters), glycerol biosynthetic process (MLOC\_40292, glycerol-3-phosphatase 1), and phytochromobilin metabolic process (MLOC\_70465, phytochromobilin oxidoreductase; Supplementary Material S8). These genes encode proteins located in chloroplasts and they were mapped in qHS3.1 hotspot which included the overlapping QTLs for sucrose content stress index, as well as for proline and maltose content under drought stress. Among CGs for this category, GO enrichment analysis led to the identification of two barley genes exclusive for drought stress, which were positioned in qHS3.2 hotspot. MLOC\_65646 is a putative ortholog of Arabidopsis gene AT5G20990, encoding molybdopterin biosynthesis protein, which is involved in molybdenum cofactor biosynthesis, whereas MLOC\_22343, encoding molybdate transporter showed enrichment in BP categories related to vacuolar transmembrane transport, molybdate ion export from vacuole, and sulfate transport. The most represented category of BPs among the CGs for osmoprotectant and hormone content stress indices was auxin mediated signaling pathway. One of the genes involved in this process was MLOC\_58506, encoding a member of the small auxin up-regulated RNA family. It was mapped on chromosome 5H in QTL confidence interval for maltose content stress index. Other highly enriched GO categories were the same as for the plant water status under drought stress: cyanidin 3-O-glucoside metabolic process, beta-glucoside metabolic process (MLOC\_36976, UDP-glucosyl transferase), and regulation of brassinosteroid mediated pathway (MLOC\_72613, L-type lectin-domain containing receptor kinase III).

#### CGs for Activity and Accumulation of Antioxidants

For the activity and accumulation of antioxidants trait category, a single QTL for α-tocopherol content under water-deficiency conditions was identified on chromosome 6H. The GO enrichment analysis among CGs for α-tocopherol led to the identification of over-represented BP categories, such as: regulation of isopentenyl diphosphate biosynthetic process, regulation of phospholipid biosynthetic process, nitric oxide metabolic process, and indole glucosinolate metabolic process (Supplementary Material S9). In-depth analysis showed MLOC\_63263, encoding a chloroplastic/mitochondrial NO-associated protein, and MLOC\_57100, encoding Omethyltransferase which is involved in flavonol biosynthesis. The analysis of the GO annotations of the positional CGs for antioxidant activity DSI revealed their involvement in: flavonol metabolic process, flavone metabolic process (MLOC\_73233, O-methyltransferase), and urea transport (MLOC\_58872, an ortholog of Arabidopsis gene AT4G01470, encoding tonoplast intrinsic protein functions as water and urea channel). Furthermore, the both selected genes for antioxidant activity stress indices were identified in qHS3.2 hotspot which overlapped 6 QTLs, among them for γ-tocotrienol and sucrose stress indices, as well as proline content under drought stress.

# DISCUSSION

The identification of genes underlying particular QTLs, followed by the elucidation of their molecular functions, is one of the paramount challenges in the modern plant genetics, as more than 2500 studies on QTL analysis in crop plants have been published so far (Kumar et al., 2017). Most studies have focused on exploiting the highly saturated molecular function maps to colocalize QTL regions with the genetic markers, based on withingene polymorphisms (Diab et al., 2008; Sehgal et al., 2012; Li W. T. et al., 2013; de Miguel et al., 2014; Mikołajczak et al., 2016). In a standard strategy, after fine mapping of a QTL, a positional cloning of the QTL has been proposed to identify a gene or genes underlying a complex trait (Kumar et al., 2017). Although, there are examples of successful positional cloning of QTLs (Salvi and Tuberosa, 2007; Collins et al., 2008; Mir et al., 2012), this procedure remains difficult and cumbersome, especially for large genomes. This could explain a small number of cloned QTLs in barley with only a few examples related to abiotic stress tolerance, such as boron (Sutton et al., 2007), and aluminum toxicity (Furukawa et al., 2007) or freezing tolerances (Francia et al., 2007).

Recently, the availability of complete genome sequences for a growing number of plant species gave the opportunity to gain an insight into the physical gene space, resulting in the straightforward searching for all possible positional candidate genes associated with QTL confidence intervals (Monclus et al., 2012; Bargsten et al., 2014; Correa et al., 2014). Therefore, an inherent bottleneck is an effective selection of the most promising CGs within QTL regions, typically including tens to hundreds of genes, most of them unrelated to the trait of interest. Various methods to prioritize CGs within the QTL regions for the target trait can be combined, including the approach based on CGs overrepresentation in the biological processes, the analysis of their differential expression at transcriptome and/or proteome levels or functional annotation of CGs based on their orthologs in related or model species (Monclus et al., 2012; Bargsten et al., 2014; Kumar et al., 2017).

Our study presents the comprehensive multistep approach for elucidation of the genetic basis of drought response mechanisms in barley. Firstly, the identification of precise chromosomal regions underlying physiological and biochemical indicators of the drought tolerance was accomplished by merging the extensive analysis of these parameters under drought stress and classical quantitative genetic approach, performed with the use of strict threshold criteria and a high-density function map. Next, the anchoring of drought response-specific QTLs to the reference barley genome sequence was performed, after the thorough evaluation of QTL data. Finally, a huge number of identified positional CGs was reduced to potentially causative genes underlying the analyzed traits, using a Gene Ontology enrichment approach, followed by additional prioritization tools. Below, the most important results of the particular steps of the study are discussed.

#### Reliability and Accuracy of QTL Mapping

QTL analysis in this study identified a total of 64 QTLs for 25 out of 40 physiological and biochemical traits measured under the optimal and drought stress conditions (absolute values), as well as for their DSI (relative values). The identified QTLs showed LOD scores ranging from 3.0 to 20.8, and explained 6.7 to 87.5% of observed variation (R 2 ). The applied LOD threshold of 3.0 to declare a presence of the QTL ensured the detection of reliable and significant QTLs involved in the formation of analyzed traits and the exclusion of small effect, tentative QTLs (Salvi and Tuberosa, 2015; Kumar et al., 2017). Using the same criterion (LOD > 3.0), Liu et al. (2017) identified 11 significant QTLs for 6 traits, out of a total of 15 stomatal and photosynthetic traits related to salinity tolerance in barley, while with a lower LOD score (2.5 < LOD < 3.0), the authors added 11 tentative QTLs for 11 traits to this list. Similarly, Wójcik-Jagła et al. (2013), in two barley mapping populations identified in total 33 QTLs (LOD > 2.5) for physiological parameters measured under shorttime drought, but only 11 QTLs were detected with a LOD score >3.0. In the present study, with the aim to search for positional candidate genes associated with QTL confidence intervals, we decided to consider only significant QTLs with a high LOD score.

 Assuming that the proportion of phenotype variation explained by a major QTL should be above 10% (Kumar et al., 2017), as much as 95% of the QTLs identified in our study can be considered as the major ones, including all 32 QTLs detected under drought stress and for DSI. Similarly, the 11 significant QTLs detected by Liu et al. (2017), mentioned above, can be classified as the major QTLs, as each of them explained at least 11.2% of variation. On the contrary, in the work of Wójcik-Jagła et al. (2013) only 4 out of 33 QTLs detected in two mapping populations could be considered as the major ones using the above criterion. Moreover, the majority of the QTLs (92.2%) of our study were precisely positioned on the linkage map within confidence intervals below 6.0 cM. It is not a common situation as usually the linkage analysis locates QTLs within the confidence intervals of 10–30 cM (Chander et al., 2008; Li et al., 2008; Capelle et al., 2010; Chen et al., 2010; Kumar et al., 2017). However, in accordance with our results, Bertholdsson et al. (2015) showed a group of 5 QTLs including 3 major QTLs for quantum yield of electron transport of PSII under waterlogging in barley, and all of them but one, were mapped within confidence intervals below 5.0 cM.

# Co-locations of Present and Known QTL Regions

Over the last two decades, QTL for a wide range of traits related to the drought tolerance, including physiological/biochemical characteristics, have been mapped to all seven barley chromosomes (Mir et al., 2012). A precise comparison among these results is not possible owing to the differences in used plant materials and maps, various traits analyzed and diverse methodology applied, however, some interesting observation can be made in regard to the results of the present study. The QTLs identified in our study were mapped to all barley chromosomes, except 4H, with the highest number of QTLs located on chromosomes 2H, 3H and 5H. Based on the overlapping confidence intervals, 11 hotspots were identified that contained together over 60% of mapped QTLs for all the traits, most on 2H (4) and 3H (3). In the previous studies with the same mapping population, aimed at the detection of QTLs underlying yield-related agronomic traits under drought and control conditions (Mikołajczak et al., 2016, 2017; Ogrodowicz et al., 2017), the important role of the regions on chromosomes 2H and 3H was also reported. Moreover, the most significant QTL hotspot on 2H, clustering QTLs for the length of main stem and yield traits (Mikołajczak et al., 2017), as well as heading date (Ogrodowicz et al., 2017), co-localized with a hotspot region (qHS2.1) of the present study. The other region of overlapping QTL hotspots between these studies was found on chromosome 5H. Similarly, Mora et al. (2016) revealed the highest number of QTLs for drought-related morphological and physiological traits on chromosomes 2H and 3H, using a distinct barley gene-pool and environmental conditions. Recently, meta-QTL analysis approach has been developed which is aimed in the integration of data from multiple QTLs studies and has a greater statistical power for the detection of, so called, meta-QTLs (MQTL), and more precise estimation of their genetic effects (Wu and Hu, 2012). Zhang et al. (2017) performed a meta-QTL analysis of drought tolerance in barley using 72 major QTLs described in several studies, and most of QTLs were located on chromosomes 2H, 3H, 5H, and 7H. As a result, MQTLs, integrating QTLs for barley drought tolerance, have been positioned on chromosomes, with some particularly important regions, common to drought and salinity tolerance on 2H (2), 3H (1), and 5H (1).

In the present study no overlapping QTLs identified for a given trait under different water regimes were found. A similar situation was observed by Wang et al. (2014) in QTL mapping for yield-related traits in barley in six diverse environments. Mora et al. (2016) also showed that 80% of the identified QTLs for analyzed traits were specific for a particular environment. According to this observation, we used only the QTLs detected under drought stress conditions and QTLs for adequate DSI, to be subjected for the co-location analysis, projection to barley genome sequence and positional CGs identification.

For the plant water status related traits, a total of 5 QTLs were detected under drought stress conditions and for adequate DSI in our study. Among them, two neighboring regions from chromosome 5H have been pointed in previous studies: QdsiWC.5H\_3 for the water content DSI and QdWL.5H\_3 for water loss under drought stress co-located with QTL region for RWC determined by Teulat et al. (2001a, 2003), in the close vicinity of dhn1 and dhn9 genes encoding dehydrins. This colocation is further confirmed as QdsiWC.5H\_3 (R 2 = 23.6%) was mapped into the position of X7136.1 gene, a functional CG which encodes a dehydrin protein. This gene was chosen for mapping based on the transcriptome analysis of parental genotypes of MCam population under drought stress.

For the photosynthetic efficiency category 13 QTLs were detected under drought stress conditions and for adequate DSI. The analysis revealed that two QTL regions from the distal part of a long arm of chromosome 2H: QdABS/RC.2H.2 for PSII light absorption flux per RC and QdET0/RC.2H (from hotspot qHS2.4) for electron transport flux in PSII per RC co-localize with two hotspots for chlorophyll a fluorescence parameters detected by Guo et al. (2008). These two chromosomal regions were also identified by Wójcik-Jagła et al. (2013) with QTLs for photochemical quantum yield of PSII (8PSII) DSI, and for photochemical quenching (qP). The region of QdET0/RC.2H seems to be especially significant as it explained the highest proportion of the variation (20.5%) and was mapped in the position of AK353596.1 gene (MLOC\_65477), a functional CG used for map construction. This gene encodes for an enzyme from the peroxidase class engaged in reactive oxygen species (ROS) detoxification. Different ROS are excessively formatted under drought stress, but they are also produced during the light phase of photosynthesis where about 10% of the pool of electrons is transferred to oxygen forming superoxide anion radical (Foyer and Noctor, 2000). The AK353596.1 gene was selected for mapping based on its differential expression between parental genotypes of MCam under drought (Janiak et al., 2018). Thus, the relationship between this gene and the identified QTL seems to be well-supported. Additionally, two other QTLs related to photosynthetic efficiency of PSII showed the co-locations with QTLs detected previously. The region of Qd(1-B)av.7H\_2 from 7H co-located with the QTL for initial fluorescence (F0) from the study of Guo et al. (2008) and QdsiDI0/RC.5H\_2 from 5H overlapped with the chromosomal region of four QTLs for DSI for PSII photosynthetic efficiency and water content (Wójcik-Jagła et al., 2013).

In the osmoprotectant and hormone content category, altogether 11 drought response-specific QTLs were detected. It should be underscored that our study for the first time identified the QTLs for glucose, maltose, sucrose and ethylene content under drought stress conditions in barley. Among them, the QTL from the short arm of 2H, QdsiGlu.2H, for glucose content DSI, overlapped with three other QTLs for WC, WL, and TR0/RC (hotspot qHS2.1, around SNP marker 5880- 2547) and was mapped near the position of another functional CG, AK356764 (MLOC\_21709, an ortholog of Arabidopsis AT3G60750), encoding transketolase. This enzyme is involved in the Calvin-Benson cycle during the light-independent phase of photosynthesis, and according to GO terms, might be also involved in salinity tolerance, gluconeogenesis and water transport. The QdsiGlu.2H region showed the co-location with QTL for WSC concentration at full turgor under drought stress (WSC100) designated by Diab et al. (2004). The most interesting result was obtained in the case of six, partially overlapping QTLs for proline, maltose and sucrose content, located in the distal part of a long arm of chromosome 3H. These QTLs explained the highest amount of phenotype variation (54.6–85.9%), were divided among three neighboring hotspots (qHS3.1, qHS3.2, qHS3.3), and coincided with the QTL regions for proline content under drought and salinity stresses identified by Sayed et al. (2012) and Fan et al. (2015), respectively.

In the activity and accumulation of antioxidants category, 3 QTLs were detected under drought stress conditions and for adequate DSI. Here again, for the first time our study determined the QTLs for α-tocopherol (QdATf.6H on 6H) and γ-tocotrienol (QdsiGTt.3H\_2) content under drought stress in barley. The region QdsiGTt.3H\_2 for γ-tocotrienol content DSI explained 82.1% of phenotype variation for this trait and was included into the hotspot qHS3.2 together with QTLs for proline and sucrose content (co-locations presented above). A third QTL region, QdsiGTf.7H\_3 for γ-tocopherol content DSI was mapped on chromosome 7H near SNP marker 6628-1302 and showed the colocation with the QTL for γ-tocopherol content under drought stress detected recently by Templer et al. (2017).

#### Prioritization of Positional CGs Related to Drought Response-Specific QTLs

Although the sequence assembly of the barley genome has been accessible since 2012 (Mayer et al., 2012), only a few studies demonstrated the genetic and physical map integration resulting in numerous putative CG identification (Fan et al., 2015; Piasecka et al., 2017; Templer et al., 2017; Zhang et al., 2017). In the present work, the QTL projection on the genome was extended by the GO enrichment analysis. This allowed an efficient selection of CGs on the basis of their relevance to biological processes related to analyzed trait categories. To the best of our knowledge, it was the first attempt to apply the Gene Ontology-based prioritization to barley QTL analysis. The resulting reduction in the number of genes was almost eight-fold: out of the 1,101 CGs corresponding to the drought response-specific QTLs, we selected 143 significant candidates predicted to be responsible for the variation in the considered traits. Similar findings have been reported by Bargsten et al. (2014) who prioritized rice QTL data leading to the ten-fold reduction in the number of candidate genes. The GO enrichment was also successfully applied to select CGs participating in the genetic determination of the grapevine cluster architecture (Correa et al., 2014), as well as the productivity, growth, and water-use efficiency in poplar (Monclus et al., 2012).

Unquestionably, the thorough examination of CGs revealed in the present study is necessary to confirm their exact role in drought stress response in barley. However, we want to point out that our findings have been supported by results of the global gene expression profiling of parental genotypes that had been used for QTL mapping in our study (Janiak et al., 2018). The differentiated expression due to drought stress was confirmed for 34 CGs from the set of 143 genes selected based on prioritization approach. In the following paragraphs, we characterize the best candidates for further experimental validation of their role in driving drought response. The genes have been carefully chosen based on one of the following criteria: (1) differentiated expression under drought stress, (2) location within overlapping QTL confidence intervals (hotspots), (3) welldocumented engagement of their orthologs in model species in drought stress response.

The majority of CGs identified within drought responsespecific QTLs for plant water status corresponded to the oxidation-reduction process. It is well-recognized that accumulation of reactive oxygen species (ROS) in the chloroplasts and mitochondria is an early response of plants to decreased water potential inside the cells during adverse stress conditions (Xiong et al., 2002; Choudhury et al., 2013). Many studies have revealed the ROS-mediated retrograde signaling pathway from chloroplasts to nucleus, which results in substantial changes in the expression of the nuclear genes that maintain chloroplasts function and other aspects of plant adaptation to environmental cues (reviewed in Chi et al., 2015). In contrast, an attention has been drawn to the role of ROS as major drivers of cellular oxidative damages (Noctor et al., 2014). Among the barley genes that correlated with the oxidation-reduction process, we highlighted MLOC\_69302 encoding ascorbate peroxidase (APX) as a key player involved in the hydrogen peroxide removal. It was reported that APX lossof-function Arabidopsis mutants accumulated more hydrogen peroxide and were significantly more sensitive to different abiotic stresses than the wild type (Koussevitzky et al., 2008). Another gene identified in our study that may be important for barley drought response, MLOC\_4579, encodes a flavonol synthase. Flavonols are members of secondary metabolites highly accumulated under various stresses which are involved in ROS scavenging leading to the enhancement of oxidative tolerance (Nakabayashi et al., 2014). Therefore, it can be suspected that the identified barley redox-linked genes contribute to regulation of ROS homeostasis under water deprivation.

The carboxylic acid biosynthesis was the second overrepresented process among CGs identified in intervals of QTLs for plant water status. Within this group, MLOC\_18300 encoding a member of 9-cis-epoxycarotenoid dioxygenase family (NCED), which catalyzes an essential regulatory step in ABA biosynthesis, was found (Seo and Koshiba, 2002). Drought stress-stimulated accumulation of ABA modulates root hydraulic conductivity and regulates shoot vs. root growth, as well as alters guard cell ion transport, promoting stomatal closure and reducing water loss by transpiration (Roychoudhury et al., 2013; Sah et al., 2016). Another identified gene, MLOC\_53947, encodes a putative plastid pyruvate dehydrogenase which contributes to transforming pyruvate to Acetyl-CoA, the precursor for fatty acid biosynthesis (Johnston et al., 1997; Mentzen et al., 2008). It was shown that prolonged drought exerts phospholipid bilayer destabilization and increases its permeability leading to the ion leakage (Bajji et al., 2002; Fang and Xiong, 2015). Therefore, an enhanced biosynthesis of fatty acids may be a crucial strategy of protecting plasma membrane integrity against the damage under drought stress, maintaining osmotic homeostasis and cell turgor pressure, as well as controlling transcellular water and ion transport.

Our analysis clearly showed a putative involvement of numerous genes encoding the heat shock proteins (HSP) in the regulation of photosynthetic efficiency under drought stress in barley. Members of different HSP families perform the molecular chaperone function by participation in the proper protein folding and assembly. Furthermore, they enable disposal of non-functional and aggregated stress-labile proteins (Santhanagopalan et al., 2015). It is well-recognized that HSP proteins play a pivotal role in the protection of photosynthetic components against stress injuries. Recent studies have revealed that stress-induced impairment of photosynthesis results in the reduction of carbon fixation and oxygen evolution, as well as disruption of the linear electron flow. One of the main effects of photoinhibition is the exposure of chloroplasts to excess excitation energy. This increases the generation of ROS by the incomplete reduction of molecular oxygen, which in turn induce remarkable damages in chloroplast proteins, lipids, and pigments (Filek et al., 2015; Pospisil, 2016). The components of photosynthetic apparatus particularly sensitive to oxidative stress are: PSII with its oxygen-evolving complex, the carbon assimilation process driven by RuBisCO, and the ATP generating system (Allakhverdiev et al., 2008). Among genes for HSP identified in this work, the most represented group was orthologs of Arabidopsis AT5G51440 encoding HSP20-like chaperone (MLOC\_568 and MLOC\_31567). HSP20 proteins are encoded by nuclear multigenic families and are located in different cellular compartments (Lopes-Caitar et al., 2013). In Arabidopsis, 19 genes encoding HSP20 were identified, whereas 36 and 23 HSP20 genes were described in poplar and rice, respectively (Scharf et al., 2001; Waters et al., 2008; Sarkar et al., 2009). An intriguing overrepresented biological process among CGs for photosynthetic efficiency under drought stress was the sulfate assimilation, represented by MLOC\_19075, an ortholog of Arabidopsis gene AT4G14680 encoding ATP sulfurylase. Many studies indicate the importance of sulfur as the component of a wide range of compounds with fundamental biological functions, including tolerance to various abiotic stresses (Prioretti et al., 2014). ATP sulfurylases catalyze the first step in the sulfate assimilation pathway. They activate sulfate (SO2− 4 ; a metabolically inert form of sulfur taken up by roots), yielding a high-energy adenosine-5 ′ -phosphosulfate (APS) that is reduced to sulfide (S2−) and incorporated into cysteine. In turn, cysteine acts as a donor of highly reactive thiol group (-SH) for numerous S-compounds (Prioretti et al., 2014; Anjum et al., 2015). One of them is glutathione, whose contribution to drought stress tolerance by ROS detoxification and stress-induced signal transduction has been extensively evidenced (Munemasa et al., 2013; Pyngrope et al., 2013; Nahar et al., 2015).

Other relevant biological processes related to the photosynthetic efficiency, over-represented among CGs analyzed in this work, were: short day photoperiodism and red or far red light signaling pathway. Within this group, we highlighted two barley genes, MLOC\_6879 and MLOC\_53845, orthologs of Arabidopsis VERNALIZATION INSENSITIVE 3-LIKE 1 (VIN3-LIKE1) and PHYTOCHROME AND FLOWERING TIME 1 (PFT1), respectively, involved in the control of flowering time. Many studies have reported a limiting effect of drought stress on plant yield, including early arrest of floral development leading to the interruption of plant reproduction. One of the mechanisms to cope with the water deficit adopted by plants is an acceleration of the flowering process to shorten plant life cycle via drought escape. In Arabidopsis, the FLOWERING LOCUS T (FT) is a master regulator of flowering promotion. It activates the expression of meristem identity genes which control the reprogramming of shoot apical meristem to form flowers (Turck et al., 2008; Kazan and Lyons, 2016). FT gene is up-regulated by a transcription factor CONSTANS (CO) and repressed by the product of FLOWERING LOCUS C (FC). Furthermore, PFT1 gene was shown to activate expression of FT through CO- dependent and independent manner. In contrary, FC is epigenetically down-regulated by the vernalization response genes (e.g., VIN3-LIKE1; Inigo et al., 2012; Kazan and Lyons, 2016). Taking these data into account, it could be suggested that drought stress occurring at the seedling stage may trigger the up-regulation of MLOC\_6879 and MLOC\_53845 genes to switch on a drought escape adaptive mechanism by promotion of flowering. This suggestion is supported by a distinctly early flowering of one of the parental genotypes used for QTL mapping, the Syrian line CamB (Ogrodowicz et al., 2017).

We identified several CGs for accumulation of osmoprotectants and hormones within the hotspot regions: qHS2.1, qHS2.3, qHS3.1, and qHS3.2, which suggests their engagement in the complex metabolic networks, probably common to different drought-responsive compounds. As an example, MLOC\_174 [a CRT (chloroquine-resistance transporter)-like transporter 3, involved in modified amino acid and glutathione transport], MLOC\_40292 (glycerol-3 phosphatase), and MLOC\_20354 (gamma-secretase subunit, involved in Notch signaling pathway) were candidate genes for the proline, sucrose, and maltose accumulation. Furthermore, a special attention was given to five genes (MLOC\_58506, MLOC\_58508, MLOC\_58507, MLOC\_65368, MLOC\_65978) encoding SAUR-like auxin-responsive proteins. Many findings suggest participation of auxin in the plant response to stresses, including drought (Du et al., 2013; Krishnamurthy and Rathinasabapathi, 2013; Rahman, 2013; Shi et al., 2014). It is well-documented, that small auxin up-regulated RNAs (SAUR) are the most numerous family among early auxin response genes, found to be tandemly duplicated in plant genomes (Wu et al., 2012). Accordingly, the SAUR genes identified in this work are clustered within single QTL chromosomal region for maltose content on chromosome 5H. Despite of progressive identification of SAUR genes in different plant species, the functions of most of them, especially in the stress response, remained elusive due to their large genetic redundancy (Ren and Gray, 2015). However, recent studies revealed the role of SAUR genes in hypocotyl elongation and leaf senescence in Arabidopsis, bacterial blight resistance in rice, and drought tolerance in wheat (Chae et al., 2012; Kai et al., 2013; Aoki et al., 2016; Guo et al., 2018). Therefore, it can be hypothesized that observed decrease in maltose content under drought stress conditions, may be regulated by auxin-mediated SAUR genes, but this requires further verification.

The most striking result of GO enrichment was the identification of two genes: MLOC\_65646 (encoding molybdopterin biosynthesis protein) and MLOC\_22343 (encoding molybdate transporter) involved in the formation of molybdenum cofactor (Moco), which acts as an active compound at the catalytic center of a wide range of oxidation-reduction enzymes (Mendel, 2013). One of these enzymes is aldehyde oxidase which converts abscisic aldehyde to ABA in the last step of the ABA biosynthesis pathway. It was shown that aldehyde oxidase requires the sulfurylated form of Moco for providing its activity (Mendel and Haensch, 2002). This reaction is catalyzed by Moco sulfurase encoded in Arabidopsis by LOS5/ABI3 gene (Xiong et al., 2001). The study of Li Y. et al. (2013) have revealed that soybean transgenic lines with overexpression of LOS5/ABI3 showed a reduced electrolyte leakage, enhanced accumulation of proline and antioxidant enzymes, as well as significant increase in ABA content and drought tolerance. The identification on the aforementioned barley genes in the QTL hotspot qHS2.3 for ABS/RC and proline content suggests an interplay between the decreasing photosynthetic efficiency and expression of genes engaged in molybdenum cofactor synthesis, which results in ABA-mediated proline accumulation.

Considering CGs related to the activity and accumulation of antioxidants, we distinguished two genes, MLOC\_57100 and MLOC\_73233, orthologs of Arabidopsis AT5G54160 (AtOMT1), encoding O-methyltransferase. Enzymes belonging to the OMT family catalyze transfer of methyl group to the oxygen atom of the wide range of low molecular weight organic compounds, including flavonoids, alkaloids, and phenylpropanoids. Their methylated products are involved in lignin biosynthesis and response to environmental cues (Lam et al., 2007). Furthermore, it has already been shown that drought stress increased OMT transcription in maize, soybean, and grape (Yamaguchi et al., 2010; Liu et al., 2013; Giordano et al., 2016). Since MLOC\_73233 was located in the hotspot qHS3.2 that includes QTLs for γ-tocotrienol, proline and sucrose content, it can be hypothesized that O-methylation of intermediate products in the metabolic pathways of these drought-response compounds plays an important role in improving barley tolerance to water deficit.

A putative candidate gene for α-tocopherol accumulation is MLOC\_63263 encoding an ortholog of Arabidopsis AT3G47450 (NITRIC OXIDE SYNTHASE1; AtNOS1). Nitric oxide (NO) is a short-lived, highly membrane permeable, inorganic free radical. It is a crucial signaling molecule mediating diverse physiological processes and defense mechanisms in plants, likewise ROS (Siddiqui et al., 2011). The increased accumulation of NO was reported in wheat and pea under osmotic stress, in parsley under water deficit conditions, and in tobacco under various types of abiotic stimuli (Gould et al., 2003; Kolbert et al., 2005). It is worth noting that NO participates in the ABAmediated control of stomatal aperture by negative feedback regulation (Wang et al., 2015). The recent study of Yang et al. (2015) showed an interplay between NO and ROS homeostasis. They demonstrated that NO positively regulates the activity of Arabidopsis ascorbate peroxidase APX1 by S-nitrosylation of cysteine residue, which increases its enzymatic activity of hydrogen peroxide detoxification.

#### Gudys et al. CGs Prioritization Within QTL Regions

#### CONCLUSION

To dissect complex drought tolerance traits and decipher genes underlying resistance mechanisms in barley we employed multidisciplinary approach integrating the latest achievements and tools of physiology, genetics and genomics. Our study proved that thorough physiological/biochemical analysis merged with precisely performed QTL mapping, followed by the use of prioritization tools to dissect genes within physical QTL confidence intervals, is an effective approach to designate the most relevant CGs underlying quantitatively inherited traits. The list of CGs selected on the basis of GO prioritization can be further narrowed down by identifying those that are differentially expressed under drought stress, located within overlapping QTL confidence intervals (hotspots), or have orthologs in model species with well-documented engagement in drought stress response. However, further detailed experimentation using mutants and/or overexpression lines is required to confirm the individual CG relation to the studied trait, before breeding application.

#### AUTHOR CONTRIBUTIONS

IS and PK conceived and designed the study; MF, JB-K, and JK designed and coordinated drought stress experiment, contributed to the analysis of the physiological data; KZ and AO performed drought stress experiment; MD, JK, BJ, KH, and KZ generated and elaborated data on physiological and biochemical traits; JG-W and IS coordinated genetic part of the study; KG, JG-W, and DG identified SSRs for mapping; AJ, KG, and JG-W performed selection of functional CGs for mapping; KG, JG-W, JS, and WU genotyped two RIL populations; KG analyzed all genotyping data and constructed two RILs maps

#### REFERENCES

and the consensus function map; AK, KM, PO, KK, and PK performed LCam population genotyping and constructed LCam genetic map; KG performed QTL mapping, integration QTLs to physical map, GO enrichment analysis, and together with IS and JG-W performed analysis and interpretation of data; KG, JG-W, and IS wrote the manuscript; All authors commented on the manuscript.

#### FUNDING

This work was supported by the European Regional Development Fund through the Innovative Economy Program for Poland 2007–2013, project WND-POIG.01.03.01-00-101/08 POLAPGEN-BD, Biotechnological tools for breeding cereals with increased resistance to drought. All information about the project can be found at www.polapgen.pl.

#### ACKNOWLEDGMENTS

The authors thank Dr. Robbie Waugh of the James Hutton Institute (JHI) for the possibility of the usage of Barley Oligonucleotide Pool Assays (BOPA1), and Dr. Joseph DeYoung of the Southern California Genotyping Consortium (SCGC) at the University of California Los Angeles (UCLA) for genotyping the RIL populations with BOPA1 and providing the SNP marker data.

#### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2018. 00769/full#supplementary-material

a complex pattern of genome-wide transcriptional regulation. Plant, J. 53, 90–101. doi: 10.1111/j.1365-313X.2007.03315.x

by cold treatment in the improvement of tissue regeneration efficiency. Acta Physiol. Plant. 34, 547–560. doi: 10.1007/s11738-011-0852-3

soybean primary root. I. Proteomic analysis reveals region-specific regulation of phenylpropanoid metabolism and control of free iron in the elongation zone. Plant Cell Environ. 33, 223–243. doi: 10.1111/j.1365-3040.2009.02073.x

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Gudys, Guzy-Wrobelska, Janiak, Dziurka, Ostrowska, Hura, Jurczyk, Zmuda, Grzybkowska, ˙ Sróbka, Urban, Biesaga-Koscielniak, Filek, ´ Koscielniak, Mikołajczak, Ogrodowicz, Krystkowiak, Kuczynska, Krajewski and ´ Szarejko. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Activation of ABA Receptors Gene GhPYL9-11A Is Positively Correlated with Cotton Drought Tolerance in Transgenic Arabidopsis

Chengzhen Liang1† , Yan Liu1† , Yanyan Li1† , Zhigang Meng1 , Rong Yan1,2, Tao Zhu1 , Yuan Wang1 , Shujing Kang1 , Muhammad Ali Abid1 , Waqas Malik1,3, Guoqing Sun1 , Sandui Guo1 \* and Rui Zhang1 \*

1 Biotechnology Research Institute, Chinese Academy of Agricultural Sciences, Beijing, China, 2 College of Agronomy and Biotechnology, Southwest University, Chongqing, China, 3 Genomics Lab, Department of Plant Breeding and Genetics, Bahauddin Zakariya University, Multan, Pakistan

#### Edited by:

Hui Xia, Shanghai Agrobiological Gene Center, China

Reviewed by:

Abid Ali, University of Agriculture, Faisalabad, Pakistan Andrew Wood, Southern Illinois University Carbondale, United States

#### \*Correspondence:

Rui Zhang<EMAIL_ADDRESS>Sandui Guo<EMAIL_ADDRESS>
†These authors have contributed equally to this work.

#### Specialty section:

This article was submitted to Crop Science and Horticulture, a section of the journal Frontiers in Plant Science

Received: 11 May 2017 Accepted: 04 August 2017 Published: 23 August 2017

#### Citation:

Liang C, Liu Y, Li Y, Meng Z, Yan R, Zhu T, Wang Y, Kang S, Ali Abid M, Malik W, Sun G, Guo S and Zhang R (2017) Activation of ABA Receptors Gene GhPYL9-11A Is Positively Correlated with Cotton Drought Tolerance in Transgenic Arabidopsis. Front. Plant Sci. 8:1453. doi: 10.3389/fpls.2017.01453 The sensitivity to abscisic acid (ABA) by its receptors, pyrabactin resistance-like proteins (PYLs), is considered a most important factor in activating the ABA signal pathway in response to abiotic stress. However, it is still unknown which PYL is the crucial ABA receptor mediating response to drought stress in cotton (Gossypium hirsutum L.). Here, we reported the identification and characterization of highly induced ABA receptor GhPYL9-11A in response to drought in cotton. It is observed that GhPYL9-11A was highly induced by ABA treatment. GhPYL9-11A binds to protein phosphatase 2Cs (PP2Cs) in an ABA-independent manner. Moreover, the GhPYL-11A-PP2C interactions are partially disrupted by mutations, proline (P84) and histidine (H111), in the gate-latch region. Transgenic Arabidopsis overexpressing GhPYL9-11A plants were hypersensitive to ABA during seed germination and early seedling stage. Further, the increased in root growth and up regulation of drought stress-related genes in transgenic Arabidopsis as compared to wild type confirmed the potential role of GhPYL9-11A in abiotic stress tolerance. Consistently, the expression level of GhPYL9-11A is on average higher in drought-tolerant cotton cultivars than in drought-sensitive cottons under drought treatment. In conclusion, the manipulation of GhPYL9-11A expression could be a useful strategy for developing drought-tolerant cotton cultivars.

Keywords: ABA, ABA receptor, GhPYL9, PP2C, seed germination, drought stress, cotton

# INTRODUCTION

Abscisic acid (ABA) is a multifunctional phytohormone that regulates multiple aspects of plant growth and development, including seed dormancy and germination (Fang et al., 2008), vegetative development (Finkelstein et al., 2002), root elongation (Xing et al., 2016), leaf senescence (Liang et al., 2014), and fruit ripening (Zhang et al., 2009). ABA has also long been known to be a key endogenous messenger in the response to abiotic and biotic stresses, such as drought, salinity, cold, and pathogen infection (Choi et al., 2000; Xiong et al., 2002; Zhu, 2002; Sheard and Zheng, 2009; Yoshida et al., 2014). ABA-regulated stress responses include changes in gene expression, increased

**71**

stomatal closure, reduced transpiration rate, protection of photosynthesis, and regulation of plant growth.

Previous studies have demonstrated that pyrabactin resistance-like (PYL), clade A protein phosphatase 2C (PP2C), and sucrose non-fermenting 1- related protein kinase 2 (SnRK2) proteins are the three core components of the ABA signaling cascade (Park et al., 2009; Sheard and Zheng, 2009; Zhu, 2016). The primary step in initiating the ABA signaling pathway is triggered by the PYL ABA receptors (Sheard and Zheng, 2009), which contain a central hydrophobic ligand-binding pocket (Iyer et al., 2001). In the absence of ABA, two β-loops, termed the gate and latch, together with several nearby structural elements form a large open pocket for ligand binding. When ABA is present, ABA binds to this pocket, leading to a series of conformational rearrangements, including gate and latch closure, which sequester ABA within the pocket. The ABAbound PYL receptor physically interacts with and sequesters the constitutive SnRK2 repressor, PP2C, allowing SnRK2 kinases autophosphorylation. ABA-activated SnRK2 kinases subsequently phosphorylate and activate downstream ABAresponsive element-binding transcription factors that induce the expression of ABA-responsive genes (Melcher et al., 2009; Miyazono et al., 2009; Nishimura et al., 2009).

The PYL proteins belong to the START/Bet V1 superfamily (Iyer et al., 2001). There are fourteen PYL members in Arabidopsis, and this family is divided into three subfamilies based on amino acid sequences similarity (Zhao et al., 2013, 2016). All PYL proteins have highly conserved amino-acid residues in key interaction domains, including ABA-binding pocket, and the two entrance loops. Despite these similar structural features, PYLs have distinct properties and different expression patterns that explain the versatility of ABA signaling in plant growth and development. For example, Arabidopsis PYR1, PYL1, and PYL2 are able to form dimers in solution, whereas PYL proteins 4–10 exist as monomers (Hao et al., 2011). PYR1, PYL1, PYL2, and PYL3 interact with several clade A PP2Cs in an ABA-dependent manner. In contrast, most other PYL proteins except for PYL7 inhibit PP2Cs independently of ABA (Bai et al., 2013; Li et al., 2013; Zhao et al., 2013). Recently, several genetic studies have elucidated the significance of these receptors in ABA response; pyr1 pyl1 pyl2 pyl4 quadruple mutants show ABA insensitivity (Park et al., 2009; Nishimura et al., 2010; Gonzalez-Guzman et al., 2012), and a pyl8 single mutant has altered ABA responses during seed germination and seedling growth, as well as decreased drought resistance (Antoni et al., 2013; Zhao et al., 2014). PYL9 enhances drought resistance by limiting transpiration water loss and regulating senescence in old leaves and growth in young tissues (Zhao et al., 2016). Moreover, both PYL9 and PYL8 act as critical regulators of lateral root formation in response to ABA (Xing et al., 2016).

Upland cotton (Gossypium hirsutum L.) is an important economic crop and is an important source of fiber and oil (Wang et al., 2012, 2017; Shan et al., 2014). Compared with the model plants Arabidopsis and rice, cotton has higher tolerance to drought and salt stress (Liang et al., 2016). However, its growth and development, as well as fiber yield and quality are also significantly affected by severe environmental conditions, especially drought stress (He et al., 2016). To date, a number of candidate genes for drought tolerance in cotton have been identified. Several stress-related genes, including GhABF2 (Liang et al., 2016), GhNAC2 (Gunapati et al., 2016), GhATAF1 (He et al., 2016), GhbHLH130 (Guang et al., 2014), GhSARP1 (Liu et al., 2016), GhMKK3 (Wang C. et al., 2016), Di19-1 and Di19-2 (Qin et al., 2016), have been shown to regulate stress response in an ABA-dependent manner. Therefore, understanding ABA signal transduction in cotton will accelerate the molecular breeding of stress-tolerant cotton cultivars.

Because ABA has an essential role in plant growth and development, understanding how each PYL protein affects plant physiology has both basic and applied agricultural significance. Although there have been a number of studies of ABA receptors in the model plant Arabidopsis and other crops, such as rice and soybean, functional analysis of PYLs in upland cotton has not been reported.  In this study, we carried out genome-wide analysis to identify all PYL genes in the two diploid progenitor species, G. raimondii and G. arboretum, and two tetraploids G. hirsutum and G. barbadense. Genome-wide expression profiling analysis revealed that GhPYL9-11A was the most highly expressed PYL gene during drought stress in upland cotton. Consistent with a function in drought stress response, overexpression of GhPYL9-11A in Arabidopsisremarkably increases ABA sensitivity and enhances drought tolerance. The future characterization of GhPYL9-11A will further our understanding of ABA-mediated stress response in cotton, and facilitate its application in breeding drought tolerant cotton cultivars.

#### MATERIALS AND METHODS

#### Plant Materials and Growth Conditions

Upland cotton (G. hirsutum) seeds of cultivar Y18R, a restores variety developed by our lab (Shi et al., 2012), were planted during the regular cotton cultivation season in 2015 and 2016 in the experimental farm at the Biotechnology Research Institute, Chinese Academy of Agricultural Sciences (Beijing, China). The recommended cultural and plant measures were adopted during the growth of cotton crop. For expression analysis, various organs including leaf, root, stem, and seed, were collected after 1 week of the onset of flowering stage. For analysis of ABA-induced GhPYL9-11A expression, cotton plants were grown in the Greenhouse with a 12-h-light (30◦C)/12-h-dark (26◦C) photoperiod with about 300 µM m−2 s −1 photon density and 45% humidity, and the leaves of 3-weeks old plants were sprayed with 100 µM ABA (Sigma– Aldrich, United States). Ten plants were collected at 0, 1, 3, and 6 h after treatment and stored at −70◦C for RNA isolation.

Arabidopsis thaliana Col-0 seed and GhPYL9-11A transgenic plants were sterilized for approximately 10 min in 10% bleach

and then rinsed in sterile deionized water three times. Sterilized seeds were sown on petri dishes containing 0.8% agar media with 1/2 MS nutrients and 1% sucrose. After stratification at 4 ◦C for 2 days, the dishes were moved to a growth chamber at 22◦C with a 16-h light/8-h dark cycle and 70% humidity. All seedlings were grown vertically before transplanting to media supplemented different concentrations of ABA and mannitol.

#### Gene Cloning and Expression Analysis

The full-length open reading frame of GhPYL9-11A was amplified from G. hirsutum L. using the Invitrogen RACE system (Invitrogen, United States). Total RNA was extracted using the EASYspin reagent (YPHBio, China). After DNase treatment, approximately 2 µg total RNA was reverse-transcribed using the ReverTra Ace qPCR RT Master Kit (Toyobo, Japan). Quantitative RT-PCR (qRT-PCR) analysis was performed using SYBR Green I PCR mix (SSoFast EvaGreen Supermix, Bio-Rad, United States) and the Bio-Rad CFX96 real time PCR system. Data were analyzed with Opticon monitor software (Bio-Rad). The cotton GhHISTON3 (GhHIS3) and Arabidopsis Actin 8 genes, both with stability expression in the different tissues, developmental stages and environmental conditions (Zhu et al., 2017), were used as internal controls. Primers used are listed in **Table S1**. Values are means ± SD of three biological repeats. Student's t-test was used for statistical analysis.

#### Vector Construction and Arabidopsis Transformation

The GhPYL9-11A full-length coding sequence was cloned into pBI121-35S to generate the pBI121-35S:GhPYL9-11A overexpression construct. The constructs were transfected into Agrobacterium tumefaciens GV3101 by electroporation, and then transformed into Arabidopsis Col-0 pants using the floral dipping method (Liu L. et al., 2015). Homozygous transgenic Arabidopsis lines were obtained, and the lines GO16, GO21, and GO22, which have high levels of GhPYL9-11A expression, were selected for further analysis. The primers used for vector construction are listed in **Table S1**.

#### Yeast Two-Hybrid Assay

The Y2H assay was performed using the MatchmakerTM Gold Yeast Two-hybrid System (Clontech, United States). Full length cDNA of GhPYL9-11A was amplified and cloned into the pGBKT7 vector, and AtABI1, AtABI2, GhPP2C1, and GhPP2C2 full length cDNAs were amplified and cloned into the pGADT7 vector. GhPYL9-11A point mutations were introduced by site-directed mutagenesis using pBI121- 35S:GhPYL9-11A as a template, and GhPYL9-11AP84S and GhPYL9-11AH111A were cloned into the pGBKT7 vector. Each pair of AD and BD constructs was co-transformed according to the manufacturer's protocol. The positive yeast liquid culture was serially diluted to OD600 = 0.6 and 3 µL of dilutions 1:10, 1:100, 1:1,000, and 1:10,000 were inoculated onto various plates [Synthetic Dropout Medium (SD) -leucine (Leu)/tryptophan (Trp), SD-Leu/Trp/histidine (His), and SD-Leu/Trp/His supplemented with 10 µM ABA]. All plates were incubated at 28◦C for 2 days and then photographed.

#### Germination Assay

Surface-sterilized seeds were sown on 1/2 MS without or with ABA (0.3 and 0.5 µM) or mannitol (150 and 300 µM). After 5 days seeds with green shoots were scored as germinated.

#### Drought Stress Treatment

To test the drought tolerance of 35S:GhPYL9-11A transgenic Arabidopsis, the homozygous F2 generation was subjected to drought stress in soil in a growth chamber. Seed from GO1, GO2, GO3, and wild type lines were imbibed at 4◦C for 2 days and directly planted in soil in the same pots (30 cm

(Phaseolus vulgaris), VA (Vigna angularis), THC (Theobroma cacao), GORA (Gossypium raimondii), GH (Gossypium hirsutum), GA (Gossypium arboreum), EG (Eucalyptus grandis), MD (Malus domestica), CICLE (Citrus clementine), PT (Populus trichocarpa), PEU (Populus euphratica). (B) Amino acid alignment of the PYL9 CL2/gate-latch region for the proteins shown in (A). The black line indicates the CL2 domain. Two amino acid residues 'VK' are marked with a red open rectangle. The residues in contact with the ABA molecule are designated by asterisk. (C–F) Comparison of the predicted 3-D structures of AtPYL9 and GhPYL9-11A. The conserved proline and histidine residues in the gate and latch region are marked with red arrows in AtPYL9 (E) and yellow arrows in GhPYL9-11A (F). (G) Expression of GhPYL9-11A in various organs determined by quantitative real time (qRT)-PCR analysis. Roots, leaves, and stems were harvested from 3-week-old plants. (H) Time-course of the response of GhPYL-11A to 100 µM ABA treatment.

in diameter). All plants were grown at room temperature (22–25◦C) under long-day conditions (16 h light/8 h dark) in the culture room. After 2 weeks of seedling emergence, drought treatment was established by withholding the water for 4 weeks. The experiment was comprised of five replications. Watering was resumed 1 day after the 4 weeks of drought

treatment, and the survival rate (SR) was calculated 3 days later.

# Chlorophyll Measurement

After measuring the fresh weight of the samples, five seedlings for each sample were used to calculated the total chlorophyll contents as the absorbance at 652 nm following (Liang et al., 2014).

# Photosynthesis Parameters and Water Loss Measurement

For measurement of photosynthesis parameters, water was withheld for 7 days beginning when plants were 14-days old. Ten independent plants were used for each GhPYL9-11A overexpression line. Photosynthesis parameters were measured as previously described by (Liang et al., 2015). For determination of water loss, the whole aboveground rosettes of 21-day-old plants were cut from the base and weighed at different time points.

# Electrolyte Leakage Measurement

The measurement of electrolyte leakage was performed according to a previously described method (Lin A. et al., 2012). Ten different leaf disks from wild-type and GhPYL9-11A overexpression plants were placed into a flask containing 10 mL distilled deionized water. After shaking at room temperature for 6 h at approximately 120 rpm, the conductivity (Ci) was measured with a conductivity meter (Leici-DDS-307A). Then, the disks were boiled for 20 min to kill the leaf tissues and shaken for 1 h to completely release the electrolytes into the solution. The conductivity (Cm) was measured again. The electrolyte leakage was calculated according to the equation (Ci/Cm) × 100%.

#### Hydrogen Peroxide Production and Measurement of Antioxidant Enzyme Activities

For hydrogen peroxide production and measurement of antioxidant enzyme activities, Arabidopsis plants were grown on soil for 14 days and leaves were collected from the plants after water was withheld for 7 days. Measurement of H2O2 production was performed according to a previously described method (Liang et al., 2015), and quantified using the Hydrogen peroxidase assay kit (Beyotime, China) according to the manufacturer's instructions. The quantification of CAT, SOD, and POD activity was performed according to a previously described method (Liang et al., 2016). Total protein was measured using Bradford protein assay kit (Sangon Biotech, China).

# Protein Sequence and Phylogenetic Analysis

For phylogenetic analysis of GhPYL9-11A homologs in plants, 20 protein sequences were obtained through National Center for Biotechnology Information database using BLASTp searches

with the GhPYL9-11A protein as a query. All sequences were aligned using Clustal W. Phylogenetic trees were constructed using MEGA 6.0 based on the neighbor-joining method. Topological robustness was assessed by bootstrap analysis with 1,000 replicates (Liu L. et al., 2015).

#### RESULTS

#### Transcriptome Analysis of Cotton PYL ABA Receptor

Pyrabactin resistance-like protein coding sequences were identified in upland cotton species, G. hirsute, G. barbadense, G. raimondii, and G. arboretum, with complete genomes or genome assemblies. We then compared the expression profiles of the PYL genes under drought stress. Twentytwo PYL proteins were identified from each of the two diploid progenitor species G. raimondii and G. arboretum, and 44 and 36 members were identified in the tetraploids G. hirsutum and G. barbadense, respectively (**Table S2**). High-resolution temporal expression profiling revealed that two GhPYL genes, Gh\_A11G08070 and Gh\_D11G1013, were most highly expressed during drought stress in upland cotton (G. hirsutum) (**Figure 1** and **Table S2**). The proteins encoded by Gh\_A11G0870 and Gh\_D11G1013 share approximately 80 and 79% amino acid identity with AtPYL9, an ABA receptor, respectively.

Both Gh\_A11G0870 and Gh\_D11G1013 belong to the PYL9 subfamily of PYL proteins. There are two members each in G. raimondii and G. arboretum, and four members each in G. hirsutum and G. barbadense. Based on phylogenetic analysis, the PYL proteins could be further divided into two distinct subgroups, PYL9-I and PLY9-II (**Figure S1**). Interestingly, all members of PYL9-I in both diploid and tetraploid cotton were highly expressed during seed germination and drought stress, whereas PLY9-II members were with low expression level (**Table S2**). These results suggest that PYL9-I proteins may play important roles in the regulation of seed germination and drought responses in cotton. Because Gh\_A11G08070 (hereafter referred to as GhPYL9-11A) belongs to the drought stressresponsive PYL9-I subfamily, its function was investigated in detail.

#### Structural Features of GhPYL9-11A

Using protein database searches, we identified a number of putative GhPYL9-11A orthologs in various higher plants. GhPYL9-11A homologs were only found in the plant kingdom, indicating that PYL9 proteins are plant-specific. The identification of PYL proteins in both monocots and dicots suggests that PYL proteins are evolutionally conserved across plant species. The existence of PYL9 in different plant species might also suggest that its function in drought stress response is also conserved. Phylogenetic analysis of the cotton GhPYL9-11A protein and its 20 homologs showed that GyhPYL9-11A is most closely related to LOC18613997 of Theobroma cacao, and shares approximately 80% amino acid sequence identity with related proteins in **Figure 2A**.

PYL9 family proteins have a pyrabactin resistance1/pyr1 like/regulatory component of ABA receptor (PYR/PYL/RCAC) like structure (amino acids 28–172) typical of ABA receptors (**Figure 2B**). We analyzed the 3D structure of this conserved motif in GhPYL9-11A using SWISSMODEL. The predicted structure is a common helix-grip structural feature consisting of four α-helixes and seven antiparallel β-sheets and is highly similar to the Arabidopsis AtPYL9 PYR/PYL/RCAC motif (**Figures 2C,D**).

#### Expression Pattern of GhPYL9-11A

Quantitative RT-PCR analysis showed that GhPYL9-11A is expressed preferentially in seeds and roots, with a low level of expression in the leaf and stem (**Figure 2E**). Moreover, time-course analysis revealed that GhPYL9-11A was rapidly induced by ABA, with expression levels increasing approximately 2.8-fold after 1 h of ABA application (**Figure 2F**). Levels of GhPYL9-11A transcription reached a maximum after 4 h of treatment, and were approximately eightfold higher than without hormone treatment. These

results clearly demonstrate that this gene is transcribed in response to ABA and thus is potentially involved in ABA response.

#### GhPYL9-11A Interacts with Clade A PP2Cs in an ABA-Independent Manner