Genetic Loci Governing Grain Yield and Root Development under Variable Rice Cultivation Conditions

Drought is the major abiotic stress to rice grain yield under unpredictable changing climatic scenarios. The widely grown, high yielding but drought susceptible rice varieties need to be improved by unraveling the genomic regions controlling traits enhancing drought tolerance. The present study was conducted with the aim to identify quantitative trait loci (QTLs) for grain yield and root development traits under irrigated non-stress and reproductive-stage drought stress in both lowland and upland situations. A mapping population consisting of 480 lines derived from a cross between Dular (drought-tolerant) and IR64-21 (drought susceptible) was used. QTL analysis revealed three major consistent-effect QTLs for grain yield (qDTY1.1, qDTY1.3, and qDTY8.1) under non-stress and reproductive-stage drought stress conditions, and 2 QTLs for root traits (qRT9.1 for root-growth angle and qRT5.1 for multiple root traits, i.e., seedling-stage root length, root dry weight and crown root number). The genetic locus qDTY1.1 was identified as hotspot for grain yield and yield-related agronomic and root traits. The study identified significant positive correlations among numbers of crown roots and mesocotyl length at the seedling stage and root length and root dry weight at depth at later stages with grain yield and yield-related traits. Under reproductive stage drought stress, the grain yield advantage of the lines with QTLs ranged from 24.1 to 108.9% under upland and 3.0–22.7% under lowland conditions over the lines without QTLs. The lines with QTL combinations qDTY1.3+qDTY8.1 showed the highest mean grain yield advantage followed by lines having qDTY1.1+qDTY8.1 and qDTY1.1+qDTY8.1+qDTY1.3, across upland/lowland reproductive-stage drought stress. The identified QTLs for root traits, mesocotyl length, grain yield and yield-related traits can be immediately deployed in marker-assisted breeding to develop drought tolerant high yielding rice varieties.


INTRODUCTION
The development of rice cultivars with improved tolerance for drought stress is important to increase the production of rainfed rice ecosystems. Efficient use of existing genetic variability in traditional cultivars by employing advanced tools and techniques is needed to combat the adverse effects of climate change on food security and agriculture sustainability. The Post-Green Revolution high-yielding, biotic stress-resistant but drought-susceptible varieties were bred for the targeted irrigated ecosystem. Exploitation of genetic variation, use of holistic modified breeding strategies, and direct selection for grain yield (Kumar et al., 2014) combining traits contributing yield advantage has been suggested as an appropriate approach to develop drought tolerant rice cultivars. Several studies involving grain yield as the main selection criterion have identified stable, consistent and large-effect QTLs for grain yield under reproductive-stage drought stress (Bernier et al., 2007;Venuprasad et al., 2009;Vikram et al., 2011;Ghimire et al., 2012;Mishra et al., 2013;Swamy and Kumar, 2013;Yadaw et al., 2013;Dixit et al., 2014a,b). Some of these identified QTLs have been deployed to develop high-yielding reproductive-stage droughttolerant rice varieties (Kumar et al., 2014). However, although the reproductive stage is the growth stage at which rice yield is most affected by drought, rice yields in farmers' fields may be reduced by drought occurring at any growth stage. Therefore, continued efforts are needed to identify QTLs and traits that confer improved rice yields under multiple types of drought stress.
Phenotypic screening involving direct selection for grain yield under reproductive-stage drought stress, non-stress (control), multiple environments (upland and lowland), locations, seasons/years has led to the identification of major and consistent effect grain yield QTLs (Kumar et al., 2014;Dixit et al., 2015;Sandhu et al., 2015).Identification of co-located genetic regions associated with grain yield under drought, root and seedling establishment traits have opened up further possibilities for improving rice yield under drought stress.
To date, most of the reports identifying major-effect QTLs for rice yield under drought have involved selection of a traditional variety as the drought tolerance donor. Many traditional upland/aus cultivars showed a high level of drought tolerance, but their drought tolerance is often linked to undesirable traits, such as low yield potential and/or tall plant height (Vikram et al., 2015). Therefore, preliminary studies confirming the efficiency of the drought donor in a breeding program is necessary (Kumar et al., 2014). One such potential donor is the traditional variety Dular, which was identified in early rice genotype screening efforts as a drought-tolerant cultivar that maintained grain yield under reduced water availability and showed deep root growth (De Datta et al., 1975). Subsequent reports highlighted Dular as a traditional variety with the ability to maintain seminal root elongation under drought, but with less lateral root formation and plasticity in response to fluctuating soil moisture (Bañoc et al., 2000), and as one of the most deep-rooted genotypes of the diverse OryzaSNP panel in solution culture, lysimeter, and in upland and lowland drought field studies (Henry et al., 2011;Gowda et al., 2012;Shrestha et al., 2013;Wade et al., 2015). Given the previous research identifying Dular as deep-rooting, we hypothesized that RILs with the highest grain yield under drought would be those with highest root growth at depth, which would be related to seedling stage crown root formation.
To harness the drought tolerance (yield and root growth) traits of Dular for use in breeding for multiple drought stress environments, this study involved recombinant inbred lines (RILs) derived from a cross of IR64-21 with Dular. The aim was to identify the genetic regions for grain yield, yield-related agronomic and root development traits under direct seeded upland and transplanted lowland reproductive-stage drought stress and non-stress conditions.

MATERIALS AND METHODS
Recombinant inbred lines (RILs-F 7 ) were developed at the T.T. Chang Genetic Resources Center at the International Rice Research Institute (IRRI), Los Baños, Philippines, using single-seed descent from F 2 progenies of the cross between IR64-21 (single plant selection from IR64, IRGC 117268) and Dular (purified line, IRGC117266). Dular is an upland adapted, drought-tolerant traditional cultivar from India (McNally et al., 2009) and has showed consistent performance in drought screening at IRRI. Dular is early maturing and has low-yield potential due to low tillering. IR64-21 is a lowland adapted, semi dwarf, medium-duration, high-yielding cultivar with good grain quality that has been widely grown across large areas of Asia and Africa. Eight field and two greenhouse experiments were conducted at IRRI (14 • 11 ′ N 121 • 15 ′ E, 21 m above sea level) during the dry seasons (DS) of 2013 and 2014 ( Table 1). The soil type in the experimental fields is a Maahas clay loam; isohyperthermic mixed typic tropudalf (Zhao et al., 2006;Venuprasad et al., 2009).

Management of Upland and Lowland Field Experiments
In this study, upland non-stress (UNS) refers to experiments in non-flooded, non-puddled, rainfed, naturally well drained fields where seeding was done when dry, whereas lowland non-stress (LNS) refers to experiments under flooded, puddled, transplanted and anaerobic conditions. Upland drought stress (URS) and lowland drought stress (LRS) refer to experiments comprising cyclic naturally imposed drought stress at the reproductive stage.
Upland experiments were conducted during the dry seasons of 2013 and 2014 with 480 F 7 RILs in α-lattice design with two replications of 2-m single-row plots with row spacing of 25 cm, except in the 2014 DS UNS experiment in which 1.5-m single row plots were used. The two parents, IR64-21 and Dular, were used as checks and were replicated five times within each replication. Seeds were dry direct-seeded in aerobic soil using a seeding rate of 5 g per row. Fertilizer with the proportion 100N:40P:40K was applied; P and K were applied as basal doses at 10 days after seeding (DAS) and N was applied in three equal splits at 10, 25, and 45 DAS. Hand weeding was done as needed and insecticide (Furadan; carbofuran, 1 kg a.i ha −1 ) was used to control insect pests when necessary. The stress experiments were sprinklerirrigated twice weekly during establishment and early vegetative growth. Tensiometers were installed at a depth of 30 cm to record the soil water tension. Drought stress was initiated at 35 DAS and plots were irrigated only when most lines were wilted and exhibited severe leaf drying and the soil water potential fell below −50 kPa. This cyclic reproductive stage drought stress allows effective reproductive-stage drought screening for broad range growth duration genotypes (Lafitte et al., 2004). Upland nonstress experiments received the same cultural practices as the stress experiments except that the irrigation was continued twice weekly up to 1 week before harvest. The total rainfall accumulated during stress period in the upland experiments (35-90 DAS) was 54.6 mm and 28.4 mm during 2013DS and 2014DS, respectively. The same 480 RILs were used in the lowland experiment conducted during 2013DS while 300 randomly selected RILs were used in the 2014DS lowland experiments. Soil in the 2014DS field study was characterized as having a water-extractable pH of 6.5, 33.3 mg kg −1 P (Olsen), 441 mg kg −1 K, 40% clay, 21% sand, and 39% silt. Seeds were sown in a raised bed nursery; 21-day-old seedlings were transplanted to the main field. In the 2013DS lowland experiments, one seedling per hill was planted in an α-lattice design with two replications of 5 m single-row plots, with hill and row spacing of 0.2 m. In the 2014DS lowland experiments, three or four seedlings per hill were planted in an α-lattice design with two replications in plots of seven 1-m rows, with hill spacing of 0.2 m and row spacing of 0.25 m. After transplanting, ∼5 cm of standing water was maintained in the lowland fields before the initiation of drought stress at 50 DAS in both 2013DS and 2014DS. Plots in the drought stress treatment were rewatered when the soil water potential dropped to −50 kPa. The total rainfall accumulated during the stress period in the lowland (50-100 DAS) was 197.4 and 14.2 mm during 2013DS and 2014DS, respectively. During the 2014DS LRS experiment, water potential dropped below−50 kPa twice (68 and 91 DAS) and was rewatered twice (79 and 102 DAS). Irrigation was maintained up to 10 days before harvest for the non-stress experiments. Furadan (carbofuran, 1 kg a.i ha −1 ) was applied to control stem borer and other insects when necessary.

Field Studies
Days to emergence was recorded in the 2014DS upland stress and lowland experiments. Seedling stage growth rate, in terms of the increase in height in centimeters per week from 28 to 56 DAS, was recorded in both stress and non-stress upland experiments of 2014DS.
Destructive sampling of three plants (shoot and root crown) was done at 42 DAS in the 2014DS upland stress experiment to evaluate crown root traits. Shoots were separated from the roots and then dried and weighed to determine shoot dry mass while the roots were stored in 50% ethanol for root counting. Numbers of crown (nodal) roots, adventitious (mesocotyl-borne) roots and seminal roots were counted manually and the sum of all axial crown root types was defined as the total root number. The mesocotyl length was measured using a centimeter scale. All crown root number and mesocotyl length values from the three root samples per plot were averaged.
In the 2014DS lowland experiments, root samples were acquired using a 4-cm diameter core sampler (fabricated at IRRI, Los Baños, Philippines) to a depth of 60 cm. One core per plot was sampled at 104-109 DAS, soil cores were divided into 15-cm segments and roots were washed by repeatedly mixing the soil with water in a container and pouring the root-water suspension over a 1-mm plastic sieve. Roots were scanned (Calibrated Color Optical Scanner STD4800, Epson) and analyzed (WinRHIZO Pro v. 2013e, Regent Instruments) to determine root length density in each depth segment. Sampled roots were dried and weighed to measure the root dry weight. Xylem sap bleeding rate, defined as the weight of xylem sap bled after detopping normalized by the shoot mass of each hill, was measured at 89 or 90 DAS as described by Morita and Jun (2002). Shoots of one hill per plot in the drought stress experiment were cut about 10 cm above the soil surface. Sap exuded from the cut stems was collected with a preweighed cotton towel wrapped in plastic. Canopy temperature was measured on six dates between 71 and 97 DAS by the tractor-based phenotyping system (Tanger et al., 2017).
In all field studies, days to 50% flowering, plant height (cm) at maturity and grain yield (kg ha −1 ) were recorded ( Table 2). Days to 50% flowering (DTF) was recorded as the number of days from seeding to the day on which 50% of the panicles had emerged. Specifically, it is at growth stage 5 (inflorescence emergence), code 55 of the BBCH-scale for rice (Lancashire et al., 1991). The height of three randomly chosen plants per plot was measured at maturity from ground level to the tip of the highest panicle and then averaged. At physiological maturity, grains were harvested 6PYT, plant height at 6 weeks after seeding (cm); 6SDW, shoot dry weight at 6 weeks after seeding (kg ha -1 ); EME, days to emergence; ARN, number of adventitious roots; CRN, number of crown roots; SRN, number of seminal roots; RAS, number of root axial sites; TRN, total number of roots; ML, mesocotyl length (mm); DR, percent deep roots; SR, percent shallow roots; LAT, percent lateral roots; RBN, root branching number; RDW, root dry weight (g); RLD, root length density (cm cm -3 ); BR, bleeding rate (sap/g); CT, canopy temperature; PU, phosphorus uptake (mg P/tiller); TWU, total water uptake (L); WUE, water use efficiency (g/L); PDW, panicle dry weight with rachis (g); RD, maximum root depth (cm).
from each plot, dried (moisture content ∼14%) and weighed to calculate grain yield in kg ha −1 . Shoot samples were collected at harvest from the 2014DS UNS, LNS, and LRS experiments, oven dried, threshed and weighed for the calculation of total biomass (kg ha −1 ) and harvest index (HI). HI was calculated as the ratio of grain weight to the total above-ground plant weight. Total numbers of panicles and tillers were counted from the whole plant samples and number of tillers and panicles per square meter were calculated from the 2014DS UNS.

Greenhouse Basket Root Experiment
In 2013WS, seedling-stage root growth of 300 RILs (the same RILs used in 2014DS lowland experiments) was measured using the basket method as described by Uga (2012). Open stainlesssteel mesh baskets (top diameter of 7.5 cm, depth of 5.0 cm and mesh size of 2 mm) were filled with sieved, dried soil from the IRRI farm. Each soil-filled basket was placed on a support made of PVC pipe (7.5-cm diameter, 15-cm height) and arranged in plastic trays (56 × 36.5 cm, 15-cm height; 28 baskets per tray). Each basket included a ring indicator dividing the basket into the upper and lower halves. Roots emerging above the ring indicator (shallow roots) represented an angle of from 0 to 50 • , while roots emerging below the ring indicator (deep roots) represented an angle of from 50 to 90 • with respect to the horizontal soil surface centered on the plant stem. The 300 RILs were grown in three replicates, with three successive germination dates for the three replicates (Jul 8, 2013, Aug 9 2013 and Sep 13 2013). Seeds were germinated in Petri dishes lined with filter paper, then transferred to the middle soil surface of each basket after 2-6 days. The nutrient solutions were monitored and recorded daily. The setup was watered using tap water for the first week, half strength Yoshida solution for the second week, and full-strength Yoshida solution from the third week onward. Average temperature and humidity were 30 • C and 60%, 30 • C and 50% and 28.5 • C and 45% in replicates 1, 2, and 3, respectively. The plants were harvested from 26 to 29 days after planting. Plant height, tiller number, root length, total number of roots and percent deep and shallow roots were recorded. Percent deep roots and shallow roots were computed as follows: % deep roots = number of deep roots total root number × 100% % shallow roots = number of shallow roots total root number × 100% Percent deep roots were used as a proxy for root growth angle (Uga, 2012). The shoots and roots were dried for 3 days in an oven at 60 • C and the root dry weight and shoot dry weight were recorded.

Greenhouse Lysimeter Study
A greenhouse lysimeter study was conducted in 2015WS to evaluate deep root growth, water uptake and phosphorus uptake under drought in the IR64-21 × Dular population. Seeds of IR64-21, Dular and the same 300 RILs that were evaluated in the 2014DS lowland field experiment were germinated for 4 days in Petri dishes and then three plants were transplanted into each lysimeter. The lysimeters consisted of PVC cylinders (height: 95 cm, diameter: 20 cm) lined with a plastic sheet and filled with 23 kg of sieved, dry upland soil to a height of 70 cm, with an additional 20 cm of lowland paddy soil added on top of the upland soil. Representative soil samples were sent to the IRRI Analytical Services Laboratory for analysis of Kjeldahl N (0.13%), available P (49.67 mg kg −1 ), Exchangeable K (2.26 meq 100 g −1 ), Mg (7.63 meq 100 g −1 ), Ca (14.87 meq 100 g −1 ), pH (6.1), clay (35.7%), sand (17%), and silt (47.3%). Each genotype was planted in three replicates in a randomized complete block design within the greenhouse and each replicate was arranged within a concrete tank within the greenhouse. Five grams of complete fertilizer (14-14-14 N-P-K) were added to each lysimeter before planting. At 10 days after planting, the plants were thinned to one seedling per lysimeter. Each lysimeter was then sealed at the top with a plastic sheet wrapped around the base of the plant to prevent water loss due to evaporation from the soil surface so that transpiration could be assessed gravimetrically. The drought stress treatment was initiated by removing the rubber stoppers from three holes at the bottom of each lysimeter to drain them at 29, 32, and 34 days after germination (DAG) in replicates 1, 2 and 3, respectively. Weekly weighing of the lysimeters and simultaneous imaging of shoots were initiated at 31, 33 and 35 DAG in replicates 1, 2, and 3, respectively using a mechanical hoist and camera suspended from a gantry crane as described by Kijoji et al. (2012) to determine water uptake and apparent leaf area. During the course of the study, temperature in the greenhouse averaged 33.35 • C; relative humidity averaged 62.74% and mid-day light levels (10-14 h) averaged 602.13 µmol m −2 s −1 . Days to flag leaf appearance and to flowering were recorded for each lysimeter. At 58 DAG, one tiller was sampled from each lysimeter to be used for phosphorus-uptake measurements. Phosphorus uptake per tiller was determined on dried, ground tissue spectrophotometrically using the method described by Murphy and Riley (1962). Harvest of each lysimeter was initiated at 85 DAG based on observations of maturity for each plant. At the time of harvest, the number of tillers was counted, plant height was measured and the shoot was cut to determine shoot dry weight (SDW) and panicle dry weight. The soil was removed from each lysimeter by pulling out the plastic liner sheet and the maximum root depth was determined. The soil segment below the depth of 60 cm was sampled and roots were carefully washed from the soil and dried to determine root dry weight (RDW) >60 cm. Total water uptake (TWU) was determined based on the difference between the initial and final weights of each lysimeter and water use efficiency (WUE) was calculated as SDW/TWU. Weekly water uptake rates were normalized by the apparent leaf area to determine normalized water uptake rates in each lysimeter.

Statistical Analysis of Phenotypic Data
Analysis of variance (ANOVA) was conducted using PBTools V 1.4.0. (PBTools, 2014). Mixed model analysis was conducted to calculate the trait means, which were later used for further analysis. The model used for ANOVA for alpha lattice design was: Where, P ijk is the measurement recorded on a plot, M is the mean over all plots and R, B, L and e refer to replications, blocks, lines and error, respectively. Data of yield experiments for computation of means, heritability and least square difference (LSD) were analyzed using PBtools 1.4.0 taking the effect of replications and block within replications. Broad sense heritability was calculated as: where H is the broad sense heritability of the experiment, σ 2 G is the genetic variance, σ 2 E is the error variance and r is the number of replications in the experiment. Correlation among traits was calculated using STAR V 2.0.1 (STAR, 2014) and was visualized graphically through multidimensional scaling (MDS) using the same program. Briefly, the MDS model accounts for each object or event as a point in a multidimensional space wherein the points are arranged within the space; therefore, the distances between pairs of points have the best fit relation to the similarities among pairs of objects (Wilkinson, 1996).
Multi-environment analysis of the 480 genotypes under eight field experiments was performed using the AMMI-1 (Additive Main Effect and Multiplicative Interaction) stability model in PBTools v.1.4.0. The AMMI model equation according to Gauch and Zobel (1996) for T genotypes and S environments was: where Y ij is the mean yield of the i th genotype in the j th environment; µ is the grand mean; g i is the i th genotypic effect; e j is the j th environment effect; λ n is the eigen value of the PCA axis n; α in and γ jn are the i th genotype j th environment PCA scores for the PCA axis n; θ ij is the residual; n' is the number of PCA axis retained in the model.
To relate seedling-stage roots traits, reproductive-stage root traits, water uptake under drought and yield, a path analysis was conducted using lavaan script in R v. 3.3 using the subset of 304 IR64-21 × Dular RILs. Parameters measured in multiple experiments or on multiple dates were averaged, including days to emergence, total root number, biomass at harvest, canopy temperature and yield under drought stress.

Genetic Analysis
Leaf Collection, DNA Extraction, PCR and PAGE The genotyping work was carried out at IRRI's Genotyping Services Lab (GSL). Eight fresh leaves from one replication of each genotype in 2014DS UNS were collected in bulk and immediately kept on ice. The leaves were collected in glassine bags with proper labeling. A portion of these leaves was placed in individual tubes, ground and stored at −20 • C for later processing. The DNA from all 480 genotypes and the two parents were extracted using a modified CTAB protocol developed by Murray and Thompson (1980). The DNA was dissolved in 200 µL of TE (Tris-EDTA) buffer. DNA solutions were stored at −20 • C. Polymerase chain reaction (PCR) was performed using 15-µL reactions. Amplifications were performed using 40 ng of DNA, 1×PCR buffer, 100 µM dNTPs, 250 µM oligonucleotide primers and 1 unit of Taq polymerase. This mix was prepared in 96-well polycarbonate plates and the thermocycling procedure followed the method described by Panaud et al. (1996). Amplified products were resolved using high-resolution 8% (v/v) polyacrylamide gel electrophoresis (PAGE) (CBS scientific, model MGV-202-33) as described by Sambrook et al. (1989). The gel was run in a 1xTBE buffer at 95 volts for 1-3 h, depending on the SSR marker product size. DNA fragments were stained with SYBER Safe TM and visualized using a UV trans-illuminator.

Parental Survey for Polymorphism and Whole Population Genotyping
A total of 400 rice SSR markers from available rice genetic sequence maps were used to perform a polymorphic survey between IR64-21 and Dular. The details of the primers such as marker position, chromosome position, expected product size and annealing temperature were taken from Gramene Genome Browser (www.gramene.org). PCR, PAGE and gel documentation of the images were performed as described earlier. The bands were scored as either 1 or 2 based on their banding pattern. A total of 115 polymorphic markers spread evenly across the genome were selected and used for generating the genotyping data of the whole population including the parents. Primers designated to each marker were used in PCR and the subsequent PAGE and gel documentation were performed as described earlier.

Linkage Map Construction and QTL Analysis
The entry means of all phenotypic traits were correlated with the genotypic data of the respective lines using two software programs; QTL Network 2.1 (Yang et al., 2008) and Windows QTL Cartographer version 2.5 (Wang et al., 2012). Two software programs were used due to the different analysis tools available in each program and also to identify consistent and stable QTLs identified by both programs. A genetic locus was confirmed only if consistent with both programs. The QTL Network software identifies putative regions within the QTLs based on onedimensional genome scan considering chosen candidate intervals as cofactors. A mixed linear model framework was used to perform the mapping procedure. F-statistics based on Henderson method III were used for hypothesis testing. The critical F-value was calculated based on 1,000 permutation tests to minimize the genome-wide type I error. A window size of 1 and a walk speed of 1.0 cM were used for the whole genome scan. A significant level of p < 0.001 was maintained for QTL detection in the whole experiment. The linkage map was constructed and analysis was then repeated using the QTL cartographer software (Wang et al., 2012). Composite interval mapping (CIM) with default parameters (300 permutation time, 5% of significance level, model 6; standard model, method 3; forward and backward method with walk speed 2 cM) was performed; LOD explained by each QTL was calculated. The flanking markers intervals and positions, peak marker and position (P) and additive effect (A) of the QTLs was calculated using QTL network and also confirmed with Windows QTL Cartographer version 2.5; the logarithm of the odds (LOD) score of the QTLs was calculated using Windows QTL Cartographer version 2.5 and h 2 (heritability) using QTL network 2.1. As the absolute additive effect changes with severity of stress , the additive effect was presented as a percentage of the population mean as: %A (percent additive effect) = A (additive effect) ×100/Population mean.

Identification of Candidate Genes and other Colocating QTLs
Functionally characterized genes within the indentified QTLs as well as previously identified QTLs colocating with the identified QTLs in this study were determined using Q-TARO (QTL Annotation Rice Online, http://qtaro.abr.affrc.go.jp/cgibin/gbrowse) database based on the physical position of both primers flanking the QTL . The chromosome number and the genomic region were entered in the genome browse panel using the format Chr#:genome start..genome end. For example, Chr1:8,000,000.10,000,000 was entered to indicate that the physical position 8,000,000-10,000,000 bp in chromosome 1 was the region inquired for candidate gene and other colocating QTL identification.

Yield and Related Trait Characterization
The mean performance of the 480 genotypes and the checks (Dular and IR64-21) is presented in

Effect of Drought Stress on Yield and Related Traits
Under upland reproductive-stage drought stress conditions, the average grain yield reduction (GYR) of Dular and IR64-21 was observed to be 73 and 99%, respectively, while the RIL population exhibited a GYR of 71-100%. A GYR of 42 and 90% were observed for Dular and IR64-21, respectively, under lowland reproductive-stage drought while a GYR of 0-99% was observed in the RIL population (Table 4). Under upland stress conditions, DTF of the RIL population was delayed by 1-36 days while it was only delayed by 0-9 days under lowland stress conditions. Plant height of the RIL population was reduced by 4-67 cm in upland stress and 9-78 cm in lowland stress compared to non-stress conditions. The biomass and harvest index of Dular decreased by 32 and 48%, respectively, while it decreased by 50 and 78%, respectively, in IR64-21. A reduction of 4-84 and 0-87% to the biomass and harvest index, respectively, was observed in the RIL population under reproductive-stage drought stress compared to non-stress.

Correlation among Traits
Correlations among traits were calculated using STAR V 2.0.1 and were visualized graphically through MDS using the same program (Figure 1). Formation of clusters indicates good correlation among the clustered traits. The closeness of the same trait measured in different environments indicates less effect of the environment on these traits. MDS analysis showed the clustering of traits; grain yield, root traits, harvest index and number of tillers in one cluster (I) while plant height, biomass, growth rate and leaf rolling were in another cluster (II). A third cluster (III) was formed by DTF, canopy temperature, bleeding rate and early vigor. In the path analysis relating traits measured at germination stage, seedling/vegetative stage, reproductive stage and harvest, correlation coefficients ranged below 0.36 (Figure 2). The strongest correlations between growth stages were between total water uptake and biomass at harvest, canopy temperature and yield, seedling stage shoot dry weight and root dry weight below 60 cm, and root growth angle and root dry weight below 60 cm. The strongest paths linking seedling stage to harvest were via (a) seedling stage shoot dry weight which was related to reproductive-stage biomass and total water uptake, which were related to yield under drought, and (b) root growth angle, which was related to root dry weight below 60 cm, which affected total water uptake and was related to yield under drought.

Identified QTLs for Grain Yield and Other Traits
Stable, consistent-effect QTLs with significant p-values were identified; three for grain yield on chromosomes 1 and 8 (Figure 3) and four for root traits on chromosomes 1, 5, 8, and 9 (Figure 4). A hotspot genetic locus was identified as qDTY 1.1 on chromosome 1 where multiple QTLs for grain yield, plant height (qPHT 1.1 ), biomass (qBIO 1.1 ), growth rate, leaf rolling (qLR 1.1 ), mesocotyl length (qRT 1.1 ), bleeding rate, water-use efficiency (qWUE 1.1 ) ( Table 5) and phosphorus uptake (qPU 1.1 ) were detected across years and environments between markers RM11943 and RM3482 (Figure 4). The novel grain yield QTLs, qDTY 1.3 and qDTY 8.1 with heritability values ranging from 3.1-5.1 and 4.4-7.6, respectively, and additive effects of 9.0-43.7% and 11.6-39.7%, respectively were identified. The grain yield QTL qDTY 1.3 on chromosome 1 between RM292 and RM11013 was consistent under both non-stress and reproductive stage drought stress under upland conditions. The grain yield QTL qDTY 8.1 on chromosome 8 between markers RM80 and RM230 showed consistent effects under both nonstress and reproductive-stage drought stress, qRT 9.1 on chromosome 9 between RM434 and RM257 was identified for root growth angle (Figure 4). QTL qRT 9.1 showed a heritability of 29%, an additive effect of 16.7 and 15.4% for % deep roots and % shallow roots, respectively. The percent deep root QTL was contributed by Dular, whereas the % shallow root QTL was contributed by IR64-21 (Table 5). A genetic region of 9.9 cM on chromosome 5 was associated with root traits including root length and root dry weight measured in the basket study and crown root number at the seedling stage measured in 2014DS URS (Figure 4) with heritability of 5.8, 7.9, and 3.5%, respectively ( Table 5).

Physiological Response of Selected Genotypes
In comparison with IR64-21 in the greenhouse lysimeter experiment, the most stable and highest yielding IR64-21 × Dular RILs ( Table 6) showed generally higher root dry weight at depth   (Supplementary Figure 2A) but similar maximum root depth (Supplementary Figure 2B) and higher WUE (Supplementary Figure 2C) but similar water uptake rates (Supplementary Figure  2D). In the field under drought, the most stable and highest yielding IR64-21 × Dular RILs showed higher percentages of the total root length as lateral roots (Figure 5B), lower sap bleeding rates ( Figure 5C) and lower canopy temperatures ( Figure 5D) than IR64-21, but similar percent deep roots ( Figure 5A).

Effect of Identified QTLs on Grain Yield
Compared to genotypes without qDTY 1.1 , the genotypes with qDTY 1.1 showed grain yield improvement (GYI) of 2.9% in FIGURE 4 | QTL likelihood curve of LOD score for root traits QTLs identified on chromosomes 1, 5, 8, and 9 using QTL cartographer version 2.5. FS, Field study (2014DS_URS); BS, basket study; LYS, lysimeter study. LOD graph generated.
2013DS and 55% in 2014DS under URS (Table 7), reflecting the higher severity of stress in 2013DS, in which observed yield reductions of the population ranged up to 99% ( Table 4).
The genotypes with qDTY 1.1 showed GYI of 70.8 and 1.6% under lowland non-stress and reproductive-stage drought stress, respectively ( Table 7). The genotypes with qDTY 1.3 , however, showed higher GYI in 2013 upland with 108.9%, than in 2014 upland with 24.1% GYI. Similar results were observed under lowland conditions (Table 7). Grain yield improvements of genotypes with qDTY 8.1 over the lines without qDTY 8.1 during 2013 and 2014, were 97.4 and 74%, respectively, in upland and 3.0 and 16.2%, respectively, in lowland conditions ( Table 7). The QTL combination qDTY 1.1 + qDTY 8.1 showed the highest GYI, followed by qDTY 1.3 +qDTY 8.1 and qDTY 1.1 +qDTY 8.1 +qDTY 1.3 across upland/lowland reproductive-stage drought stress. qDTY 1.3 and qDTY 8.1, individually and combined, showed a high GYI even under severe drought stress, while qDTY 1.1, alone, resulted in high GYI under non-stress and medium stress but not under severe upland stress.

DISCUSSION
With the root and drought response traits measured in this study from a large population (>300 RILs) at times ranging from seedling stage to harvest, we aimed to (1) integrate the relationships among traits at different growth stages in terms of their effects on yield and (2) identify the genetic loci controlling the traits that have an effect on yield under drought. We hypothesized that, although traits measured at seedling and vegetative stage are chronologically distant from the harvest; some key intermediate traits measured over the growing seasons could link seedling-stage traits with grain yield under drought. Given the wide differences between IR64-21 and Dular in terms of their phenotype and their response to reproductive-stage drought stress, we had the opportunity to explore and observe which traits have an effect on yield under drought and which combinations of QTLs were most effective under drought in the RIL population. Delayed flowering, decreased plant height and reduced grain yield are some of the effects of drought stress (Lafitte et al., 2004;Zhao et al., 2010;Vikram et al., 2011;Kumar et al., 2014). The reduction in grain yield upon imposition of stress on the populations was 71-100% in the upland experiments and 0-99% in the lowland experiments indicating the severity of stress ( Table 4). Different levels of stress among experiments are always desirable to eliminate the effect of yield potential and to select better drought-resistant germplasm (Bernier et al., 2007). Under reproductive-stage drought stress, Dular outperformed IR64-21 supporting the suitability of breeding efforts involving Dular as a source of genetic loci that enhance grain yield under drought stress.
Clustering of root traits with grain yield in the MDS analysis signifies the role of root traits in maintaining water or nutrient uptake under reproductive-stage drought stress conditions in this population. According to the path analysis, the strongest links between early-and late-season traits were from seeding stage shoot dry weight or seedling stage root growth angle (Figure 2). The identification of a relationship between seedlingand vegetative-stage root traits from the greenhouse experiments (root growth angle measured in the basket experiment and root GY, grain yield; GR, growth rate; PHT, plant height at maturity; LR, leaf rolling; ML, mesocotyl length; DR, % deep roots; SR, % shallow roots; RL, root length; RDW, root dry weight; CRN, crown root number; WEU, water use efficiency; PU, phosphorus uptake; DS, dry season; UNS, upland non-stress; URS, upland reproductive-stage drought stress; LNS, lowland non-stress; LRS, lowland reproductive-stage drought stress; GH, green house experiment; P, position of the peak of the QTL; %A, additive effect of the QTL as percent of the population mean. dry weight below 60 cm in the lysimeter experiment) rather than from the field may have been due to the different root growth components measured in the two types of environments (crown root numbers vs. angles) as well as the high degree of variability inherent to field root studies. This study identified QTL qRT 9.1 for root growth angle, which is adjacent to the previously reported QTL qGY 9.1 , and qGY 9.2 for grain yield and qEVV 9.1 for early vegetative vigor in the Aus276/3 * IR64 population under direct-seeded conditions (Sandhu et al., 2015), for root length (Price et al., 1999(Price et al., , 2002Uga et al., 2010;Sandhu et al., 2013), root thickness (Steele et al., 2006(Steele et al., , 2013, and root number (Li et al., 2005). Furthermore, meta analysis from a 675-root QTL database from 12 populations (Courtois et al., 2009) showed the presence of root architecturerelated QTLs colocated or adjacent to qRT 9.1 . Although the effect of the identified genetic regions requires a proper validation in different genetic backgrounds and environments, under the present climate change scenario, successful introgression following marker-assisted breeding may be used to improve the root system architecture of the widely cultivable popular rice varieties. The contribution of deep-and shallow-root enhancing alleles by Dular and IR64-21, respectively, signifies the role of both parents in maintaining the nutrient-water balance in the promising lines, as IR64-21 is a lowland-adapted parent with shallower root growth and Dular is rainfed-adapted parent with a deeper root system. The genomic locus qDTY 1.1 stood out as a hotspot for grain yield, agronomic and physiological traits. The major, consistent and large effect of qDTY 1.1 on grain yield under non-stress, reproductive-stage drought stress, and direct seeded conditions has been reported previously in different backgrounds, Swarna (with additive effect of 13%), IR64 (with additive effect of 8, 11, 15, and 24%) and MTU1010 (using different donors such as N22, Dhagaddeshi, Aus276 and Kali Aus (Vikram et al., 2011;Sandhu et al., 2014Sandhu et al., , 2015Sandhu et al., , 2016. This region has also been observed to be associated with plastic responses to drought stress, including enhancing deep root growth and shoot growth regulation (Vikram et al., 2015;Wade et al., 2015). qDTY 8.1 from this study was observed to be located adjacent to qGY 8.1 , identified in an Aus276/3 * IR64 population under direct-seeded conditions, validating the effect of drought-tolerant loci in the IR64 background across variable environmental and cultivation conditions. Similarly, the genetic locus qRT 5.1 (with additive effect of 3.6, 8.9, and 4.5%), associated with root traits on chromosome 5 in the present study, was located in the same region previously associated with QTLs for nutrient uptake (with additive effect of 12 and 23%) and root traits (with additive effect of 8 and 6%) under direct-seeded conditions (Sandhu et al., 2015) in Aus276/3 * IR64 population. This indicates the possibility of presence of conserved allelic region in different donors and help to harness the potential of multiple donors enhancing water-nutrient uptake under stress conditions. This colocation signifies the relationship of root traits with nutrient uptake in improving grain yield under drought and direct-seeded conditions. The identification of qDTY 1.1 , qDTY 8.1 , qRT 5.1 , and qRT 9.1 in earlier studies clearly demonstrates the stability of these QTLs across different genetic backgrounds and the identification of a novel QTL, qDTY 1.3 , presents an alternate resource for introgression of drought tolerance into susceptible varieties.
The better performance of yield-stable genotypes across variable growing environments appears to be largely attributed by drought tolerance-related root traits ( Table 6). Longer mesocotyls were expected to result in earlier time to emergence and a greater number of crown roots at the seedling stage, however, only the latter was confirmed in this study. The observation that the most stable and highest yielding IR64-21 × Dular RILs differed from IR64-21 in some (particularly % lateral roots) but not all deep-rooting traits (such as maximum root depth) suggests that the differences in root architecture may be more related to lateral root growth rather than nodal root elongation. Furthermore, the differences in traits related to water uptake between IR64-21 and the most stable and highest yielding IR64-21 × Dular RILs indicate that the drought response of these lines is conferred by a combination of root architecture and root functional traits. This observation is in agreement with previous reports that Dular contrasts with IR64-21 in terms of root hydraulic parameters under drought, including xylem cell diameter, root hydraulic conductivity, sap bleeding rate and aquaporin expression (Henry et al., 2011). These results are also supported by the MDS analysis in which root traits clustered with grain yield under drought stress.
Previously identified genes and other QTLs in the QTL regions identified in the study may give insights as to why these QTLs confer grain yield enhancement and stability under drought stress conditions (Supplementary Tables 1,  2). Three functionally characterized genes were identified within qDTY 1.3 , namely S-Adenosyl-l-methionine synthetase3 (OsSAMS3), S-Adenosyl-l-methionine synthetase2 (OsSAMS2) and photoassimilate defective1 (phd1) (Li C. et al., 2011;Li W. et al., 2011). OsSAMS3 and OsSAMS2 affect fertility, germination rate and flowering time in rice while phd1 affects photosynthetic activity and biomass and grain production. QTLs identified in this study colocate with previously identified QTLs associated with root and grain yield-related agronomic traits such as root branching index and number of filled grains per panicle (Zhuang et al., 2001;Li et al., 2005;Horii et al., 2006) Figure 3). Within the QTL region qDTY 8.1 on chromosome 8 identified in the present study, the gene rice authentic Hiscontaining phosphotransfer 1, associated with lateral and crown root development, osmotic adjustment and seed setting under drought, was identified in knockdown experiments (Sun et al., 2014). Other genes, such as OsMADS7 that controls floral organ formation in rice and proton gradient regulation 5 that controls the plant's photosynthetic capacity (Cui et al., 2010;Nishikawa et al., 2012), were also identified to be within qDTY 8.1 . Drought tolerance genes, such as dehydration-responsive element-binding transcription factor 1G, heat shock factor class B 2b (Chen et al., 2008;Xiang et al., 2013) and QTLs for root-to-shoot ratio, maximum root length, relative germination vigor and filled grain weight per plant, have also been reported within qDTY 8.1 (Zhuang et al., 1997;Price et al., 2002;You et al., 2006).

CONCLUSION
The identified genetic regions associated with grain yield, yield attributing and root development traits under reproductivestage drought stress and the identified promising lines in this study may facilitate future marker-assisted QTL pyramiding/introgression breeding programs to improve rice yield under drought conditions. Fine mapping, precise introgression of genetic regions with positive interactions and evaluation of consistent performance across environments and genetic backgrounds may further improve grain yield and may provide insight into the physiological and molecular characterization of genetic loci.

AVAILABILITY OF DATA AND MATERIALS
The data sets supporting the results of this article are included within the article.

AUTHOR CONTRIBUTIONS
MC was involved in conducting the experiment, recording observations and drafting the article; NS and SD was involved in experimental analysis, interpretation of data and revising the manuscript; NAAS was involved in conducting the experiment and recording observations; MN and KM developed the RILs; AH was involved in designing the greenhouse experiment and the critical revision of the manuscript; MD and KM were involved in revision of the manuscript; and AK conceived the study and contributed to the critical revision of the manuscript. All authors approved the final version of the manuscript.