Transcripts of wheat at a target locus on chromosome 6B associated with increased yield, leaf mass and chlorophyll index under combined drought and heat stress

Drought and heat stress constrain wheat (Triticum aestivum L.) yields globally. To identify putative mechanisms and candidate genes associated with combined drought and heat stress tolerance, we developed bread wheat near-isogenic lines (NILs) targeting a quantitative trait locus (QTL) on chromosome 6B which was previously associated with combined drought and heat stress tolerance in a diverse panel of wheats. Genotyping-by-sequencing was used to identify additional regions that segregated in allelic pairs between the recurrent and the introduced exotic parent, genome-wide. NILs were phenotyped in a gravimetric platform with precision irrigation and exposed to either drought or to combined drought and heat stress from three days after anthesis. An increase in grain weight in NILs carrying the exotic allele at 6B locus was associated with thicker, greener leaves, higher photosynthetic capacity and increased water use index after re-watering. RNA sequencing of developing grains at early and later stages of treatment revealed 75 genes that were differentially expressed between NILs across both treatments and timepoints. Differentially expressed genes coincided with the targeted QTL on chromosome 6B and regions of genetic segregation on chromosomes 1B and 7A. Pathway enrichment analysis showed the involvement of these genes in cell and gene regulation, metabolism of amino acids and transport of carbohydrates. The majority of these genes have not been characterized previously under drought or heat stress and they might serve as candidate genes for improved abiotic stress tolerance.


Introduction
Climate change is a threat to future food security. Prolonged drought periods and heatwaves are predicted to become more common by the end of the century having a major impact on economically important crops such as wheat [1,2]. recorded [3][4][5]. In 2016, yield losses of between 17 and 45% of total production (corresponding to~8 million tonnes) were recorded in France, one of the biggest wheat producers in Europe. Losses were associated with high Autumn temperatures and the compound effects of precipitation changes that are predicted to recur more frequently in the future [6]. In Australia, nine of the ten warmest years on record have occurred since 2005 with rainfall and wheat yields below average, so that the climate was both hot and dry [7,8]. To minimize yield losses and to keep up with future food demand, the development of climate resilient wheat varieties is required. One way to develop more resilient crops is the identification and integration of quantitative trait loci (QTL) and the underlying genes associated with abiotic stress tolerance. QTL have been identified for yield in low-yielding growth environments experiencing drought, heat or combined drought and heat stress (reviewed in [9]). However, phenotyping for grain yield on its own might not be enough to contribute significantly to cultivar improvement given the complexity of the genetic control of stress tolerance (i.e., multigenic, low heritability with high genotype by environment interactions) [10,11]. The dissection of these QTL into their component physiological traits, which then can serve as target traits for breeding in dry and hot climates, is of similar importance. Potential key traits that have been suggested for drought and heat stress tolerance are the regulation of the water use in plants and the adaption of photosynthetic assimilation to improve radiation use efficiency [9,12,13].
Differentially expressed candidate genes under combined drought and heat stress have been identified in tetraploid durum wheat. Combined drought and heat stress triggered the expression of genes encoding trans-membrane proteins, as well as proteins involved in fatty acid β-oxidation [14,15]. In bread wheat, significant genetic marker-trait associations were identified under combined drought and heat stress on group 7 chromosomes and were associated with phytoene synthase 1, integral membrane glycoproteins and a protein conferring rust resistance [16].
We developed lines (NILs) that targeted a QTL on the long arm of chromosome 6B of hexaploid bread wheat (Triticum aestivum L.). The QTL (QYld.aww-6B.1) had previously been identified in three independent studies. In semi-controlled conditions [17], the allele predominantly occurring in Asian accessions contributed to higher total grain weight per plant, single grain weight and harvest index for the heat response under drought, while decreasing screenings (% small grain weight) under drought. However, a different effect was observed in hot conditions in the field with the Asian allele reducing harvest index [18]. In controlled conditions [19], the QTL was associated with single grain weight, leaf chlorophyll content following heat stress and for the heat susceptibility index. To find candidate genes associated with the QTL, we performed a gene expression analysis of developing grains collected during early and late drought as well as with drought and heat stress during the grain-filling period. Further, by studying physiological traits such as water use and photosynthesis, we aimed to identify important tolerance mechanisms associated with drought and heat stress tolerance at this locus, which could potentially be used as target traits in crop breeding.

Plant material
NILs targeting the QYld.aww-6B. 1 were developed from an existing nested association mapping population. The nested-association mapping population consisted of 73 diverse ('exotic') donors which were crossed with two recurrent ('non-exotic') Australian varieties (cvs. Gladius and Scout), back-crossed with the corresponding recurrent parent and selfed over three generations. A total of 772 recombinant inbred lines of 34 families at BC1F4 were genotyped with the SNP marker "BobWhite_c27364_296" (S1 Table). The development of molecular markers and genotyping was performed using Kompetitive Allele Specific Polymerase Chain Reaction (KASP TM ) technology. Subsequently, the selected SNP marker was used to genotype 663 single seed descendants from heterozygous recombinant inbred lines (BC1F5), identifying pairs of NILs carrying the allele from either the non-exotic or exotic parent at the target region. Additional molecular markers developed in-house (S1 Table) revealed a NIL pair, resulting from a cross between the Australian variety 'Gladius' and the Algerian variety 'Taferstat'. Seeds (BC1F6) for each line, derived from a single plant, were used for phenotyping and genotyping. DNA of two seedlings of each line was isolated using the DNeasy 96 Plant Kit (Qiagen, Germany) according to the manufacturer's instructions, of which one sample carrying the nonexotic allele failed (S2 Table, Fig 1). DNA of two additional samples from plants carrying the non-exotic allele was isolated (as before) of which one was similar (sample ID: 22) and one different in its phenotype (sample ID: B23, reduced biomass and plant height) in comparison to the other replicates. A targeted genotyping by sequencing (tGBS) analysis based on data from the 90k SNP Illumina array was carried out by Agriculture Victoria Research (AgriBio, Australia), resulting in 9,424 markers which contained no missing values and could be mapped to a unique position within the bread wheat 'Chinese Spring' reference genome sequence, RefSeq v1.0 [20].

Plant growth conditions
Growth chamber and glasshouse conditions. Plants were vernalized according to [21]. For vernalization, 32 seeds were initially germinated for 24 hours in Jiffypots™ at room temperature (20 o C) in a reach-in chamber, followed by 4 weeks at 4-8 o C with a 2-hour photoperiod (2h/22h day/night) and well-watered conditions. At the end of vernalization, sixteen seedlings of each line were transferred to plastic pots together with their Jiffypots™ and placed on balances on a gravimetric platform (Droughtspotter, Phenospex, Netherlands) which recorded pot weight at 30 min intervals, with precision irrigation. Plastic pots (240 mm high x 165 mm diameter) were filled with 3.5 kg dry weight of a soil mix (1: 1: 1, coco peat-based potting mix: clay loam paddock soil: sand), supplemented with a slow-release fertilizer and covered with a double layer of foam mesh to minimize evaporation. A metal frame was placed around each plant for support. The combinations of sampling dates, treatments and lines were randomized to pots in the glasshouse according to a triple-split-unit design with four biological replicates. Each replicate block was divided into two areas, each of which was, in turn, divided into two subareas. The two sampling dates (i.e. during treatment and at maturity) were randomized to the areas within replicates and the two treatments (i.e. drought or combined drought and heat) were randomized to the two subareas within each area. The lines were randomized to pairs of pots within an area so that NIL contrasting allele pairs occurred together. Pots were manually irrigated and adjusted to their target weight corresponding to 20% soil water content and -0.44 MPa (S1 Fig) during the first two weeks and then automatically irrigated whenever the pot weight dropped 0.5% below target weight. Temperatures in the glasshouse were set to 22/15 o C day/night with a 16 hours day length period and supplemental LED lighting providing approximately 400 μmoles.m 2 .s -1 at the canopy.
Drought and heat treatment. Plants were grown under well-watered conditions (20% soil moisture) except for a 9-day drought period commencing three days after anthesis. All plants were scored every day, with treatment starting three days after the anthesis of each individual. Drought was imposed by lowering the target weight to 12% soil water content, corresponding to -0.72 MPa (S1 Fig). Half of the plants were additionally subjected to heat stress (combined drought and heat) by transferring them to a neighbouring glasshouse with 35/25 o C day/night and similar light settings during the last three days of the drought period. During the heat treatment, plants were irrigated manually to the target weight four times a day. After nine days, plants were moved back to their spot on the gravimetric platform and kept under well-watered conditions until reaching physiological maturity. Transfer of pots and changes in treatments were carried out daily at 8:00 am to maintain a consistent treatment period. Water use, temperature and relative humidity were recorded at 30 min intervals in both gravimetric and neighbouring glasshouses throughout (S2 Fig). thickness), chlorophyll index, photosynthetic capacity (normalized to 25˚C), electron transport capacity and mitochondrial respiration rate (normalized to leaf dry mass).
At physiological maturity, plant height of the primary tiller of the four replicates of each NIL was measured as the distance between the base of the stem and the top of the spike excluding awns. Spike length was measured from the base of the first spikelet to the tip of the last spikelet of the primary tiller. The number of spikes per plant was counted. Samples were dried at 37 o C for 10 days and total above-ground biomass per plant was determined including stems, leaves and spikes of all tillers. Grain traits per primary tiller and for the remainder of the spikes were analysed separately. To differentiate between small (< 2.0 mm) and large grains (> 2.0 mm) a wheat grain sieve (2.0 mm, Graintec, Australia) was used. Grain weight and grain number were determined for grains > 2.0 mm. Single grain weight was measured as the average of grain weight divided by the number of grains. Screenings was defined as the difference in percentage of small grains (< 2.0 mm) compared to total grain weight.

Statistical data analysis of phenotypic traits
Linear mixed models analyses of the phenotyping data were conducted using the R packages ASReml [25] and asremlPlus [26], with the data for the screenings trait being transformed so that they were normally distributed. The models included random terms for the subareas and pairs of pots that were part of the experimental design and allowed for unequal variances between pots for the combinations of treatments and families. Statistical tests were conducted to assess whether these models could be simplified. Also, for each trait, the differences between the combinations of the two NILs of interest here and the treatments were investigated by choosing a model to describe them based on the P-values calculated for both the NIL and treatment main effects and for their interaction. Predicted means (BLUEs) and least significant differences for a 5% level of significance [LSD (5%)] were obtained under the chosen model for the four NIL-treatment combinations. The average daily water consumption was estimated over three consecutive multiple-days periods: during the first six days of drought (3-8 DAA), during the last three days of drought or combined drought and heat stress (9-11 DAA) and during recovery . The daily water use index was calculated by dividing the aboveground biomass measured at maturity by the average daily water consumption during these periods. We used the term 'water use index' rather than 'water use efficiency' because, although we reduced the evaporation during the experiment to a minimum, it is not directly measured using the gravimetric platform. Physiological traits were assessed across the same three time periods as the water consumption and, in addition, for the pre-treatment period (0-2 DAA).

Gene expression analysis by RNA sequencing
Sampling and RNA isolation. The first two spikelets of the primary tiller with extruding anthers of four independent plants per genotype were marked at anthesis. For RNA sequencing, two developing grains from the two spikelets were collected on the fifth day of drought treatment (i.e. eight days after anthesis (8 DAA), spikelet number one) and on the last day of either drought or combined drought and heat treatment (i.e. 11 days after anthesis (11 DAA), spikelet number two), snap frozen using liquid nitrogen and stored at -80 o C. All samples were collected between 10:00 and 11:00 am. A total of four replicates per line, treatment and time point were collected. RNA was isolated using the Spectrum Plant Total RNA kit (Sigma-Aldrich, United States) including alpha-amylase (E-BLAAM 54.0 U/mg, Megazyme, United States) and DNAse (On-Column DNase I Digestion Set, Sigma-Aldrich, United States) treatments. Samples were further purified using the RNA Clean & Concentrator-5 kit (Zymo Research, United States) according to the manufacturer's instructions. Samples were subsequently sent for quality check (Agilent Bioanalyzer 2100, Agilent Technologies, United States) and concentration measurement (Qubit 2.0 fluorometer, Thermo Fisher Scientific, United States) to the ACRF Cancer Genomics Facility (Adelaide, Australia). 15 samples with RNA integrity number (RIN) of at least 8.9 collected at 11 DAA (four replicates per treatment and line, except for one sample where RNA extraction failed) and eight samples collected 8 DAA were sent for sequencing and data analysis to NovoGene (Beijing, China). One additional sample collected 8 DAA (ID: B23, sample name: D8_AA4) was included because it differed in its phenotype from the other replicates (i.e. reduced plant height and biomass).
RNA sequencing. Following quality check of the RNAs, mRNA was enriched using oligo beads, followed by a random fragmentation and a single and then double stranded complementary deoxyribonucleic acid (cDNA) synthesis. Further steps including purification, terminal repair, poly-A-tailing, ligation of sequencing adapters, fragments selection and polymerase chain reaction (PCR) enrichment were performed to obtain the final cDNA libraries. Library concentration was first quantified using a Qubit 2.0 fluorometer and then diluted to 1 ng/μl before checking insert size on an Agilent Bioanalyzer 2100 and quantifying to greater accuracy by quantitative PCR. RNA sequences were obtained through the Illumina sequencing platform (Illumina, United States).
Data analysis. Raw data output from Illumina were transformed to sequence reads by base calling and recorded in a FASTQ file. Initial quality checks included an estimation of error rates for each base along reads and guanine-cytosine content distribution, followed by the removal of reads containing adapters or which were of low quality. After the quality checks, sequences were mapped against the IWGSC reference genome version 2 (https://wheat-urgi. versailles.inra.fr/Seq-Repository/Assemblies) using the hierarchical indexing for spliced alignment of transcripts (HISAT) software [27]. Total gene expression was estimated by counting the reads that mapped to genes or exons. Thereby, read count was not only proportional to the actual gene expression level, but also to the gene length and sequencing depth. In order to compare gene expression levels of different genes, the fragments per kilobase of transcript sequence per millions base pairs sequenced (FPKM) were calculated, taking into account the effects of both sequencing depth and gene length [28]. Gene expression levels were analysed using the HTSeq software [29]. Pearson correlations were calculated to reveal differences in gene expression between samples. Differences in gene expression between treatments, timepoints and alleles including all biological replicates were estimated using the software DESeq [30,31]. The resulting P-values were adjusted using the Benjamini and Hochberg's approach for controlling the false discovery rate. Genes with an adjusted P-value <0.05 were assigned as differentially expressed.
Genes with similar expression patterns were clustered together and represented in a heat plot. Gene ontology (http://www.geneontology.org/) enrichment analysis of differentially expressed genes was carried out using the GOseq software [32] based on the Wallenius noncentral hyper-geometric distribution. To assign the differentially expressed genes to their putative biological function and pathway, a Kyoto Encyclopedia of Genes and Genomes (KEGG, https://www.genome.jp/kegg/kegg2.html) pathway enrichment analysis was conducted using Oryza sativa japonica as a reference genome. We used KOBAS software to test the statistical enrichment of differential expression genes in KEGG pathways. Genes which were identified as differentially expressed between NILs across treatments and timepoints were aligned to the International Wheat Genome Sequencing Consortium (IWGSC) Chinese Spring (CS) RefSeq v1.0 [20] using BLASTN with an e-value cutoff of 10 −40 in order to find their putative locations in the wheat genome.

Genotyping by sequencing of near-isogenic lines
Near-isogenic lines (NILs) were developed targeting a QTL interval of~7 Mbp on the long arm of chromosome 6B. tGBS data of the five NIL pairs at BC1F6 (Fig 1, S2 Table) indicated 93.7% similarity between NILs and a~17 Mbp interval, including the target region, segregating between NIL pairs on chromosome 6B (RefSeq v1.0). Further segregating regions were detected on chromosomes 1B and 7A. NILs carrying the non-exotic allele (Gladius) were 95.7% similar to each other and segregated on chromosomes 1B and 4A. NILs carrying the exotic allele from Taferstat were heterozygous for regions on chromosomes 1B, 3B and 7A with a genotypic similarity of 96.1%. As expected, NILs segregating for the target region were genotypically more different than replicates of the same NIL. Genetic differences between replicates (3.7-4.5%) could be due to residual heterozygosity of the plants which were selected and separately propagated from BC1F5.

Phenotypic data
Yield-related traits. Flowering time occurred over a time frame of 13 days with the first plant flowering 80 days after sowing. NILs carrying the exotic allele at the target region on chromosome 6B flowered, on average, two days earlier (p = 0.011) in comparison to NILs carrying the non-exotic allele (Fig 2). Similarly, plant height and biomass were promoted by the exotic allele under drought and combined drought and heat (p � 0.022) with average increases of 7.5 cm in length and 0.6 g in weight. Biomass, harvest index and grain weight were reduced under combined drought and heat stress in comparison to drought in both NILs (p � 0.05). Single grain weight was reduced in the primary tiller under combined drought and heat stress compared to drought, but similar under both stresses when measured per plant. In contrast, grain number was reduced under drought and heat stress per plant but not per primary tiller. Spike length, spike number, plant height and screenings were similar in both treatments.
Significant differences between alleles at the 6B QTL were found in six of the eleven yield component traits (Fig 2). Grain weight of the primary tiller and whole plant was increased by the allele derived from Taferstat in both treatments (p � 0.011) and, in addition, screenings per plant were reduced in NILs carrying this allele (p � 0.03) in both treatments. Grain number per primary tiller was increased in NILs carrying the exotic allele in both treatments, a potential result of the longer spikes in NILs carrying the exotic allele. Grain number per plant showed the same trend but differences were not significant. Similarly, single grain weight was significantly different under drought between both NILs per primary tiller, but not per plant.
Water use. The average daily water use did not differ significantly between the treatments during the drought treatment period and was 84.7 ml per day (Fig 3). For the drought treatment the average daily water use did not increase as the treatment prolonged (i.e., 9-11 DAA), being 82.1 ml per day, but it increased by 46.8 ml per day during the recovery period. On the other hand, the average daily water use increased markedly during the combined drought and heat treatment, when it was 117.1 ml per day; it then remained little changed during the recovery period at 104.64 ml per day. On average, plants under drought stress decreased their water use 52 DAA, i.e. 41 days after treatment, whereas plants under combined drought and heat stress reduced their water use 36 days after treatment. The daily water use was, however, similar in both NILs with NILs carrying the exotic allele tending to have a slighter higher water use throughout the experiment. Daily water use index ranged between 0.57 and 0.62 g biomass ml -1 day -1 during the drought treatment (3-11 days, Fig 3) but was reduced in plants subject to combined drought and heat treatment and during recovery (0.32-0.50 and 0.35-0.43 g biomass ml -1 day -1 , respectively). During recovery (12-42 days), daily water use index was increased in NILs carrying the exotic allele (p = 0.002) compared with NIL carrying the non-exotic allele. Similar to the water consumption, no significant differences between NILs prior to and during the treatment were observed, although NILs carrying the exotic allele tended to have a higher water use index.
Photosynthesis-related traits. Photosynthesis-related traits were relatively stable throughout the whole measurement period in both NILs and treatments, except for an increase (p � 0.005) in leaf dry mass, photosynthetic capacity (Vcmax) and electron transport capacity (J) and a decline (p = 0.013) in the respiration rate under combined drought and heat compared to drought alone at 9-11 DAA (Fig 4). NILs carrying the exotic allele had increased leaf dry mass per area (i.e., thicker leaves) and chlorophyll index in comparison to NILs carrying the non-exotic allele at all times and in both treatments (p � 0.05). Photosynthetic and electron transport capacity showed a similar trend but were only significantly higher in lines carrying the exotic allele at 9-11 DAA (p = 0.038) and 0-2 DAA (p = 0.045), respectively. Mitochondrial respiration rate was similar over the measurement period with no significant differences between NILs.

Gene expression analysis of developing grains
RNA sequencing data were of good quality with 94.7-99.0% of clean data in each sample. The percentage of mappable reads for all samples was above 70% (S5 Table). Most reads could be mapped to exons (57.6-76.2%), followed by intergenic regions (22.5-41.5%). The smallest proportion was mapped to intron regions (0.8-2.0%). Some of the reads were mapped to more than one exon (17.5-26.3%) potentially due to repetitive DNA within a chromosome or a gene copy on one of the other two homeologous chromosomes. Read densities were similar for the positive and negative strands of each chromosome.
Overall, all samples showed similar gene expression levels with the majority of genes poorly expressed (FPKM <1 for �55.9% of the total number of genes) (S6 Table). Genes with a medium to high expression, i.e., FPKM between 1-3, 3-15, 15-60, accounted for 11.3-13.2%, 19.3-21.7% and 6.1-8.0% of the genes, respectively. Genes with a very high expression (FPKM > 60) accounted for 1.8-2.3% of the genes. Biological replicates were 90-99% similar in their gene expression (S7 Table), except for sample "D8_AA4". The sample "D8_AA4" was included as additional sample because of its phenotypic difference (i.e., reduced plant height and biomass) in comparison to the other replicates. Pearson correlations (R 2 ) between "D8_AA4" and the other replicates collected at 8 DAA carrying the non-exotic alleles ranged between 0.79 and 0.81, at the limit of the suggested threshold of 0.80 for reliable replicates. tGBS data did not reveal any difference in the genome of plants with reduced biomass and plant height in comparison to plants showing the common phenotype. The markers used for tGBS might not have covered the region associated with this phenotype, or an epigenetic modification might have caused this difference.
A total of 42,393 genes were similarly expressed in both NILs and treatments. Genes differentially expressed between the two treatments at 11 DAA (Fig 5 - Table) accounted for 5,507 and 28,371 of the genes, respectively. Of these, 2,278 and 10,579 were differentially expressed between treatment or timepoint in both NILs. Genes differing between timepoints were principally involved in plant and grain development such as cell number regulation, DNA repair mechanisms and transport and hydrolysis of sugars (S9 and S10 Tables).
NILs carrying the exotic allele differed by 3-11% in their number of expressed genes to NILs carrying the non-exotic allele (S7 Table) with a total of 2,082, 358 and 164 differentially expressed genes under drought at 8 DAA, under drought at 11 DAA and under drought and heat at 11 DAA, respectively. Differentially expressed genes at 8 DAA were mainly located on the long arms of chromosomes 1B (64 genes), 4B (1,037), 6B (37) and 7A (18) and on the short arm of chromosome 4B (763 genes) (S8 Table). Most of these genes encode proteins which are cellular components of the cytoplasm and the endomembrane system (S9 Table). Upregulated genes on chromosome 4B in NILs carrying the exotic allele were mostly involved in alanine, aspartate and glutamate metabolism (Fig 6A, S10 Table). Genes differently expressed between NILs at 11 DAA under drought were similar to those detected at 8 DAA and mostly, but not significantly, associated with cell components (Fig 6B). Differentially expressed genes under combined drought and heat at 11 DAA were dominantly located on the long arms of chromosomes 1B (51 genes), 6B (48 genes) and 7A (16 genes); only one of the genes was located on chromosome 4B (S8 Table). Most of these genes encode binding proteins which are involved in pathways such as glutathione metabolism, plant-pathogen interactions and RNA transport (not significant) (Fig 6C).
A total of 67 high confidence and 8 low confidence genes were differentially expressed between NILs carrying the opposite allele across treatments and timepoints, with 27 Table 1). Thirty five genes were upregulated and 40 downregulated in NILs carrying the exotic allele. Both, upregulated and downregulated genes, are involved in cell and gene regulation,  protein binding, disease resistance, carbohydrate transport and metabolic pathways. In addition, two of the upregulated genes are associated with the arginine and proline metabolism (TraesCS6B02G456400) as well as with the development of anatomical structures (TraesCS1B02G269500), whereas one of the downregulated genes (TraesCS1B02G288900) is associated with the Golgi vesicle transport and three encode proteins located in the chloroplast.

Discussion
The identification of QTL and the underlying genes associated with grain yield and yield stability following abiotic stress can be valuable for the development of new, high yielding varieties. Identified QTL and their associated molecular markers can be used for marker-assisted selection, a method which enables the selection of genotypes in large populations with reduced need for costly and time-consuming phenotyping in the field [33]. Knowing the genes and their function, on the other hand, can provide information on key mechanisms associated with stress tolerance and can be used for the direct modification of current cultivars by methods such as CRISPR [34]. Of similar importance is the physiological dissection of these QTL into their components, which can then serve as target traits for breeding in dry and hot climates. Using NILs, we studied a target QTL on the long arm of chromosome 6B and its effect on yield components, water use and photosynthesis-related traits. The QTL was previously identified in three independent studies [17-19], contributing to grain weight, single grain weight, harvest index and leaf chlorophyll content, but with opposite allelic effects under semicontrolled compared with field conditions. NILs carrying the allele from the exotic parent at chromosome 6B had an increase in grain weight, grain number and single grain weight as well as a decrease in screenings. Results were therefore consistent with results observed in semi-controlled conditions [17]. Differences in grain number, single grain weight and screenings were significant per primary tiller but not per plant. Anthesis was recorded only for the primary tiller. Other tillers may have differed in their development compared to the primary tiller, varying the timepoint and impact of the stress, buffering thus the differences between NILs. Differences in the development of tillers could also explain the observed variation in yield components within the same plant. Single grain weight was more severely affected by the combined drought and heat stress at the primary tiller compared to drought alone, whereas grain number was more affected in the whole plant. Reductions in single grain weight are commonly associated with post-anthesis stress, whereas a decline in grain number is often observed following pre-anthesis stress [19,35]. In our experiment, the stress treatments began three days after anthesis of the primary tiller, a post-anthesis stress. Final grain weights both per tiller and per plant were increased in plants carrying the exotic allele in both treatments. This indicated that the beneficial effect of the QTL was not dependent solely on post-anthesis growth and the development of tillers, but more likely related to either whole life cycle plant physiology or a change in grain filling or senescence, post-treatments.
Photosynthesis-related traits, including leaf dry mass, chlorophyll index, photosynthetic and electron transport capacity as well as mitochondrial respiration rate were fairly stable throughout the whole measurement period in both NILs, except for a peak during the combined drought and heat treatment. To avoid photoinhibition and allow acclimation, plants can optimise light absorption through leaf and chlorophyll movement and anthocyanin accumulation, as well as the energy balance of their photosystem, through modifications of CO 2 fixation, photo-and mitochondrial respiration and cyclic electron flow [36,37]. In our case, photosynthetic and electron transport capacity were reversibly modified in response to the occurrence and absence of the heat stress. This suggests that all NILs were able to acclimate to higher temperatures.
NILs carrying the exotic allele maintained their photosynthetic and electron transport capacity above NILs carrying the non-exotic allele with significant differences observed at 0-2 and 9-11 DAA. NILs carrying the exotic allele also had an increased leaf dry mass (i.e. thicker leaves) compared to NILs carrying the non-exotic allele before, during and after treatments, and increased chlorophyll index. A higher photosynthetic capacity would suggest a higher contribution to grain yield [38,39], which was the case in our study.
Water use and photosynthesis-related traits have previously been suggested to be important traits to increase wheat yield potential under drought and heat stress [9,13]. Particularly, an increased water use efficiency has often been hypothesized in literature to be associated with a higher stress tolerance [40][41][42][43], while others argue that the effective use of water (i.e. maximal soil moisture capture for transpiration) and not water use efficiency is important for crop improvement [12,44].
In our case, water use was similar, although tendentially higher in NILs carrying the exotic allele, throughout the experiment. Both NILs adapted their water use equally to the changing water availabilities. Water use was increased in both NILs during the combined drought and heat treatment compared to drought alone, potentially due to an increase in transpiration for cooling. Pots were constantly re-irrigated to a certain target weight throughout the experiment instead of withholding irrigation completely, enabling a constant delivery of water. A higher water use and transpiration rate, therefore, had no negative effect in this experiment but might be a problem in a situation in which water availability was more restricted. Water use was also increased during the recovery phase following both treatments, likely due to the higher water availability, and resulting in an overall lower water use index during this phase. Although there were no differences between NILs for water use, NILs carrying the exotic allele had an overall higher water use index with a significant increase during recovery. Increased water use index in these NILs was due to increased biomass and this post-treatment difference suggested an accelerated stress-induced senescence in Gladius allele NILs compared with the exotic. The QTL at this locus was first identified for grain-filling heat stress tolerance for the 'stay-green' trait, i.e., delayed chlorophyll loss of the flag leaf [19]. Our results are consistent with these and underscore the importance of this trait. We additionally found increased leaf dry mass and leaf chlorophyll index throughout, including pre-stress, and in drought as well as combined drought and heat in the exotic NIL. The identification of the underlying physiological traits contributed by the beneficial alleles is particularly relevant to refine trait phenotyping priorities for growing environments where drought and heat stress occur simultaneously.
Segregation regions between NILs were observed not only at our target region on chromosomes 6B, but also on 1B and 7A. The region on chromosome 1B co-located on the Ref Seq v1.0 [20] with a previously identified QTL for anther extrusion in wheat [45]. The segregation region on chromosome 7A co-located with QTL identified for thousand kernel weight and spikelet number per spike [46][47][48]. One of our differentially expressed genes (TraesC-S7A02G479800, Table 1), encoding a putative Eukaryotic translation initiation factor 5B, was located within the exact same region on the Ref Seq v1.0 as the previously identified QTL and was upregulated in NILs carrying the exotic allele in all treatments and timepoints. A second gene (TraesCS7A02G484800), encoding a G-type lectin S-receptor-like serine/threonine-protein kinase B120 was 1.3 Mb distant from the QTL and downregulated in NILs carrying the exotic allele. The increase in grain number and single grain weight in NILs carrying the exotic allele might have therefore been the result of the interaction of both QTL -6B and 7A -or of one alone.
Most genes differentially expressed between NILs were observed at an early drought stage (i.e. 8 DAA), whereas this number decreased with an increase of stress duration (i.e. at 11 DAA) and intensity as shown under combined drought and heat treatment. Genes differentially expressed between NILs carrying the opposite allele were located on chromosomes 1B, 4B, 6B, 6D and 7A, i.e. mostly in regions of genotypic differences observed between NILs (Fig  1). Genotyping by sequencing data, however, did not suggest regions of segregation on chromosomes 4B and 6D. Genotypic differences might have had a trans-regulatory effect. Several of the differentially expressed genes between NILs, in particular those under drought and heat, were associated with regulation of gene expression and RNA processing supporting the hypothesis of a trans-acting control.
The majority of differentially expressed genes under drought were located on chromosome 4B with upregulated genes in NILs carrying the exotic allele being involved in alanine, aspartate and glutamate metabolism. All three amino acids have been observed to increase in developing grains of drought-resistant wheat plants when subjected to drought [49] and form part of the photorespiratory cycle [50,51]. Differentially expressed genes under combined drought and heat stress were dominantly located on chromosomes 1B, 6B and 7A and were associated with the metabolism of glutathione, a component of the antioxidative system in plants, which is synthetised from glycine, a by-product of photorespiration [51]. Photorespiration, the fixation of O 2 by ribulose-1,5-bisphosphate, is usually considered a negative trait in crop yield, degrading sugars which have been produced in energy-consuming reactions during photosynthesis, releasing CO 2 and producing reactive oxygen species which can harm the cell. The ascorbate/ glutathione cycle is an important component in the regeneration of antioxidant scavengers and glutathione reductase and peroxidase have been previously observed to be specifically induced by drought and heat stress in tobacco [52].
A total of 36 differentially expressed genes across treatments and timepoints could be mapped to our target region on the long arm of chromosome 6B ( Table 1). Several of these genes were associated with a wide range of functions such as carbohydrate transport, gene regulation, protein binding, disease resistance and various enzymatic activities, whereas others were of unknown function. To our knowledge, none of these genes has been previously characterized under drought or heat stress in wheat, except for one (TraesCS6B02G456400) whose predicted protein functions in hydroxylation of the well-studied amino acid proline. Proline has been shown to accumulate in plants in response to drought, heat and combined drought and heat stress [53,54]. Its accumulation has been associated with tolerance mechanisms such as reactive-oxygen species scavenging, osmotic adjustment, signalling and the stabilisation of proteins [55,56] as well as with an increase in grain yields [54,57].
Seed weight was increased in NILs carrying the exotic allele by 24-32% under drought and by 6-23% under combined drought and heat stress. An increase of only 6% could mean, in theory, a contribution of 1 million tonnes per annum in Australia alone. Marker-assisted selection and introduction of such favourable alleles to breed higher yielding cultivars for regions affected by drought and heat stress is a promising approach. First steps have been achieved by introgressing the favourable allele into an Australian background in a nested-association mapping population [17]. However, with the selection for a single allele, such significant yield increases are less likely, as traits such as yield are under multigenic control and allelic effects are highly dependent on the environment. Nonetheless, this work is significant progress towards the identification of genes for drought and heat tolerance which can enable the discovery of new regulatory pathways in wheat and, possibly, the development of novel alleles via genome editing technologies.
We could not draw a conclusion about the opposite allele effects observed in the field by Garcia et al.
[18] based on two trials in South Australia. However, performances in field trials are influenced by season to season environmental variations, additional and interacting biotic and abiotic stresses and, therefore, a number of field trials are required to confirm results. The higher harvest index promoted by the Australian allele in the field might have resulted from one of these factors rather than from actual drought or heat stress. On the other hand, a positive effect of the non-Australian allele was always found in pot experiments. Pot systems often differ in their water relations, the structure and temperature of the soil as well as the available root space and competition between plants compared to field conditions, all of which strongly influence root architecture and physiology as well as the interactions in the rhizosphere [58]. These factors could play a role in whether the exotic allele has a positive or negative effect on yield.

Conclusions
Allelic effects of QYld.aww-6B.1 on grain weight, single grain weight, grain number and screenings under drought and combined drought and heat stress were consistent with results from the genome-wide association study [17]. An increase in yield was also associated with thicker leaves, higher leaf chlorophyll index, a higher photosynthetic capacity and a higher water use index. Using gene expression analysis, we could narrow down our target region on chromosome 6B to 36 potential candidate genes with a further 39 genes of interest differentially expressed across treatments and timepoints in the NIL. Genes from four other regions within the genome (i.e., on chromosomes 1B, 4B, 6D and 7A) were differentially expressed between both NILs. Two of the four regions could be mapped to previously detected QTL, of which one has been identified for thousand kernel weight. Candidate differentially expressed genes were usually associated with genetic segregation, illustrating the value of the nested-association mapping population for the rapid incorporation and validation of the beneficial exotic alleles. The majority of these genes have not previously been associated with drought or heat stress tolerance in wheat and might serve as candidate genes for crop improvement in dry and hot climates. Further analysis regarding their involvement in the observed changes in physiology and yield components is required.
Supporting information S1 Fig. Soil water potential-soil water content curve of the drought soil mix. Eight pots of the same size and filled with the same substrate mix as used in the experiment were first watered and then drained until reaching~5% soil water content. The water content and water potential of the soil were measured daily by coring a of sample of 15 mm diameter and 200 mm length. The water content in soil was determined by calculating the weight difference between fresh soil sample and the oven-dried soil sample (at 65 o C for 72 hours), divided by the oven-dried sample. The water potential of the fresh soil sample was measured using a water potential meter (WP4C, Meter Group, United States) in continuous mode until the value maintained stable. Dots, raw values of the eight pots. Blue line, logarithmic trendline.  Table. KASP marker assisted selection for the development of near-isogenic lines. Grey indicates the allele derived from the non-exotic parent (Gladius) in the target regions of chromosome 6B. Green indicates the allele derived from the exotic parent (Taferstat). (XLSX) S2 Table. Targeted genotyping by sequencing data of near-isogenic lines for the target QTL on chromosome 6B. The target region (T) and segregating regions (S) between NILs of the same pair as well as between replicates are marked in yellow. Exotic, allele from the Algerian parent (Taferstat); non-exotic, allele from the Australian parent (Gladius).  Table. Quality assessment of the RNA sequencing data. Raw reads refer to the total number of reads, clean reads to the number of filtered reads. GC content, guanine-cytosine content; Q20, percentage of bases whose correct base recognition rates are greater than 99% in total bases. Q30, percentage of bases whose correct base recognition rates are greater than 99.9% in total bases. (XLSX) S5 Table. Mapping of RNA sequences to the wheat genome. Total reads, total number of filtered reads (clean data); total mapped, total number of reads that could be mapped to the reference genome CS RefSeq v1; uniquely mapped, number of reads that were uniquely mapped to the reference genome; multiple mapped, number of reads that were mapped to multiple sites in the reference genome; reads map to '+', number of reads that map to the positive strand; reads map to '-', number of reads that map to the negative strand; exons, percentage of reads mapped to exons; introns, percentage of reads mapped to introns; intergenic, percentage of reads mapped to intergenic regions; non-splice reads, reads that were mapped entirely to a single exon; splice reads, reads that were mapped to two exons.