Identification and use of high oil yielding genotypes is important for increasing oil content in safflower cultivars. Breeding efforts in America led to the development of cultivars with increased seed oil content ranging from 45 to 55% (Bergman et al., 1985; Rubis et al., 2001). However, such improvements are lacking among Indian cultivars which have oil content ranging from 27 to 35%. Evaluation of oil quantity in the 531 accessions by NIRS identified 15 accessions with high oil content (>40%). These would serve as important breeding material in safflower. All the high oil yielding accessions (>40%) had low 100-seed weight (3–4 gm) in our study. This observation is in consonance with earlier reports, which suggest that increased hull thickness enhances seed weight but reduces oil content (Ranga Rao et al., 1977; Dajue and Mündel, 1996). Safflower oil has a desirable fatty acid composition. High linoleic lines of safflower are favored for animal feed and in the paint and varnish industry (Knowles, 1989; Bergman et al., 2001) while high oleic lines are nutritionally desirable because of its hypo-cholesterolemic effect and greater oxidative stability (Fuller et al., 1967). A high oleic line of Indian origin (Knowles and Bill, 1964) was effectively utilized in various safflower breeding programs in the USA (Mündel and Bergman, 2009). The safflower collection used in this study contained 17 accessions with high oleic acid content and 15 accessions with high linoleic acid content (**Table 1**). Accessions with desirable traits identified in the present study could be incorporated in breeding programs for crop improvement.

### Assessment of Regional Gene Pools Based on UPGMA Analysis of Phenotypic Data

Accessions from the Indian Subcontinent and American region were distributed in all the five clusters (**Figure 1**) suggesting that they harbor maximum phenotypic diversity for the studied

traits. Knowles (1969), based on morphological analysis of accessions from the Indian subcontinent, reported them as a uniform assemblage resulting from a single introduction. In contrast, our assessment indicates that accessions from the Indian subcontinent are phenotypically diverse. Morphological diversity among accessions from the Indian subcontinent was reported in earlier studies (Kupsow, 1932; Chavan, 1961; Hanelt, 1961). Indian accessions have also been shown to harbor significant genetic diversity (Kumar et al., 2015). The American germplasm was found to be phenotypically diverse in the current study but was genetically conserved (Kumar et al., 2015). Near East and Iran-Afghanistan accessions clustered together based on phenotypic data, supporting our earlier proposal of considering them as a single gene pool (Kumar et al., 2015). Accessions from European region were distributed in several clusters based on phenotypic data similar to the observation obtained through molecular data analysis. Interestingly, accessions from Far East, Turkey and Egyptian region were present in all the clusters although they exhibited low genetic diversity (Kumar et al., 2015). These results indicate that UPGMA analysis based on phenotypic data alone is unable to accurately define the genetic relationships among safflower accessions.

#### Composite Core Collections Effectively Capture the Global Molecular, Phenotypic, and Geographical Variability of the Crop

In recent years, increased availability of molecular resources has enabled their utilization in development of core collections in crop species (Belaj et al., 2012; El Bakkali et al., 2013) but until now, no such attempts have been made in safflower. Use of molecular markers for development of core collections is advantageous as they reflect diversity at the DNA level as opposed to morphological markers wherein different genotypes might show similar phenotypic traits due to environmental effects. Additionally, molecular markers are more effective in identifying and minimizing redundancy. Several studies have emphasized on use of maximization (M) strategy for development of highly robust core collections (Bataillon et al., 1996; McKhann et al., 2004). The M strategy retains maximum number of alleles at each locus and is considered as the most powerful approach for maintaining diverse alleles (Schoen and Brown, 1993). MSTRAT and POWERCORE programs have been successfully used for construction of core collection in various plant species such as grapes, olive and sesame (Le Cunff et al., 2008; Belaj et al., 2012; Zhang et al., 2012). A combination of molecular markers and maximization (M) strategy has been utilized for the first time in our study for construction of a core collection in safflower.

Earlier studies reported significant GE interactions in safflower and emphasized on multi-location and multi-seasonal trials to evaluate heritability of characters for their effective utilization in breeding programs (Singh et al., 2004; Mahasi et al., 2006). In our study, MANOVA analysis indicated prominent GE interactions (**Table 4**). Therefore, seasonal datasets were treated independently for developing core collections. The six core collections thus generated, efficiently captured the entire range of trait variability but failed to include complete genetic diversity represented in 19 clusters derived earlier (Kumar et al., 2015) through Bayesian analysis. Additionally, many accessions were common between different core collections. For example, in core collections developed using POWERCORE, 10 accessions were common between CC1 (marker-based) and CC2/CC3 (phenotype-based). Only 4 accessions were unique to CC1 while 16 and 17 accessions were unique to CC2 and CC3, respectively. In MSTRAT-derived core collections, 19 accessions were common between CC4 (marker-based) and CC5/CC6 (phenotype-based). The number of accessions unique to CC4, CC5, and CC6 were 7, 28, and 35, respectively. The presence of common accessions between core collections derived using different types of data indicates an overlap in genetic and phenotypic components of the studied accessions. These accessions represent a subset of genotypes that are highly diverse at both molecular and phenotypic level.

The core collections developed using each program were merged to derive a more robust and non-redundant composite core collection (CartC1 by POWERCORE and CartC2 by MSTRAT). The vast phenotypic diversity of the initial collection was retained in both collections. Accessions with desirable agronomic traits and extreme phenotypes, which were present in very low numbers in the entire collection and scattered in the initial core collections were captured in the composite core collections (**Table 8**). Both the composite core collections provided comprehensive coverage of allelic diversity and had representation from all the 19 BAPS clusters identified earlier for safflower (Kumar et al., 2015; **Table 7**). Evaluation indices (MD%, VD%, VR%, CR%, I, H) for CartC1 and CartC2 were

comparable and reflect their effectiveness in capturing diversity of the crop (**Table 6**). Our approach of deriving independent core collections from molecular and phenotypic data and their subsequent merger to create composite core collections avoided trade-off between the diversity captured using the molecular and phenotypic data sets.

Geographical distribution influences the extent of genetic variability of a species. The effect is more prominently seen in case of in-breeding species (Rao and Hodgkin, 2002). Geographical patterning is evident in safflower which is highly self-pollinating in nature and is grown in different agro-climatic regions across the world (Knowles, 1969; Ashri, 1975; Chapman et al., 2010). The two composite core collections showed minor variations in representation of the 10 regional gene pools. CartC1 included 8 regional pools excluding Sudan and Kenya while CartC2 contained representation from all the 10 regional gene pools (**Table 5**). Similar to the entire collection, both CartC1 and CartC2 showed predominance of accessions from Indian subcontinent and America accounting for ∼50% of the total entries. In contrast, the earlier core collection developed by Johnson et al. (1993) had a major proportion of accessions (∼46%) from the Mediterranean region and South-West Asia while the core collection derived by Dwivedi et al. (2005) consisted of ∼78% accessions from South and South-East Asia.

The number of accessions in a core collection is an important factor determining its effective utilization (Brown and Spillane, 1999). The core collections developed earlier for safflower consisted of 210 accessions (Johnson et al., 1993) and 570 accessions (Dwivedi et al., 2005) while the composite core collections developed in the present study are comparatively smaller with 57 (CartC1) and 106 (CartC2) accessions. The larger size of the core collections developed in earlier studies could be due to the larger number of accessions in their initial germplasm collection. However, the advantage of the present study is that the initial collection used for development of composite core collections has been characterized extensively for both molecular and phenotypic diversity and the generated core collections have therefore effectively captured the global genetic and phenotypic diversity of the crop. Additionally, CartC1 has better utility value in comparison to CartC2 due to its smaller size and comparable diversity.

The present study is the first attempt where molecular diversity data has been used in conjunction with phenotypic data and geographical distribution to develop core collections in safflower. The small size of the composite core collections would be advantageous for field studies and association mapping. These collections will provide access to genetically diverse and agronomically important germplasm that would be useful in widening the genetic base of the crop and facilitate characterization of genetic determinants of trait variability. This information can be used to design more effective breeding programs to increase the global utility of safflower as an oilseed crop.

#### AUTHOR CONTRIBUTIONS

AJ, SG conceived and designed the experiments. SK, HA performed the experiments, analyzed the data. MV helped in

#### REFERENCES

collection of phenotypic data. AR provided necessary support in the statistical analysis. SK, HA, SG, AJ wrote the manuscript. MA, AK, AJ, SG provided facilities for completion of experiments and reviewed the manuscript.

#### ACKNOWLEDGMENTS

This work was supported by the DST-PURSE grant of Department of Science and Technology, Government of India (grant no. Dean(R)/2009/868) provided to the University of Delhi. HA was aided by a research fellowship from University Grants Commission, India. We sincerely thank the reviewers for their critical comments and suggestions to improve the manuscript.

#### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fpls.2016. 01554

Llimensee, and W. J. Peacock (Cambridge: Cambridge University Press), 161–170.

rice (Oryza rufipogon Griff.) populations in China. PLoS ONE 10:e0145990. doi: 10.1371/journal.pone.0145990

for diversified utilization of germplasm. Crop Sci. 49, 1769–1780. doi: 10.2135/cropsci2009.01.0014

**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 © 2016 Kumar, Ambreen, Variath, Rao, Agarwal, Kumar, Goel and Jagannath. 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) or licensor 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.

# Trichome-Related Mutants Provide a New Perspective on Multicellular Trichome Initiation and Development in Cucumber (*Cucumis sativus* L)

#### Xingwang Liu1, 2, Ezra Bartholomew1, 2 , Yanling Cai 2 and Huazhong Ren1, 2 \*

1 College of Horticulture, China Agricultural University, Beijing, China, 2 Beijing Key Laboratory of Growth and Developmental Regulation for Protected Vegetable Crops, China Agricultural University, Beijing, China

Trichomes are specialized epidermal cells located in aerial parts of plants that function in plant defense against biotic and abiotic stresses. The simple unicellular trichomes of Arabidopsis serve as an excellent model to study the molecular mechanism of cell differentiation and pattern formation in plants. Loss-of-function mutations in Arabidopsis thaliana have suggested that the core genes GL1 (which encodes a MYB transcription factor) and TTG1 (which encodes a WD40 repeat-containing protein) are important for the initiation and spacing of leaf trichomes, while for normal trichome initiation, the genes GL3, and EGL3 (which encode a bHLH protein) are needed. However, the positive regulatory genes involved in multicellular trichrome development in cucumber remain unclear. This review focuses on the phenotype of mutants (csgl3, tril, tbh, mict, and csgl1) with disturbed trichomes in cucumber and then infers which gene(s) play key roles in trichome initiation and development in those mutants. Evidence indicates that MICT, TBH, and CsGL1 are allelic with alternative splicing. CsGL3 and TRIL are allelic and override the effect of TBH, MICT, and CsGL1 on the regulation of multicellular trichome development; and affect trichome initiation. CsGL3, TRIL, MICT, TBH, and CsGL1 encode HD-Zip proteins with different subfamilies. Genetic and molecular analyses have revealed that CsGL3, TRIL, MICT, TBH, and CsGL1 are responsible for the differentiation of epidermal cells and the development of trichomes. Based on current knowledge, a positive regulator pathway model for trichome development in cucumber was proposed and compared to a model in Arabidopsis. These data suggest that trichome development in cucumber may differ from that in Arabidopsis.

Keywords: unicellular, *arabidopsis*, multicellular, cucumber, mutants, trichome-related genes, regulator pathway

#### INTRODUCTION

Plant trichomes are highly specialized epidermal protrusions that are located on the surfaces of leaves, stems, petioles, sepals, seed coats, and other aerial organs. Their diversity is almost as great as the number of species on which they are found. Morphologically, they can be unicellular or multicellular as well as secretory glandular or non-glandular (Hülskamp et al., 1998; Hülskamp, 2004; Tissier, 2012; Chen et al., 2014). In the model plant Arabidopsis thaliana,

#### *Edited by:*

John Doonan, Aberystwyth University, UK

#### *Reviewed by:*

Henrik Buschmann, University of Osnabrück, Germany Amy T. Hark, Muhlenberg College, USA

> *\*Correspondence:* Huazhong Ren<EMAIL_ADDRESS>
#### *Specialty section:*

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

> *Received:* 20 May 2016 *Accepted:* 22 July 2016 *Published:* 10 August 2016

#### *Citation:*

Liu X, Bartholomew E, Cai Y and Ren H (2016) Trichome-Related Mutants Provide a New Perspective on Multicellular Trichome Initiation and Development in Cucumber (Cucumis sativus L). Front. Plant Sci. 7:1187. doi: 10.3389/fpls.2016.01187 extensive studies have been performed on unicellular trichome development, especially on leaves (Hülskamp, 2004; Ishida et al., 2008; Pesch and Hülskamp, 2009). Classical molecular genetic approaches have identified several regulators that work in distinct developmental processes, such as trichome initiation/formation, endo-reduplication, and branch construction, and growth orientation (Schwab et al., 2000; Szymanski et al., 2000; Chen et al., 2014). The regulatory pathways of unicellular trichome development comprise both positive (mutants develop fewer trichomes) and negative (mutants develop more and/or clusters of trichomes) transcription factors (Ishida et al., 2008; Balkunde et al., 2010; Grebe, 2012).  They bind together to form an MYB-bHLH-WD40 (MBW) trimeric complex that activates the downstream gene GL2 (GLABRA2), which initiates trichome differentiation (Pesch and Hülskamp, 2009).  They compete with the R2R3 MYB protein GL1 and bind to bHLH proteins, including GL3/EGL3, to suppress trichome initiation in adjacent cells (Wang et al., 2014).

Cucumber (Cucumis sativus L.), as one of the most important vegetable crops, is also covered with trichomes ranging from the stems and leaves to the flowers, branches, fruits, and tendrils. During early fruit development, deep ridges along the length of the fruit cover the fruit surface, and densely spaced fruit trichomes are randomly scattered relative to the ridges (Liu et al., 2011; Chen et al., 2014). Trichomes on cucumber fruit are called spines. In cucumber, the fruit spine combines with tubercles to form the warty fruit trait, which is a very important fruit quality trait (Chen et al., 2014; Li et al., 2015). Compared to warty fruit, smooth fruit, with no spines, or tubercles, is very important for the breeding of freshly eaten cucumber types, as they are easier to clean, package, transport and store (Zhang et al., 2010; Yang et al., 2014). Moreover, smooth fruit is becoming increasingly popular due to its attractive and distinctly shiny appearance (Li et al., 2015; Pan et al., 2015; Cui et al., 2016). Despite the importance of fruit spines in cucumber breeding for external quality, there are limited reports of the regulation of fruit spine development and few detailed characterizations of cucumber genes with disturbed trichome development. To understand trichome development in cucumber, we must determine the crucial regulators for its initiation and formation based on cucumber trichome-related mutants.

We begin this review by focusing on cucumber trichomerelated mutants and the role that they play in trichome formation in plants with multicellular trichomes. Several mutants, such as trichome-less (tril) (Wang et al., 2016), glabrous 3 (csgl3) (Pan et al., 2015; Cui et al., 2016), tiny branched hair (tbh) (Chen et al., 2014), micro-trichome (mict) (Zhao et al., 2015a), and glabrous 1 (csgl1) (Li et al., 2015), inhibit trichome development via different mechanisms, but all of these mutants cause a reduction in one type of spine. The possible mechanisms involved in modulating trichome development in cucumber mutants are summarized herein.

## MUTANTS WITH DISTURBED TRICHOME IN CUCUMBER

Five trichome-related mutants have been reported in cucumber (Chen et al., 2014; Li et al., 2015; Pan et al., 2015; Zhao et al., 2015a; Cui et al., 2016; Wang et al., 2016). Wild-type cucumber fruits have two types of trichomes, both of which are multicellular (**Figures 1**, **2**) (Chen et al., 2014). Type I trichomes are small, glandular trichomes, with a 3- to 5-cell base topped with a 4- to 8-cell head. Type II trichomes, which predominate, are much larger, non-glandular trichomes that are composed of a base and stalk (**Figure 2**). Compared to wild-type cucumber plants, all mutants appeared to be glabrous with no noticeable trichomes on the leaves, stems, tendrils, floral organs, or fruits (**Table 1**, **Figures 1**, **2**). Each mutant has a trichome type except for the mutants tril and csgl3 (**Figure 2**). Trichomes have been reclassified into three morphologies according to their shape on leaves and fruit surfaces. Type I trichomes, which are found in the tbh, mict, and csgl1 mutants, have a small papillar-shaped head (**Figures 2A,D,F**). Type II trichomes, which exist in the tbh and mict mutants, consist of one to five rounded cells without the pyramid-shaped head or pieshaped base (**Figures 2B,E**). However, type III trichomes occur only in tbh mutants; this type of trichome was not previously described. Here, we classified these trichomes into the new type III according to the branching at the top of the trichome (**Figure 2C**). Trichome phenotypes from tbh, mict and csgl1 indicated that these three mutants are all involved in trichome development and not in trichome initiation. Interestingly, the tril and csgl3 mutants showed a completely glabrous phenotype on the shoot and fruit epidermis (**Figures 2G,H**), suggesting that TRIL and CsGL3 may be upstream positive regulators of TBH, MICT, and CsGL1 for the regulation of multicellular trichome development and may affect epidermal cell initiation.

**Abbreviations:** TRIL, trichome-less; TBH, tiny branched hair; MICT, microtrichome; TTG1, transparent testa glabra1; bHLH, basic helix-loop-helix; GL3, glabra3; EGL3, enhancer of glabra3; MBW, MYB-bHLH-WD40; TRY, triptychon; CPC, caprice; ETC1, enhancer of tryand caprice1; HD-Zip, homeodomain-leucine zipper; PDF2, protodermal factor2; ATML1, Arabidopsis thaliana meristem layer1; MADS, MCMI-AGAMOUS-DEFICIENS-SRF4.

Compared to the multicellular trichomes found on aerial organs, such as leaves and fruits, the trichomes found on underground organs, such as root hairs, are characterized as single-celled, unbranched, elongated, and soft-structured with small tumors attached. Based on the published results, there was no difference between the wild-type and the tbh, mict, and csgl1 mutants (Chen et al., 2014; Li et al., 2015; Zhao et al., 2015a), indicating that the TBH, MICT, and CsGL1 mutants may not be involved in root hair formation. In contrast, the root length and number of branches of the tril mutant increased, suggesting that root hair formation in cucumber might be regulated by TRIL (Wang et al., 2016).

Other remarkable differences (except for the common phenotype) may exist among each mutant, such as dwarfism, branching, leaf curvature, and petal opening rates. In mict and tbh, rounded-head trichomes were found on the hypocotyls,

(Li et al., 2015), and Picture (G) is from (Wang et al., 2016).

TABLE 1 | Previous morphological studies of trichome in *tbh*, *csgl1*, *mict*, *tril and csgl3*.

but none were found in csgl1, suggesting trichome distribution is specific (Chen et al., 2014; Li et al., 2015; Zhao et al., 2015a).

#### *TBH, MICT*, AND *CSGL1* ARE THE ALLELIC WITH ALTERNATIVE SPLICING

To decipher the molecular defects in the mict, csgl1, and tbh mutants, a map-based cloning approach was undertaken by three independent groups to isolate these genes. The results showed that TBH, MICT, and CsGL1 are allelic and that there was a 2649 bp fragment deletion from -189 to 2460 bp of the start codon in Csa3M748220 in the mict, csgl1, and tbh mutants (Li et al., 2015; Zhao et al., 2015a) (**Figure 3A**). According to the gene ID that provided in the references, we extracted their proteins from the cucumber genome database (http://www.icugi.org/cgi-bin/ICuGI/genome/search.cgi) and analyzed their protein domains online (https://blast.ncbi.nlm. nih.gov/Blast.cgi). These genes were predicted to encode a class I homeodomain-leucine zipper (HD-Zip) protein consisting of different amino acid residues. For example, CsGL1 and TBH contain 240 amino acids with the conserved HD (65AA-120AA) and Zip domains (121AA-164AA), but the MICT protein consists of 242 amino acid residues with a HALZ motif. Notably, there were two isoforms of this gene in the cucumber genome database (Csa3M748820.1 and Csa3M748820.2); therefore, we infer that the different protein lengths encoded by the same gene suggest that Csa3M748220 may exist by alternative splicing in cucumber. Subcellular localization showed that the Csa3M748220 coding sequence was fused to GFP (35S) in the nuclei of tobacco and onion epidermal cells (Zhao et al., 2015a,b). Moreover, a transcriptional activation activity assay in yeast found that Csa3M748220 had weak activity as a transcriptional activator (Zhao et al., 2015a). Therefore, based on the above results, Csa3M748220 has the typical features of a transcription factor. HD-Zip proteins are unique to the plant kingdom and can be classified into four groups, I-IV, based on their distinctive traits of DNA-binding specificities, gene structures, and common motifs (Abe et al., 2003; Hülskamp et al., 2005). HD-Zip I genes have been demonstrated to be involved largely in biological processes, such as abiotic stress responses, meristem regulation,and trichome development (Hanson et al., 2001; Himmelbach et al., 2002; Hjellström et al., 2003; Saddicl et al., 2006; Zhao et al., 2015a,b).

#### *TRIL*, *CSGL3* ARE THE ALLELIC WHICH OVERRIDE THE EFFECT OF *TBH/MICT/ CSGL1* ON THE REGULATION OF MULTICELLULAR TRICHOME INITIATION AND DEVELOPMENT

Thereare few reports on the regulatory genes that control multicellular trichome development in cucumber. This presents a good opportunity to use these mutants to analyze their regulatory mechanisms in trichomes of cucumber. As mentioned above, tril and csgl3 show a completely different glabrous morphology from that of the other three trichome-developmentrelated mutants, and only the epidermal cells, including stomata and encircling guard cells, were visible (**Figure 2G**). This indicates that TRIL and CsGL3 function in trichome cell fate determination. A transcriptome profiling analysis among wildtype, tril, tbh, mict, and csgl1 revealed that TBH, MICT, and CsGL1 were not expressed in the tril mutant (baseMean 0) but were highly expressed in the wild-type (baseMean 518.25) (Zhao et al., 2015a). TRIL was mapped to Cas6M514870, a member of the class IV HD-Zip family that shares 66.7% identity with PROTODERMAL FACTOR2 (PDF2, At4g04890.1), a shoot epidermal cell differentiation-related gene, and 35% identity with GL2 (At1g79840.2), a gene that initiates trichome differentiation in Arabidopsis. In the tril mutant, there is a 5008 bp insertion fragment after the first exon (**Figure 3B**). CsGL3 was also mapped to Csa6M514870 by Pan and Cui (Pan et al., 2015; Cui et al., 2016). In Pan's study, the loss of function of the CsGL3 was due to the insertion of a 5-kb-long terminal repeat (LTR) retrotransposon in the 4th exon of CsGL3 (**Figure 3B**). The insertion location is different from that of tril. Cui used three markers (InDel-19, dCAPs-2, and dCAPs-11) designed from the sequence of Csa6M514870, which co-segregated with the trait. In addition, Csa6M514870 was found to harbor 3 single-base substitutions: T→C (611 bp), G→A (820 bp), and G→A (865 bp) at the fourth exon, resulting in a change in the amino acid sequence. We still do not know how the CsGL3 changed in its corresponding mutant csgl3.

In the tbh, mict, and csgl1 mutants, TRIL/CsGL3 was overexpressed at different levels due to the different stages at which the samples were collected (Chen et al., 2014; Li et al., 2015; Zhao et al., 2015a). Moreover, the trichome phenotype analyzed in the F2 population between the tril and mict mutants suggests that the TRIL gene has a significant influence on trichome initiation and in determining the fate of epidermal cells,

whereas the MICT/TBH/CsGL1 gene only influences trichome development at the shoot. Genetically, based on the doublemutant phenotype in csgl3csgl1, the TRIL/GsGL3 gene is assumed to act upstream of the MICT/TBH/CsGL1 gene and control the expression of MICT/TBH/ CsGL1 (Zhao et al., 2015a; Wang et al., 2016).

#### MODELS FOR TRICHOME PATTERNING IN CUCUMBER MAY DIFFER FROM MODELS IN *ARABIDOPSIS*

In the model plant Arabidopsis, the activator-inhibitor model has guided intuitive modeling and experimental design for a long time because it offers a reasonable explanation for the apparently paradoxical situation in that trichome-promoting and trichomeinhibiting genes are both expressed strongly in trichomes (Pesch and Hülskamp, 2009). Here, we focus on a positive regulatory model for trichome development between Arabidopsis and cucumber. The main reason is that all of the mutants in cucumber reveal fewer or no trichomes visually, suggesting that they encode a positive regulator of trichome development. Solid evidence of the genetic basis of trichome initiation has identified genes that (a) control the entry into the trichome pathway and (b) control the spacing of initiation events. Lossof-function mutations in Arabidopsis thaliana have suggested that the genes GL1 (which encodes an MYB transcription factor) and TTG1 (which encodes a WD40 repeat-containing protein) are important for the initiation and spacing of leaf trichomes (Galway et al., 1994; Walker et al., 1999), while for normal trichome initiation, the genes GL3, and EGL3 (Payne et al., 2000; Zhang et al., 2003), which encode helix-loop-helix (bHLH) proteins, are needed (**Figure 4A**).

In cucumber, the TRIL/CsGL3 gene encodes a HD-Zip IV protein. The loss-of-function mutant showed that tril/csgl3 has a trichome-less phenotype, indicating that this gene is a positive regulator of trichome development. Its homologous genes PDF2 and ARABIDOPSIS THALIANA MERISTEM LAYER1 (ATML1) in Arabidopsis also encode HD-Zip IV protein family members (Nakamura et al., 2006). The pdf2atml1 double mutant displays striking defects in shoot epidermal cells (Abe et al., 2003). In tomato, another typical multicellular trichome plant, Wo RNAi transgenic tomato plants showed similar shoot and root epidermal cell formation as that of the tril mutant. Wo shares an amino acid sequence identical to that of TRIL/CsGL3 and belongs to the same HD-Zip IV protein family. These results indicate that the role of TRIL/CsGL3 in the initiation and control of multicellular trichrome pathways. Based on previous results, a proposed model explaining trichome patterning in the positive regulator pathway in cucumber was built (**Figure 4B**). In this model, HD-Zip transcription factors may bind to other types of transcription factors to generate a specific complex to control cucumber multicellular trichome formation.

Based on transcriptional data from all of the mutants, several candidate genes should be focused on extensively (Chen et al., 2014; Zhao et al., 2015a,b). For example, the CsMYB6, CsWIN1, and CsGL2 genes were down-regulated not only in tril/csgl3 but also in tbh, mict, and csgl, indicating that those genes were involved in multicellular trichome development in cucumber.

#### DISCUSSIONS

In the past decade, substantial progress has been made in delineating the genes that control trichome development in

developmental progression.

a complex to positively regulate trichome development. The core complex in MBW (A) may differ from the HD-Zip complex (B). The bold arrows in the boxes indicate

cucumber. Several groups provide evidence to suggest that a number of transcriptional activators, such as HD-Zip I and HD-Zip IV, play a role in fine tuning the spatial and temporal distribution of trichomes. Researchers have tried to demonstrate the possible mechanisms of these transcriptional activators.

However, much remains unknown and needs to be elucidated in future research. The functions of the HD-Zip IV and HD-Zip I genes require further investigation through genetic transformation in cucumber plants. Our current knowledge about the gene regulatory networks is largely limited to the unicellular trichomes in the model plant Arabidopsis. However, little is known about the regulatory network that controls the development of multicellular trichomes in cucumber. Evidence indicates that the common regulatory mechanisms of unicellular trichomes in Arabidopsis or of multicellular trichomes in cotton involve plant-specific genes that function distinctively. An analysis of the differential expression data generated by RNA-seq can offer new information for identifying putative key multicellular development transcription factors in cucumber. This is why we focused on the critical transcription factor genes. A new set of more than 42 transcription factor genes, including Homeodomain, MCMI-AGAMOUS-DEFICIENS-SRF4 (MADS), and WRKY domains, has been identified in a transcriptome analysis of all cucumber trichomerelated mutants (Zhao et al., 2015a). These transcription factor genes are unique to plants and are involved in a range of activities; many are associated with multicellular trichrome development and other species-specific development processes. The study of the interaction among those transcription factors is an emerging area of research because these factors probably share many biological functions for trichome development. All of this information warrants further investigation.

The aim of this review is to provide readers with a summary of the progress of cucumber trichome development and to

#### REFERENCES

encourage plant scientists to further investigate the mechanism of trichrome initiation and development and their regulatory network. We have summarized the mutants related to fruit trichomes, the key genes that control trichome initiation and development, gene relationships and a possible model for fruit trichome positive regulatory mechanisms that is different from Arabidopsis in core transcriptional factor numbers. All of this may help us to better understand the advances in the study of cucumber trichomes. At the same time, more investigations must be conducted.

Because cucumber is a horticultural crop of worldwide importance and fruit spines directly affect its commercial quality, an extensive characterization of cucumber trichomes will not only help us to understand the underlying molecular mechanisms involved in multicellular trichome development but will also pave the way for creating new cucumber varieties with desired trichome growth and density. Moreover, cucumber trichomes may serve as a model system for studying the development of multicellular trichomes.

#### AUTHOR CONTRIBUTIONS

HR and XL organized the review, wrote the first draft and generated **Figures 1**–**4**. EB and YC contributed to a second draft. All of the authors revised the manuscript multiple times. HR and XL performed the final revision of the manuscript, which was read and approved by all authors.

#### ACKNOWLEDGMENTS

We thank the members of the Ren lab for helpful discussions for this review. This work was supported by National Research and Development Program (2016YFD0101705), Beijing Agricultural Innovation Consortium (BAIC01-2016) and Beijing Agricultural Scientific and Technological Program (20160415) to HR.

TRANSPARENT TESTA GLABRA1 (TTG1) gene. Plant Mol. Biol. 57, 67–81. doi: 10.1007/s11103-004-6768-1

**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 © 2016 Liu, Bartholomew, Cai and Ren. 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) or licensor 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.

# Interactive Effects of Elevated [CO2] and Water Stress on Physiological Traits and Gene Expression during Vegetative Growth in Four Durum Wheat Genotypes

#### Susan Medina1, 2, Rubén Vicente1 \*, Amaya Amador 3 and José Luis Araus 1

1 Integrative Crop Ecophysiology Group, Plant Physiology Section, Faculty of Biology, University of Barcelona, Barcelona, Spain, 2 Crop Physiology Laboratory, International Crops Research Institute for Semi-Arid Tropics, Patancheru, India, 3 Unitat de Genòmica, Centres Científics i Tecnològics, Universitat de Barcelona, Barcelona, Spain

#### Edited by:

Paul Christiaan Struik, Wageningen University and Research Centre, Netherlands

#### Reviewed by:

Iker Aranjuelo, Instituto de Agrobiotecnología (CSIC-UPNA), Spain Fulai Liu, University of Copenhagen, Denmark Salvatore Ceccarelli, Rete Semi Rurali, Italy

\*Correspondence:

Rubén Vicente<EMAIL_ADDRESS>
#### Specialty section:

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

Received: 06 June 2016 Accepted: 04 November 2016 Published: 22 November 2016

#### Citation:

Medina S, Vicente R, Amador A and Araus JL (2016) Interactive Effects of Elevated [CO2] and Water Stress on Physiological Traits and Gene Expression during Vegetative Growth in Four Durum Wheat Genotypes. Front. Plant Sci. 7:1738. doi: 10.3389/fpls.2016.01738 The interaction of elevated [CO2] and water stress will have an effect on the adaptation of durum wheat to future climate scenarios. For the Mediterranean basin these scenarios include the rising occurrence of water stress during the first part of the crop cycle. In this study, we evaluated the interactive effects of elevated [CO2] and moderate to severe water stress during the first part of the growth cycle on physiological traits and gene expression in four modern durum wheat genotypes. Physiological data showed that elevated [CO2] promoted plant growth but reduced N content. This was related to a down-regulation of Rubisco and N assimilation genes and up-regulation of genes that take part in C-N remobilization, which might suggest a higher N efficiency. Water restriction limited the stimulation of plant biomass under elevated [CO2], especially at severe water stress, while stomatal conductance and carbon isotope signature revealed a water saving strategy. Transcript profiles under water stress suggested an inhibition of primary C fixation and N assimilation. Nevertheless, the interactive effects of elevated [CO2] and water stress depended on the genotype and the severity of the water stress, especially for the expression of drought stress-responsive genes such as dehydrins, catalase, and superoxide dismutase. The network analysis of physiological traits and transcript levels showed coordinated shifts between both categories of parameters and between C and N metabolism at the transcript level, indicating potential genes and traits that could be used as markers for early vigor in durum wheat under future climate change scenarios. Overall the results showed that greater plant growth was linked to an increase in N content and expression of N metabolism-related genes and down-regulation of genes related to the antioxidant system. The combination of elevated [CO2] and severe water stress was highly dependent on the genotypic variability, suggesting specific genotypic adaptation strategies to environmental conditions.

Keywords: climate change, durum wheat, elevated [CO2 ], genotypic variability, stable isotopes, transcript levels, vegetative growth, water stress

## INTRODUCTION

Food security is facing new challenges nowadays due to the increase in the world population and the impacts of climate change on agriculture and food supply. Wheat is a very important crop for the human diet, ranking in fourth position in terms of the world's most important crops by production quantity after sugarcane, maize, and rice (FAO, 2013). Although bread wheat dominates global wheat production, durum wheat is an economically and culturally important staple crop in the Mediterranean region, used for the production of pasta, bread, burghul, couscous, and freekeh (Habash et al., 2009). In the second half of the twentieth century, local durum wheat landraces were replaced by improved semi-dwarf cultivars, which showed higher yield and harvest index (Soriano et al., 2016). In the early 1970s, introduction of germplasm from CIMMYT (International Maize and Wheat Improvement Centre) increased grain yield (Sanchez-Garcia et al., 2013). Improvement in wheat yield per unit area constitutes one of the largest challenges to be addressed by breeding programs, covering numerous research areas (McKersie, 2015). Projections of wheat production assume that the growth rate will be lower than the historical growth rates reported in the second half of the twentieth century (Bort et al., 2014; Nakhforoosh et al., 2015), with insignificantly higher yields in modern wheat genotypes released in recent years (Sanchez-Garcia et al., 2013). It is unlikely that any improvements will support the increase in world population or mitigate against future extreme weather events (Araus et al., 2002; Alexandratos and Bruinsma, 2012; Trnka et al., 2014).

Observations of the climate system confirm that Earth's mean surface temperature is increasing rapidly as a consequence of the anthropogenic emissions of CO2 and other greenhouse gases (IPCC, 2013). The atmospheric concentration of CO2 ([CO2]) has increased by more than 40% since the beginning of the industrial revolution and is expected to double by the end of this century (IPCC, 2013). As atmospheric [CO2] is currently a limiting factor for C3 photosynthesis, the primary effect of a short-term exposure to elevated [CO2] includes an initial stimulation of photosynthesis due to both enrichment of substrate for ribulose bisphosphate carboxylase oxygenase (Rubisco) carboxylation and inhibition of competitive Rubisco oxygenation which may eventually contribute to a higher biomass (Stitt and Krapp, 1999; Long et al., 2006). High [CO2] also induces a stomatal closure leading to a better leaf water status. However, growth over the long-term under elevated [CO2] leads to a down-regulation of photosynthetic capacity, which has been related to a decline in Rubisco protein content and activity, together with a higher carbohydrate accumulation and a decline in N concentration and protein content in wheat (Aranjuelo et al., 2011, 2013; Vicente et al., 2015a,b). This phenomenon suggests that regulatory mechanisms may occur in the plant, e.g., end-product inhibition, carbon sink limitation, biomass dilution effects, or a decline in nutrient uptake and/or assimilation (Stitt and Krapp, 1999; Vicente et al., 2015a). Moreover, elevated [CO2] leads to an altered expression pattern of genes involved in the photosynthetic apparatus, the distribution of C, respiration, and N metabolism in durum wheat (Vicente et al., 2015b).

Increasing greenhouse gas emissions may cause further warming together with rainfall reduction in the next decades, which will increase the frequency and intensity of drought in the Mediterranean basin (Habash et al., 2009; IPCC, 2013; McKersie, 2015). For the Iberian Peninsula it is predicted that drought stress can occur at any growth stage of wheat (Russo et al., 2015), with the grain-filling phase being the most studied. However, the number of studies focusing on drought stress during early growth is limited. Although rainfall has been traditionally most abundant and evapotranspiration the lowest during winter, the occurrence of drought in winter months during the early stages of the crop cycle has been reported in recent times (Russo et al., 2015). This can further constrain wheat growth and thus final grain yield, mostly through a decrease in the ear density and number of kernels per unit crop area (Araus et al., 2008; Rebolledo et al., 2013). In addition, a constitutive (i.e., in absence of water stress) rapid development of wheat plants (early vigor) could be a positive trait and relevant for further avoiding drought stressrelated consequences at both early and late growth stages. Early vigor could benefit plant growth and yield by increasing resource acquisition, shading the soil, preventing evaporation from it, and suppressing weeds (Maydup et al., 2012; Bort et al., 2014; Pang et al., 2014). As a consequence, differences in early growth (tillering and further stem elongation) will affect the number of fertile stems (and thus the ear density) and the size of the ears (and thus the potential number of grains per ear), which are the main contributors determining grain yield (Guo et al., 2016).

Plant responses to water stress define a complex and sophisticated regulatory network comprising physiological, biochemical, and molecular mechanisms. In wheat, some of these responses include inhibition of plant growth and photosynthetic capacity, together with a wide range of physiological responses, including changes in stomatal closure and decreases in transpiration, Rubisco efficiency, and chlorophyll content as well as an increase in oxidative stress among other responses (Budak et al., 2013; Nezhadahmadi et al., 2013). Such responses are modulated by stress severity. Cessation of watering showed a progressive reduction in leaf relative water content, water potential and photosynthesis in durum wheat (Habash et al., 2014). Liu et al. (2016) reported a progressive inhibition of photosynthetic activity as water stress is more severe in fieldgrown bread wheat, probably due to non-stomatal limitations, which led to lower grain yields even at moderate water stress. Furthermore, water stress in wheat leads to complex changes in the expression of some genes, including those involved in photosynthesis, respiration, N metabolism, lipid metabolism, transcription factors, signal transducers, and synthesis of protective proteins (Habash et al., 2009, 2014; Budak et al., 2013; Yousfi et al., 2016). These changes in gene expression occurred mainly in the early phases of the stress (Habash et al., 2014).

Plant responses to elevated [CO2] or water stress are influenced by the duration and level of the environmental factor, the growth stage, and the genetic variability. Studies carried out with different durum wheat genotypes demonstrated that the responsiveness to elevated [CO2] (Aranjuelo et al., 2013), water stress (De Leonardis et al., 2007; Aprile et al., 2013; Habash et al., 2014), and the combination of both (Erice et al., 2014) is genotype specific. Moreover, the growth stage greatly influences the response of durum wheat to elevated [CO2] (Aranjuelo et al., 2011; Vicente et al., 2015a) and drought (Liu et al., 2016). In addition, the interactive effects of environmental conditions and genotypic variability cannot be anticipated from the individual effects of these treatments (Ceccarelli et al., 1991). Some studies have shown positive effects of elevated [CO2] on water stress tolerance of different bread wheat varieties (Harnos et al., 2002; Wall et al., 2006; Robredo et al., 2011; Bencze et al., 2014). A positive synergistic effect of elevated [CO2] and water stress has been reported to decrease gs , and thus leads to an improvement in water use efficiency at the stomatal and whole plant level (Bencze et al., 2014; Pazzagli et al., 2016). The decrease in photosynthesis under water stress is often mitigated by elevated [CO2] (Bencze et al., 2014), resulting in increased levels of carbohydrates for the development of new tissues or filling grain (Wall et al., 2006). However, such positive effects of elevated [CO2] in improving stress tolerance are not always achieved (Hudak et al., 1999; Pleijel et al., 2000). Bencze et al. (2014) reported that drought at elevated [CO2] led to a stimulation of the antioxidant enzyme system in bread wheat, which suggests a high level of oxidative stress. Erice et al. (2014) showed that the stimulation of plant growth by elevated [CO2] was only found in durum wheat genotypes with high harvest indices and optimal water supply. Therefore, additional efforts are still necessary to deepen our understanding of the interactive effect of [CO2] and water regime in durum wheat.

The aim of this work was to determine the physiological and molecular mechanisms involved in the adaptive response of four semi-dwarf (i.e., post-Green Revolution) durum wheat cultivars to different [CO2] and water regimes. Durum wheat genotypes were grown under controlled conditions at ambient and elevated [CO2] and two different water regimes (fully irrigated and moderate/severe water stress). We assessed plant growth, physiological traits, stable C and N isotopic signatures, and transcript levels for stress-responsive genes that could be good indicators of durum wheat's adaptation to future climate conditions at vegetative growth stages. The genes selected corresponded to key enzymes in the metabolism of C (the Rubisco large and small subunits, RBCL and RBCS, respectively, and phosphoenolpyruvate carboxylase, PEPC) and N (the cytosolic and plastidial glutamine synthetases, GS1 and GS2, respectively), as well as proteins involved in stress responses (dehydrins 11, DHN11, and 16, DHN16, catalase, CAT, and superoxide dismutase, SOD). Rubisco is the key enzyme for photosynthetic CO2 assimilation, and its activity is highly responsive to atmospheric [CO2] (Vicente et al., 2011; Carmo-Silva et al., 2015). PEPC is a cytosolic enzyme that catalyzes the βcarboxylation of phosphoenolpyruvate to produce oxaloacetate, which is involved in anaplerotic functions. GS1 and GS2 play a central role in N metabolism: the former is thought to be involved in the primary assimilation of ammonium from nitrate reduction and photorespiration, while the latter is mainly involved in the transport of N through the plant and N recycling from catabolic processes. The function of the dehydrin family is not completely understood, but these proteins are involved in conferring stress tolerance (Kosová et al., 2014). Catalases and superoxide dismutases are primary antioxidant enzymes involved in the elimination of reactive oxygen species (ROS) such as the cytotoxic H2O2 produced by photorespiration (Luna et al., 2005) and the superoxide generated during photosynthetic electron transport (Xu et al., 2010; Huseynova et al., 2014). Thus, our study combines the effects of genotypic variability and future environmental conditions, integrating plant performance with gene expression, and aims to identify traits associated with better performance during vegetative growth.

#### MATERIALS AND METHODS

#### Plant Material and Growth Conditions

The experiment was conducted with four semi-dwarf durum wheat [Triticum turgidum L. ssp. durum (Desf.)] genotypes: Mexa (year of commercial release: 1977), Regallo (1988), Burgos (1997) and Ramirez (2006). These cultivars represent high-yield genotypes released in the last forty years that are (or were) widely cultivated in the Mediterranean regions of Spain. The study of these genotypes could provide information about the adaptation of modern cultivars to climate change and whether there are differences between them associated with the year they were released. The experiment was conducted from May to July 2015 in two controlled environment chambers (Conviron E15; Controlled Environments, Winnipeg, MB, Canada) in the Experimental Facilities of the Faculty of Biology at the University of Barcelona. A total of 96 durum wheat plants (24 for each genotype) were sown in 2 L pots containing a mixture of standard substrate:perlite (1:1, v/v) and were grown with a long light period of 16 h, a photosynthetic photon flux density (PPFD) of 350 µmol m−2 s −1 , a day/night temperature of 23/17◦C and a relative humidity of 60%. During the entire experiment, half of the pots were cultivated under atmospheric [CO2] (400 µmol mol−1 ) in one chamber, while the other half grew under elevated [CO2] (790 µmol mol−1 ) in the other chamber with injection of CO2 from an external bottle (Carburos Metálicos S.A., Barcelona, Spain). The temperature, relative humidity and [CO2] within each chamber were continuously monitored by Conviron series controllers (CMP3243 Controlled Environments Ltd., Winnipeg, MB, Canada). The technical staff of the Experimental Facilities of the Faculty of Biology tested the growth conditions of each chamber periodically with external sensors: an HMP75 humidity and temperature probe and a GMP222 CO2 probe for use with an MI70 series hand-held indicator (Vaisala, Vantaa, Finland). Similarly, the PPFD was periodically verified with an LI-188B quantum/radiometer/photometer (LI-COR Inc., Lincoln, NB, USA).

The plants were uniformly irrigated every 2 days with 50% Hoagland's nutrient solution over a 25 day period. After that (Zadoks 21), the water stress was imposed; one half of the plants of each genotype and [CO2] were maintained under wellwatered conditions (100% pot capacity, PC) until the end of the experiment, while the other half were subjected to water stress conditions. The maximum soil volumetric water content of each pot was evaluated at the beginning of the experiment as the difference between pot weight after watering with the excess water drained and the pot dry weight. Thus, pots were watered by direct measurements of the pot weight and the water supply was adjusted to the pot water conditions established for each water regime. In the water-stressed plants the watering was progressively restricted by 10% PC every 2 days. First, after 8 days the water-stressed plants received a 60% PC (moderate water stress) and this irrigation regime was strictly maintained for 10 days (see **Figure 1** for a schematic representation of the experimental design). At the end of this period (Zadoks 26), equal numbers (48) of well-watered and water-stressed plants were sampled. The youngest fully expanded leaf was collected, immediately frozen in liquid nitrogen and stored at −80◦C for gene expression, C and N content and stable isotope analyses. After that, the whole plant was harvested and dried in an oven at 60◦C for 72 h for biomass analysis. Second, in the remaining half of the plants (48), the progressive water limitation continued for 8 more days until water-stressed plants received a 30% PC (severe water stress). As in the moderate water stress, the irrigation conditions in well-watered and water-stressed plants were maintained for 10 days. Later, these 51-day-old plants (Zadoks 28–32) were collected following the procedure described above. The moderate and severe water stresses were defined in this experiment based on similar reductions in irrigation and stomatal conductance used in other studies (Galmes et al., 2007; Liu et al., 2016). The pots were rotated three times a week to avoid edge effects in the growth chambers over the course of the experiment. We used a rotatory randomized complete block design with three replicates (one plant per pot) per factor combination ([CO2], water level and genotype) at each sampling.

#### Physiological Traits

Prior to harvest a hand-held portable spectroradiometer (GreenSeeker, NTech Industries, Ukiah, CA, USA) was used to estimate the normalized difference vegetation index (NDVI) of each plant (only at the second sampling date). Relative chlorophyll content was measured with a Minolta SPAD-502 chlorophyll meter (Spectrum Technologies, Plainfield, IL, USA). Stomatal conductance (gs) was measured using a Decagon SC-1 Leaf Porometer (Decagon Device, Inc., Pullman, WA, USA). Both chlorophyll content and stomatal conductance of the adaxial surface were recorded in the central segment of the same youngest fully expanded leaf between 3 and 5 h after the start of the photoperiod. In addition, plants were collected to determine the leaf, shoot, root, and plant dry weights as indicated above, while the roots were washed in tap water until all substrate was removed. The number of tillers and the root to shoot dry weight ratio (root/shoot) were then determined.

## C and N Content and Stable Isotope Signatures

A fraction of the youngest fully expanded leaf was finely powdered and then 1 mg of this leaf material was used for the measurements of total C and N content (as a percentage of

leaf dry weight) and the stable C (13C/12C) and N (15N/14N) isotope ratios. Measurements were carried out using an elemental analyzer (Flash 1112 EA; ThermoFinnigan, Bremen, Germany) coupled with an isotope ratio mass spectrometer (Delta C IRMS; ThermoFinnigan), operating in continuous flow mode, at the Scientific Facilities of the University of Barcelona. As has been described previously (Bort et al., 2014; Yousfi et al., 2016), the 13C/12C ratio was expressed in δ notation: δ 13C (h) = [(13C/12C)sample/(13C/12C)standard − 1] × 1000. The standard refers to international secondary standards of known 13C/12C ratios (IAEA CH7 polyethylene foil, IAEA CH6 sucrose, and USGS 40 L-glutamic acid) calibrated against Vienna Pee Dee Belemnite calcium carbonate. The same δ notation was used for the 15N/14N ratio (δ 15N) using N2 in air as standard.

#### Quantitative Reverse Transcriptase PCR Amplification

Frozen leaf samples were ground with liquid nitrogen and subsequently RNA was isolated from 100 mg of this material with Ribozol RNA Extraction Reagents (Amresco, Solon, OH, USA) according to the manufacturer's instructions. RNA quantity and quality was measured using a NanoDrop ND-1000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). RNA integrity was checked by 1.5% (w/v) agarose gel electrophoresis. Total RNA (1 µg) was treated with PerfeCTa DNase I RNase-free (Quanta Biosciences, Gaithersburg, MD, USA) to eliminate residual genomic DNA. cDNA was synthesized using a qScript cDNA Synthesis Kit (Quanta Biosciences) following the manufacturer's instructions. The qRT-PCR assays were performed in optical 384-well-plates with the LightCycler 480 System (Roche Applied Science, Penzberg, Germany) in the Centres Científics i Tecnològics de la Universitat de Barcelona (CCiTUB), in a reaction volume of 10 µL: 5 µL of PerfeCTa SYBR Green FastMix (Quanta Biosciences), 200 nM of each gene-specific primer and 1 µL of diluted cDNA (1:10). The thermal profile was as follows: initial denaturation for 30 s at 95◦C, PCR cycling (45 cycles) for 5 s at 95◦C, 15 s at 60◦C, and 10 s at 72◦C, and a final step of 95◦C for 5 s and 60◦C for 60 s to obtain the dissociation curve. Two technical replicates were analyzed per biological replicate. Specific primers for genes encoding the Rubisco large subunit (NC\_021762), phosphoenolpyruvate carboxylase (Y15897), plastidial glutamine synthetase (DQ124212), dehydrin 11 (AJ890140), and superoxide dismutase (KP696754) were designed in Primer-BLAST (http:// www.ncbi.nlm.nih.gov/tools/primer-blast/) using the following criteria: Tm = 60 ± 1 ◦C, primer length of 18–25 bases, GC content of 30–70% and product size of 60–150 bases. The specificity of PCR amplification was confirmed by the presence of unique amplicons of the expected length on 3.5% (w/v) agarose gels. The genes encoding the ADP-ribosylation factor and the RNase L inhibitor-like protein, previously identified as potential reference genes (Vicente et al., 2015b), were used to normalize qRT-PCR data after the evaluation of their expression stability in this study. All primers used for gene expression analysis and their symbols are listed in Supplementary Table S1. The values of the cycle threshold (Ct) were calculated using the LightCycler

1.5 software (Roche Applied Science). The quantification of the relative gene expression was analyzed using the comparative Ct method 2−11Ct (Schmittgen and Livak, 2008), and the data were presented as the log2 fold change.

#### Data Analysis

The effects of [CO2] (ambient and elevated), water regime (wellwatered and water stressed), genotype (Mexa, Regallo, Burgos, and Ramirez), and their interaction on plant growth, chlorophyll content, gs , and C and N contents and isotope composition were determined through a three-factor (2 CO2 × 2 water regimes × 4 genotypes) analysis of variance (ANOVA) for each sampling date (moderate and severe waters stress; see Supplementary Table S2) with GenStat 6.2 (VSN International Ltd, Hemel Hempstead, UK). Further, and given the implicit complexity of the design, each genotype was analyzed through a two-factor ANOVA (2 CO2 × 2 water regimes) for both sampling dates. All factors were treated as fixed independent variables. When the F-ratio was significant (P < 0.05), the least significant difference (LSD) test was used to assess differences between treatment means. Clustered heat maps of relative gene expression were built in the R statistics environment (R Development Core Team, 2008) to study the effects of elevated [CO2] and water stress on transcript levels. A correlation matrix was generated in R for evaluating the relationships between all parameters analyzed. Visualization of significant correlations was performed using Cytoscape software (Shannon et al., 2003).

#### RESULTS

#### Effect of [CO2], Water Regime, and Genotype on Plant Growth

Total biomass of the plant and its different fractions (leaves, shoot and root), the root/shoot ratio, and the number of tillers were analyzed through two-factor ANOVA ([CO2] × water regime) for each genotype (**Tables 1**, **2**). Moderate and severe water stress were established with reductions of 40 and 70% in the water supplied to the pots and average decreases of 34 and 57% in gs , respectively, compared to well-watered plants (data not shown). Growth under elevated [CO2] led to significant increases in biomass compared to ambient [CO2] (**Tables 1, 2**). At the first sampling date, elevated [CO2] increased root biomass in Mexa, Regallo and Ramirez (and also in Burgos, P = 0.074), but only increased plant biomass in Regallo. The root/shoot ratio also increased in Regallo and Burgos under elevated [CO2]. At the second sampling date, elevated [CO2] increased plant biomass due to higher shoot and root biomass compared to ambient [CO2], with larger increases in plant biomass in Mexa and Regallo under well-watered conditions in comparison to water stressed conditions. As a consequence of the increases in both shoot and root dry weights by elevated [CO2], the root/shoot ratio was not altered, except in Regallo. Moreover, in this genotype an increase in the tillers per plant was also observed under elevated [CO2] but only in well-watered conditions.

Moderate water stress did not lead to statistical differences in biomass, the root/shoot ratio, or the number of tillers between

TABLE 1 | Total leaf (LDW), shoot (SDW), root (RDW) and plant (PDW) dry weight, root/shoot ratio, and number of tillers per plant in four durum wheat genotypes grown under ambient or elevated [CO2 ] and well-watered or moderate water stress conditions (100 vs. 60% pot capacity).

Significant effects for elevated [CO2] (C), water stress (W) and their interaction (C × W) were determined by two-factor ANOVA (P). Values with the same letter are not significantly different for the interaction [CO2] × water level. Significant P values are marked in bold (P < 0.05).

well-watered and water-stressed plants (**Table 1**). However, severe water stress led to significant changes in these parameters, while NDVI was also affected (**Table 2**). Plant biomass generally decreased under severe water stress compared to well-watered conditions and was associated with decreases in leaf, shoot, and root dry weights. Water restriction decreased the number of tillers per plant in Burgos under severe water stress, while this reduction was not significant in the other genotypes. Additionally, the NDVI values were lower in water-stressed plants compared to well-watered plants, irrespective of the [CO2] and the genotype (**Table 2**). In general, at severe water stress the interaction [CO2] × water regime × genotype showed that the root/shoot ratio strongly increased in Regallo, especially under ambient [CO2] and well-watered conditions (Supplementary Table S2). The Burgos and Mexa cultivars had higher shoot dry weight than Ramirez and Regallo, while root dry weight was higher in Regallo (Supplementary Table S2). Furthermore, significant [CO2] × genotype interaction showed that Burgos and Regallo under elevated [CO2] increased tiller production, whereas Ramirez and Regallo plants under ambient [CO2] had lower tillering (Supplementary Table S2).

#### Effect of [CO2], Water Regime, and Genotype on Chlorophyll Content, gs, C and N Content and C and N Isotope Composition

The interactive effects of [CO2] and water regime on chlorophyll content, gs , and C and N contents and isotope composition were analyzed in the youngest fully expanded leaf through twofactor ANOVA for each genotype during vegetative growth under moderate (**Table 3**) and severe water stress (**Table 4**). At moderate water stress, elevated [CO2] compared to ambient [CO2] decreased N content in Mexa and Regallo, and δ 13C regardless of the genotype, while it increased chlorophyll content in Mexa, gs , and C content in Regallo, and δ 15N in all genotypes except in Regallo. Water stress reduced gs and increased

TABLE 2 | Total leaf (LDW), shoot (SDW), root (RDW) and plant (PDW) dry weight, root/shoot ratio, number of tillers per plant, and normalized difference vegetation index (NDVI) in four durum wheat genotypes grown under ambient or elevated [CO2 ] and well-watered or severe water stress conditions (100 vs. 30% pot capacity).

Significant effects for elevated [CO2] (C), water stress (W) and their interaction (C × W) were determined by two-factor ANOVA (P). Values with the same letter are not significantly different for the interaction [CO2] × water level. Significant P values are marked in bold (P < 0.05).

δ 15N in Regallo, and increased δ 13C in Burgos and Ramirez. Three-factor ANOVA showed significant interactions for δ 15N (Supplementary Table S2). The [CO2] × genotype interaction mainly showed that δ 15N was higher in Ramirez and Burgos at elevated [CO2] and in Regallo at both [CO2], whereas the lowest values were observed in Ramirez at ambient [CO2]. The water regime × genotype interaction indicated that δ 15N was higher in Mexa and Regallo under water stress than in the other genotypes.

At severe water stress, elevated [CO2] relative to ambient [CO2] decreased the N content in Regallo and Burgos and δ 13C regardless of the genotype, while it increased δ 15N in Burgos (**Table 4**). Furthermore, gs in Burgos and N content in Regallo decreased under severe water stress compared to well-watered conditions. In addition, under well-watered conditions gs was higher in Burgos than in other genotypes (Supplementary Table S2). Chlorophyll content was lower in Regallo, whereas δ 13C was higher in Burgos, as compared to other genotypes (Supplementary Table S2).

#### Effect of [CO2] and Water Regime on Gene Expression for Each Durum Wheat Genotype

Treatment effects on transcript levels were evaluated for each genotype using nine genes that encode enzymes of primary C and N metabolism and stress-responsive proteins (Supplementary Table S1). Elevated [CO2] and water stress led to changes in

TABLE 3 | Chlorophyll content, stomatal conductance (gs), N and C content, and N and C isotope composition (δ 15N and δ 13C, respectively) in four durum wheat genotypes grown under ambient or elevated [CO2 ] and well-watered or moderate water stress conditions (100 vs. 60% pot capacity).

Significant effects for elevated [CO2] (C), water stress (W) and their interaction (C × W) were determined by two-factor ANOVA (P). Significant P values are marked in bold (P < 0.05).

gene expression depending on genotype and the level of water restriction (**Table 5**; Supplementary Figure S1). At moderate water stress, elevated [CO2] decreased transcript levels of RBCL, RBCS, and GS2 relative to control conditions (ambient [CO2] and well-watered conditions), particularly in the Mexa and Regallo genotypes. Under water stress the transcript levels for these enzymes markedly increased in Ramirez. Elevated [CO2] caused a generalized increase in the transcript levels of PEPC and GS1, particularly when it was combined with moderate water stress. Transcript abundances of the dehydrins, DHN11 and DHN16, were generally higher under elevated [CO2] and well-watered conditions, but lower under ambient [CO2] and water stress relative to control conditions. However, DHN11 and DHN16 showed opposite expression patterns under elevated [CO2] and water stress. [CO2] enrichment and moderate water stress decreased transcript levels of CAT and SOD in Mexa, Regallo, and Burgos compared with control conditions, whereas in Ramirez they did not change significantly.

Gene expression analysis indicated greater genotype-specific differences under severe water stress than under moderate water stress (**Table 5**; Supplementary Figure S1). In Mexa under elevated [CO2] and well-watered conditions there were higher transcript levels of RBCL, RBCS, PEPC, GS1, GS2, and CAT and lower levels of DHN16, relative to control conditions. Severe water stress did not substantially alter gene expression. In Regallo most of the transcripts studied were lower in all treatment combinations than in control conditions. However, DHN16 and SOD transcripts increased under ambient [CO2] and water stress, and these together with CAT and GS1 also increased under elevated [CO2] and water stress. In the case of Burgos, elevated [CO2], water stress and their combination strongly reduced transcript levels in comparison to control conditions, especially for GS1 and DHN16, while SOD transcripts increased under water stress and elevated [CO2] × water stress as observed in Regallo. In Ramirez elevated [CO2] led to a reduction in the transcript levels of RBCL, RBCS, and SOD and an increase in PEPC and GS1 compared to control conditions. Water stress increased PEPC and DHN16 transcript levels relative to control conditions, while under the combination of elevated [CO2] and water stress

TABLE 4 | Chlorophyll content, stomatal conductance (gs), N and C content, and N and C isotope composition (δ 15N and δ 13C, respectively) in four durum wheat genotypes grown under ambient or elevated [CO2 ] and well-watered or severe water stress conditions (100 vs. 30% pot capacity).

Significant effects for elevated [CO2] (C), water stress (W) and their interaction (C × W) were determined by two-factor ANOVA (P). Significant P values are marked in bold (P < 0.05).

greater transcript abundances were observed for most of the genes.

### Correlation Network of Physiological Traits and Gene Expression

A Pearson correlation matrix was generated using the mean values for each treatment combination, genotype and sampling date (n = 32) of the physiological traits and transcript levels (Supplementary Table S3), excluding NDVI, which was only measured at severe water stress, and δ 13C, which was influenced by C composition of the CO2 bottles used in the elevated [CO2] chamber (Aljazairi et al., 2015). Of the 190 correlations between parameters, there were 28 positive and 19 negative significant correlations (P < 0.05) that are represented in an association network (**Figure 2**). Most of the significant correlations were observed between physiological traits and transcript levels independently. Positive correlations were found among leaf, shoot, and plant dry weights, between the leaf and shoot dry weights with the number of tillers, and between root and plant dry weights. The root/shoot ratio was positively correlated with root dry weight and negatively correlated with leaf and shoot dry weights, the number of tillers and N content. Furthermore, δ 15N was also negatively correlated with N content and the number of tillers. Chlorophyll content was correlated positively with leaf, shoot, root, and plant dry weights, and negatively with gs . On the other hand, positive correlations were found between N content with leaf and shoot dry weights and the number of tillers, and negative correlations between N content with root dry weight, and between gs with root and plant dry weights. In the case of transcript levels, RBCL was correlated with RBCS, GS1, GS2, and PEPC, whereas RBCS correlated with GS2 and DHN11, GS2 with DHN11, and PEPC with GS1. Furthermore, some relationships were found between physiological traits and gene expression (**Figure 2**). Positive correlations appeared between DHN16 with plant biomass (leaf, shoot, root, and plant dry weights), CAT with gs and SOD with C content. Moreover, negative correlations were found between chlorophyll content with RBCL, RBCS, GS2, and CAT, also between root dry weight with RBCL, RBCS, GS2, and DHN11, and finally plant dry weight with RBCL.

TABLE 5 | Transcript changes in four durum wheat genotypes grown under ambient or elevated [CO2 ] and well-watered or water stressed conditions: (A) moderate and (B) severe water stress.

White indicates no change, blue up-regulation, and red down-regulation in each treatment relative to the treatment under ambient [CO2] and optimal water supply for each genotype, as shown in the color bar for a log2 scale. RBCL, Rubisco large subunit; RBCS, Rubisco small subunit; PEPC, phosphoenolpyruvate carboxylase; GS1, cytosolic glutamine synthetase; GS2, plastidial glutamine synthetase; DHN11, dehydrin 11; DHN16, dehydrin 16; CAT, catalase; SOD, superoxide dismutase.

−3 0 3

## DISCUSSION

Although, substantial efforts have been made in recent years to identify traits associated with wheat performance during early growth (Maydup et al., 2012; Rebolledo et al., 2013; Bort et al., 2014; Pang et al., 2014; Wilson et al., 2015), little attention has been paid to the effect of interactions between elevated [CO2] and water stress in durum wheat. The effects of water restriction on crop growth have been mostly studied with the view of improving drought impacts at late growth stages in Mediterranean environments. However, projections of future climate change in the Iberian Peninsula predict major rainfall limitations and higher evapotranspiration during winter months (Russo et al., 2015) and therefore early-season drought is a matter of concern. In this context, we describe the effects of elevated [CO2] and water stress during the first part of the growth cycle in four durum wheat genotypes on physiological traits and expression of nine genes that respond to changes in [CO2] and water levels (Ali-Benali et al., 2005; Budak et al., 2013; Vicente et al., 2015b; Yousfi et al., 2016). The coordination of these parameters under the different combinations of factors is discussed.

## Changes in Physiological Traits of Durum Wheat Genotypes under Different Water Regimes and [CO2] Levels

A moderate water stress in 43-day-old plants did not significantly alter plant growth (**Table 1**). Long-term exposure to elevated [CO2] led to higher root biomass relative to ambient [CO2] independently of genotypic variability, in concordance with reports from other crop species (Madhu and Hatfield, 2013). This increment was associated with higher plant growth in Regallo and higher root/shoot ratios in Regallo, Burgos, and Mexa (**Table 1**). In fact, under elevated [CO2] root growth is often more stimulated than the aerial part of the plant, although it depends on genotype × environment variation (Stitt and Krapp, 1999; Madhu and Hatfield, 2013). A severe water stress

in 51-day-old plants showed greater effects on plant growth than moderate water stress (**Table 2**). [CO2] enrichment generally led to an increase in plant biomass by increasing root and shoot biomass and tillering, particularly under optimal water supply. This could be due to the effects of [CO2] fertilization on the net photosynthetic rate (Long et al., 2006; Vicente et al., 2015b), especially in genotypes with large harvest indices such as post-Green Revolution cultivars (Aranjuelo et al., 2013). It could also be caused by carbohydrate accumulation, which may lead to increases in the number of tillers (Stitt and Krapp, 1999). On the other hand, severe water stress constrained plant growth (dry matter and NDVI), in agreement with earlier studies in durum wheat (Erice et al., 2014; Nakhforoosh et al., 2015; Yousfi et al., 2016), with Ramirez and Burgos being the genotypes most affected. According to Marti et al. (2007), we suggest that progressive water restriction during the vegetative stage constrained the photosynthetic area, which may cause negative effects on final biomass and yield.

Chlorophyll content, gs , and N and C contents and isotope compositions at moderate and severe water stress did not reveal statistical significance for the interactions [CO2] × water regime and [CO2] × water regime × genotype (**Tables 3, 4**; Supplementary Table S2). Stomatal conductance (gs) generally decreases under elevated [CO2] and drought stress due to an increase in internal [CO2] and as a water saving strategy, respectively (Long et al., 2006; Nakhforoosh et al., 2015; Vicente et al., 2015b; Pazzagli et al., 2016). The average gs values decreased under water restriction at moderate and severe water stress, but it was only significantly decreased in some genotypes (**Tables 3, 4**). On the other hand, elevated [CO2] did not alter gs at this growth stage, except for an increase in gs under moderate water stress in Regallo, which could favor CO2 assimilation and consequently biomass accumulation under this water regime (**Tables 1, 3**). Earlier studies have shown a decrease in gs under water stress (Peremarti et al., 2014; Pazzagli et al., 2016), while negligible changes have been reported under elevated [CO2] in tomato and durum wheat, and increases have even been recorded for Regallo (Vicente et al., 2015a; Pazzagli et al., 2016). Therefore, the growth stage and the severity of the water stress influenced stomatal closure, while elevated [CO2] had minor effects on gs during vegetative growth.

Elevated [CO2] generally decreased N content in the present study (**Tables 3, 4**), which has been observed in C3 plants through shifts in N uptake and/or assimilation (which agrees with the changes in transcript levels of N-metabolism enzymes; see below) together with other uncertain mechanisms, e.g., the biomass dilution effect, increased N loss, and sink limitation (Stitt and Krapp, 1999; Aranjuelo et al., 2011; Vicente et al., 2015a,b). N content was also diminished by severe water stress in Regallo, in agreement with previous studies in durum wheat (Yousfi et al., 2012, 2016). Chlorophyll content only increased under elevated [CO2] in Mexa at the first sampling date, but the effect disappeared at the second sampling (**Table 3**). [CO2] enrichment and water stress did not modify C content in leaves, suggesting that the decrease in N content was not simply due to N dilution caused by rapid growth (Taub and Wang, 2008). Overall, our data showed that the decrease in N content in plants grown under elevated [CO2] and water stress during vegetative growth is genotypically dependent.

The δ 13C and δ 15N have been used as potential physiological tracers in plants under elevated [CO2] and water limitation (Aranjuelo et al., 2011; Yousfi et al., 2012, 2016; Araus et al., 2013; Bort et al., 2014). Elevated [CO2] and water stress caused an increase in δ 15N, although these effects depended on the genotype and were attenuated or disappeared in severe water stress relative to the moderate stress treatment (**Tables 3, 4**; Supplementary Table S2). Variations in δ 15N in response to the growth conditions, together with N content, could indicate shifts in N metabolism (Bort et al., 2014), although δ 15N is determined by many processes that are not completely understood (Ariz et al., 2015). Nevertheless, the higher δ 15N could suggest lower N availability, because N absorption and assimilation cannot fractionate between the 14N and 15N isotopologues under such environmental factors (Lopes and Araus, 2006; Tcherkez, 2011). Additionally, this could reflect a decrease in N translocation from the root to the shoot (Lopes and Araus, 2006). Moreover, δ 13C increased in some genotypes under moderate water stress, regardless of the [CO2] considered, but this increment, also observed under severe water stress, did not reach statistical significance (**Tables 3, 4**). Elazab et al. (2012) and Bort et al. (2014) also showed a δ 13C increase in flag leaves of different durum wheat genotypes under water stress at later growth stages, which could be associated with higher water-use efficiency (Araus et al., 2008, 2013; Tardieu, 2013; Bort et al., 2014). A stronger water stress does not always lead to larger changes in δ 13C, particularly when analyzed in dry matter, as noted in previous studies in rice (Kano-Nakata et al., 2014) and Pinus tabuliformis (Ma et al., 2014). In addition, δ 13C was strongly reduced at high [CO2] because of the very negative δ 13C of the CO2 used to increase the [CO2] within the growth chamber (Aljazairi et al., 2015).

## Expression of Stress-Responsive Genes in Durum Wheat Genotypes under Different Water Regimes and [CO2] Levels

Strong differences in gene expression were observed between treatments and among the different genotypes studied (**Table 5**; Supplementary Figure S1). In our study, RBCL and RBCS showed a common expression pattern (**Table 5**), confirming the coordinated expression of both subunits necessary for the assembly of the Rubisco holoenzyme (Suzuki and Makino, 2012). At the first sampling date, gene expression of RBCL and RBCS was down-regulated in response to elevated [CO2] no matter which water regime was considered, in agreement with other wheat studies (Aranjuelo et al., 2013; Habash et al., 2014; Vicente et al., 2015b). This down-regulation was associated with lower N content and higher δ 15N in a genotype-dependent manner. The former could be explained by non-selective decreases in N or reallocation of N within the plant under elevated [CO2] (Aranjuelo et al., 2011; Vicente et al., 2015a). The latter was probably associated with changes in N uptake, assimilation or redistribution within the plant (Araus et al., 2013). At the second sampling date, elevated [CO2] decreased the N content in Regallo and Burgos, which was related to down-regulation of transcript levels of Rubisco subunits and N-assimilation enzymes (GS1 and GS2), and higher root and plant biomass. These shifts could indicate that plant biomass might increase under elevated [CO2] in a genotype-dependent manner even when transcript levels of Rubisco subunits decrease during vegetative growth. This could be due to the remobilization of an N over-investment in Rubisco to reuse it in developing new tissues (Richards, 2000; Vicente et al., 2011; Carmo-Silva et al., 2015). However, the decrease in Rubisco transcript levels under water stress did not indicate the greater photosynthetic efficiency that was hypothesized under elevated [CO2]. Instead it was associated with lower plant biomass, which might suggest an inhibition of CO2 assimilation and plant growth in concordance with previous studies (Hayano-Kanashiro et al., 2009; Peremarti et al., 2014).

PEPC is a multifaceted key enzyme that in C3 plants is linked to the provision of Krebs cycle intermediates, and its overexpression in transgenic wheat improved drought tolerance and grain yield (Qin et al., 2015). PEPC expression has not been widely studied during early growth in durum wheat plants. In the current work it was induced under the combination of elevated [CO2] and moderate water stress in most genotypes, whereas at severe water stress genotypic variation determined its expression pattern (**Table 5**). The induction could be related to its major role in providing C skeletons for amino acid and lipid biosynthesis (González et al., 2003). This may be due to an increase in the enzyme's substrates, such as carbohydrates, typically found under elevated [CO2] and water stress (Khoshro et al., 2013; Vicente et al., 2015b). These results indicate that further work is necessary to broaden our understanding of the biological role of PEPC and its implication in plant growth, especially in genotypes (i.e., Ramirez) with an up-regulation of gene expression under stress conditions.

At moderate and severe water stress, GS1 and PEPC expression was significantly coordinated, as were the expressions of the GS2 and Rubisco genes (**Table 5**; Supplementary Table S3). Under severe water stress, GS1 and GS2 expression was more influenced by genotypic variability than environmental conditions. Yousfi et al. (2016) also reported genotypic differences in the expression of these genes under drought stress, with a general down-regulation under stress conditions. Lower N contents and transcript abundances for RBCL and RBCS under water stress and especially under elevated [CO2] were associated with higher repression of the GS2 gene, indicating a coregulation of primary C and N metabolism (Stitt and Krapp, 1999; Vicente et al., 2015b, 2016). In some treatments, mainly at the first sampling date, opposing gene expression patterns were observed between GS1 and GS2. This fact, together with the coordination of GS1 with PEPC, might indicate a predominant remobilization of C and N compounds and an inhibition of primary N assimilation under water stress and elevated [CO2]. Thus, the results support a significant coordination between C and N metabolism at the transcript level under conditions of elevated [CO2] and water stress. In addition, the pattern of gene expression for GS1 and GS2 supports the use of these genes as indicators of N metabolism under

water stress conditions, as reported previously (Nagy et al., 2013).

DHN11 and DHN16 encode for two dehydrins that belong to group 2 of late embryogenesis abundant (LEA) proteins (Ali-Benali et al., 2005). The up-regulation of dehydrin genes under water restriction is often associated with stress tolerance, although their specific role as osmotically active compounds is still unknown (Kosová et al., 2014). Moderate and severe water stress reduced DHN11 gene expression regardless of the [CO2] level compared with control conditions. In the case of the DHN16 gene, moderate water stress mostly up-regulated its expression, whereas under severe water stress the opposite occurred (**Table 5**). Elevated [CO2] at the first sampling date mostly enhanced DHN11 and DHN16 gene expression, while at the second sampling date its combination with severe water stress led to a wide range of changes in transcript levels in a genotypedependent manner. Our results showed that the pattern of gene expression could differ between dehydrins, in concordance with previous studies (Ali-Benali et al., 2005; Melloul et al., 2013; Kosová et al., 2014). Additionally, the severity of the water stress, [CO2] enrichment and the genotype influenced dehydrin transcript levels.

CAT and SOD enzymes form part of the system responsible for lowering ROS and avoiding oxidative stress. In general, gene expression of CAT and SOD was repressed under moderate water stress regardless of [CO2] (**Table 5**). Such repression was only maintained for CAT at severe water stress in Burgos, while their expression was up-regulated in the other genotypes under elevated [CO2] × severe water stress. This could suggest a higher demand for ROS control, which would indicate a limitation to the transfer of electrons through photosystems to drive C assimilation (Martins et al., 2016). Enzyme activity and CAT gene expression have been reported to decrease under elevated [CO2] in wheat, possibly due to the inhibition of photorespiration, while they increased only in response to severe drought (Luna et al., 2005; Xu et al., 2010; Vicente et al., 2015b). The available studies reporting changes in SOD gene expression and protein content under such conditions are contradictory, reporting different pattern of changes (Kim et al., 2006; Li et al., 2008; Caruso et al., 2009; Xu et al., 2010). Our results highlighted that water regime and genotype were key factors influencing the expression of genes involved in the antioxidant system, indicating a greater need for protection against oxidative damage under severe water stress.

#### Coordination between Physiological Traits and Transcript Levels in Durum Wheat Grown under Different Environmental Conditions during Vegetative Growth

The different changes in plant growth parameters indicate that the responsiveness to elevated [CO2] and water stress during early growth depends on (i) the duration of the treatment, because [CO2] enrichment results in greater increases in plant biomass in older plants; (ii) the severity of the water stress, which is more pronounced under severe water stress; (iii) and the genotypic variability. In general, elevated [CO2] stimulated plant growth and reduced N content, which at the transcript level was related to a down-regulation of Rubisco and N assimilation genes and up-regulation of genes that take part in C-N remobilization. Moderate water stress did not lead to gross changes in physiological traits, but severe water stress restricted plant growth and N content, while changes in gs and δ 13C suggested a water-saving strategy relative to wellwatered conditions. The transcript profile suggested an inhibition of primary C fixation and N assimilation, differences between dehydrins and a genotypic variation in gene expression under severe water stress, with an induction of genes involved in antioxidant machinery. The stimulation of plant biomass under elevated [CO2] did not compensate for plant growth limitation under water restriction. Lastly, we observed different genotypic responses to environmental factors, as also reported in barley (Ceccarelli et al., 1991). Regallo showed the lowest plant biomass and chlorophyll and N contents, which was related to a repression of genes for N assimilation and induction for dehydrins, SOD and CAT, while the opposite results were recorded for Burgos (data not shown). Therefore, increased plant growth was linked to up-regulation of N assimilation and down-regulation of stressresponsive genes, suggesting lower oxidative damage.

Considering different environmental conditions predicted for the future climate scenario and genotypic variations, network analysis was used to identify physiological traits, and transcript levels that are correlated during vegetative growth in durum wheat (**Figure 2**). Early growth is a positive trait for improving plant tolerance in water-limited environments that has the potential for larger final plant biomass and yield (Wilson et al., 2015). Plant growth parameters were positively correlated with each other in most cases, suggesting that early plant growth is driven by all plant fractions and tiller production, as reported in other studies (Rebolledo et al., 2013; Wilson et al., 2015). Regardless of genotype, the positive correlation between root and plant biomass was mainly due to the stimulation of root biomass under elevated [CO2], in agreement with previous reports (Madhu and Hatfield, 2013; and citations therein). In contrast, water restriction (mainly severe water stress) limited both root and shoot biomass, which are often diminished under severe drought conditions (Nezhadahmadi et al., 2013). Positive effects of elevated [CO2] on root biomass could mitigate drought effects on plant growth by allowing better exploitation of water and nutrients from deep soil layers (Madhu and Hatfield, 2013).

N content was correlated negatively with the root/shoot ratio and positively with the tillers per plant and shoot biomass, and this was probably due to the typically higher N content observed in shoots relative to roots (Vicente et al., 2015a). Hence, greater vegetative growth in durum wheat requires high amounts of N, which in turn will be conditioned by N availability. δ 15N has been proposed as an indicator of responses to stress, such as water stress, N starvation and salinity (Yousfi et al., 2012, 2016; Bort et al., 2014), although it has had little attention for studies of elevated [CO2] (Ariz et al., 2015). Here we observed a negative correlation of δ 15N with N content and tillers per plant, with elevated [CO2] being the main factor that increased δ 15N in our experiment. Nevertheless, the fractionating processes of N metabolism affecting δ 15N under elevated [CO2] and water stress are not fully understood (Tcherkez, 2011).

Leaf chlorophyll content has been extensively used as an indicator of different physiological and agronomical components, particularly at later growth stages (Araus et al., 2008). The network analysis confirmed that chlorophyll content is a positive trait for vegetative growth in durum wheat, and this can be easily implemented in most of studies because this measurement is simple, quick, and non-destructive with modern portable devices. Effects of elevated [CO2] and water stress on gs have been widely studied (Long et al., 2006; Pazzagli et al., 2016), including the proposal of gs as a trait indicator of drought stress tolerance (Nagy et al., 2013). In our study gs was negatively correlated with chlorophyll content and root and plant biomass. This could highlight that increased vegetative growth was related to stomatal closure, maybe as a water saving strategy or as a direct response to elevated [CO2].

The positive correlations among the transcript levels of the genes encoding RBCL, RBCS, GS1, GS2, and PEPC supported a balanced coordination between C and N metabolism under elevated [CO2] and water stress. On the other hand, our results underlined the key role of Rubisco and GS in plant responses to environmental conditions (Nagy et al., 2013; Carmo-Silva et al., 2015; Vicente et al., 2015b; Yousfi et al., 2016). We showed negative associations between transcript levels of Rubisco subunits and GS2 with chlorophyll content and plant biomass. This fact could indicate that a stimulation of plant growth may be associated with a lower investment of resources (mainly N) in Rubisco protein, especially under elevated [CO2], thus leading to a higher nitrogen efficiency (Pang et al., 2014; Carmo-Silva et al., 2015). The negative correlation between transcript levels of CAT and chlorophyll content highlighted that the up-regulation of CAT expression was a response to the high H2O2 levels generated under stress conditions (Luna et al., 2005), which could promote chlorophyll degradation (Upadhyaya et al., 2007). Interestingly, transcript levels of CAT were positively correlated with gs , although a negative correlation should be expected since greater gs leads to lower photorespiration rates and consequently lower H2O2 generation (Luna et al., 2005). We found a positive relationship between C content and transcript accumulation for SOD, not previously reported to our knowledge. Higher SOD expression might suggest a better ROS control that triggers an efficient electron transfer and C fixation. In our study, DHN11 transcript accumulation was negatively associated with root biomass, while transcripts for DHN16 were positively linked with plant biomass. These results suggest promising functions for DHN16 in stress tolerance during vegetative growth, as Kosová et al. (2014) proposed in a study examining wheat seed development.

In summary, parameters such as chlorophyll and N content, gs and δ 15N, and the expression of RBCL, RBCS, GS2, DHN11, and DHN16 genes were identified as good indicators for the selection of genotypes with better performance during early plant growth under elevated CO2 and water stress. Additionally, network analysis underlined the relevance of N metabolismtraits such as N content, δ 15N, GS1, and GS2, in the genotypic response of durum wheat to future environmental scenarios in the Mediterranean basin.

## CONCLUSION

We conclude that [CO2] effects on plant growth had greater impacts than moderate or severe water stress during vegetative growth of durum wheat. Whereas, elevated [CO2] generally led to increases in plant growth, water stress had a negative effect, preferentially as the water stress develops over time. In addition, the interactive effects of both [CO2] and water regime depends on genotypic variability. Gene expression profiles at moderate water stress were mainly affected by environmental conditions among the different genotypes. However, with further water restriction, genotype-specific differences were found to affect gene expression more than environmental conditions. These facts reflect a wide range of adaptation mechanisms in durum wheat under elevated [CO2] and water stress during vegetative growth, probably due to the complex regulatory network that takes place with both factors. Moreover, our study did not show a clear trend concerning the genetic advance in response to future climate change scenarios. Our results evidenced for durum wheat the need to take into account the genotypic variability for a greater understanding of plant adaptation to climate change. Moreover, the correlation network demonstrated that the combination of phenotyping and gene expression analysis is a useful approach to identify phenotype-genotype relationships and their behavior in response to different environments during vegetative stages.

## AUTHOR CONTRIBUTIONS

SM and JLA conceived and designed the experiments. SM, RV, and AA contributed to the experimental work. SM, RV, and JLA analyzed the data and interpreted the results. RV wrote the paper under the supervision of JLA, and SM and AA revised the manuscript. All authors have read and approved the final manuscript.

## FUNDING

This study was supported by the Spanish National Programme for Research Aimed at the Challenges of Society of the Ministry of Economy and Competitiveness (grants No. AGL2013-44147- R and AGL2016-76527-R). SM was the recipient of a fellowship "Presidente de la República PRONABEC-III" from Peruvian Government.

## ACKNOWLEDGMENTS

We thank the Unitat de Genòmica of the CCiTUB, Josep Matas (Servei de Camps Experimentals), Adrián Gracia of the University of Barcelona, and Marco Betti of the University of Seville for technical assistance.

## SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fpls.2016. 01738/full#supplementary-material

### REFERENCES

on drought variability in the Iberian Peninsula. Front. Environ. Sci. 3:1. doi: 10. 3389/fenvs.2015.00001

**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 © 2016 Medina, Vicente, Amador and Araus. 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) or licensor 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.

# SNP-Based QTL Mapping of 15 Complex Traits in Barley under Rain-Fed and Well-Watered Conditions by a Mixed Modeling Approach

Freddy Mora1 , Yerko A. Quitral 2 , Ivan Matus 3 , Joanne Russell 4 , Robbie Waugh4 and Alejandro del Pozo2 \*

#### Edited by:

Soren K. Rasmussen, University of Copenhagen, Denmark

#### Reviewed by:

Xiaoquan Qi, Institute of Botany-The Chinese Academy of Sciences, China Evelyne Costes, Institut National de la Recherche Agronomique, France

> \*Correspondence: Alejandro del Pozo<EMAIL_ADDRESS>
#### Specialty section:

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

Received: 24 January 2016 Accepted: 08 June 2016 Published: 27 June 2016

#### Citation:

Mora F, Quitral YA, Matus I, Russell J, Waugh R and del Pozo A (2016) SNP-Based QTL Mapping of 15 Complex Traits in Barley under Rain-Fed and Well-Watered Conditions by a Mixed Modeling Approach. Front. Plant Sci. 7:909. doi: 10.3389/fpls.2016.00909 1 Instituto de Ciencias Biológicas, Área de Biología Molecular y Biotecnología, Universidad de Talca, Talca, Chile, 2 Centro de Mejoramiento Genético y Fenómica Vegetal, Facultad de Ciencias Agrarias, PIEI Adaptación de la Agricultura al Cambio Climático (A2C2), Universidad de Talca, Talca, Chile, 3 Centro Regional de Investigación Quilamapu, Instituto de Investigaciones Agropecuarias, Chillán, Chile, 4 The James Hutton Institute, Dundee, Scotland

This study identified single nucleotide polymorphism (SNP) markers associated with 15 complex traits in a breeding population of barley (Hordeum vulgare L.) consisting of 137 recombinant chromosome substitution lines (RCSL), evaluated under contrasting water availability conditions in the Mediterranean climatic region of central Chile. Given that markers showed a very strong segregation distortion, a quantitative trait locus/loci (QTL) mapping mixed model was used to account for the heterogeneity in genetic relatedness between genotypes. Fifty-seven QTL were detected under rain-fed conditions, which accounted for 5–22% of the phenotypic variation. In full irrigation conditions, 84 SNPs were significantly associated with the traits studied, explaining 5–35% of phenotypic variation. Most of the QTL were co-localized on chromosomes 2H and 3H. Environment-specific genomic regions were detected for 12 of the 15 traits scored. Although most QTL-trait associations were environment and trait specific, some important and stable associations were also detected. In full irrigation conditions, a relatively major genomic region was found underlying hectoliter weight (HW), on chromosome 1H, which explained between 27% (SNP 2711-234) and 35% (SNP 1923-265) of the phenotypic variation. Interestingly, the locus 1923-265 was also detected for grain yield at both environmental conditions, accounting for 9 and 18%, in the rain-fed and irrigation conditions, respectively. Analysis of QTL in this breeding population identified significant genomic regions that can be used for marker-assisted selection (MAS) of barley in areas where drought is a significant constraint.

Keywords: drought, marker segregation distortion, RCSL, physiological trait, kinship

## INTRODUCTION

Abiotic stresses can significantly reduce crop yields and restrict the latitudes and soils on which commercially important species can be cultivated (Lobell and Field, 2007; Jacobsen et al., 2012). The development of drought-tolerant genotypes as well as genotypes with higher water-use efficiency is of global interest as populations continue to increase and water availability decreases (Araus et al., 2002, 2008; Cattivelli et al., 2008). In barley (Hordeum vulgare L.), a large number of morphological and physiological traits are linked to drought tolerance (Chen et al., 2010; del Pozo et al., 2012) which exhibit strong environmental interactions (Tondelli et al., 2006). Increasing tolerance to drought stress has become a major goal for barley breeding programs particularly in light of prolonged drought periods as a result of climate change (Wehner et al., 2015).

To determine the genetic basis of complex traits, important genetic, and genomic resources have been developed in a wide range of species (Kota et al., 2003; Mora et al., 2015), including barley (Kota et al., 2003; Close et al., 2009; Zhou et al., 2015). Single nucleotide polymorphism (SNP) markers are the most abundant sequence variations encountered in eukaryotic genomes (Griffin and Smith, 2000), and in barley offer the potential for generating very high-density genetic maps (Close et al., 2009), providing a useful tool for quantitative trait locus/loci (QTL) mapping for marker-assisted selection (MAS; Sato and Takeda, 2009; Szucs et al., ˝ 2009).

Conventionally, quantitative trait locus/loci (QTL) mapping is carried out using markers that follow a Mendelian segregation ratio, which depends on the population under investigation (Xu, 2008). The phenomenon that alleles at a locus deviate from the Mendelian expectation has been defined as segregation distortion (Zhang et al., 2010), which has been encountered in many commercially important species, such as maize (Lu et al., 2002), rice (Xu et al., 1997), tomato (Paterson et al., 1988), and barley (Malosetti et al., 2011). Biologically, segregation distortion can be due to chromosome loss, genetic isolation mechanisms, and the presence of viability genes. Non-biological factors such as scoring errors and sampling errors also can contribute to segregation distortion (Alheit et al., 2011).

Construction of genetic maps and QTL analysis using distorted markers is risky because the basic assumption of Mendelian segregation is violated, and, consequently most of these markers are removed from the subsequent QTL analyses (Luo et al., 2005; Xu, 2008). Moreover, several authors found that distorted markers influence the estimation of genetic intervals and the markers order on a same chromosome (Lorieux et al., 1995a,b; Zhu et al., 2007). Because these markers are routinely removed, valuable information around these regions is lost (Luo et al., 2005). However, if distorted markers are handled properly, their effects on genetic map construction and QTL identification and mapping can be significantly improved (Xu, 2008; Alheit et al., 2011; Hashemi et al., 2015).

From an analytical standpoint, many QTL mapping procedures have been developed for Mendelian populations, but few are available for markers that do not segregate in a typical Mendelian fashion (Zhu et al., 2007; Xu, 2008). Recently, Malosetti et al. (2011) found severe allele frequency distortions in many chromosomal regions in a RIL population of barley. The authors stated that violation of the basic assumptions implies that the genetic covariance between genotypes (i.e., genetic relatedness) in the population is not homogeneous, and analogous to association mapping studies, a QTL mapping mixed model that account for this heterogeneity should be used to avoid false QTL detection. In addition, mixed models have been used to investigate QTL-by-environment interaction (Malosetti et al., 2004; Boer et al., 2007) and to map QTL for several traits simultaneously (Malosetti et al., 2008).

In barley, a large number of mapping populations have been developed to map QTL. Further, advanced mapping populations, including near-isogenic lines (NILs), chromosome segment substitution lines (CSSLs), and recombinant chromosome substitution lines (RCSLs, Schmalenbach et al., 2008, 2011; Sato and Takeda, 2009; Naz et al., 2012), have also been developed to facilitate the genetic dissection of complex traits. As a consequence, many QTL controlling complex traits including agronomic, and morphological traits, yield component, disease resistance, tolerance to abiotic stress, and malting quality have been identified (Pillen et al., 2004; Talame et al., 2004; Li et al., 2005; Gyenis et al., 2007). However, studies that have identified QTL for drought-related morphological and physiological traits are still scarce in barley (Chen et al., 2010; Mir et al., 2012; Sayed et al., 2012; Kalladan et al., 2013; Li et al., 2013; Wójcik-Jagła et al., 2013; Honsdorf et al., 2014; Mansour et al., 2014; Naz et al., 2014). Thus, a breeding program has been developed in the Mediterranean area of central Chile, in areas where drought is a significant constraint to yield. This study identifies SNP markers associated with 15 complex traits (including physiological and morphological traits) in a breeding population consisting of 137 RCSLs of barley (Matus et al., 2003), evaluated under contrasting water availability. Given that markers showed a very strong segregation distortion, a QTL mapping mixed model was employed to account for the heterogeneity in genetic relatedness between genotypes.

## MATERIALS AND METHODS

#### Plant Material and Field Evaluation

A set of 137 RCSL were evaluated in field conditions in Santa Rosa (35◦ 78′ S, 72◦ 17′ W), in the Mediterranean climatic region of central Chile, under rain-fed and fully irrigated conditions in 2008–2009. The average annual temperature in this region is 13◦C, the minimum average is 3◦C (July), and the maximum 28.6◦C (January). Monthly maximum and minimum temperatures and precipitation during 2008 are in Table S1. The annual precipitation in 2008 was 992 mm but the amount during the growing season was 120 mm. The soil is a sandy loam, Humic Haploxerand, Andisol. For fully irrigated conditions, four furrow irrigation of 50 mm water was applied from heading to maturity.

The RCSL population was developed using the advanced backcross strategy of Tanksley and Nelson (1996). The accession of Hordeum spontaneum (Caeserea 26-24 from Israel) was the donor parent and Hordeum vulgare subsp. vulgare "Harrington" (a North American malting quality standard) was the recurrent parent (Matus et al., 2003). The recurrent parent was used as the female and the donor as the male to obtain the F1 generation. Finally, the lines were obtained using two backcrosses with the recurrent parent and six generations of self-pollination (BC2F6).

The field trial was arranged in a 14 × 10 alpha-lattice design with two replications ("Harrington" cultivar was replicated more times for arrangement v = 14 × 10 = 140). Fertilizer and field management practices recommended for optimum barley production were used (Inostroza et al., 2009; del Pozo et al., 2012). The 15 different morphological, agronomic, and physiological traits measured during this experiment are described in **Table 1**.

#### Phenotypic Data Analysis

A mixed linear modeling approach was employed for phenotypic data analysis using the MIXED procedure in SAS. Field data were analyzed on the basis of the statistical model (Stich et al., 2008):

$$\mathcal{y}\_{ijno} = \mu + \mathfrak{g}\_i + l\_j + (\mathfrak{gl})\_{ij} + r\_{nj} + b\_{onj} + \varepsilon\_{ijno}$$

where yijno is the observed phenotype for the ith RCSL at the jth location in the oth incomplete block of the nth replicate, µ is an intercept term (overall population mean), gi is the genotypic effect of the ith entry, lj is the effect of the jth location, rnj is the effect of the nth replicate of the jth location, bonj is the effect of the oth incomplete block of the nth replication of the jth location, and εijno is the residual or within-plot error. To analyze the effect of environment (rain-fed and fully irrigated) on the RCSL population, environment (contrasting environmental condition), genotype, and genotype-environment interaction were considered as fixed effects while replication and block were analyzed as random effects. Error variances were assumed to be heterogeneous among locations according to Stich

TABLE 1 | List of 15 morphologic, agronomic, and physiological traits measured in a RCSL population of barley under two contrasting environmental conditions in southern Chile.

et al. (2008). Adjusted entry means Mi were calculated for each RCSL as:

$$M\_i = \hat{\mu} + \hat{\mathfrak{g}}\_i$$

where µˆ and gˆi are the generalized least-squares estimates of µ and gi , respectively. Adjusted entry means, calculated by environment and in a combined analysis, were latter used in the analysis of QTL. The CORR procedure of SAS was used to estimate Pearson correlations (r, n = 137) between pairs of traits.

#### DNA Extraction and Genotyping

Leaf tissues were harvested from each plant and crushed with liquid nitrogen. Genomic DNA was extracted from 200 to 300 mg of leaf sample using Qiagen DNeasy Plant mini kit (QIAGEN Co.). DNA concentration was determined by nanodrop (Thermo Fisher Scientific) and adjusted to 100 ng/µl. DNA samples were sent to the Southern California Genotyping Consortium, Illumina BeadLab at the University of California, Los Angeles (UCLA) for the OPA-SNP assay with the 1536 plex detection platform of barley OPA 1 (BOPA1) developed by Close et al. (2009). SNP markers are distributed across the seven barley chromosomes. The SNP loci were designated by HarvEST:Barley unigene assembly #32 numbers (http:// harvest.ucr.edu/). Genotyping of the 137 RCSLs and their parents was performed by Illumina GoldenGate assay. The order of polymorphic markers from BOPA1 was performed using MEGA5 software (Tamura et al., 2011). Moreover, the chromosome segments introgressed into H. vulgare cv. Harrington from Hordeum spontaneum were estimated from the graphical haplotypes (Van Berloo et al., 2008) in each of the recombinant lines selected.

#### QTL Mapping by a Mixed Modeling Approach

For the segregation data of each SNP marker, deviations from the Mendelian ratios (1:1 ratio for RSCL populations) were tested using the Chi-square test. Given that SNP markers did not segregate in an expected Mendelian ratio (details in the Results Section) a mixed model (Stich et al., 2008) was employed for QTL analysis of a designed cross using a structured variance–covariance matrix, where the structure was induced by selection (Malosetti et al., 2011). This QTL mapping mixed model accounts for the heterogeneity in genetic relatedness between genotypes. The hypothesis of association of SNP markers with the 15 target traits was tested using the following mixed model implemented in the program TASSEL 3.0 (Bradbury et al., 2007):

$$M\_i = \mu\_1 + \varkappa\_i \alpha + \mu\_i + e\_i \tag{1}$$

where Mi is the adjusted entry mean of the ith RCSL, µ1 is an intercept term, xi is the SNP genotype of ith RCSL, α is the additive allele substitution effect (SNP effect), ui is the residual genetic background effect of the ith entry, and ei is the random residual effect. ui is assumed to follow a normal distribution with variance-covariance matrix G = σ 2 u · 2· K, where K is the coefficient of co-ancestry matrix between entries. The threshold used for declaring an association significant was P < 0.01. According to Malosetti et al. (2011) this mixed modeling approach for QTL detection can accommodate the extra genetic covariance by embedding kinship information in the model, leading to appropriate tests, and minimizing the rate of false QTL or gene detection. Variance components were estimated using the Restricted Maximum Likelihood (REML) method. Additionally, false discovery rate (FDR)-adjusted p-values were calculated using PROC MULTTEST in SAS software. The Bayesian information criterion (BIC) was used to compare the simple model that ignore kinship information and the model with a structured variance–covariance matrix (kinship information).

#### RESULTS

#### Phenotypic Data Analysis

The statistical analysis of fixed effects for the 15 complex traits under study are summarized in **Table 2**. According to the F-values (type III tests of fixed effects), the 137 RCSL presented significant differences at P < 0.01 in most traits under study. Significant differences were observed between both contrasting water regimes for most traits, including biological yield (BY) and grain yield (GY). In fact, the average GY under rain-fed was reduced by 81% in relation to fully irrigated condition, indicating that the barley RCSL lines were exposed to a severe water stress. Also, there was a strong reduction in BY (−59.9%), dry weight at tillering (DWT; −39.9%), tiller number (TN; −48.1%) and thousandkernel weight (TKW; −28.1%) (**Table 2**). Peduncle length (PL), peduncle extrusion (PE), plant height (PH), and the physiological trait IPAR also evidenced significant environmental effect. There were significant interactions between RCSL and environment for PE, PH, TN, hectoliter weight (HW), kernels per spike (KS), and GY.

Estimates of phenotypic correlations among traits are shown in **Table 3**. Grain yield was positively and significantly correlated with TN, BY, HW, HI, and KS, and negatively correlated with PL and PH (P < 0.01, **Table 3**). The physiological trait IPAR was significantly and positively correlated with biological yield (r = 0.6) and grain yield (r = 0.41). Relative water content (RWC), PE, and chlorophyll fluorescence (Fv/Fm) were not significantly correlated (P > 0.05) with most of the variables studied.

#### Kinship Matrix to Avoid False-Positives

A large number of spurious QTLs were detected when the genetic covariance matrix is ignored in the mixed model. In fact, 444, 460, and 516 false-positives were detected for all traits studied in the rain-fed and fully irrigated conditions, and combined analysis, respectively, using a structured variance–covariance matrix in a mixed model, where the structure was induced by selection. This QTL mapping mixed model accounts for the heterogeneity in genetic relatedness between genotypes. **Figure 1** shows the results of segregation distortion analysis based on the P-values, obtained by the chi-square test (P-values were plotted

TABLE 2 | Summary of statistical analysis of fixed effects for 15 complex traits measured in a RCSL population of barley under rain-fed and well-watered conditions.

Data are presented as phenotypic means with minimum and maximum values in parentheses. GEI, genotype-environment interactions; NS, not significant.

\*Significant at the 0.05 probability level.

\*\*Significant at the 0.01 probability level.

TABLE 3 | Pearson correlation coefficients among mean variables (adjusted entry means) measured in a RCSL population of barley under rain-fed and well-watered conditions.

\*Significant at the 0.05 probability level.

\*\*Significant at the 0.01 probability level.

on a −Log10 scale). According to the Bayesian information criterion (BIC) there was significant evidence against the simple models that ignore kinship information in most of the traits studied (Table S2); between the two competing models, the best model is the one that has the smallest BIC-value. In the combined analysis, BIC was smaller for the mixed model that includes the kinship information (K model), for most of the traits. Therefore, the use of a QTL detection model that accounts for the heterogeneous genetic relatedness between RCSLs lines, caused by the uneven sharing of genetic background, is clearly necessary.

## SNP-Based QTL Mapping of Complex Traits

The number of QTLs detected for all the traits under study, including the chromosome number and the percentage of the total phenotypic variation explained by the QTLs, are summarized in **Table 4** (details of all significant SNP–trait associations are given in Table S3). The Manhattan plots of the genome-wide QTL study showed a lower number of significant QTLs under rain-fed (**Figure 2**) in comparison with fully irrigated conditions (**Figure 3**). In fact, 57 QTLs were detected in rain-fed, which accounted for 5–22% of the phenotypic variation (KS and PH, respectively). In the fully irrigated condition, a total of 84 SNPs were significantly associated with the traits studied, which explained from 5 to 35% of phenotypic variation (TN and HW, respectively). In the combined analysis, 92 QTLs were detected to be highly significant, which accounted for 5–21% of the phenotypic variation.

Most of the QTLs were co-localized on chromosomes 2H and 3H; i.e., about 53% (30/57) and 62% (52/84) in the rain-fed and irrigation conditions, respectively. The chromosomes with lowest number of QTLs found were 4H and 6H, with two QTLs each, which accounted for 5.1–8.3% of the phenotypic variation. Only two QTLs were detected for grain yield in rain-fed condition, SNPs 2711-234 (1H) and 1923-265 (1H), which explained 7 and 8.6% of phenotypic variation, respectively, but importantly, they were stables across the two contrasting water regimes. Under full irrigation, six QTLs were detected for grain yield, SNPs 2711- 234 (1H), 1923-265 (1H), ConsensusGBS0598-3 (3H), 9018-522 (3H), 4105-1417 (3H), 7045-950 (3H), and explained from 5.2 to 18%. Different to grain yield, the putative QTLs underlying biological yield were detected on chromosomes 5H and 7H, in the rain-fed condition, and on chromosomes 2H, 3H, and 5H in the fully irrigated condition.

PV (%), percentage of the phenotypic variation explained by SNP markers; NQ, number of significant QTLs; ChN, chromosome number.

All QTLs detected for BY, Fv/Fm, DWT, HI, IPAR, KS, PE, PH, SL, RWC, TKW, and TN were environment-specific. Although most QTL-trait associations were environment specific, some stable associations were also detected for HW, SNPs 2711-234 (1H) and 1923-265 (1H), and PL, SNPs 9282-205 (3H), 3965-353 (3H), 2335-1614 (3H), and 5212-1409 (7H).

Fifteen and ten SNPs, respectively, were associated with more than one trait in the fully irrigated and rain-fed conditions. Importantly, the SNPs 4105-1417, 7045-950, 9018- 522, ConsensusGBS0598-3 (all on chromosome 3H) were concomitantly associated with grain yield, harvest index, and kernel per spike, in the fully irrigated condition, which accounted for 5.2–20% of the phenotypic variation. The mentioned SNPs ConsensusGBS0598-3, 7045-950, and 4105-1417 were also associated to biological yield, and the SNP 9018-522 with spike length. In both environmental conditions, grain yield and HW shared QTLs linked to the SNPs 1923-265 and 2711-234 (on chromosome 1H); both QTLs correspond to stable associations across environments.

QTLs with >15% of the phenotypic variation explained by SNP markers are shown in **Table 5**; these are relatively moderate (to major) QTLs detected in this RSCL population of barley. The results confirmed the presence of more QTLs with >15% in the fully irrigated (9) than the rain-fed condition (2), and the majority of these QTLs were environment-specific. In the fully irrigated condition, a relatively major QTL was detected underlying HW, i.e., SNP marker 1923-265 on chromosome 1H at 140 cM, which accounted for 35.3% of the total phenotypic variation. In the same genomic region, the SNP marker 2711-234 (on 1H at 139 cM) explained 27% of the phenotypic variation for HW. These QTLs were stable across the contrasting environmental conditions. Interestingly, the major locus 1923-265 (1H) associated with HW, was also detected for grain yield in both environmental conditions, which explained between 8.6 and 18.1% of the phenotypic variation, in the rain-fed and fully irrigated conditions, respectively. Finally, one genomic region on chromosome 4H, was moderately associated with plant height (i.e., SNPs ABC08009-1-2-304 and 954-1377) and was environment-specific in the rain-fed condition.

#### DISCUSSION

In this study, a RCSL population consisting of 137 lines was evaluated for 15 complex traits, including plant height, grain yield, and yield-related traits in two contrasting environment conditions. Severe allele frequency distortions was evidenced along the seven linkage groups in the RCSL population. By incorporating kinship into the model, as proposed by Malosetti et al. (2011) for populations that have undergone some selection resulting in a departure from Mendelian segregation ratios, we identified QTLs associated with key traits in two contrasting field conditions.

Environment-specific genomic regions were detected for the majority of the traits (12/15). These findings are consistent with the study conducted by Wang et al. (2014), in which most of the QTL for different traits varied between environments in a doubled haploid population of barley.

Most of the economically important traits in barley are inherited quantitatively. Plant height (PH), for instance, is under polygenic control, and represents one of the most important agronomic traits for barley (Wang et al., 2014; Zhou et al., 2015). In this study, the highest number of QTL was detected for plant

height (20 QTL in both conditions) which were observed on chromosomes 2H, 3H, 4H, and 6H. Similarly, Honsdorf et al. (2014) found the highest number of associations for this trait, which were located on all chromosomes except 5H. Inostroza et al. (2009) found SSR-trait associations for PH on chromosomes 1H, 2H, 4H, 5H, 6H, and 7H, evidencing a genome-wide distribution. Interestingly, in the present study, a significant genomic region (explained 22% of phenotypic variance) on chromosome 4H, which comprises the SNPs ABC08009-1-2-304

and 954-1377, controls PH in barley under rain-fed conditions. The majority of the QTLs were detected on chromosomes 2H (4/7) under fully irrigated and 4H (5/7) under rain-fed condition. This result is partially consistent with Pasam et al. (2012), who found 32 associations with plant height with the majority located on chromosomes 2H and 3H. Malosetti et al. (2011) found significant SNP associations with PH on chromosomes 2H, 3H, 5H, and 7H, with two important QTLs on chromosomes 3H (127.1 cM) and 5H (69.3 cM), which are known to carry semi-dwarfing genes in barley (Malosetti et al., 2011; Wang et al., 2014). The use of semi-dwarf genes has greatly improved barley TABLE 5 | Relatively moderate (or major) QTLs detected in a RSCL population of barley (only QTLs with >15% of the phenotypic variation explained by SNP markers), evaluated under well-watered and rain-fed conditions by a mixed modeling approach.

PV(%), percentage of the phenotypic variation explained by SNP markers; ChN, chromosome number; Pos, SNP position in cM.

yields with controlled plant height being used to reduce yield loss arising from lodging and to increase the harvest index (Wang et al., 2014). Zhou et al. (2015) found a major QTL for plant height mapped at 105.5 cM on chromosome 3H, which had a LOD score of 13.01 and explained 44.5% of phenotypic variation. In this study, one QTL was detected at 197 cM on chromosome 3H (SNP 9610-1195), in the rain-fed condition, which explained 7.3% of phenotypic variation.

In the fully irrigated condition, a relatively major QTL was found for HW on chromosome 1H at 140 cM (SNP 1923-265), which explained 35.3% of phenotypic variation. This result is consistent with the study of Rode et al. (2012) who found three QTLs for HW on chromosome 1H, including the SNP 1923-265. In rain-fed condition, this SNP explained 8.6% of phenotypic variation for grain yield, and under fully irrigated condition, it explained 5.7 and 18.1% for RWC and GY, respectively. In contrast, Rode et al. (2012) did not find any SNP controlling HW associated with another related-trait. As expected, the SNP marker 2711-234 (in the same genomic region of SNP 1923- 265 on chromosome 1H, at 139 cM) that explained 27% of the phenotypic variation for HW, was associated with grain yield in both environmental conditions, explaining 10 and 7% of phenotypic variation. The correlation coefficient between both traits was positively correlated (r = 0.48; P < 0.01; **Table 3**).