Mapping consistent additive and epistatic QTLs for plant production traits under drought in target populations of environment using locally adapted landrace in rice (Oryza sativa L.)

ABSTRACT A total of 36 colocated QTLs associated with plant production traits were mapped using IR62266/Norungan recombinant inbred lines under drought and non-stress conditions in target populations of environment. Multi-model composite interval mapping and multi-environment analysis detected seven consistent QTLs for plant production traits on chromosomes 1, 2, 3, 6, 10 and 11. QTLs for plant height, number of tillers, biomass, spikelet fertility and grain yield under drought identified on chr. 1 (RM5389-RM11943-RM5794), chr. 2 (RM324-RM5390 and RM12460-RM423-RM5345), chr. 3 (RM6329-RM5475-RM16030), chr. 6 (RM508-RM585-RM217), chr. 10 (RM7361-RM8207) and chr. 11 (RM209-RM6499) in this study have been reported as meta-QTLs for yield under drought. Consistent positive allelic contribution from the landrace, Norungan was shown at the two QTLs on chr. 2, while IR62266 showed consistent positive allelic effect at the QTL on chr. 10. At other four regions the effect was confounded with the stress condition, for instance the region on chr. 3 identified for straw yield showed additive effect from Norungan in all stress trials, while in the non-stress trials the positive allele came from IR62266. Epistasis for grain yield was detected in the regions of non-additive effects and its relative contribution was small. Relatively small percentage of grain yield variation was explained by additive x environment interaction (0.28 - 1.77%), but was larger than that explained by additive effect QTL across environments (0.05 - 0.35%). These consistent QTLs may have genes evolutionarily conserved in response to drought and could be ideal candidates for rice yield improvement in water-limited environments. GRAPHICAL ABSTRACT


Introduction
Rice (Oryza sativa L.) is a staple food crop for more than half of world's population, cultivated over 167 million ha globally (Food Agriculture Organization (FAO), 2018). Rainfed lowland and upland rice occupy 31 and 12% of the global rice-growing area, respectively (Kamoshita et al., 2008), where yields are seriously affected by water deficit. Drought is thus the major environmental constraint to rice production and yield stability in rainfed areas. Asia is the most drought-affected region from 2003 to 2013, with 40% of total crop and livestock production losses amounting to 28 USD billion (Food Agriculture Organization, 2015). In Asia nearly, 46 Mha rainfed lowland and 10 Mha upland rice suffer from drought stress (Pandey et al., 2007). In India alone, 14.4 and 6.0 Mha of rainfed lowland and upland rice are affected by drought, respectively (Mahajan et al., 2017).
Drought, in general reduces the accumulation of biomass and shortens the life cycle of the plants (Blum, 2005) causing yield reduction ranging 13-35% (Jongdee et al., 2006). The recent predictions in global climate change suggest a further increase in water deficits in coming decades (Wassmann et al., 2009), leading to an increased frequencies and intensities of drought. It has been reported that by 2025, the irrigated wetland rice (13Mha) in Asia may experience physical water scarcity, whereas the irrigated dry-season rice (22Mha) may suffer from economic water scarcity (Tuong & Bouman, 2003).
Rainfed rice yields an average of 1.5 t/ha compared with 3.1 t/ha in irrigated rice (Rosegrant et al., 2002). Low yields, large yield gap, risk of crop failure, and insecurity of livelihoods are predominant in rainfed rice ecosystems (Hossain & Narciso, 2004), which in turn increases the economic vulnerability of poor rice producers in less-favored rainfed environments, where more than 30% of the population is already considered extremely poor, living on less than USD 1.90 per day (Mottaleb et al., 2017). The development of rice cultivars with inbuilt drought resistance is thus important in increasing productivity and yield stability as well as alleviating poverty in communities' dependent on rainfed production (Venuprasad et al., 2009). Sufficient progress in breeding for drought tolerance has been made in recent years, however, is still limiting with a lack of reliable and efficient screening techniques (Mohanty et al., 2016). Direct selection for yield under stress in managed stress environments (MSEs) (Babu et al., 2003;Venuprasad et al., 2007) and target populations of environment (TPEs) (Kumar et al., 2008) is considered promising to improve drought tolerance in rice. A TPE is the set of all environments, farms, and future seasons in which an improved variety will be grown (Fischer et al., 2012). However, direct selection for yield under drought in MSE and TPEs is difficult because of differences in the timing and severity of drought over seasons (Prince et al., 2015). Alternatively, quantitative trait loci (QTL) mapping followed by marker-assisted breeding (MAB) could be an effective approach to identify genomic regions linked to crop performance in stressful environments and pyramiding desirable alleles to improve drought resistance in crops (Ashraf, 2010). Large number of QTLs associated with several droughtresistance traits (Kamoshita et al., 2008) and yield under drought have been identified in rice, but only few have been successfully deployed in MAB till date (Sandhu & Kumar, 2017), due to large intervals and lack of repeatability across environments and genetic backgrounds (Bernier et al., 2007). Identifying large consistent effect meta-QTLs for yield and secondary traits under drought will hasten MAB for drought resistance.
India is one of the centers of diversity of rice (McNally et al., 2009), harboring drought-tolerant landraces adapted to local ecosystems. These landraces endowed with genes of agronomic value such as drought tolerance offer great scope in increasing genetic variability in breeding program (Vikram et al., 2016). Populations involving local landrace of indica ecotypes and indica accessions will results in superior indica lines with adaptation to local target rainfed environment. Novel QTLs for yield and yield components can be identified by utilizing genetic variation among indica ecotypes (Suji et al., 2012). Recently, fewer QTLs for drought resistance have been detected using populations derived from landraces in managed stress environments (Sandhu & Kumar, 2017). However, identifying consistent meta-QTLs for drought resistance using locally adapted landraces under drought stress predominant in TPE is essential  due to heterogeneity in soil physico-chemical properties and drought scenarios across different TPEs. Identifying consistent QTLs for drought resistance, especially using populations derived from locally adapted rice lines (Atlin et al., 2006) will improve the efficacy of MAB.
Previously, using same mapping population we have identified genomic regions that can be called as QTL hot spots and fine mapped one large effect QTL for yield under drought (Suji et al., 2012;Muthukumar et al., 2015). However, the average distances between marker were fairly large with several chromosome breakpoints rendering them less practical for any downstream marker trait association (Pootakham et al., 2015) and usually the hotspot regions undergo more recombination than the other regions (Smith, 1994). The utility of the linkage map information is often limited to the genetic background of the mapping population (Galeano et al., 2011). Further to find genes of interest and marker-assisted selection through fine mapping, map-based cloning and comparative genomics study requires good genetic map with greater genome coverage (Galeano et al., 2011) providing validation of marker order.
Further, the consistency of the identified genomic regions, epistasis interaction and QTL x Environment interaction under different stress are important in dissecting complex quantitative trait like drought resistance beside elucidating the contribution of single gene/QTL. The genetic architecture of a complex trait consists of the actions of genes at a single locus as well as the interlocus interactions and gene x environment interactions (Qi et al., 2016).
Thus, the present study was conducted with an objective of saturating the genetic linkage map and identifying consistent large-effect QTLs, to detect additive and epistatic effect QTLs and their interaction with environment for plant production traits under drought using rice lines derived from local landrace, Norungan adapted to TPE. The landrace, Norungan is traditionally grown in rainfed TPE over considerable area by subsistence farmers because of its inherent adaptation to drought and flooding. Identifying consistent QTLs for drought resistance in this landrace will help in rainfed rice improvement.

Mapping population
Norungan (NG), an indica landrace well adapted to rainfed ecosystem from Tamil Nadu, India is drought resistant with its deep and thick roots with superior hardpan penetration ability (Babu et al., 2001) IR62266-42-6-2 (IR62266), a lowland indica ecotype, efficient in extracting water from depth due to capacity for osmotic adjustment (Babu et al., 1999;Siopongco et al., 2009) and has a thin root system (Babu et al., 2001). These two lines were used as parents in developing a recombinant inbred population. A subset of 132, out of the total 232 recombinant inbred (RI) lines from the cross IR62266 x Norungan was used in the present study for mapping consistent QTLs for phenology, plant production and yield traits under drought stress in TPE and non-stress environment.

Field trials
Field trials were conducted under upland conditions in experimental fields of Tamil Nadu Agricultural University, India at Agricultural Research Station, Paramakudi, India during rainfed seasons (September-January) 2013-14 (DS2013) and 2014-15 (DS2014) located in rainfed TPE. The site, soil, and drought stress characteristics, and rainfall of the experimental locations are summarized in Table 1 and Supplementary Figure 1, respectively.

Drought stress experiments in TPE (DS2013, DS2014)
A subset of 143 F 10 RI lines along with parents was evaluated for drought response. The rice lines were planted in plots of 2.0 × 0.4 m 2 with three replications in randomized complete block design. Briefly, seeds were hand dibbled into the dry soil before monsoon with a seed rate of 80 kg ha −1 . NPK fertilizers were applied at the rate of 50:25:25 kg ha −1 . Insect and weed control measures were applied periodically, as required. Stress plots were completely rainfed from sowing to harvest.
The DS2013 trial received a total rainfall of 217 mm during the crop growth as against long-term average rainfall of 475 mm for this site. There was no rainfall right from emergence up to harvest and the rice lines experienced very severe drought stress. Due to very acute drought, the bore wells dried and no irrigation could be given even to save the plants. During DS2014 the rainfall was normal with total rainfall of 483 mm, evenly distributed during the cropping season and was considered normal rainfed year with no visual symptoms of drought stress.
Data on days to 50% flowering (DFF) were recorded as the number of days from emergence to flowering in 50% plants in a plot for each RI line. Data on plant height (PH), number of total (NTT) and productive tillers (NPT), panicle length (PL), and spikelet fertility (SF) (ratio of number of filled grains/total number of grains (filled + unfilled (chaff grains)) per panicle expressed as percentage and yield per plant were recorded from three randomly selected hills per RI line in each replication. All the plants in each plot were harvested at maturity and grain yield (GY) and straw yield (SY) were recorded after sun drying. Total above-ground biomass (TB) was computed by summing grain and straw yields. Harvest index (HI) was calculated as the ratio of grain weight to total above ground biomass for each RI line.
Drought score (DS) was recorded during peak stress in DS2013 trial using a scale of 0-9 based on the intensity of leaf drying, leaf rolling standardized for rice (IRRI (International Rice Research Institute), 1996). The soil water table depth was measured at two days' interval using eight piezometers (15 x 100 cm) made of polyvinyl carbonate that were installed across the field in DS2014. Due to acute drought in DS2013, piezometers could not be installed in dry and hardened soil.

Statistical analysis of phenotypic data
Analysis of variance for each trait in different trials was done as mixed models using PROC MIXED procedure of SAS V. 9.4 (SAS Institute Inc, 1990). The broad sense heritability (H) of different traits was estimated within a year for each treatment. The required variance components for calculating heritability were obtained as explained by Fehr (1987). The SAS programme PROC VARCOMP with REML method was used for broad sense heritability calculation. The linear relationship, correlation, and summary statistics for each trial were calculated using PROC CORR and PROC SUMM programme in SAS.

Genotyping
Fresh leaf samples of each RI lines and parents were collected from field-grown young seedlings and freeze-dried. DNA was extracted from the leaf tissues using cetyl trimethyl ammonium borate buffer (Gawel & Jarret, 1991).
The quantity and quality of DNA was assessed in 0.8% agarose gel and quantified by using NanoDrop® ND-1000 Spectrophotometer. Sample concentration values in ng/µl based on absorbance at 260 nm were determined. The DNA concentration was adjusted to 50 ng/µl from the spectrophotometer absorbance values. Polymerase chain reactions (PCR) were performed in a volume of 15 µl in thermal cycler (BIO-RAD, USA and Eppendorf Master Cycler Gradient, Eppendorf, Germany). The reaction mixture contained each primer of 1 µM concentration (Sigma Aldrich, USA and Eurofins Genomics, Bangalore, India), 100 µM deoxy nucleotide, 1X Taq buffer, 0.02 U Taq polymerase (Bangalore Genei, India) and 50 ng of template DNA. After 5 min of initial denaturation at 94°C, the PCR involved 36 cycles of amplification, each cycle comprising 1 min denaturation at 94°C, 1 min annealing at 57°C, 1 min extension at 72°C, and the final extension step of 5 min at 72°C. The amplified products were separated on 3% Metaphor agarose gels and 6% PAGE (Bio-Whittaker Molecular Applications, Vallensbaek Strand, Denmark) (Sambrook & Russell, 2006).

Parental polymorphism and Genetic Linkage map construction
A framework genetic map comprising 114 markers was constructed previously in this laboratory using the total 232 RI lines of this mapping population (Muthukumar et al., 2015;Suji et al., 2012). In the present study, the parents, IR62266 and Norungan, were further screened using 388 additional rice microsatellite markers and 96 were found to be polymorphic between the parents. These 96 polymorphic markers were used in saturating linkage map and in genotyping the subset of 132 RI lines used in this study. The genotypic data generated with 132 RI lines were tested for χ2 goodness of fit against a 1:1 segregation ratio and was used in reconstruction of the linkage map with additional markers using their physical position in reference to rice variety, Nipponbare. The linkage map was reconstructed using Map Manager QTX software (Manly et al., 2001) using the Haldane mapping function. The reconstructed map covers a total length of 1408.02 cM with an average marker distance of 6.70 cM (Supplementary Figure 2).

QTL analysis
QTL analysis for five trials was carried out with the line mean of each trait using the composite interval mapping (CIM) approach in WinQTLCART v2.5 (Wang et al., 2005). Genome-wide significant threshold at p < 0.05 was calculated by running 1000 permutation from which logarithm of the odds ratio LOD score of ≤2.5 was observed for trait analysis. Model six was selected with control marker numbers (Cofactors) of 5 and window size of 10 cM. The forward regression method was used for selecting the control markers. QTLs were searched with walk speed of 1 cM. The additive effect and phenotypic variance explained by each QTL (R 2 ) at maximum LOD score is reported. QTLs with LOD score of 2.5 and above are reported. The data from (DS2011, DS2012 and NS2012) Muthukumar et al. (2015) was reanalyzed for QTL detection with possibility of new QTLs detection using improved genetic map.

Multi environment, epistasis and QTL interaction analysis
Multiple environment trial (MET) analysis was carried out using all five field trials data, three field trials (two drought stress (DS2011, DS2012) and one non-stress (NS2012)) from previous study (Muthukumar et al., 2015) and two field trials (DS2013, DS2014) from current study. Multi-model composite interval mapping (MCIM) of QTLNetwork2.1 (Yang et al., 2008) was performed and stable QTLs across environments were identified. Epistatic QTL and QTL interaction was performed based on methodology summarized by Yang et al. (2007) following 2D genome scan and running 1000 permutations.

Phenotypic variation in plant phenology and production traits
Significant variation in phenotypic response was observed among the RI lines and parents for various traits under drought stress in TPE and non-stress (irrigated) conditions. The phenotypic mean and range for various traits along with standard deviation and broadsense heritability for the field trials are presented in Table 2. The traits observed showed a normal distribution among the lines in all the trials. Transgressive segregation was observed among RI lines over the parents for most traits (data not shown). Norungan flowered earlier than IR62266 in DS2013 and DS2014 trials. The reduced plant height and an average drought score of 4.6 in DS2013 indicated the intense drought severity as compared to DS2014. Mean grain yield across RI lines was 966 and 3964 kg ha -1 in DS2013 and DS2014, respectively. Mean grain yield under drought stress was reduced by 75% in DS2013 compared to DS2014.

Phenotypic correlation among traits
Phenotypic correlations among traits for each trial are presented in Supplementary Table 1. In all the trials, grain yield has significant positive correlations with plant height, number of total tillers and productive tillers, panicle length, spikelet fertility, straw yield, total biomass, and harvest index. Grain yield was strongly correlated with harvest index (0.70 in DS2013) and total biomass (0.71 in DS2013) under drought stress. However, grain yield was negatively correlated (−0.26) with days to 50% flowering under drought stress in DS2013.

Genetic map construction and QTLs mapped for plant production traits
To saturate the existing genetic linkage map constructed with 114 SSR markers (Muthukumar et al., 2015), a set of 338 SSR markers spread throughout the genome were additionally used in parental polymorphism screening in this study. Out of 388 SSR markers, 96 were polymorphic (24.8% polymorphism). Hence, a total of 210 markers, distributed over 12 chromosomes, were used in reconstruction of the linkage map in this study.
A total of 68 genomic regions were identified for various traits under drought stress and non-stress conditions. Results of CIM and putative QTLs identified under stress and non-stress conditions with peak or nearest marker(s) linked to various traits are given in Table 3. The QTLs explained phenotypic variance ranging from 6.05% to 54.66% in drought stress trials and 4.72 to 53.33% in nonstress trial across various traits measured.

QTLs for DFF
Eleven genomic regions were identified on chromosome 2, 3, 4, 6, 9, and 11 by using CIM (Table 3). These QTLs were detected in DS trials with large phenotypic variance ranged Table 2. Mean, range, and broad sense heritability of 132 RI lines and parents for phenology and plant production traits under field trials conducted in rainfed target populations of environment in Paramakudi (DS2013 and DS2014), India.  Table 3. QTLs identified for phenology and production traits in IR62266/Norungan RI lines field trials conducted in rainfed target populations of environment in Paramakudi (DS2011, DS2012, DS2013, and DS2014) and irrigated non-stress environment in Coimbatore (NS2012), India.

QTLs for PH
Eleven genomic regions were linked to plant height on chromosome 1, 5, 6, 7, and 10 in both stress and nonstress conditions ( Table 3). The phenotypic variance was ranged from 4.72 to 53.33%. The large phenotypic value (53.33% in NS2012, 51.8% in DS2014) was at RM11943 on chromosome 1 with positive allele from IR62266. Four regions, RM5389 (DS2011, DS2012), RM414 (DS21014, NS2012) RM11943 (DS2014, NS2012) on chromosome 1 and RM18828 (DS2014, NS2012) on chromosome 5 were consistent for two trials. The positive allele at peak marker RM5389, RM414, and RM18828 was coming from Norungan. The region with peak marker RM5389 was identified in drought stress while RM414 and RM18828 were identified in a non-stress environment shows they are environment responsive.

QTLs for SF
A total of 10 major effect regions (R 2 > 10) were identified to be linked with spikelet fertility on Chromosome 1, 2, 6, and 7 (Table 3). Four QTL regions on chromosome 2 and 6, and one each on chromosome 1 and 7. All the four QTL regions on chromosome 2, RM5345, RM318 in DS2012 and RM3501, RM324 in DS2014 showed allelic effect from IR62266 while all the four regions on chromosome 6, RM6917, RM8270, and RM5957 in DS2014 and RM585 in DS2013 on chromosome 6 showed allelic effect from Norungan.

Co-located Multi environment and additive effect QTLs across trials
Multi-environment trial (MET) analysis using multi-model composite interval mapping (MCIM) detected a total of 23 additive effect QTLs for seven traits, viz., days to 50% flowering (1), plant height (4), number of total tillers (4), number of productive tillers (2), spikelet fertility (1), grain yield (4), straw yield (3), total biomass (4) ( Table 4) across five seasons. The regions identified using MCIM also showed to be linked to several traits in single trial or single trait in different trials, which can be referred to as pleiotropic effect of genomic region. For instance, the region RM5389-RM5794 on chromosome 1 was identified for several traits (Table 3). The markers RM11943-RM1361-RM3520-RM414-RM5794 between co-localized region RM5389-RM5794 on chromosome 1 has been detected for additive effect QTL using MCIM for plant height, number of productive tillers, straw yield, and total biomass ( Table 4) showing stability of the genomic regions for plant production traits with allelic contribution from landrace except at RM11943-RM1361 for plant height. MCIM analysis identified RM423-RM5345 on chromosome 2 for grain yield with significant large additive effect (210.55 kg/ha -1 ) from donor parent but with low phenotypic variance (R a 2 = 0.29%) ( Table 4). The large significant additive effect from donor parent was observed in nonstress trial (Data not shown) while in all other rainfed trial at TPE the allelic effect was from IR62266 showing environment-specific allelic expression. The loci with RM423 has showed both additive and additive x environment interaction effect (both A and AE) (Figure 1). Another region, RM324-RM5390 on chromosome 2 was identified to be linked with number of total tillers with additive effect from Norungan and significant expression only in DS2011 which also shows environment specificity (Table 4). The genomic region RM16030-RM4404 on Chromosome 3 is colocated for yield traits with additive effect from donor parent (Table 3). MCIM analysis also identified the distal region, RM6329-RM5475-RM16030 associated with spikelet fertility and straw yield ( Table 4). The region RM6329-RM5475 linked to spikelet fertility showed an additive effect from Norungan. While the region RM5475-RM16030 linked to straw yield showed an additive effect from IR62266.
Interestingly for straw yield all the drought stress trial (DS2011, DS2012, DS2013), the allelic effect is from Norungan while for non-stress (NS2012, DS2014) the allelic effect is from IR62266 (Data not shown). MCIM analysis also identified another region RM2334-RM7097 on  chromosome 3 to be linked with grain yield with additive allelic effect (28%) from IR62266. The higher phenotypic variance explained by additive x environment (Rae 2 = 1.56%) than that of phenotypic variance across environment (Ra 2 = 0.26%) shows environment interaction at RM2334-RM7097 region (Table 4). The region on chromosome 6 between RM585-RM217 was previously identified to be linked to grain yield under drought found to have consistent large effect (R 2 31.2 to 37.9%) and hence was fine mapped to 94 kb in same TPE using these RI lines in this laboratory. In this study, one or more markers located within this QTL region are detected to be linked to days to 50% flowering, spikelet fertility, grain yield, straw yield, and harvest index under drought (Table  3). MCIM analysis identified the region RM4173-RM6194-RM276 to be linked with days to 50% flowering and number of total tillers (Table 4, Supplementary  Figure 3). For flowering, the positive allele came from landrace Norungan while for number of total tillers positive allele came from IR62266.
MCIM analysis identified RM7361-RM8207-RM304 on chromosome 10 as an additive effect QTL for grain yield and total biomass with additive effect from recipient parent, IR62266. The region showed environment interaction for both the traits and explained higher additive x environment interaction variance (Rae 2 = 1.77% (GY), 0.99% (TB)) than that of additive effect across environments (Ra 2 = 0.35% (GY), 0.08% (TB)) ( Table 4). For both the traits, IR62266 contributed positive allele in nonstress condition only while Norungan contributed positive allele for all the stress trials in TPE (Data not shown).
The region RM209-RM6499 was identified for plant height, grain yield, and total biomass using MCIM analysis (Table 4). The additive allelic contribution for grain yield was from IR62266 while Norungan contributed for plant height and total biomass. Though the phenotypic variation was small, the phenotypic variation for grain yield and total biomass explained by additive x environment interaction (Rae 2 = 1.43 (GY), 1.11 (TB)) was larger than that explained by additive effect QTL across environments (Ra 2 = 0.05 for both traits) showing the interaction of complex traits with environment.

Epistatic QTLs, epistatic x environment interaction
A total of 16 pairs of significant epistatic QTLs (P< 0.01) were identified by the combined analysis of the multienvironment phenotypic values of stress (DS) and nonstress (NS) trials (Table 5). The phenotypic variance explained by each QTL across environment (Raa 2 ) was relatively smaller and ranged from 0.04 to 1.03%, than the phenotypic variation explained by each QTL x environment interaction (Raae 2 ) ranged from 0.08 to 2.00% indicating importance of QTL x environment interaction for drought tolerance.
2D genome scan revealed epistatic interaction of different plant production traits. The peak marker RM4472 positioned at 88.6 cM on chromosome 2 has epistatic (I) and epistatic x environment effect (IE) interaction with the RM5957 on chromosome 6 positioned at 90.8 cM for days to 50% flowering (Supplementary Figure 3).
Similarly, significant epistatic interaction was observed for plant height (Supplementary Figure 4). The peak marker RM11943 on chromosome 1 (140 cM) has both additive and additive x environment interaction effect (both A and AE) and has an epistasis x environment interaction effect (IE) with region on chromosome 11 (77.5 cM). Another peak marker RM10916 on chromosome 1 positioned at 55.7 cM has epistatic x environment interaction effect (IE) with region near to RM5746on chromosome 12 (18.2 cM), while only epistatic main effect (I) with other region between RM28048 and RM28076 positioned at50.6 cM on same chromosome 12 (Supplementary Figure 4).
We also observed positive interaction between major consistent QTLs (Figure 1). A total of five pairs of epistatic QTLs were detected for grain yield, among them only four pairs of epistatic QTLs showed significant aa interaction (P< 0.01). The additive epistatic effect was from donor parent, landrace Norungan at one epistatic QTL interaction locus while at other three locus effect was from recipient parent, IR62266 (Table 5). The donor parent epistatic locus, RM424-RM324 on chromosome 2 positioned at 32.3 cM has both epistatic (I) and epistatic x environment interaction effect (IE) with region RM286-RM6623 on chromosome 11 (8 cM). Similarly, two epistatic QTLs on chromosome 3, RM6832-RM3513 positioned at 91.5 cM and RM4404-RM3586 positioned at 131 cM also has both epistatic (I) and epistatic x environment interaction effect (IE) with region RM6194-RM276 on chromosome 6 (22.4 cM) and RM18758-RM18780 on chromosome 5 (78 cM), respectively. Another epistatic QTL region RM225-RM6917 on chromosome 6 (14.7 cM) has both epistatic (I) and epistatic x environment interaction effect (IE) with region RM6075-RM3840 on chromosome 8 (102 cM). The region on chromosome 1 positioned at 44.4 cM have epistatic x environment interaction effect (IE) with region on chromosome 2 positioned at 47.5 cM (Figure 1). The epistatic locus where positive allele came from donor Norungan showed lower aae (epistasis x environment effect) phenotypic variation compared to aa (epistasis effect) phenotypic variation and it was reverse for locus where positive allele came from IR62266. The locus where positive allele came from IR62266 has greater phenotypic variance explained by epistasis x environment than that of explained by epistasis across environment indicating epistatic effects were affected by environment (Table 5).

Discussion
The expression of genetic variability for droughttolerance traits requires a mapping population to be exposed to natural drought stress predominant in TPE (Muthukumar et al., 2015). Identifying major and consistent QTLs for drought tolerance using populations derived from locally adapted rice lines under drought predominant in TPE and the use of those QTLs in markerassisted breeding could be an effective strategy.

Variation for phenology and plant production traits
In this study, the landrace, Norungan, was used as a donor to identify consistent large-effect QTLs for plant production traits under drought stress in TPE. The drought tolerant donor performed better compared to IR62266 in terms of grain yield in all the trials indicating spill over effect of elite line suggesting potential yield as a strong determinant of yield even under stress. DS2013 considered as a severe drought year with the 75% mean yield reduction in DS2013 compared to DS2014 which received normal and evenly distributed rainfall with no visible symptoms of drought in the rice lines.
The landrace, Norungan flowered earlier than IR62266 in severe drought stress trial DS2013. Rice is very sensitive to drought especially during reproductive stage. Early flowering is a drought escape mechanism generally noticed in rice landraces (Ghimire et al., 2012) and drought tolerant donors.

Genetic map reconstruction
The low level of polymorphism (24.8%) despite a relatively large number of SSR markers (210) used, resulted in poor genome coverage and large gaps in the linkage map. This low polymorphism is due to the fact that both parents are indica ecotypes. Low level of polymorphism had already been noticed between the parents in intra-sub-specific (Ali et al., 2000) and even in inter-sub-specific crosses of rice (McCouch et al., 1988).
The discontinuity in linkage groups was reduced with saturation of linkage map and resulted in continuous linkage group, thus helped to reduce unlinked markers and improve the resolution of the linkage map and validating marker order.

Consistent QTLs for plant production traits
For MAB to be effective, it is desirable to identify consistent large effect QTLs for plant production traits across seasons and environments. The region RM5389-RM5794 on chromosome 1 consistently identified for plant production traits. The region including marker RM5389, RM11943, and RM5794, on ch.1 is a QTL hotspot for tolerance against different abiotic stresses in rice (Kamoshita et al., 2008; Vikram et al., 2011). This region is also consistently associated with yield under drought (qDTY1.1) in different populations of rice (Ghimire et al., 2012;Prince et al., 2015;Rajurkar et al., 2019;Venuprasad et al., 2012;Vikram et al., 2011). Meta-QTLs (MQTL1.3, MQTL1.4) for yield under drought have been mapped in this region  and drought responsive genes underlying this genomic region have been identified using rice lines derived from drought tolerant landrace, Dhagaddeshi and drought susceptible IR20 (Borah et al., 2017).
The region RM423-RM5345 on chromosome 2 identified using MCIM analysis with large significant additive effect from donor parent for grain yield. The marker RM12460 just above to this region was also identified for HI and DFF in DS2014. This region RM12460-RM423-RM5345 on chromosome 2 was reported for grain yield under drought (qDTY2.2) from drought tolerant Aday sel previously (Swamy et al., 2013). High and consistent additive effect of qDTY2.2 under severe lowland drought stress and aerobic conditions on GY was reported by Dixit et al. (2012). This region (qDTY2.2) was consistently linked to grain yield under drought (Shamsudin et al., 2016). The qDTY2.2 region has been reported previously to have an effect on other drought-related traits apart from GY, including drought tolerance index, canopy temperature, osmotic adjustment, and leaf water content (Robin et al., 2003;Yue et al., 2005). Similarly, QTLs for root traits viz. root length (MacMillan et al., 2006) and root thickness (Champoux et al., 1995) have been reported within and near to qDTY2.2 region, respectively. Colocation of QTLs for yield and secondary traits suggests importance of this region in rice yield improvement under drought (Swamy et al., 2013).
Another region RM324-RM5390 on chromosome 2 identified to be linked with number of total tillers with additive effect from Norungan with significant expression only in moderate stress showing environment specificity. The marker RM324 was also identified for spikelet fertility with separate year QTL analysis in DS2014. The distal marker RM7426 to region RM324-RM5390 was identified for grain yield, straw yield, and biomass in DS2013. Using separate year analysis, for all the traits, at marker RM324 and RM7426 allelic contribution came from IR62266. The region RM324-RM7426 on chromosome 2 was earlier reported as one of the major QTLs (qDTY2.1) for drought tolerance under lowland and upland environments using Apo and Swarna recombinant lines (Venuprasad et al., 2009). Venuprasad et al. (2009) also reported allelic effect of both the parents in stress and non-stress condition indicating requirement of opposite alleles to optimize performance under stress and non-stress condition. Further, Dixit et al. (2012) fine mapped qDTY2.1 region to 1.6 cM and showed segments of donor and recipient parents are present in the qDTY2.1 region with differential expression under varying level of stress condition. The study by Dixit et al. (2012) also showed varying additive effect under moderate and severe stress at qDTY2.1 with peak marker RM324, showing environment specificity of this region with peak marker RM324. The different allelic effect also indicates different physiological mechanism, such as high osmotic adjustment from IR62266 with response to severity of stress, conferring effect on yield and biomass.
The region RM6329-RM5475-RM16030 on chromosome 3 was identified for spikelet fertility and straw yield while another region RM2334-RM7097 was identified for grain yield using MCIM analysis. The regions, RM6329-RM5475 for spikelet fertility and RM5475-RM16030 for straw yield and RM2334-RM7097 for grain yield showed shifting allelic effect from either of parent under drought stress and nonstress. The shifting allelic effect from parent suggest the plasticity of genomic region to environment and their adaptation to respective environment. CIM analysis also reported, distal region RM16030-RM4404 on Chromosome 3 colocated for yield traits with additive effect from donor parent. The consistent region with RM16030 for GY was earlier identified as qDTY3.1 with large effect (R 2 = 31%) on yield in both lowland drought and aerobic environments and reported shifting allelic effect from drought tolerant Apo and high yielding Swarna under stress and non-stress condition, respectively (Venuprasad et al., 2009). Swamy et al. (2011) reported meta QTL (MQTL3) with large effect (R 2 = 20%) on grain yield in this region. Shamsudin et al. (2016) pyramided qDTY2.1 and qDTY3.1 into elite Malaysian cultivars MR219 and MRQ74 using MAB and reported yield advantage under drought in MAB derived lines.
One or more markers located within previously fine mapped region RM585-RM217 on chromosome 6 were reported to be linked to several traits in this study. MCIM analysis also identified the region RM4173-RM6194-RM276 to be linked with days to 50% flowering and number of total tillers. The markers, RM6194 and RM4173 were also found to be colocated for yield and plant production traits using CIM. Interestingly, allelic effect for days to 50% flowering came from donor parent Norungan while all other plant production traits and grain yield the allelic effect is from IR62266. The allelic effect from IR62266 at RM4173 for number of total tillers complements each other in both, CIM and MCIM analysis. The positive allelic effect for flowering shows drought escape mechanism which is also evidenced from early flowering in Norungan than IR62266 in both rainfed trials. While the allelic effect from IR62266 for grain yield, harvest index, plant height, number of total tillers and number of productive tillers shows drought adaptation with higher osmotic adjustment in IR62266. Selectively combining early flowering, drought adaptation traits and high yield potential through molecular breeding may help increase rice production under water limited environment.
The region on chromosome 10 between RM7361-RM8207-RM304 identified for grain yield and total biomass as additive effect QTL using MCIM analysis. This region was reported for grain yield under drought (qDTY10.1) in earlier study (Swamy et al., 2013;Vikram et al., 2011) and also reported to be a meta QTL for yield under drought (MQTL10.1) . The additive effect for grain yield and biomass is from recipient parent, IR62266. Recipient parent, IR62266 showed an allelic effect in non-stress trial while all rainfed trials in TPE the additive effect is from Norungan, which suggest environment specificity and allelic plasticity of this region. Allelic plasticity in terms of qDTY10.1 is reported earlier, the recipient parent (MTU1010) contributing allele in stress and non-stress condition  while allelic effect from donor parent (Aday sel) in stress condition (Swamy et al., 2013). The marker RM7361 was linked to plant height and number of productive tillers in DS2014 using CIM analysis and allelic effect was from Norungan, confirming the environment specificity and allele plasticity of the region for plant production traits.
MCIM analysis identified the region RM209-RM6499 on chromosome 11 colocated for plant height, grain yield, and total biomass. The allelic effect for grain yield was from IR62266 while for plant height and biomass allelic contribution was from Norungan. CIM analysis also identified RM209 to be colocated with 50% flowering, grain yield, straw yield, and biomass. The allelic effect at RM209 for grain yield in DS2011 was from donor parent using CIM analysis, MCIM analysis complements the allelic effect from donor parent, Norungan for grain yield and total biomass in all rainfed trials in TPE, suggesting environment-specific allelic contribution of the region. The region with marker RM209 on chromosome 11 was earlier reported as a QTL region affecting grain yield under drought stress by affecting net photosynthetic rate, stomatal conductance, transpiration rate (Zhao et al., 2008). The consistency of RM209 for plant production traits with positive allele contribution from drought tolerant Norungan in rainfed environment and its association with photosynthetic traits shows the significance of this region in drought adaptation in rainfed rice.
The region on chromosome 12 with marker RM28089 colocated for grain yield and harvest index in DS2011 with positive allele contribution from donor, drought tolerant parent, Norungan. The marker RM28089 was also reported for a major drought grain yield QTL on chromosome 12, (qDTY12.1), with positive allele for qDTY12.1 contribution by the tolerant parent IR74371-46-1-1. (Bernier et al., 2007).

Additive and epistatic QTLs, and environment interaction
The higher phenotypic variation explained by QTL x environment interaction indicated QEI as an important component of drought tolerance. The phenotypic variance explained was 0.04-1.03%, while variance explained by QTL x environment interaction was 0.08-2.00%. Similar low phenotypic variance was observed for yield-related traits under low temperature in rice (Yang et al., 2017). In the present study, the phenotypic variance of epistatic QTLs (R aa 2 ) ranged from 0.04% to 1.03% while phenotypic variance of additive QTLs (R a 2 ) ranged from 0.05% to 2.50%. The relative lower phenotypic variance of epistatic QTLs than that of the additive QTLs indicate the minor role of epistasis in plant traits, similar minor role of epistasis was reported in rice yield traits under low temperature (Yang et al., 2017).
Epistasis for grain yield was detected in the regions of non-additive effects and its relative contribution was small. Relative small percentage of grain yield variation was explained by additive x environment interaction (Rae 2 = 0.28-1.77%), but it was larger than that explained by additive effect QTL across environments (Ra 2 = 0.05-0.35%) showing environment interaction with complex trait like grain yield. Epistatic interaction for traits related to drought tolerance and major grain yield QTLs under drought was observed previously (Sandhu et al., 2014). The lower R ae 2 (0.28%) or R aae 2 (0.27%) compared to R a 2 (0.29%) or R aa 2 (0.34%) at region RM423-RM5345 and RM424-RM324, respectively, which linked with grain yield on chromosome 2, and where donor drought tolerant parent Norungan contributed positive allele, suggest donor parent adaptation to environment with lower phenotypic variance explained by additive/epistasis x environment interaction. The interaction between different traits and non-additive QTLs could be useful in improvement of more than one trait simultaneously and will help genomic dissection of complex traits.

Conclusion
Six genomic regions linked to plant production traits under drought in TPE and non-stress conditions have been detected with combined analysis using rice lines derived from locally adapted landrace in this study. The region RM5389-RM5794 on chromosome 1 linked with straw yield and total biomass seemed to be a QTL hotspot colocated with plant height, and number of productive tillers. The region RM423-RM5345 on chromosome 2, linked to grain yield with positive allele from landrace Norungan. Similarly, RM6329-RM5475-RM16030 on chromosome 3 linked with spikelet fertility and straw yield harbors MQTLs (MQTL3) for grain yield. The region on chromosome 6, PSM202-RM508 linked with total biomass. QTLs for several secondary traits under drought have been reported earlier at this genomic region in rice. The region RM7361-RM8207 on chromosome 10 associated with yield and colocated with plant height, number of productive tillers. Further, this QTL region was identified as yield QTL (qDTY10.1) in rice. The region RM209-RM6499 on chromosome 11 was also associated with yield and affect yield through their linkages with DFF, number of productive tillers, and biomass. These six regions identified using combined analysis are earlier identified as yield QTLs under drought and their detection using this unique mapping population derived from locally adapted landrace further validates these QTLs. These QTLs may be robust candidates for yield improvement in rice under water-limited environments through map-based cloning of underlying genes and MAB strategies.