The K/Z RIL population was developed at the USDA Dale Bumpers National Rice Research Center, Stuttgart, Arkansas, USA by crossing diverse parental genotypes from different subspecies, Kaybonnet (tropical japonica) and ZHE733 (indica) by SSD method to create segregating progenies with high genetic variability for selection desirable genes. Root growth and development vary among the rice genotypes studied here. Greub (2015) identified that the root length of Kaybonnet showed the longest root compared to diverse rice genotypes, such as Bengal, Sipirasikkam (GSOR 310428), O. glaberrima, IR64, Nagina-22 (N22), and Vandana. Root length is probably the most important architectural trait for drought avoidance, by enabling the roots to reach deeper water levels in the ground. Additionally, Greub (2015) also determined that Kaybonnet had higher metaxylem number than Bengal, O. glaberrima, IR64, Nagina-22 (N22), Vandana, and parental line ZHE733. The number and size of the metaxylem vessels are associated with their conductivity of water. Among these lines, Bengal and Nagina-22 (N22) are examples of drought-resistant rice genotypes, while IR64 and Nippobare as the reference for drought-sensitive genotypes.
The differences of the root anatomy and root morphology between Kaybonnet (drought resistant) and ZHE733 (drought sensitive) displayed in figure 1 and supplementary table 2. Kaybonnet exhibited longer root than ZHE733, and also Kaybonnet has two metaxylem number while ZHE733 only has one metaxylem in the stele area. Additionally, total aerenchyma area in Kaybonnet (0.34 mm2) is wider than ZHE733 (0.29 mm2). According to Zhu et al. (2010) aerenchyma improves drought resistance by decreasing root metabolic costs and greater water acquisition from dehydrated soil. Thus, Kaybonnet exhibited more drought resistance compared to ZHE733.
Based on the screening of grain yield under DS conditions at reproductive stage (R3), among the two parental genotypes, Kaybonnet is drought resistant while ZHE733 displays a drought sensitive phenotype. Under DS conditions, the number of filled grains per panicle in Kaybonnet reduced 20%, while ZHE733 reduced 50% (Fig. 2). The progeny of a cross between drought resistant and sensitive genotypes are useful to study the inheritance of drought resistance, and identification of important QTLs for variation in grain yield under DS conditions (Kumar et al., 2018; Islam et al., 2012). A total of 198 K/Z RILs was studied for filled grains per panicle number, and out of which, 13.13% were found to be highly drought resistant lines, 11.11% moderately drought resistant lines, and 75.75% drought sensitive lines, suggesting there are multiple factors involved in inheritance of drought resistant and sensitive phenotypes in the population.
Variation in morphological traits of RILs under drought stress conditions
Rice is more sensitive to DS conditions compared to the other cereal crops such as wheat, rye, and barley (Huang et al., 2014). The two parental lines are diverse in the morphological traits such as plant height (PH), tiller number (TN), heading day (HD), flag leaf width (FLW), and leaf rolling score (DLR). Kaybonnet (derived from variety Katty and Newbonnet) compared to ZHE733 is the donor parent, with high stature, low TN, late HD, wider FLW, and low DLR while ZHE733, the recurrent parent, has short stature, higher TN, early HD, narrow FLW, and high DLR under WW conditions. In addition to the two parental lines having variation for morphological traits when subjected to DS conditions, Kaybonnet showed overall superior performance to ZHE733 because Kaybonnet is drought resistance while ZHE733 is sensitive to DS conditions, thus exhibiting less resistance. Under DS conditions, PH of Kaybonnet reduced to 14.29%, while ZHE733 reduced to 39.22%. Meanwhile, TN of Kaybonnet reduced to 20%, and ZHE733 reduced to 14.3%. These results are in agreement with previous studies showing that water deficit conditions have a negative effect on plant development and growth due to a loss of turgor (Hsiao et al., 1970; Specht et al., 2001; Farooq et al., 2011; Todaka et al., 2015; Kumar et al., 2019). The reduction of the plant height and productive tiller number under DS conditions are associated with the reduction of the cell cycle processes, cell expansion and elongation (Mantovani and Iglesias, 2008).
In the RIL population, morphological traits such as PH, TN, HD, FLW, and DLR showed normal frequency distribution (Figs. 3 and 4) signifying quantitative inheritance, and thus the morphological traits were suitable for QTL analysis (Fang et al., 2019). Within the RIL population, there was a wide range of morphological traits that showed variation across WW and DS conditions for PH and TN, while HD, FLW, and DLR measured in DS conditions exhibited transgressive segregation pattern.
Within the RIL population, all the morphological traits under DS showed a significant reduction compared to WW conditions (Supplementary Table 1). Based on the PH and TN characteristics in WW conditions, the RIL population is more skewed toward parent ZHE733, plant height greater than 50 cm and tiller number greater than 8. The average PH in the RIL population under DS conditions showed 50% reduction, which is more than either of the parents, while the average of TN reduced to 13.90% which is less than either of the parents. Based on the PH reduction in the DS conditions compared to WW conditions, 54.04% of lines in the RIL population showed more than 50% reduction suggesting more than 50% of the lines are sensitive to DS conditions. Meanwhile, based on the TN reduction, 71.72% of lines in the RIL population showed drought resistant with less than 30% reduction suggesting that majority of the lines have reduced the tiller number and have characters of ZHE733. Among the RIL population, 28 lines showed less reduction (≤ 30%) for both plant height and tiller number traits and these lines could be potential genetic resource to develop drought resistant varieties.
Genetic variation for grain yield components under reproductive stage drought stress
The effect of DS conditions on rice are different at every growth stage. Rice is very sensitive to DS conditions, especially at the reproductive stage, and grain yield is dramatically reduced even under slight DS conditions (Kamoshita et al., 2008; Palanog et al., 2014; Dwiningsih, 2020b). The effective way to reduce grain yield loss under DS conditions is to develop drought-resistant rice varieties, that is however very challenging due to the complexity of the drought resistant trait associated with grain yield components. The identification of drought QTLs associated with the grain yield components under water defisit conditions, and their utilization in molecular breeding, is an alternative method of increasing breeding effectiveness. Moreover, identified QTLs can be incorporated in a marker assisted breeding (MAB) strategy (Venuprasad et al., 2011).
The mean values of the grain yield components showed a significant difference in performance of both parents and population under DS conditions compared to WW, including biological yield (BY), spikelet per panicle number (SP), filled grain per panicle number (FG), panicle length (PL), and primary panicle branch number (PPB) (Table 1). ZHE733 showed a greater reduction in overall grain yield components compared to Kaybonnet. Under DS conditions, ZHE733 showed significant reduction (P < 0.05) in BY, SP, FG, PL, and PPB; 51.49%, 38.89%, 62.96%, 16.13%, and 12.5%, respectively. In addition, Kaybonnet showed lower reduction compared to ZHE733 in BY, SP, FG, and PL; 4.73%, 28.85%, 24.24% and 11.06%, respectively. The trait PPB in Kaybonnet showed greater reduction 31.58% compared to ZHE733.
In the RIL population, grain yield components showed a continous frequency distribution under both WW and DS conditions, indicating polygenic control of the traits as shown in Fig. 5A (BY), Fig. 5B (SP), Fig. 5C (FG), Fig. 5D (PL), and Fig. 5E (PPB). Transgressive segregation was observed for the grain yield components under WW and DS conditions. Moreover, there were significant differences (P < 0.05) among the RILs for all the grain yield components under both WW and DS conditions. Under DS conditions, the average of BY, SP, FG, PL, and PPB in the RIL population reduced 12.5%, 38.75%, 47.62%, 15.82%, and 16.92, respectively compared to WW conditions. Based on the BY, SP, PL, and PPB reduction in the DS conditions compared to WW conditions, most of the lines in the RIL population is skewed toward drought resistant phenotypes with reduction less than 30%. Meanwhile, based on the FG, most of the lines belong to drought sensitive with reduction more than 50%. There are 12 lines that showed less reduction (≤ 30%) for all the grain yield componets and these RIL lines could be a potential genetic resource to develop drought resistant rice varieties.
Table 1 The average and range values of grain yield components of K/Z RIL population under WW and DS conditions
|
Traits
|
Treatments
|
Kaybonnet
|
ZHE733
|
K/Z RIL Population
|
|
Average
|
Range
|
|
Biological yield
|
WW
|
14.80
|
13.40
|
20
|
4 - 49
|
|
DS
|
14.10
|
6.50
|
17.50
|
2 - 47
|
|
Spikelet per panicle number
|
WW
|
104
|
90
|
133.35
|
45 - 328
|
|
DS
|
74
|
55
|
81.68
|
26.8 - 188.4
|
|
Filled grain per panicle number
|
WW
|
66
|
43.2
|
43.07
|
2 - 176
|
|
DS
|
50
|
16
|
22.56
|
0 - 97
|
|
Panicle length (cm)
|
WW
|
21.7
|
18.6
|
21.24
|
11.90 - 33.30
|
|
DS
|
19.3
|
15.6
|
17.88
|
11.64 - 26.62
|
|
Primary panicle branch number
|
WW
|
19
|
8
|
11.17
|
3 - 22
|
|
DS
|
13
|
4
|
9.28
|
5.4 - 15.6
|
Variation in root architectural traits under ABA conditions
The parents, Kaybonnet and ZHE733 had contrasting responses under ABA conditions, where Kaybonnet as drought-resistant parent showed more sensitivity to ABA compared to ZHE733 as the drought sensitive parent. Both parents experienced a reduction in root length (RL), total root number (TRN), shallow root number (SRN), deep root number (DRN), and root fresh weight (RFW) under ABA conditions compared to control conditions. Kaybonnet showed a greater reduction in RL, TRN, SRN, DRN, and RFW compared to ZHE733 (Supplementary Table 2). Under ABA conditions (3 µM), Kaybonnet showed high reduction in RL, TRN, SRN, DRN, and RFW; 61.94%, 62.12%, 100%, 66.22%, and 50%, respectively. However, ZHE733 showed less reduction under ABA conditions (3 µM) in RL, TRN, SRN, DRN, and RFW; 17.59%, 5.43%, 4.68%, 5.79%, and 6.45%, respectively. Furthermore, under ABA conditions (5 µM), ZHE733 also displayed less reduction in RL, TRN, SRN, DRN, and RFW compared to Kaybonnet.
The K/Z RIL population showed variation in their root architectural traits in response to ABA (0 µM, 3 µM, and 5 µM), such as reduction in RL (Fig. 6A), increase in RSR (Fig. 6B), reduction in TRN (Fig. 6C), decrease of SRN (Fig. 6D), shortening in DRN (Fig. 6E), and lessening RFW (Fig. 6F). Furthermore, the distribution of root architectural traits studied under ABA conditions showed a near normal distribution. The range of RL, RSR, TRN, SRN, DRN, and RFW of the K/Z RIL population under ABA conditions also exhibited a wide variation (Supplementary Table 2). Additionally, the average of root architectural traits of the RILs lie between the parents in ABA conditions. The average of RL, TRN, SRN, DRN, and RFW in the RIL population under ABA conditions 3 µM also showed reduction; 18.69%, 16.46%, 33.78%, 9.26%, and 2.94%, respectively. In addition, under ABA conditions 5 µM; the average of RL, TRN, SRN, DRN, and RFW in the RIL population also showed reduction; 21.16%, 21.07%, 36.51%, 12.61%, and 23.53%. Although both parents and the RIL population exhibited a reduction in RL, TRN, SRN, DRN, and RFW under ABA conditions, the reduction was greater in Kaybonnet. Root to shoot ratio (RSR) increased 38% and 16% for Kaybonnet under ABA conditions 3 µM and 5 µM, while for RIL population under ABA conditions 3 µM and 5 µM increased 22.92% and 22.69%. Previous studies have demonstrated that ABA sensitivity is associated with drought stress resistance through its effect on the stomatal movement (Lim et al., 2015; Duan et al., 2008; Todorov et al., 1998). Lim et al. (2015) also indicated that drought-sensitive rice plants grown in media with 2 or 5 µM ABA had significantly longer roots and shoots compared to the plants in control media. These data suggested that the drought-sensitive rice plants were insensitive to ABA and exhibited the ABA-dependent pathway in response to drought stress.
Correlation of morphological traits and grain yield components under WW and DS conditions with root architectural traits under ABA conditions
Correlation analysis increases an understanding of the overall contribution of various rice plant traits to each other (Gibert et al., 2016). A Pearson’s correlation coefficient analysis was carried out on morphological traits and grain yield components under WW and DS conditions, and also on root architectural traits under ABA conditions, to analyze the correlations among them. Significant correlations were observed among all of the traits studied. Furthermore, FG-DS as the major trait among the grain yield components under DS conditions showed significant positive correlations with most of the morphological traits, including HD, TN-WW, PH-WW, PH-DS, FLW, and DLR. FG-DS also has positive correlations with other grain yield components, such as PL-WW, PL-DS, PPB-WW, SP-WW, SP-DS, and FG-WW (Supplememtary Table 5). Additionally, FG-DS exhibited significant positive correlations with most of the root architectural traits under ABA conditions such as RL-ABA0, RL-ABA3, RL-ABA5, RSR-ABA0, RSR-ABA3, RSR-ABA5, TRN-ABA0, TRN-ABA3, TRN-ABA5, DRN-ABA0, SRN-ABA3, SRN-ABA5, DRN-ABA3, DRN-ABA5, RFW-ABA0, RFW-ABA3, and RFW-ABA5 (Supplementary Table 6). However, significant negative correlations were also observed between FG-DS and the morphological trait like TN-DS, and also with the grain yield components such as BY-WW, BY-DS, and PPB-DS. Several root architectural traits also showed significant negative correlations with FG-DS, including SRN-ABA0.
In this study, a positive correlation found between FG-DS with most of the morphological traits, the other grain yield components, and the major root architectural traits under ABA conditions, indicate that the rice drought-resistant plants maintain their grain yield under DS conditions through development of cell elongation, maintenance of cellular membrane integrity, and regulation of osmotic stress tolerance via ABA-mediated cell signaling (Kanbar et al., 2009; Ramegowda et al., 2015; Ding et al., 2016; Basu et al., 2016; Catalos et al., 2017; Nada et al., 2018; Hassaoni et al., 2018; Li et al., 2019; Kim et al., 2020). Furthermore, the negative correlation between FG-DS with BY under WW and DS conditions indicate that there is more assimilate distribution into the grains compared to the other parts of the plant. Moreover, ABA sensitivity was found correlated with the drought resistance of rice plants (Lim et al., 2015; Duan et al., 2008; Todorov et al., 1998).
Based on the correlation analysis, three lines (100007, 100036, and 100135) of the RIL population showed drought resistant, based on the morphological traits (PH and TN) and grain yield components (BY, SP, FG, PL, and PPB) and ABA sensitivity response characteristics similar to Kaybonnet as drought resistant parent. These three lines showed less than 30% reduction under DS conditions for PH, TN, BY, SP, FG, PL, and PPB. Meanwhile, for root architectural traits (RL, TRN, SRN, DRN, and RFW), these three lines showed more than 50% reduction. Furthermore, based on the morphological traits and grain yield components under WW conditions, these three lines displayed short PH and high number in TN, BY, SP, PL, FG, and PPB which become a good genetic resource to develop drought resistant varieties. Under WW conditions, the range of PH, TN, BY, SP, PL, FG, and PPB were 42 – 51 cm, 5 – 9, 17 – 24 g, 104 – 130, 17.5 – 23.7 cm, 100 – 103, and 10 – 12, respectively.
Relationships of root anatomy and drought resistance
Many studies have been reported that various root anatomical phenes such as metaxylem, aerenchyma, and cortical cells properties influence plant performance and productivity under DS conditions in cereal crops including maize (Chimungu et al., 2014a, b; Chimungu et al., 2015a, b), rice and wheat (Kadam et al., 2015); and also in legume crops such as soybean (Prince et al., 2017), common bean (Peña-Valdivia et al., 2010), chickpea, groundnut, pigeonpea, and cowpea (Purushothaman et al., 2013). Under DS conditions, plants allocate more energy and carbon resources to root growth rather than to shoot growth, which can increase water acquisition (Lynch and Ho, 2005; Palta and Gregory, 1997; Sharp and Davies, 1979). Aerenchyma and cortical cells contribute to increase root growth for soil exploration in the deeper soil layer, consequently, greater water acquisition, and improved productivity by reducing the root metabolic costs leading more carbon resources allocation to improve root growth (Jaramillo et al., 2013; Lynch, 2013). Understanding root adaptive mechanisms is important to maintain higher productivity under DS conditions.
Metaxylem vessels play an important role in water and nutrients uptake from the soil, and subsequent transport within the plant cells (Steudle and Peterson, 1998). A study by Kadam et al. (2015), metaxylem vessel properties are important traits to increase water-use efficiency under DS conditions. Parental Kaybonnet and drought-resistant lines showed higher number of metaxylem vessels compared to ZHE733 and drought-sensitive lines (Fig. 7). Kaybonnet has two metaxylem vessels (Fig. 1) and in an average, drought-resistant lines have three metaxylem vessels while ZHE733 and drought-sensitive lines only have one metaxylem vessels. Increasing of the metaxylem number leading to increase root hydraulic conductivity that decrease the metabolic costs for exploring water in deeper soil layer allowing water uptake efficiency and higher yield under DS conditions (Prince et al., 2017). The increased root hydraulic conductivity also leads to increase in performance of shoot physiological processes under DS conditions. The increased number of metaxylem vessels, improves water uptake efficiency under DS conditions thereby causing increased stomatal conductance and internal capture of carbon dioxide and thus overall maintains photosynthesis activity under DS conditions. In contrast, drought-resistant lines displayed smaller metaxylem size compared to drought-sensitive lines. The average of the metaxylem size of the drought-sensitive lines are 0.0008 mm2 while drought-resistant lines are only 0.0005 mm2. According to Richards and Passioura (1989), the decrease of metaxylem size allows greater yield under DS conditions, because decreased metaxylem size is correlated to greater hydraulic conductivity, allowing greater water acquisition.
Aerenchyma is the enlarged air space in the root cortex resulting from programmed cell death that affects root respiration (Evans, 2004). Aerenchyma showed a positive correlation with drought resistance in rice. Kaybonnet and drought-resistant lines also showed higher percentage of aerenchyma than ZHE733 and drought-sensitive lines. The average coverage of aerenchyma in drought-resistant lines is 60% of the total root cross section while drought-sensitive lines is only 40%. Higher aerenchyma is associated with decreased root respiration costs, leading to more carbon allocation and improved root growth, thereby increasing water acquisition from the deep soil layer of drying soil and higher grain yield under DS conditions (Zhu et al., 2010).
Another root anatomical phene, important for improved drought resistance is cortical cell. Cortical cell area is the total transverse root cortex minus aerenchyma area. There are average of eleven cortical cell file number in drought-sensitive lines that is higher compared to five in the drought-resistant lines. Cortical cell file numbers are positively linked to higher root respiration (Jaramillo et al., 2013; Lynch, 2013) therefore the drought-sensitive lines have higher respiration compared to drought-resistant lines, causing reduced root growth, lower water acquisition and decreased productivity under DS conditions. According to Rieger and Litvin (1999), reduced cortical cell file number also showed a positive correlation to a smaller radial path that affect root hydraulic conductivity, thereby increasing water acquisition under DS conditions.
High-density genetic linkage map with GBS markers
A genetic linkage map is an important tool to explore the plant genome, and to obtain information of allele introgression during plant breeding efforts (De Soursa et al., 2015). By using high-density genetic linkage map leads to narrow down the location of QTLs into a specific region and predict more accurate candidate genes for gene cloning, then validation with reverse genetics approaches (Hattori et al., 2009). Based on the GBS analysis, 28,598 SNP markers were obtained from 200 samples (2 parental lines and 198 RILs) with heterozygosity level of 1.3% and non-parental alleles at 0.4%. The filtering process of the SNP markers was done based on the missing data (≤90%), minor allele frequency (MAF < 1%), polymorphic markers between parents, recombinant frequency, and percentage of heterozygosity. Furthermore, 4133 filtered SNP markers were obtained, and were used in the development of the high-density genetic linkage map by using QTL IciMapping software version 4.2.53 with the Kosambi mapping function.
The number of SNP markers mapped to each chromosome varied from 182 SNP markers found on chromosome 12 to 562 SNP markers on chromosome 1, with an average of 344.42 SNP markers per chromosome. The total length of the genetic linkage map was 6063.12 cM (varied from 343.72 cM on chromosome 10 to 676.52 cM on chromosome 2), with the average of 505.26 cM per chromosome. A calculated average genetic distance between two SNP markers across the chromosome was 1.58 cM (ranged from 0.92 cM on chromosome 6 to 2.89 cM on chromosome 12). The density of the genetic linkage map was 0.69 SNP markers per cM or an average of 1 SNP marker every 1.5 cM. This high-density genetic linkage map covered 373 Mb of the rice genome and can be used to identify QTLs with higher resolution and reliability for application in rice breeding under DS conditions. Moreover, the high-density genetic linkage map lead to identification and selection of candidate genes within the QTL regions that are involved in improving drought resistance in the rice plants, and towards developing drought-resistant rice varieties. Many previous studies have also used high-density linkage maps for QTL mapping (Bhattarai and Subudhi, 2018; Sabar et al., 2019; Barik et al., 2019; Melandri et al., 2020; Barik et al., 2020; Dwiningsih et al., 2021b).
QTL mapping of morphological and yield traits under reproductive stage drought stress conditions and root architectural traits under ABA conditions
The identification of QTLs for morphological traits, grain yield components, and root architectural traits is important to understand the genetic complexity of the drought-related traits. The genetic variation of drought-related traits is regulated by numerous genes that have a large effect on the traits (Baisakh et al., 2020). In this research, 41 QTLs were identified for morphological traits and grain yield components under WW and DS treatments, and also root architectural traits under ABA treatments (Supplementary Table 3, Fig. 8). The identified QTLs varied under WW, DS, and ABA treatments. Additionally, the number of QTLs varied for morphological traits, grain yield components, and root architectural traits. A total of 12 QTLs were identified for multiple morphological traits, grain yield components, and root architectural traits suggesting that these QTLs have pleiotropic effect. For this study, the LOD values more than 2.5 was selected in order to filter unassosiated candidate regions. The maximum values of LOD was 28.8 and PVE was 13.809% (Supplementary Table 3).
Most of the QTL identified in this study were mapped to approximately the same locations as previous reports (Mu et al., 2003; Courtois et al., 2003; Lanceras et al., 2004; Bernier et al., 2007; Bernier et al., 2009; Venuprasad et al., 2009; Mishra et al., 2013; Yadaw et al., 2013; Wang et al., 2014; Palanog et al., 2014; Saikumar et al., 2014; Prince et al., 2015). Chromosome 1 harbored the highest number of QTLs for PH (Monna et al., 2000; Lanceras et al., 2004; Zhou et al., 2016; Jiang-xu et al., 2016; Yadav et al., 2019a; Zeng et al., 2019; Xu et al., 2020). The most QTLs associated with grain yield were located on chromosome 5 and 6 (Solis et al., 2018; Yadav et al., 2019b; Baisakh et al., 2020; Takahashi et al., 2001; Bernier et al., 2007; Xu et al., 2020). Moreover, chromosome 10 had the highest number of QTLs for root architectural traits (Xu et al., 2001; Lou et al., 2015; Kitomi et al., 2015; Gimhani et al., 2018). Eventhough there were overlapping genomic regions with the previous studies, the results of this study predicted QTLs/genomic regions in close association with genes of interest.
Candidate genes underlying QTL regions
Identification of candidate genes within a QTL region is useful for marker-assisted pyramiding to develop drought-resistant rice varieties (Bhattarai et al., 2018), and is also important for developing transgenic rice with enhanced drought resistance (Varshney et al., 2011). In this research, we identified candidate genes involved in many biological processes, molecular functions, cell components, and drought response (Supplementary Table 3). The QTL clusters contain candidate genes with large pleiotropic effects. The total number of 184 candidate genes with an average of 5 genes per QTL were identified in 41 QTLs for morphological traits, grain yield components under WW and DS conditions, and root architectural traits under ABA conditions (Supplementary Table 3). These candidate genes were distributed unevenly on different chromosome. There were 15 candidate genes with high LOD and PVE score that have been previously studied for drought stress. In total 11 novel candidate genes with high LOD and PVE score were discovered in this analysis with unknown annotations. Many known genes present within these QTL regions include genes with homology to APETALA2 (AP2)/ETHYLENE RESPONSE FACTOR (ERF) transcription factor, malate dehydrogenase protein, photosystem II oxygen evolving complex protein PsbQ family protein, WRKY transcription factor, MYB transcription factor, Zinc Finger (ZFN) protein, endoplasmic reticulum protein, DEAD-box RNA helicase, glycosyl transferase protein, Late Embryogenesis Abundant (LEA) protein, and no apical meristem protein (NAC).
A transcription factor family well-known for drought response like APETALA2 (AP2)/ETHYLENE RESPONSE FACTOR (ERF) (Lata and Prasad, 2011; Mizoi et al., 2012; Licausi et al., 2013; Phukan et al., 2017), has family members present in several of the QTL regions such as QTLs for PH-WW on chromosome 1 (LOC_Os01g66270); RL-ABA3, RL-ABA5 on chromosome 2 (LOC_Os02g54160); and FG-WW, FG-DS, RL-ABA3, and RL-ABA5 on chromosome 5 (LOC_Os05g49010). Furthermore, the QTL regions for PH-WW on chromosome 1 is adjacent to the semi-dwarfing gene sd1 locus (38.3 Mb). A strong linkage has been found previously between sd1 and QTL for drought-related traits (Vikram et al., 2015). The sd1 locus is also associated with underground and above ground traits in rice, such as plant height and root architectural traits (Yadav et al., 1997). Reduction in plant height under DS conditions, is the adaptation of rice plants to DS. An important QTL for grain yield components under DS, FG located on chromosome 5, overlaps with 12 candidate genes. There is also an overlap between QTL for root architectural traits, RL under ABA conditions and FG under DS, suggesting that ABA is involved in the drought stress resistance mechanism. Under DS conditions, ethylene biosynthesis is increased and interacts with AP2/ERF, and finally a response to water deficit (Abiri et al., 2017; Nakano et al., 2006; Ma et al., 2014). Additionally, AP2/ERF responds to ABA in order to help activate ABA dependent and independent stress responsive genes. Transgenic rice with over-expression of an AP2/ERF showed an increase in drought resistance (Pan et al., 2012). An understanding of the AP2/ERF gene functions in the drought resistance mechanisms in rice, may provide valuable information to facilitate the improved adaptation of rice to DS conditions.
An important candidate gene, LOC_Os03g56280, known to regulate malate dehydrogenase in response to drought stress was found in the QTL regions for BY-DS, RL-ABA3, and RL-ABA5, on chromosome 3. Another candidate gene that is responsible for carbohydrate metabolism was detected on chromosome 9 (LOC_Os09g08120) in the QTL regions for BY-DS and RL-ABA5. Malate dehydrogenase is an enzyme that catalyzes the oxidation of malate to oxaloacetate by using NAD(H)/NADP(H) as a cofactor. Additionally, this enzyme can be expressed in different parts of the rice plants, such as root, leaf, panicle, and stem and was induced in the presence of water deficit (Nan et al., 2020). Transgenic plants over-expressing malate dehydrogenase exhibited increasing drought-resistance compared to wild-type. Malate dehydrogenase was also identified as a drought responsive protein by Agrawal et al. (2016). Under DS conditions, drought-resistant genotypes accumulate a higher level of malate dehydrogenase, that protects membranes from damage by reactive oxygen species (ROS) (Guo et al., 2018). By elucidating the function of malate dehydrogenase, the drought response in rice can be better understood.
A gene encoding photosynthesis function (LOC_Os04g44190) is present in the QTL regions for BY-WW, RL-ABA3, and RL-ABA5 on chromosome 4. This gene is involved in the light reaction of photosystem II (PSII) and is known to control stomatal closure, and protect plants from dehydration (Sasi et al., 2018). Under DS conditions, the photosynthetic rate is decreased due to reduction in photosynthetic electron transport and carbon assimilation, resulting in reduction of grain yield. Moreover, PSII is a pigment-protein complex in thylakoid membranes that is responsible for oxygen evolution, water splitting, and plastoquinone reduction (Lu, 2016). LOC_Os04g44190 encoding a PsbQ family protein that belongs to the class of PSII extrinsic proteins, and under DS conditions this protein changes expression due to change in PSII efficiency (Sasi et al., 2018). Therefore, PsbQ protein plays an important role in drought stress resistance.
A WRKY transcription factor that is involved in drought stress response and plant development was detected on chromosome 5 (LOC_Os05g49210) of the QTL regions for FGC and FGD. Shen et al. (2012) reported that transgenic rice with over-expression of OsWRKY30 showed improved drought resistance. Likewise, silencing WRKY genes in transgenic rice demonstrated increased drought sensitivity. In addition, expression of the WRKY transcription factor induced ABA accumulation under DS conditions, leading to stomatal closure and reduction in water loss (Chen et al., 2010; Schroeder et al., 2001). Yan et al. (2015) also reported that the expression of a WRKY transcription factor was increased by ABA treatment.
Genes present on chromosome 5 include LOC_Os05g49240 in the QTL regions for SP-DS and PPB-DS; and also on chromosome 10 (LOC_Os10g41460) in the QTL regions for BY-DS and RL-ABA3, were identified MYB transcription factor, a known transcription factor in drought response (Tang et al., 2019; Baldoni et al., 2015; Li et al., 2015; Dai et al., 2007; Ma et al., 2009; Yang et al., 2012; Xiong et al., 2014; Quan et al., 2010). Transgenic rice with over-expression of OsMYB6 exhibited increased resistance to drought compared to wild-type and contained higher proline catalase (CAT) and superoxide dismutase (SOD) activity. Additionally, OsMYB6 transgenic rice plants also showed higher expression of abiotic stress-responsive genes under DS conditions (Tang et al., 2019). Katiyar et al. (2012) also reported that the expression of MYB genes is controlled by drought. Increasing drought resistance correlated with over-expression of MYB genes and ABA accumulation (Xiong et al., 2014).
A stress-responsive transcription factor, Zinc Finger (ZFN) protein (LOC_Os06g49080), is located on chromosome 6 in the QTL regions for FG-WW, BY-WW, RL-ABA3, and RL-ABA5. The ZFN protein was reported to improve drought resistance in plants, suggesting that the ZFN protein contribute to the higher yield under DS conditions via control of stomatal closure (Huang et al., 2009; Ciftci-Yilmaz et al., 2000; Mukhopadhyay et al., 2004; Sakamoto et al., 2004).
A gene involved in sugar metabolism, OsSAC1 was found on chromosome 7 (LOC_Os07g02520) in the QTL regions for SP-DS. OsSAC1 regulates sugar partitioning in the carbon metabolism of young leaves and developing leaf sheaths (Zhu et al., 2018). OsSAC1 encodes an endoplasmic reticulum protein that causes sugar accumulation in the rice leaves and can be used to produce energy and construct carbon skeletons.
The genomic region on chromosome 8 (LOC_Os08g06344) in the QTLs for FG-DS, RL-ABA3, and RL-ABA5 encodes a DEAD-box RNA helicase which was reported to improve drought resistance in rice (Nawaz and Kang, 2019; Vashisht et al., 2006; Macovi et al., 2012). Over-expression of OsRH58, a chloroplast DEAD-box RNA helicase, in transgenic rice showed improved drought resistance, displayed by better survival rate than the wild-type under DS conditions (Nawaz and Kang, 2019). Furthermore, gene expression of the OsRH58 was increased under drought.
A gene regulating late embryogenesis abundant (LEA) protein was detected on chromosome 12 (LOC_Os12g02700) underlying the QTL regions for BY-DS. LEA protein has a major role in drought resistance in plants (Xiao et al., 2007; Duan and Cai, 2012; Magwanga et al., 2018; Liang et al., 2019; Chen et al., 2019; Kamarudin et al., 2019). Under DS conditions, LEA genes showed higher expression in the drought resistant plants compared to drought sensitive. In support of LEA functions, Xiao et al. (2007) reported that transgenic rice with over-expression a LEA protein gene OsLEA3-1 exhibited higher grain yield compared to wild-type under DS conditions.
Within the QTL region for SP-DS on chromosome 12, LOC_Os12g29330, (OsNAC139) was identified. OsNAC139 is a member of the NAC transcription factor family that is known to control plant response to drought (Kikuchi et al., 2003; Nakashima et al., 2007; Takasaki et al., 2010; Shim et al., 2018) by producing no apical meristem (NAM)/NAC protein. Rice contains 151 NAC genes (Puranik et al., 2013), from which several studies have reported that over-expression of OsNAC genes leads to improved drought resistance (Nakashima et al., 2007; Hu et al., 2008; Yokotani et al., 2009; Zheng et al., 2009; Jeong et al., 2010).
Among all the candidate genes identified within the QTL regions, various transcriptomes correlated with drought stress resistance were detected. Drought resistance of the rice plants can be either associated with metabolic regulation or osmoregulation. In conclusion, this study has revealed a number of potential candidate genes for developing drought-resistant rice varieties.
RT-qPCR validation of the key functional genes identified within the QTL regions regulating drought-related traits and ABA sensitivity
Information about the differences in the expression of drought resistance genes between drought-resistant and sensitive genotypes and their relationship to morphological characterisctics, grain yield components, and root architectural traits under DS conditions has been limited. However, the study of the differences in gene expression between resistance and sensitive genotypes could improve the efficiency and opportunities of developing drought resistant varieties.
The identified candidate genes within the QTL regions regulate morphological traits and grain yield components under WW and DS, and also root architectural traits under ABA treatments, were further analyzed to exemplify their roles in increasing drought stress resistance in rice by identifying their expression under different conditions. In this research, QTLs controlling morphological and yield traits under reproductive stage drought stress conditions were identified and compared to those identified under WW conditions. The QTLs under DS conditions were different from those identified under WW conditions. This might be due to the difference in expression of the genes for morphological and grain yield traits under DS and WW conditions. Plant respond and adapt to the drought stress through several processes, such as physiological, biochemical, and molecular processes that are regulated by transcriptional regulators. When rice plants are exposed to drought stress, certain genes are activated or repressed. Proteins, as the products of the activated genes, will protect the plants from the damage of drought stress (Dai et al., 2007). The known genes and transcription factor families that have been proven to regulate the plant response to drought stress are AP2/ERF, WRKY, MYB, NAC, NAP, and bZIP (Wu et al., 2017; Mao et al., 2017; Tang et al., 2017; Butt et al., 2017; Zhu et al., 2018; Sun et al., 2018). Understanding the regulation of gene expression in response to drought stress is important to develop drought-resistant rice varieties. About 26 out of 184 candidate genes obtained from the QTL regions were selected for identifying their expression between drought-resistant parent (Kaybonnet) and drought sensitive parent (ZHE733) under DS conditions. These 26 candidate drought resistance genes comprise of 15 candidate genes with known annotations to be responsive to drought stress, and 11 candidate loci comprising of genes from the traits with high LOD and PVE with unknown annotations (Supplementary Table 3).
Out of 15 annotated candidate genes, LOC_Os01g66270 (ERF/Ethylene Response Factor) and LOC_Os10g41460 (MYB protein) showed up-regulation in Kaybonnet relative to ZHE733 under DS conditions when compared to WW conditions (Fig. 9C). All other 13 candidate genes, LOC_Os02g54160 (APETALA2/ERF transcription factor), LOC_Os03g56280 (Malate dehydrogenase), LOC_Os04g44190 (Light reaction photosystem II), LOC_Os05g49010 (APETALA2/ERF transcription factor), LOC_Os05g49210 (WRKY transcription factor), LOC_Os05g49240 (MYB protein), LOC_Os06g49080 (Brassinosteroid/BR signaling), LOC_Os07g02520 (Regulation of sugar partitioning in carbon-demanding young leaves and developing leaf sheaths), LOC_Os08g06344 (DEAD-box RNA helicase), LOC_Os09g08120 (Carbohydrate metabolic process), LOC_Os11g30760 (Galactosyltransferase activity), LOC_Os12g02700 (Late embryogenesis abundant/LEA protein), and LOC_Os12g29330 (No apical meristem/NAM protein domain) eventhough showed down-regulation in Kaybonnet relative to ZHE733 under DS conditions when compared to WW conditions do not support the phenotypic traits associated with each loci. These genes are inherently up-regulated in WW (Fig. 9A) and DS conditions (Fig. 9B) compared to ZHE733, so eventhough the relative amount was down-regulated in DS compared to WW (Fig. 9C), the intrinsic gene expression values of that genes are higher in Kaybonnet compared to ZHE733 suggesting a cause for better DS phenotype of Kaybonnet compared to ZHE733.
In order to verify the results from QTL mapping for two candidate regions for the traits, the loci high LOD and PVE score were studied for gene-expression between Kaybonnet and ZHE733. One of the region for the trait PHC has LOD score 24.42 and PVE 13.81%. The window region of 25 Kb upstream and downstream was selected where the polymorphic SNP was detected on chromosome 2. This region has 6 candidate genes, LOC_Os02g44590, LOC_Os02g44599, LOC_Os02g44610, LOC_Os02g44620, LOC_Os02g44630, and LOC_Os02g44642 (Supplementary Table 3) with unknown annotations. The second region selected for gene-expression study also had high LOD score 21.82 and PVE 8.49% and is associated with traits BYD; and the region spanning the polymorphic marker selected included 5 genes: LOC_Os10g07030, LOC_Os10g07040, LOC_Os10g07050, LOC_Os10g07060, and LOC_Os10g07080 with unknown annotations on chromosome 10.
Four candidate genes (LOC_Os02g44590, LOC_Os02g44610, LOC_Os10g07040, and LOC_Os10g07080) out of 11 candidate genes with unknown annotations showed up-regulated in Kaybonnet relative to ZHE733 under DS conditions when compared to WW conditions (Fig. 8C). Seven other candidate genes (LOC_Os02g44599, LOC_Os02g44620, LOC_Os02g44630, LOC_Os02g44642, LOC_Os10g07030, LOC_Os10g07050, and LOC_Os10g07060) even though showed down-regulation in Kaybonnet relative to ZHE733 under DS conditions when compared to WW conditions (Fig. 8C), but the relative gene expression under DS (Fig. 8B) and WW conditions (Fig. 8A) showed higher expression in Kaybonnet and ZHE733, suggesting that higher intrinsic values of candidate genes in Kaybonnet compared to ZHE733 are probably enough for the traits. LOC_Os10g07040 and LOC_Os10g07080 up-regulated in Kaybonnet relative to ZHE733 under DS conditions when compared to WW conditions (Fig. 8C) and also showed higher relative gene expression under DS (Fig. 8B) and WW conditions (Fig. 8A) in Kaybonnet.
Among the candidate drought resistance genes not annotated to drought stress function, LOC_Os10g07040 showed high up-regulation in Kaybonnet compared to ZHE733 under DS and WW conditions (Figs. 9A, 9B, 9C), correlated with chalcone synthase, according to the MSU rice reference genome annotation release 7.0 that is involved in the drought stress response in rice (Hu et al., 2017), Arabidopsis (Nakabayashi et al., 2013), and tobacco (Hu et al., 2019). The other candidate drought resistance gene that also showed high up-regulated in Kaybonnet compared to ZHE733 under DS and WW conditions (Figs. 9A, 9B, 9C) was LOC_Os10g07080 is related to myosin (Jiang et al., 2004) and transposon protein (Cho et al., 2017) that regulates cell growth and developmental processes in rice. These candidate genes could be functioning in a cumulative manner in order to show a measurable positive effect on improving drought resistance in rice and the effect of genes can further be exploited to develop drought resistant cultivar.
A large number of genes were up-regulated in Kaybonnet (drought resistant parent), indicating that the drought resistant cultivar had higher capability to modulate drought resistance genes when exposed to DS conditions, thereby enhancing its resistance level compared to drought sensitive parent (ZHE733). Modulation of a higher number of up-regulated expressed genes with different transcription factor gene families is a crucial characteristics of drought resistant genotypes. Similar results were also obtained by Hayano-Kanashiro et al. (2009) who showed that drought-resistant maize genotypes inducted more genes compared to the sensitive genotypes under DS conditions. All 15 candidate drought resistance genes identified within QTL regions have been strongly associated with direct roles in drought stress resistance. For example, transcription factors MYB, NAP, NAC, ZIP, and APETALA2/ERF are responsive to dehydration induced by water deficit conditions (Tang et al., 2019; Baldoni et al., 2015; Li et al., 2015; Dai et al., 2007; Ma et al., 2009; Yang et al., 2012; Xiong et al., 2014; Quan et al., 2010; Lata and Prasad, 2011; Mizoi et al., 2012; Licausi et al., 2013; Phukan et al., 2017; Nakashima et al., 2007; Hu et al., 2008; Yokotani et al., 2009; Zheng et al., 2009; Jeong et al., 2010). These results provide strong evidence for genes expressed under DS conditions being involved in various physiological, biochemical, and molecular processes within the rice, in order to reduce the effects of drought stress, thereby enhancing their ability to resist the drought stress and maintain their grain yield production under DS conditions. Therefore, the up-regulation of the drought genes in Kaybonnet compared to ZHE733 provide important information to characterize the function of candidate drought resistance genes and to understand the drought stress mechanisms in rice.
Candidate genes within QTL regions involved in regulatory response to drought include a large family of genes expressed under DS conditions. Proteins expressed by known and candidate drought resistance genes played important roles in (1) cellular protection, including structural adaptation and osmotic adjustment, and (2) drought responses by interaction with other proteins and transcription factors, such as MYB, NAP, NAC, bZIP, and APETALA2/ERF. Under DS conditions, in drought resistant genotype Kaybonnet, exogenous ABA significantly improved the expression of ABA biosynthetic genes suggesting Kaybonnet genotype must be maintaining the water potential and cellular activity of the cell by closing the stomata.
Based on the RT-qPCR results, it may also be suggested that there is a correlation between gene expression, transcriptional regulation, and resistance to drought across resistance and sensitive genotypes. Therefore, the up-regulation of the drought genes and novel candidate genes in Kaybonnet compared to ZHE733 provided an important information to characterize the function of candidate drought resistance genes. All in all, these results enhance our understanding of the role of candidate drought resistance genes in the regulation of drought stress response, and this research has also revealed a number of potential candidate drought resistance genes that could be used to develop rice cultivars with greater drought resistance.
Genetic diversity in 26 loci across K/Z RIL population
The 26 drought loci described earlier are distributed across all 12 chromosomes with two distinct clusters on chromosome 2 (six loci) and chromosome 10 (five loci). The up- and down-stream region of the 26 drought loci was used for dissecting the genetic diversity in- and surrounding these loci in the diversity panel. Depending upon the proximity to the nearest neighbor gene, the region of interest was defined between 400 bp to 22Kbp for upstream, and 800 bp to 18Kbp for the downstream region.
We identified a total of 7,475 SNPs across these 26 loci in the 200-genotype of K/Z RIL population. The diversity panel represents six distinct subpopulations (Ind, indica; Aus, aus; TeJ, temperate japonica; TrJ, tropical japonica, Adm, admixture and Aro, aromatic). Each individual SNP has multiple effects on individual loci depending upon the ‘impact’, ‘functional class’, ‘type’ and ‘region’ where it is present as shown (Supplemental Table 7). The SNPs were predominantly found in the upstream/promoter and intergenic regions. Five SNPs resulted in a stop codon gained in three loci (one stop gained in LOC_Os09g08120.1 and two stops gained in LOC_Os10g07030.1 and LOC_Os10g07050.1 each). SNP annotation showed that although most variants are moderate/modifying in effect, there are seven high impact and 200 low impact SNP effects that are potential candidates for further testing and validation. Moreover, the two QTLs on chromosomes 2 and 10 containing six and five loci, respectively show a distinct pattern of variation as compared to the overall variation in 26 loci.
To determine how the variation present in these 26 loci differentiates the 200 genotypes, we performed a principal component analysis (PCA) on the vcf files containing the SNP datasets. We performed three separate runs of PCA for i) all 26 loci ii) six-loci cluster on chromosome 2 and iii) five-loci cluster on chromosome 10 containing 7,475 SNPs, 565 SNPs and 4,375 SNPs, respectively (Supplemental Figure 8A-C). In the PCA for 26 loci, the first principal component (PC1) explains 34% of the variation in the data, and clearly separates a subset of indica subpopulation from the rest of the genotypes (Supplemental Figure 8.A). This is in contrast to the PCA performed for 6.5 million SNPs representing the whole genome where the indica subpopulation segregates as a single group (Gill et al. unpublished data). SNPs on chromosome 10 cluster also show a similar pattern of separation of indica subpopulation into two distinct groups with PC1 explaining 49% of the variation in the data (Supplemental Figure 8.C). However, the SNPs on chromosome 2 cluster show an interesting pattern where PC1 clearly separates the japonicas (both temperate and tropical) from the indica and aus subpopulations, and explains 84% of the variation in the data in doing so (Supplemental Figure 8.B). This indicates the presence of novel variation in the six drought loci on chromosome 2 (LOC_Os02g44590, LOC_Os02g44599, LOC_Os02g44610, LOC_Os02g44620, LOC_Os02g44630, LOC_Os02g44642) that distinguish the japonicas from the indica subpolulations.