Linkage Mapping Revealed Non Association of the Green Revolution Genes (Rht1 and Rht2) with Drought Tolerance in Wheat (Triticum aestivum)

1. Introduction

The "Green Revolution" in the 1960s revolutionized global agriculture, particularly through enhancing the wheat production, through innovations including the introduction of dwarfing genes for the high yielding variety (HYV) development. Led by the International Maize and Wheat Improvement Center (CIMMYT), the movement significantly boosted yields in highly populated countries like India, Pakistan, and Mexico [1,2]. Green revolution played a key role in ensuring food security in past five decades through the introduction of the short stature, lodging resistant and fertilizer responsive high yielding wheat varieties. This major paradigm shift was mainly achieved by introducing semi-dwarfing Rht genes into the locally adapted low yielding, low fertilizer responsive and lodging prone tall traditional varieties. The reduced plant height gene in wheat located on two homeologs (Rht1 B1b and Rht1D1d) cause a reduced response to the gibberellin hormones which are well known for increasing plant height in plant species [3]. These semi-dwarf plants, with higher Harvest Index (HI), produced more robust spikes filled with heavy grains compared to their taller counterparts [4,5,6,7,8]. The green revolution varieties became popular among farmers in both irrigated and rain-fed ecologies, and led to a gradual replacement of the landraces and traditional varieties. Impact of the green revolution was more prominently observed in the irrigated areas but many semi-dwarf wheat varieties are usually grown on millions of hectares area in rainfed regions due to their high grain yields [9]. However, drought stress has been considered as a major limiting factor in the expression of Rht genes in yield enhancement [10,11]. Some researchers opined that drought, in particular, has a strong impact on yield in semi-dwarf and dwarf wheats compared to the taller ones [12]. The plant’s response to drought varies depending on its duration and severity, as well as the developmental stage of the plant [13,14]. The taller wheat landraces are comparatively less susceptible to drought unlike semi-dwarf high yielding varieties [15].

Interestingly, genomic regions harboring these Rht genes were found to be associated with drought tolerance [16]. As an example, in wheat a drought QTL, qDSI.4B.1 harboring Rht 1 B1b gene has been identified to be associated with drought tolerance [17,18,19]. However, linkage verses pleiotropy relation of the Rht genes with drought tolerance is still in debate. These facts clearly signify the importance of adjoining genomic regions of green revolution genes for drought tolerance. Unlike rice, there were two genes in wheat conferring semi dwarfness. Therefore, it becomes imperative to investigate the impact of drought stress on major effect Rht genes i.e. Rht1 B1b and Rht1D1d.

The high density genomics provides a powerful alternative approach for an in-depth genomic and genetic analyses in a way to solve complexities in trait understanding. The emphasis is beginning to move from genotyping-by-assay to genotyping-by-sequencing (GBS) thanks to next-generation DNA sequencing technology. The DArT-seq GBS technology, which enables users to select genome fractions that largely correspond to active genes, was developed by Diversity Arrays Technology (DArT), based in Canberra, Australia. Restrictions enzymes are used in this method to keep the low copy sequences from coming into contact with the repetitive DNA. These low copy sequences are helpful for the identification of markers. Representative snippets are then sequenced on Next Generation Sequencing (NGS) methods [20,21]. DArTseq GBS uses a combination of restriction enzymes to produce high-density SNPs and PAV (presence and absence variations) markers at a reasonable cost [22]. A typical DArT experiment involves sequencing roughly 200,000 genomic fragments ten times on average, with approximately 2,000,000 tags per sample. Most of the samples in each experiment are processed in duplicate, which enables the tight selection of markers based on score repeatability and eliminates any sequence variants that are not true SNP markers. Additional metadata generated by the analytical pipeline (DArTsoftS) makes it easier to select and filter markers further. According to [23], this enables clients to choose specific sets of markers that are most appropriate for their needs.

In this study it was hypothesized that the genetic linkage map of respective population would precisely look in to the probable collocation of drought QTLs and Rht genes in wheat population subjected to analysis. In this study we have used three RIL populations, PBW343 × Kingbird, PBW343 × Kenya Swara and PBW343 × Muu which were originally used for the construction of first DArT-seq genetic linkage map [24]. In addition to these three, another population was used for this study which segregated for both the Rht genes simultaneously.

2. Materials and Methods

2.1. Plant Material

There were four recombinant inbred line (RIL) populations were used for the study. One population was developed with the cross of a CIMMYT’s breeding line, Kachu/Kiritati/Kachu with a Mexican landrace ‘Jal 95.4.3’. The other three populations were derived from popular Indian wheat variety, ‘PBW343’ which were previously used for disease evaluation [24]. A total of 270 F_4:5 and F_4:6 recombinant inbred lines (RIL) of Jal 95.4.3 × Kachu/Kiritati/Kachu were used for the study. The F7 and F8 RILs of PBW343 derived populations were subjected to the evaluation. The three populations namely, PBW343 × KINGBIRD #1, PBW343 × Kenya Swara, PBW343 × MUU comprised 155, 191 and 124 genotypes respectively.

2.2. Phenotyping for Different Irrigation Regimes

Phenotypic evaluation for the key agronomic traits was carried out at two locations: (1) The three PBW343 derived were evaluated in crop seasons of 2015/2016 and 2016/2017 at CIMMYT’s experimental station near Ciudad Obregon (27 20°N, 109 54°W) and (2) 270 RILs of Jal 95.4.3 × Kachu/Kiritati/Kachu were evaluated in crop seasons of 2016/2017 at Indian Agricultural Research institute, Pus New Delhi, India (28.080°N 77.120°E). All experiments used similar agronomic management strategy, with the exception of irrigation. The experimental designs used for the field experimentation was alpha-lattices with two or three replications and plot sizes of 2.0 and 4.8 m² in both seasons as well as locations. All experiments were sown in the second fortnight of November. The well-watered experiment received about 600 mm of water, whereas the drought experiment had about 200 mm of total soil moisture over the crop season. Irrigation was provided in early crop growth stages (before tillering). The population was phenotyped for the agronomically relevant traits including- of days-to-heading (DTH), plant height (HT), days-to-maturity (DTM), grain yield per plot (GY) and NDVI. Data recording was performed as explained by [25].

2.3. DNA Extraction and DArT Sequencing

Genomic DNA of the Jal 95.4.3 × Kachu/Kiritati/Kachu RIL population was extracted from fresh leaf tissue of 2-week-old seedlings following the modified CTAB method standardized in CIMMYT followed by DNA quality & quantity check through using NanoDrop 2000 spectrophotometer (ND2000 V3.5, NanoDrop Technologies, Inc.). Further, DNA samples were shipped to Diversity Arrays Technology Pty Ltd, Canberra, Australia for destructive DNA analysis and genotyping using the DArTseq protocol using 38,611 silico DArTs. DArT loci with unknown chromosome positions were omitted from analysis followed by filtering of markers with more than 5% missing data, a total of 1170 markers distributed across the 21 chromosomes were used for analysis. The polymorphic information content (PIC) values of the silico DArTs that were used ranged from 0.02 to 0.50, the repeatability values were 1, the mean call rate was 0.93 with a range of 0.84 to 1, and the read mean depth was 14.92 with a range of 5 to 399. The DArT seq derived SNPs of three populations (PBW343 × KINGBIRD #1, PBW343 × Kenya Swara, PBW343 × MUU) were obtained from previous report of [24].

2.4. Diagnostic Marker Analysis for Phenological Traits

The Sequence tagged site (STS) markers associated with plant height, photoperiod and vernalization in wheat that are reported on MASWheat database (http://maswheat.ucdavis.edu/protocols/index.htm) i.e. Rht-B1, Rht-D1, Ppd-D1 and VrnA1 were used for genic characterization in the study. These gene-based markers were genotyped using PCR protocols and gel electrophoresis procedures described in this database.

2.5. Statistical Analysis

The Analysis of variance (ANOVA) was conducted for all traits separately for estimating variance components for evaluation of the significance of genotype, treatment and trial effects and their interactions in the three RIL populations PBW343 × Kingbird, PBW343 × Kenya Swara and PBW343 × Muu.

The statistical analysis was performed with the help R-project version 3.1.1[26]. The following linear mixed model for the analysis of variance was used to estimate the phenotypic means of the entries:

$${\mathrm{P}}_{\mathrm{ijk}}=\mathrm{M}+{\mathrm{R}}_{\mathrm{i}}+{\mathrm{B}}_{\mathrm{j}}\left({\mathrm{R}}_{\mathrm{i}}\right)+{\mathrm{L}}_{\mathrm{k}}+e_{\mathrm{ijk}}$$

The measurement on a plot was P_ijk, the replications, blocks, lines, and errors were represented by R, B, L, and e, respectively; the entry mean is represented by M. Replications and blocks within replicates were taken as random and entries were considered as fixed variables while estimating the entry means. Season effects were similarly treated as random when evaluating the entry means across years. Pearson’s coefficient was followed for correlation analysis among different traits. The broad sense heritability was estimated while considering all factors including genotypes as random. It (H) was calculated using formulae:

$$H={\sigma }_g^2/{\sigma }_p^2\mathrm{and}{\sigma }_p^2={\sigma }_g^2+\left({\sigma }_e^2/r\right)$$

Where,

σ²_p = phenotypic variance,

σ²_g = genotypic variance,

σ²_e = error variance and

r = replications in each season.

3. Linkage Map Construction and QTL Mapping:

To test the segregation of markers with a 1:1 segregation ratio, chi squared tests were used. Recombination frequencies were converted into centi Morgan (cM) values using the Kosambi function [27]. The order of SNP markers and distances between adjacent ones were determined using Join Map 3.0 software [28]. The QTL analysis was performed with ICIM 2.0. For every trial, a QTL analysis was performed for every trait. Analysis was performed following the forward and backward stepwise regression with a window size of 10 cM and a walk speed of 2.0 cM, QTLs were found using Composite Interval Mapping (CIM). Each trait's LOD score for QTL significance varied based on 1,000 permutations (Alexander et al. 2012). Using single-factor analysis from a general linear model approach, the coefficient of determination (R2) was used to assess the percentage of phenotypic variation (PV) explained by a QTL.

4. Results

4.1. Phenotypic Evaluation of RILs Under Drought and Well-Watered Conditions

All four populations revealed phenotypic variations for the traits under investigated. In the PBW × Kingbird population, DTH, PHT and GY ranged 60 - 89 days, 65 - 98 cm and 522 - 6705 Kg/ha respectively across drought-stress experi ments of two years. Similarly, under well-watered experiments of two seasons DTH, PHT and GY ranged 65 - 93 days, 73 - 95 cm and 1193 - 8181 Kg/ha respectively. DTH, PHT, and GY ranged from 66 - 88 days, 44 - 117 cm, and 960 - 6498 Kg/ha, respectively, throughout two years of drought-stress experiments in the PBW × Muu population. In the same way, throughout two seasons of well-watered tests, the ranges for DTH, PHT, and GY were 62 - 90 days, 45 - 117 cm, and 1072 - 9179 kg/ha, respectively. During two years of drought-stress tests in the PBW × Kenyaswara population, the ranges of DTH, PHT, and GY were 58 - 86 days, 70 - 118 cm, and 333 - 5644 Kg/ha, respectively. Similarly, during two well-watered test seasons, the ranges for DTH, PHT, and GY were, respectively, 58 - 87 days, 74 - 116 cm, and 601 - 8573 kg/ha. The Jal 95.4.3 × Kachu/Kiritati/Kachu was screened under drought and well-watered conditions of crop season 2016-17. DTH, PHT, and GY ranged from 60 - 84 days, 71 - 117 cm, and 800 - 4520 Kg/ha, in drought stress, and, 65 - 86 days, 80 - 113 cm, and 700 - 9950 Kg/ha in well-watered situations respectively (Supplementary Tables S2).

The grain yield reduction in drought stress experiments of two seasons was maximum in PBW × Kingbird (32.0, 10.4%) followed by PBW × Kenyaswara (31.1, 9.6%) and PBW × Muu (26.7, 8.0%). Therefore the grain yield reduction in first season was around three times as compared to the second. All three populations depicted similar trend in terms of drought severity. In the ‘Jal 95.4.3 × Kachu/Kiritati/Kachu’ this reduction was 32.1% in drought stress experiment of 2016-17. The traits investigated in this study in all four populations showed normal/ nearly normal distribution pattern (Supplementary Figures S1 a-d). The means, variances and correlations have been presented in Supplementary Tables S2 & S3).

4.2. Linkage Map of Jal 95.4.3 × Kachu/Kiritati/Kachu Population

In total, 1170 markers were mapped in ‘Jal 95.4.3 × Kachu/Kiritati/Kachu’ population with 941 unique positions. The ratio of unique positions on the linkage map varied from 67 to 97 %. A total of 39.8%, 45.6% and 14.6% markers were mapped on the A, B and D genomes, respectively. In the linkage map, A, B and D genomes covered distances of 945, 1122 and 803 cM respectively. The total genetic length of the consensus map was 2870 cM, and average marker distance was 3.04 cM, reached by calculating the average distance between two adjacent unique positions (Figure 1). The number of markers varied from 12 (4D) to 144 (2B), and unique positions 10 (1A) to 115 (2B) respectively (Figure 2, Table 1).

4.3. QTL Mapping

The QTL analysis in ‘Jal 95.4.3 × Kachu/Kiritati/Kachu’ population revealed a total of five genomic regions. The QTL, qDWW-5A was associated with the days to heading explaining phenotypic variances of up to 12.63%. Another QTL on chromosome 5A (qww-5A.1) showed significant association with days to heading. For the plant height, three QTLs were identified on chromosome 4B, 4D and 7A explaining phenotypic variance in the range of 2.63 -13.19% (Table 2).

The plant height QTLs qDWW-4B.2 and qDWW-4D.1 co-segregated with Rht-B1 and Rht-D1 genes respectively. The qDWW-4B.2 and qDWW-4D.1 QTL flanks showed allele similarity with the Rht-B1 and Rht-D1 genes up to 97 and 94% (Figure 3A). Similarly, qDWW-5A.1 and qWW-5A.2 shared similar allele pattern with Vrn-1 gene markers up to 95 and 97% respectively (Figure 3B).

A total of five consistent QTLs were identified on the five different chromosomes in three PBW343 derived populations. qDWW-1B.1 showed significant effect on grain yield under drought as well as well-watered conditions of years 2016-17 and 2015-16 respectively in ‘PBW/Muu’ population. Phenotypic variation explained by this QTLs was above 10% in both environments. Similarly, qDWW-2B.1 was found significantly associated with DTH and DTM in ‘PBW/Muu’ population and explained phenotypic variation up to 12.66 %. In the same population qDWW-4B.1 explained phenotypic variation of 11.2-13.8 % for plant height and for the NDVI, 7.8% (Table 3).

The QTL, qD-2D.1 explained phenotypic variation up to 13.56% in the ‘PBW/Kenyaswara’ population for NDVI. On chromosome 5A, qDWW-5A.1 depicted a clear consistent effect over two seasons for DTH under both drought and well-watered conditions. The QTL explained a phenotypic variation in range 10.55 – 19.19 % (Table 3).

None of the QTLs were found to show consistent effect across the three populations. Also, we could not detect any QTL showing significant effect in the ‘PBW/Kingbird’ population.

4.4. Rht Gene Segregation Pattern in RIL Populations

The PBW × Kingbird population segregated neither for Rht-B1 nor for the Rht-D1, whereas, PBW × Kenyaswara segregated for Rht-B1 gene. The PBW × Muu and ‘Jal 95.4.3 × Kachu/ Kirtati/ Kachu’ populations segregated for both the Rht genes i.e. Rht-B1 and Rht-D1.

The 'TT' allele of Rht-B1 gene showed advantages of 12.8 and 30.1 % over 'CC' allele in 2015-16 well watered and drought conditions respectively in PBW × Kenyaswara. Similarly, the 'TT' allele of this gene showed advantage over 'CC' allele of 1.2 and 3.95 % in 2016-17 well watered and drought conditions respectively in the same population. The 'CC' allele of Rht-B1 gene depicted yield advantage of 14.2 and 0.1 % over 'TT' allele in 2015-16 well watered and drought conditions respectively in PBW × Muu population. In the 2016-17 season, 'CC' allele had 0.6% yield enhancement over 'TT' under drought stress situation, whereas, under well-watered condition, the 'TT' allele resulted yield increase of 8.7% over 'CC' allele of Rht-B1 gene. In the ‘Jal 95.4.3 × Kachu/ Kirtati/ Kachu’ population, the 'CC' allele of Rht-B1 gene showed slight yield enhancements of 0.6 and 3.7 % under well-watered and drought stress situations respectively (Table 4).

The 'GG' allele of Rht-D1 gene yielded 4.5 and 11.0 % higher than 'TT' under well-watered conditions of 2015-16 and 2016-17 respectively in the PBW343 × Muu population. Under drought stress environments of 2015-16 and 2016-17, 'GG' allele yielded 1.3 and 12.7% higher than the 'TT' allele respectively in the same population. Contrarily, in the ‘Jal 95.4.3 × Kachu/ Kirtati/ Kachu’ population, 'TT' allele showed a mild yield increase of 1.4% over the 'GG' allele under well-watered situation, whereas, under drought stress environment, 'GG' allele depicted 3.6 % yield enhancement over 'TT' allele in the same season i.e. 2015-16 (Table 4).

5. Discussion

In two-season drought stress trials, the PBW × Kingbird population experienced the greatest drop in grain yield, followed by the PBW × Kenyaswara and PBW × Muu populations. Consequently, compared to the second season (8 - 10.4%), the first season's decline (26.7 - 32%) in grain yield was around three times greater. These results indicated that temperature effect might have impacted the grain yield reduction along-with drought stress in the first season. It was reported that drought stress alone may cause up to 50% grain yield loss in wheat [29]. According to few reports, moisture stress alone may cause a 30% decrease in wheat grain yield in Europe [30]. Recent research on drought-related meta-analyses of wheat showed a yield drop of 20–25% with a 40% water deficit [31,32]. In C306, under combined drought and heat stress treatment, there was a 30% decrease in grain weight per spike [33]. Reduction in yield in the RIL populations investigated in our study could be due to reduced growth and development in plants under stress condition as depicted by reduced plant height in population lines. It is also reported that due to significant decreases in plant development and shoot production, drought stress can lower wheat yields [34]. Further, previous reports also suggest that during the reproductive stage, drought stress results in sterility and premature floret abortion [35,36]. The drought stress related results in our experiments corroborate with the published research reports including meta-analysis studies. Similar trend depicted by three populations in terms of grain yield reduction under drought in two consecutive years, normal distribution of traits investigated, means and variances provide validatory evidences toward importance of the data for genetic analysis (Supplementary figures S1 a-d, Supplementary Table S2).

The RIL populations derived from PBW343 and Jal 95.4.3 × Kachu/Kiritati/Kachu were genotyped by DArT-seq approach. Linkage distances in PBW343 populations reported were used in this study for identifying the marker-trait associations [37]. The ‘Jal 95.4.3 × Kachu/Kiritati/Kachu’ population revealed a total genetic distance of 2870 cM with an average density of 3.04 cM in this study (Table 1). Similar to our report, in year 2020 a group of scientist genotyped a RIL population with the DArT-seq markers and developed a linkage map with 4439 markers covering a genetic distance of 2851 cM, with average marker density of 1.6 cM [38]. In 2023 group of scientist led by Rathan reported genotyping of a RIL population with 909 DArTseq markers that comprised a total genetic length of 4665 cM [39]. Our analysis with different RIL populations clearly validated that an adequately segregating population can be efficiently utilized for genetic mapping using DArT-seq platform.

The response of Rht-B1 alleles for grain yield was not consistent across populations in our study (Table 4). In 2005 Butler et al., also reported similar inconsistent results [40]. The allele which increased grain yield in PBW × Kenyaswara population in 2015-16 season decreased in ‘PBW × Muu’ and ‘Jal 95.4.3 × Kachu/Kiritati/Kachu’ populations. Similarly, the grain yield increasing Rht-D1 allele in ‘PBW × Muu’ population decreased grain yields in ‘Jal 95.4.3 × Kachu/Kiritati/Kachu’ population. Results clearly indicate the non-association of Rht-B1b and Rht-D1b genes with grain yield under reproductive stage drought stress in wheat. In 2016, Perez et al., group reported somewhat contrasting results after evaluating 158 RILs under field trials in Canada in rainfed conditions, 2007-2008, 2012 [41]. In this study, the class with wild-type Rht-D1a allele showed 436-322 kg/ha more grain yield as compared to its counterpart RIL class with semi-dwarf allele Rht-D1b. Most of the studies carried out in relation to Rht alleles are focused on early drought stages or uncharacterized stress environments [41]. Unlike previous reports, drought stress in our experiments were was not imposed in early growth stages which probably, nullified the impact of association of tall allele (Rht-B1a/Rht-D1a) with coleoptile length. In 1997, Flintham et al., suggested that the advantage of semi-dwarfing alleles (Rht-B1b or the Rht-D1b) is reduced in drought or heat stressed environments, due to decreased seedling emergence caused by reduced coleoptile length [42]. Most modern semi-dwarf wheat varieties harbouring Rht-B1b or Rht-D1b have short coleoptiles and low yields under drought stress relative to tall plants [43,44]. Further, a group of scientists provided evidences of association of Rht-B1b allele with reduction in coleoptile length in wheat [45]. The coleoptile length plays an important role in mitigating losses under early drought stress environments, so, in comparison to semi-dwarf genotypes (Rht-B1b/Rht-D1b), the taller cultivars with Rht-B1a/Rht-D1a allele combination show a higher emergence rate and ultimately the grain yield [46,47].

The QTL analysis results in the ‘Jal 95.4.3 × Kachu/Kiritati/Kachu’ population clearly indicates toward significance of the Vrn-1 genes for drought adaptation. This is a flowering related gene and escape is practically most important adaptation mechanism for mitigating losses due to drought stress situations. The association of qWW-5A.2 with NDVI further validated its significance for drought mitigation (data not presented). The Rht-B1 and Rht-D1 genes showed phenotypic variances of up to 4.11 and 13.19 % in drought and up to 6.24 and 12.75 % in the well-watered conditions. These results clearly reveal the non-association of Rht-B1 and Rht-D1 genes with drought adaptation in wheat ‘Jal 95.4.3 × Kachu/Kiritati/Kachu’ population. A novel QTL for plant height was identified on chromosome 7A. Even though this QTL explained phenotypic variances of 3.18 and 4.18 under drought and well-watered conditions respectively, their LOD significance values were quite high i.e. 6.68 and 7.87 which suggests toward its in-depth analysis in future studies (Table 2).

Results presented in this study can be efficiently utilized in formulating a molecular breeding strategy focusing the accumulation of positive alleles in the pre-breeding and breeding germplasm pools that can serve as foundation of the next generation wheat varieties for climate resilience.