Unraveling the genetic complexity underlying sorghum response to water availability

Understanding the adaptation mechanisms of sorghum to drought and the underlying genetic architecture may help to improve its production in a wide range of environments. By crossing a high yielding parent (HYP) and a drought tolerant parent (DTP), we obtained 140 recombinant inbred lines (RILs), which were genotyped with 120 DArT and SSR markers covering 14 linkage groups (LGs). A subset of 100 RILs was evaluated three times in control and drought treatments to genetically dissect their response to water availability. Plants with early heading date (HD) in the drought treatment maintained yield (YLD) level by reducing seed number SN and increasing hundred seed weight (HSW). In contrast, early HD in the control treatment increased SN, HSW and YLD. In total, 133 significant QTL associated with the measured traits were detected in ten hotspot regions. Antagonistic, pleiotropic effects of a QTL cluster mapped on LG-6 may explain the observed trade-offs between SN and HSW: Alleles from DTP reduced SN and the alleles from HYP increased HSW under drought stress, but not in the control treatment. Our results illustrate the importance of considering genetic and environmental factors in QTL mapping to better understand plant responses to drought and to improve breeding programs.


Introduction
Sorghum bicolor L. Moench. is native to arid and semi-arid tropical environments and a drought-tolerant cereal. In general, sorghum growing seasons in Sub-Saharan Africa are characterized by initial rainfalls with subsequent periods of drought. Sorghum plants with high vigor and fast growth rates during early developmental stages may be advantageous in regions affected by drought early in the season [1]. A plant's response to drought can be categorized into three adaptive strategies, i.e., drought escape (e.g by early flowering), drought tolerance (e.g. by improving water-use efficiency), and drought avoidance (e.g. by increasing water uptake and reducing water loss, [2][3][4]. Evaluating the natural variation of these responses by testing large numbers of genotypes in several environments improves the understanding of PLOS  genotype by environment interactions (G×E) ( [5], which in turn allows to select breeding lines with improved yield stability and helps to identify superior alleles across different environments [6]. Dissecting the genetic components underlying G×E can be achieved via mapping quantitative trait loci (QTL) and their effects in different environments, i.e., by estimating QTL by environment interaction (Q×E) effects [7][8][9].
The common approach to explore QxE is to use multi-environment analysis [7,8]. However, to develop indirect selection strategies for yield via its genetically correlated components, the multi-trait multi-environment (MTME) approach [10] is superior over independent multi-environment or multi-trait approaches [11]. The power of MTME is basically due to its ability to map QTL with different effects, such as QTL with synergistic pleiotropic effects, where one allele affects two or more traits in the same direction. If the allele is antagonistically affecting two different traits, one allele improves one trait while the other allele enhances other traits. Another possibility is the conditional neutrality of a QTL, i.e. the QTL has an effect on a trait in one environment but has no effect on the same trait in other environments. Additionally, the same allele can have unequal effects on the same trait in two environments, i.e. the effect is strong in one environment but weak in other environments [9,11]. Considering such effects in breeding programs is crucial since selecting for QTL, which are mapped in single environments, might lead to undesired responses in other environments [9].
In the present study, we focused on investigating the relation between vegetative growth of shoots and roots and yield components as well as yield. To genetically dissect these relationships, we used 100 recombinant inbred lines (RILs) genotyped with 120 DArT and SSR markers.
We used the raw data from three experiments in the present manuscript, although a transformed form of parts of the data was used earlier in a manuscript, which was submitted, published and later retracted from the Journal of Agricultural Sciences [12]. The retracted paper was published without permissions of all coauthors. Retraction was mainly done due to the mentioned fact and due to several minor errors and mistakes throughout the manuscript: Even though the data included yield components and yields, the title suggested an analysis of only vegetative traits. The discussion of the manuscript submitted earlier did not include comparisons with earlier results, which are covered in the present paper. The retracted paper did not account for interactions between genotypes and water availability. Consequently, the QTL mapping approach did not distinguish between QTL with main, conditionally neutral and antagonistic pleiotropic effects. In the present study, we overcame those limitations by using the more powerful MTME-QTL-mapping-approach.

Developing and genotyping the RIL population
The RIL population was developed at the Grain Crops Institute (GCI), Potchefstroom, South Africa from a cross between a high yielding parent (HYP) with superior grain quality under normal conditions and strongly reduced yields under drought stress and a breeding line with intermediate yielding abilities, which was described as drought tolerant (DTP). From the F 1 , 140 RILs were advanced to F 4 by selfing.
DNA was extracted from leaf tips of F 4 seedlings using the cetyl-trimethylammonium bromide (CTAB) method. Genotyping was carried out at Diversity Arrays Technology Pty. Ltd. (DArT), Yarralumla, Australia using 184 polymorphic DArT markers. In addition, nine informative expressed sequence tag (EST) derived simple sequence repeat (SSR) markers [13] were used. Polymerase chain reaction (PCR) was carried out using a T-Gradient PCR machine (Biometra, Göttingen, Germany). The PCR protocol had a denaturation temperature of 94˚C, an annealing temperature of 52˚C, and a polymerization temperature of 72˚C. The first 25 cycles with 30 s for each step were followed by eight cycles with an annealing time extended to 45 s and a polymerization time of 60 s. We used DY-682 labeled M13 primers in the PCR reactions (Eurofins MWG, Ebersberg, Germany). Amplification products were separated by polyacrylamide gel electrophoresis using an LI 4200 sequencer (Licor Inc. Lincoln, USA). The genetic map was constructed with the 193 markers using JoinMap 4 (Kyazma B.V., Wageningen, The Netherlands, www.kyazma.nl) and the multipoint maximum likelihood mapping function [14]. Low informative markers such as monomorphic markers, markers with a high number of missing scores and those with more than 75% allele skewedness towards either A or B alleles were removed.

Experimental setup and plant phenotyping
Three experiments were carried out under controlled greenhouse conditions using a subset of 100 RILs, which were grown in a control and a drought treatment with two replications arranged in a completely randomized block design. The first two experiments were conducted in order to estimate pre-flowering drought stress effects on plant growth during the vegetative phase. The third experiment was carried out in order to estimate drought stress effects on yield and yield components. Plants grown in the control treatment were watered every second day, whereas plants grown under drought stress were watered until most plants were in the fourthleaf stage. The greenhouse conditions were kept at average day/night temperatures of 25.8/ 15.9˚C and mean day/night relative humidity was 37.4/65.2%. The lengths of the stress cycles were 18, 21 or 43 days, respectively, in the three experiments.
In the first experiment, two seeds were sown in each of the 12.5 x 50 cm polyvinyl chloride pots filled with 9.4 kg dry sandy soil and supplemented with 1100 ml nutrient solution, which corresponds to 80% of the maximum soil water holding capacity (WHC). After emergence, plants were thinned to a single plant per pot and fertilized with 0.15% Scotts Universal Orange (Scotts Marysville, Ohio, USA) solution (N:P:K 16:6:26) twice a week. Evaporation was minimized by covering the soil with 200 g of gravel. The second experiment was basically a replication of the first experiment, while both focused on growth and development of vegetative plant parts before flowering. The third experiment was to analyze effects of drought on grain yield and yield components.
After harvesting the first and the second experiments, leaf area (LA) was measured using a LI-3100 area meter (Licor Inc., Lincoln, NE, USA). Harvested plants were stored in plastic bags until LA measurement to minimize errors due to transpiration losses and senescence. Roots were washed carefully, placed in a water bath and scanned with a flatbed scanner. Total root lengths (TRL) were measured using WinRhizo (Regent Instruments Inc., Quebec, Canada). Dry weights of leaves, stems (SDW), and roots were measured after drying plant parts at 105˚C until weight constancy. Above ground dry matter (AGDM) was calculated as the sum of LDW and SDW. Specific leaf area (SLA) was calculated as ratio between LA and LDW. For the third experiment, heading date (HD) was determined, seed number (SN) per plant was counted, hundredseed weight (HSW) was measured and (YLD) was calculated.

Statistical analysis
Statistical analysis was performed using SPSS 2. We used the following mixed linear model with fixed environment and random genotype effects: where μ is the overall mean, α i is the effect of the ith environment, ß j is the effect of the jth genotype, γ ij is the interaction effect of the jth genotype with the ith environment and ε ijk is a random error." Broad-sense heritability for each trait was estimated as the ratio between the genetic variance Vg, i.e. the variance among all lines, and the total phenotypic variance Vt, i.e. Vt = Vg + Ve, where Ve is the error variance, i.e. the variance among replications.

QTL mapping using the multi-trait multi-environment approach
Means of the 18 measured traits (S1 Table) were used to map QTL with Genstat 16 (VSN International, Hemel Hempstead, UK) using the multi-trait multi-environment analysis (MTME) approach [10,11]. The whole genome was scanned by simple interval mapping (SIM) with a distance of 30 cM to separate selected QTL. To estimate the allelic effect and the explained phenotypic variance of each QTL per trait and treatment, backward selection on the significant cofactors was used.

Linkage map construction and population phenotyping
We constructed a genetic map with 14 linkage groups using a total of 120 markers, i.e. 112 DArTs and 8 SSRs markers. The map covered 1212 cM with an average marker distance of 10 cM (Fig 1 and S2 Table).
Means of the 18 measured traits of the parental lines and the RIL population are shown in Table 1. LA, TRL and RDW of parental lines showed antagonistic response to water availability, while HYP showed larger values in the control treatment, DTP showed larger values in the drought treatment. A similar situation was observed in case of HD since HYP was earlier in the control treatment while DTP was earlier under drought stress. HYP had higher SN, HSW and YLD than DTP in both treatments. Transgression beyond the two parents was observed for all traits. Heritability was moderate to high for all measured traits and ranged between 0.51 and 0.91. ANOVA showed significant GxE for LDW, SDW, AGDM, TRL and RDW in the first but not in the second experiment. HD, SN, HSW and YLD showed significant GxE in the third experiment.

QTL mapping using the MTME approach
In total, 133 significant QTL (p < 0.05) were detected in ten hotspot regions and mapped mainly to LGs 2, 2a, 3, 5, 6, 8, 9, 9a, and 10a (Fig 3, Tables 5 and 6). QTL of a hotspot on LG-3 had positive effects from the HYP allele for LA_C1, LA_C2, LA_D2, LDW_C1, LDW_C2, LDW_D2, SDW_D1, AGDM_C1, AGDM_C2, AGDM_D2, RL_D2, RDW_C1, and RDW_D2. The hotspot showed conditional neutrality for some traits since no QTL were mapped under drought stress in the first experiment for LA_D1, LDW_D1 and AGDM_D1. The main effect QTL cluster on top of LG-5 had a positive effect from DTP for LA_C1, LA_C2, LDW_C1, LDW_C2, SDW_C2, AGDM_C1, AGDM_C2, TRL_C1, TRL_C2, RDW_C1, RDW_C2, and HD_D3. A conditional neutrality for several traits was observed here as well since no QTL were found for LA, LDW, AGDM, TRL, and RDW in the drought treatment. In another cluster on top of LG-6, the DTP allele had positive effects on SDW_C2, SDW_D1, AGDM_C2, AGDM_D1, SLA_D2, TRL_C2, HD_C3, HD_D3, and SN_D3, while the HYP allele had positive effects on SLA_C1, HSW_C3, HSW_D3, SN_C3, YLD_C3, and YLD_D3. Within this cluster, an antagonistic effect was observed for SN since the positive effect was from the DTP allele in the drought and from the HYP allele in the control treatment. Considering the opposite allelic effects on different traits revealed antagonistic pleiotropic effects between HD and both HSW and YLD since the DTP allele had a positive effect on HD and the HYP allele increased HSW and YLD in both environments.

Discussion
The population used here was genotyped with 120 DArT and SSR markers covering 14 LG and a total length of 1212 cM, which is comparable with the length of the sorghum consensus map, which had a size of 1355.4 cM [15]. DArT markers were used because they are affordable and represent a powerful high-throughput marker system suitable for QTL mapping. However, we are aware that the use of additional SNPs would be necessary to provide equal genome  coverage and to allow direct comparisons with recent or future studies, in which SNPs are used [16][17][18].
Observing the performance of the two parents under drought stress revealed that both parents avoided drought by increasing their TRL and reducing their LA, for more water uptake The map represents the 14 linkage groups in columns and shows significant QTL across all traitenvironment combinations using the multi-trait-multi-environment (MTME) approach. Light to dark blue indicates a significant positive effect from the DTP allele and yellow to red indicates a significant positive effect from the HYP allele (the darker the color, the higher the significance). C and D refer to control and drought treatments, respectively. The numbers 1, 2, and 3 refer to the three experiments. LA refers to leaf area, LDW to leaf dry weight, SDW = stem dry weight, AGDM to above ground dry matter, SLA is specific leaf area, TRL refers to total root length, RDW to root dry weight, HD to heading date, HSW is hundred seed weight, SN refers to seed number, and YLD to yield. https://doi.org/10.1371/journal.pone.0215515.g003 Genetics of sorghum response to drought and reduced transpiration. However, both parents did not escape drought by earliness. LA of the whole RIL population showed responses similar to the parental lines in both treatments of the first two experiments. In contrast to the two parents, the RIL population on average exhibited early heading in response to the drought treatment. Drought stress effects on vegetative  plant growth, i.e. LDW, SDW and AGDM, were more severe during experiment 1, which is probably resulting from higher temperatures and radiation, since the experiment was carried out during spring/summer in Germany, while the third experiment was conducted in autumn. However, radiation and temperature were not explicitly measured during the experiments, so that we are not able to draw clear conclusions from these results. Significant correlations between the measured traits allow to better understand response patterns to water availability. For example, plants heading early under drought stress reduced SN and increased HSW to maintain high YLD, a commonly observed mechanism and fitness tradeoff if stress occurs before or at anthesis. In contrast, plants heading early in the control treatment showed increased SN, HSW, and YLD, which are all desirable traits for selecting QTL were mapped using the multi-trait multi-environment analysis.
LG refers to linkage group and pos refers to marker position in cM. R 2 is the percentage of total phenotypic variance explained by each QTL. Effects with positive values represent a positive contribution of the DTP allele to the trait, while negative values represent a positive contribution of the HYP allele. C and D refer to the control and drought treatments, respectively. Numbers 1, 2, and 3 refer to the three experiments. LA = leaf area, LDW = leaf dry weight, SDW = stem dry weight, AGDM = above ground dry matter, SLA = specific leaf area, TRL = total root length, RDW = root dry weight, HD = heading date, HSW = hundred seed weights, SN = seed number and YLD = yield. https://doi.org/10.1371/journal.pone.0215515.t006 Genetics of sorghum response to drought breeding lines. Such adaptive responses can be genetically dissected and explained by the antagonestic pleiotroic effects of the QTL clusters mapped to LG-2, LG-2a, and LG-6, which were identified since GxE was incorporated into the QTL model. For example, the DTP allele on LG-6 had positive effects on HD in both treatments and on SN in the drought treatment, whereas the HYP allele showed positive effects on HSW and YLD in both treatments and on SN in the control treatment.
Our results revealed significant G×E effects for all traits measured in the two treatments (control and drought) in the first and the third experiments but not in the second one. However, significant environment and genotype effects were observed for all measured traits in the second experiment. This was reflected by significant Q×E for all traits and enabled us to distinguish the QTL effects. For example, three QTL mapped for LA on top of LG-3, LG-5 and LG-8, showed conditional neutrality, they were associated with water availability in the first and second experiment.
We expected that the 10 QTL hotspots of the present study would overlap with previously mapped QTL for similar traits. To facilitate a comparison to previous studies, marker positions were compared based on the sorghum consensus map [15,19]. An earlier study that used DArT markers [20] mapped a QTL for grain yield that was associated with sPb-3361 at 140.7 cM on chromosome 2. The QTL cluster on top of LG-2a was associated with the marker sPb-2229, which was mapped at 142.9 cM in the consensus map. The QTL cluster on LG-2a included QTL for several traits including HSW in the drought treatment. Another study detected a QTL for stay green and panicle length on Chromosome 2 [13]. Another QTL for chlorophyll fluorescence was mapped at DArT markers sPbn-2229 on LG-2a [21]. The QTL cluster on LG8 was associated with the marker sPb-7889 which mapped at 74 cM in the consensus map. A stay green QTL was associated with sPb-1661 which was mapped at 73.9 cM in the consensus map [20]. The QTL cluster on LG-9 overlapped with an earlier detected QTL for maturity [22] and another QTL with pleiotropic effects on flowering time and HSW [23]. We mapped QTL for several traits on LG9 at 52.5 and 140 cM according to the sorghum consensus map. These results indicate possible pleiotropic effects of the QTL on morphological traits and yield components, which were already proven in recent studies [24,25]. The main QTL cluster on LG-6 was associated with the majority of traits measured in both water regimes. Earlier studies on sorghum reported a number of significant QTL on chromosome 6 that were associated with several traits measured under drought [26] as well as other environmental constraints such as thermal [27] and cold [18] stress and sorghum ergot [28]. Altogether, these results indicate the major role of chromosome 6 on sorghum growth and development under various environmental conditions making it an interesting target for future breeding programs.
Since the pot height used in the present study was much smaller than maximum rooting depth of sorghum plants, it was expected that investing additional energy in developing longer roots did not improve water uptake or increase YLD. Therefore, we assume that the negative correlation observed between TRL and both HSW and YLD in the drought treatment is an artifact of pot size.

Conclusion
Understanding crop response to drought and the underlying QTL is essential to increase crop productivity under drought conditions which is the ultimate goal for breeding programs.
In that respect, mapping the QTL cluster on LG 6 with the observed antagonistic pleiotropic effects is very important to genetically dissect the significant antagonistic response by reducing SN and increasing HSW and YLD under drought conditions as an adaptive mechanism to cope with drought stress. In total, we detected 14 QTL clusters mapped on 11 LGs for the measured traits as a first step towards identifying genes governing those traits.
Supporting information S1 Table. Phenotypic data used for QTL mapping. Trait appreviations are as following; Leaf area (LA), total root lengths (TRL), leaves dry weight (LDW), stems dry weight (SDW), roots dry weight (RDW), above ground dry matter (AGDM) was calculated as the sum of LDW and SDW, Specific leaf area (SLA) was calculated as ratio between LA and LDW, heading date (HD), seed number (SN), hundred seed weight (HSW) and yield (YLD).