Four parental lines were derived through three cycles of selection for F_v/F_m during heat stress from originally 1274 spring wheat cultivars of diverse geographical origin (Sharma et al., 2012). Out of these four cultivars, three cultivars named 810 (IPK-2845, origin Afghanistan), 1039 (IPK-8183, origin Pakistan), and 1313 (IPK-28703, origin Pakistan) were selected as heat tolerant parents owing to their consistent high F_v/F_m during all three heat stress screenings (Sharma et al., 2012) as well as in the subsequent physiological validation for maintaining higher overall photosynthesis (Sharma et al., 2015). Similarly, the fourth cultivar named 1110 (IPK-9705, Kloka WM1353, Germany) was selected as the most heat susceptible cultivar owing to its consistent low F_v/F_m and net photosynthesis under heat stress in our previous studies (Sharma et al., 2012, 2015). Bi-parental crossings were made with each of these three heat tolerant cultivars (810, 1039, and 1313) as male parent and the heat susceptible cultivar (1110) as the common female parent. Seeds of a single F₁ plant from each cross were then selfed to produce F₂ progeny. A total of 140 F₂ individuals in each of the three bi-parental cross combinations 1110 × 810, 1110 × 1039, and 1110 × 1313 constituted the three mapping populations.

A total of 420 F₂ plants across three populations together with 10 plants each of the four parental cultivars (i.e., 460 plants in total) were completely randomized and divided into five blocks with 92 plants each. To be able to heat stress all the plants around anthesis seeds were sown 1 week apart between the blocks. Seeds were sown individually in plastic pots (11 cm diameter; 0.59 L) TEKU^®, VCD series, Pöppelmann, GmbH & Co. KG, Kunststoffwerk, Werkzeugbau, Germany) with peat substrate (Pindstrup 2, Pindstrup Mosebrug A/S, Denmark) under greenhouse conditions at 15 ± 3°C day and 12°C night temperatures, 50–70% relative humidity and ambient CO₂ concentration with regular irrigation and standard nutrient solutions consisting of 185 ppm (w/v) N, 27 ppm (w/v) P, 171 ppm (w/v) K, 20 ppm (w/v) Mg, and full micro nutrients, electric conductivity 2.0 mS m^-1 and pH 5.8.

Phenotypic evaluation for heat tolerance was performed by following the previously developed protocol, which was standardized while screening the parental cultivars (Sharma et al., 2012). Plants in each block were transferred to a growth chamber (MB-teknik, Brøndby, Denmark) with 16 h of photoperiod at 250–300 μmol m^-2 s^-1 photosynthetic photon flux density (PPFD), 400 ppm CO₂ and 60–70% RH. The intention behind the relatively low PPFD was to avoid the confounding effect of photo-inhibition on the heat stress treatment. Plants were acclimated at 20°C for 3 days before the temperature was increased 10°C per hour to 40°C (day/night) and kept constant for three consecutive days as heat stress treatment. The phenological stage of each plant was registered according to the BBCH-scale for cereals (Lancashire et al., 1991) and the heat stress treatment was given when the majority of the plants reached the growth stage in the range between heading and beginning of anthesis (BBCH scale 51–61) on day 0 of the treatment. The plants were watered frequently to avoid drought stress. During the acclimatization period, two penultimate leaves on each plant were fixed as the sampling leaves for measurements. Two measurements per plant were taken on day 0 (before transferring to growth chambers as control), day 1, day 2, and day 3 of heat stress. For each measurement, a leaf segment (5–8 cm long) was detached from each sampling leaf followed by clipping with a dark adaptation leaf clip (Hansatech Instrument, King’s Lynn, England). The clipped leaf samples were arranged on a plastic tray lined with moist tissue paper and enclosed in a plastic bag in order to avoid evaporative water loss during the 30 min of dark adaptation at room temperature (∼22°C). F_v/F_m was measured on the adaxial leaf surface with a saturating flash of 3000 μmol m^-2 s^-1 for the duration of 1 s using a Plant Efficiency Analyzer, Handy PEA (Hansatech Instrument, King’s Lynn, England). The F_v/F_m data were power transformed [(Y ^λ - 1)/λ] in order to obtain approximate uniform normal distribution of residuals across all the values (Box and Cox, 1964). The λ value was 11, which was initially determined from residual plots created using SAS ver. 9.2 (SAS Institute Inc., Cary, NC, United States). The two measurements from each day were averaged after the transformation and the data from the five blocks were pooled. For each of the parental lines (n = 10) and mapping populations (n = 140) the effect of heat stress over the duration of 3 days was calculated with paired t-test at p < 0.05. Since there were no replicates, the interaction between the block and genotype could not be tested. However, when the same growth chambers and experimental setups were used in our previous three-tiered phenotyping experiments to identify the parental lines no significant block effect was found in any of the three experiments and thus, the data could be pooled (Sharma et al., 2012).

Leaf samples from 1-month-old seedlings of each of the 460 plants used for phenotypic evaluation were freeze dried for 72 h using a freeze dryer (Christ Alpha 1-4, GmbH, Germany). The genomic DNA was isolated using the method described by Orabi et al. (2014) with slight modifications where, DNA was fished out using an inoculation loop followed by washing with 75% ethanol. Then pellets were air dried and suspended in 100 μl of TE buffer pH 8 (10 mM Tris.HCl and 1 mM EDTA) and allowed to dissolve at 4°C for a week. Concentration of the extracted DNA was measured with a nanodrop 2000 UV-Vis spectrophotometer (Thermo Fisher Scientific, Wilmington, United States) and the purity and integrity of the DNA was tested using 1% agarose gel electrophoresis.

The DNA samples (50 ng μl^-1 in TE buffer pH 8) from all the 460 plants were outsourced to the Diversity Arrays Technology Pty. Ltd. (DArT P/L, Yarralumla, ACT, Australia)¹. The samples were genotyped using wheat DArTseq, a genotyping by sequencing platform that provides high-density single nucleotide polymorphism (SNP) as well as presence/absence variation (PAVs also called SilicoDArTs) markers. Additionally, this platform also provided the 69-nucleotide DNA sequence of the representative fragments derived from the genome complexity reduction process achieved through the use of a combination of restriction enzymes (Kilian et al., 2012). This high throughput genotyping on our mapping populations provided a total of 34955 DArTseq markers out of which, 30178 were PAVs and 4777 were SNPs.

SNPs scores were re-coded to reveal parental origin (“A” referring to the common female parent 1110) for each of the three populations for co-dominant (A, B, H) and dominant (either A/C or B/D) by a custom Java-based personal program (developed by SA). From the original 4777 SNPs, monomorphic as well as SNPs showing highly distorted segregation were eliminated. The linkage groups for the three mapping populations, 1110 × 810 (with 1042 SNPs), 1110 × 1039 (with 790 SNPs), and 1110 × 1313 (with 877 SNPs) were constructed using JoinMap^® 3 (Van-Ooijen and Voorrips, 2001). In the absence of information on the chromosomal location of these SNPs during this first round of analysis (in 2013), maps were constructed from linkage groups containing 2–100 markers at increasing log of odds (LOD) scores. A total of 39, 34, and 54 linkage groups were obtained in the 1110 × 810, 1110 × 1039, and 1110 × 1313 populations, respectively. These linkage maps were used for the QTL mapping of F_v/F_m in order to check the presence of any significant QTLs using a Java-based software, SuperQTL (developed by SA). This software was based on multi-QTL interval mapping (Jansen and Stam, 1994) and was also used in a previous study (Kuzina et al., 2011). The LOD threshold for declaring a QTL significant was estimated using 1000 datasets with permuted trait values at p < 0.05 and was found to be 4.1 in the 1110 × 810, 4.4 in the 1110 × 1039, and 4.3 in the 1110 × 1313 populations (Supplementary Table S1).

The sequence information of the important SNPs around the QTL region (from the first round) was used to retrieve the sequence of corresponding wheat contigs from the CerealsDB 2.0 database (Wilkinson et al., 2012). These sequences were subsequently mapped in silico to the bread wheat chromosome-based survey sequences in URGI², to local databases of the wheat A (Ling et al., 2013) and D (Jia et al., 2013) draft genome sequences established in CLC main workbench³ (last date of access to the various databases: July 16, 2013). This bioinformatics approach provided us a hint that the identified QTL regions might be in the short arms of chromosome group 3, and group 1.

To confirm this finding, primer information for a total of 50 simple sequence repeats (SSRs) from the chromosome group 3 and 21 SSRs from the chromosome group 1 were searched from the database, Graingenes 2.0⁴ (Roder et al., 1998; Somers et al., 2004). Polymerase chain reaction (PCR) for SSR analysis was done with the three-primer approach for fluorescent labeling of the PCR products (Schuelke, 2000). PCR amplification was carried out as described by Orabi et al. (2014) and the SSR fragments were analyzed by capillary electrophoresis using a 3130xl Genetic Analyzer (AB/Hitachi, Thermo Fisher Scientific Inc., MA, United States). Out of a total of 71 SSRs tested, 12 SSRs from chromosome 3 and 15 SSRs from chromosome 1, were found polymorphic and hence, used for genotyping the three mapping populations.

In the meantime, a consensuses map on the genetic location of DArTseq markers (Li et al., 2015) and updates on the URGI wheat genome assembly (see text footnote 2) and Ensembl Plants wheat⁵ became available. The 69-nucleotide DNA sequence of the 34,955 DArTseq markers was used for BLAST search against a local depository of wheat sequences based on the Triticum_aestivum.IWGSC1+popseq.30 (downloaded from URGI and Ensembl Plants wheat, updated until August 2015). Only those markers that had a defined chromosomal position, 100% sequence identity and alignment, no gaps and mismatches, lowest e-value and highest bit scores were filtered in order to assign the chromosomal location. Further, the linkage group and position of these markers (when available) in the wheat DArTseq map (Li et al., 2015) was also assigned. This information was then used in differentiating linkage groups to assign chromosomes during map construction. The markers having missing scores of more than 11% and/or showing distorted segregation (p < 0.001) were removed before any further analysis. Since all the PAVs, some SNPs and a few SSRs were dominant markers, a linkage map was first calculated based on only the co-dominant SNPs using the regression mapping algorithm with Kosambi mapping function in JoinMap^® 4.1 (Van-Ooijen, 2006). The grouping information obtained from this co-dominant SNPs map in each population was then referred while grouping the rest of the markers during final mapping.

The final mapping of the linkage groups as well as the QTL analysis was done by QTL IciMapping software V4.0 (Meng et al., 2015). With the BIN function the redundant markers were removed by the missing rate of 11% that left the markers with fewest/no missing scores as representative marker of the bin. Before binning, there were a total of 4154, 4301, and 3819 polymorphic markers genotyped on 140 F₂ individuals each of the 1110 × 810, 1110 × 1039, and 1110 × 1313 populations, respectively. The output from the BINNING was then used to create linkage groups by LOD threshold of 20 using the MAP functionality. Ordering of markers within the group was done by nnTwoOpt algorithm, followed by rippling step with SAD (sum of adjacent distances) with a window size of 10. Rippling fine-tunes the order of the markers within the linkage groups by permutation of the number of markers specified in the window size at a time.

The resulting final linkage maps of each population were then used to scan for QTLs following the bi-parental population BIP functionality in the same QTL IciMapping software v4.0, which is based on the inclusive composite interval mapping (ICIM). The QTL analysis was done for the trait F_v/F_m separately for day 0 (control), day 1, day 2, and day 3 of heat stress. The set parameters were “deletion” for missing phenotypes, mapping steps at 1 cM and probability in stepwise regression (PIN 0.001) with ICIM-ADD mapping method and LOD threshold at 3.

After confirming that the chromosomal location of the QTLs was 3B and 1D, the physical position of the markers (when available) was searched from Ensembl Plants wheat and URGI databases. The QTL region between the flanking markers was estimated based on the one LOD drop off interval both on 3B and 1D chromosomes. The marker sequences were blasted against the Ensembl Plants wheat (see text footnote 5) to find the Traes numbers of the genes present around the QTL regions. The 3B region was updated through wheat_Jbrowse⁶ (URGI wheat_Jbrowse released September 14, 2015) hosted at URGI. These Traes numbers were searched in the UniProt in TrEMBL⁷ (release 2015_09 of September 16, 2015 of UniProtKB/TrEMBL) to obtain more information including protein domain, family, molecular and biological functions of the potential candidate genes (last assessed of various databases December 2015). However, only those genes with known function and/or related to stress and photosynthesis were counted as potential candidate genes.

The figures on phenotypic evaluation and linkage maps were drawn by Sigmaplot v11 (Systat Software, San Jose, CA, United States) and MapChart v2.3 (Voorrips, 2002), respectively.

The average phenological growth stage as per BBCH-scale for cereals (Lancashire et al., 1991) of the plants at the time of phenotypic evaluation (day 0) was considered (Figure 1). In the population 1110 × 810, the two parents respectively had a growth stage of 53 and 59, while the F₂ population maintained stage 56, showing a non-significant variation between the two parents as well as between the F₂ plants within the population. In the population 1110 × 1039, the heat tolerant parent 1039 and the F₂ population were on growth stage of 56 and 53, respectively with a non-significant variation. In the third population 1110 × 1313, the growth stage was 53 in both the parents and 54 in the F₂ population showing a non-significant variation (Figure 1). Overall, the growth stage ranged between 53 and 59 in the parental lines and 53 and 57 in the F₂ populations (Figure 1), where the BBCH scale between 51 and 61 corresponds to the stage between inflorescence emergence (heading) and beginning of anthesis. All these non-significant variations demonstrate that the four parental cultivars were fairly uniform in their developmental rates and therefore, there was no clear segregation for phenological growth stage in their respective F₂ population. However, comparison of the three populations shows that the 1110 × 810 population in average was slightly earlier than the 1110 × 1039 to reach anthesis (Figure 1).

All the plants maintained a fairly constant F_v/F_m around 0.84 at day 0 (20°C, control condition) while the 3 days heat stress treatment at 40°C led to a decrease in F_v/F_m (Figure 2A). However, the heat stress induced decrease in F_v/F_m was modest in the three tolerant parents, becoming significant only on day 3 in 810 (with 1.2% decrease) and 1313 (with 2.5% decrease) while being non-significant in 1039 (Figure 2A). On the other hand, there was a successive reduction in F_v/F_m in the heat susceptible parent 1110 as the duration of heat stress progressed, with overall reduction of 6.3% on day 3 as compared to its day 0 value (Figure 2A). The three mapping populations showed an intermediate response, indicating a clear segregation for the trait, with the overall reduction ranging between 4 and 6% on day 3 as compared to their respective day 0 values (Figure 2B).

Out of a total of 34955 DArTseq markers (including 4777 SNPs and 30178 PAVs) plus 27 polymorphic SSRs resulting from the genotyping of three mapping population, the final linkage map constructed using QTL IciMapping software V4.0 consisted of 1651 markers in 1110 × 810, 1752 markers in 1110 × 1039, and 1672 markers in 1110 × 1313 populations (Table 1 and Supplementary Table S2). In all the populations, there were a total of 27 linkage groups covering all the 21 wheat chromosomes, confirmed by the information obtained from the genetic location of these markers in the wheat consensus map (Li et al., 2015) as well as from BLAST against the Triticum_aestivum.IWGSC1+popseq.30 genome assembly. However, in population 1110 × 1313 only two markers were mapped on chromosome 3D (Table 1). The chromosomes 3B, 6A, and 7A each formed two linkage groups. A small fraction (1–3%) of the total markers located to three short linkage groups with unknown chromosomal location termed not assigned (Table 1). The total length of genome across the three linkage maps was around 5000 cM with a marker distribution of 36–42% to the A genome, 47–51% to the B genome and 9–12% to the D genome.

QTL analysis using ICIM resulted in the identification of one significant QTL for F_v/F_m in each population (Table 2). The QTLs in population 1110 × 810 (QHst.cph-3B.1 and QHst.cph-3B.2) and 1110 × 1313 (QHst.cph-3B.3) were located on chromosome 3B while that in the population 1110 × 1039 (QHst.cph-1D) was located on chromosome 1D. For each of the three mapping populations the QTL region identified corresponded to the one found during the first round of mapping as shown by the presence of common flanking SNPs (Supplementary Table S2); with a slightly higher LOD score and percentage phenotypic variation explained but the same additive and dominance effects (Table 2). The LOD score of each of the three QTLs increased with increased duration of the heat treatment, making it significant only on day 3. There were no significant QTLs in other regions of the genome in any of the three populations (similar to the first round), except that the QTL peak on day 2 (QHst.cph-3B.1) shifted slightly upstream of the day 3 QTL (QHst.cph-3B.2) in the 1110 × 810 population (Table 2). In this population, the QHst.cph-3B.2 peak was positioned at 20 cM flanked by markers 1218388s (SNP) and Xgwm389 (SSR), with the LOD score of 5.7 while the QHst.cph-3B.1 peak with the LOD score of 6.4 was positioned at 5 cM flanked by markers Xgpw8020 (SSR) and 1061426s (SNP) (Table 2). The QHst.cph-3B.1 and QHst.cph-3B.2, respectively explained about 22.1 and 25.4% of the phenotypic variation in the population with an additive effect of 0.002 and 0.003. The heat tolerant parent 810 donates the positive allele for F_v/F_m in both the QTLs with the same dominance effect of 0.002 (Table 2). Using one-LOD drop off interval as an approximate confidence interval for the QTL position, the identified QHst.cph-3B.1 and QHst.cph-3B.2 QTLs spanned between 0 and 9 cM and 17 and 22 cM, respectively, where the two flanking markers were the only markers mapped at these particular intervals (Figure 3). However, there were more markers mapped in the vicinity (including the binned ones, a total of 28 markers mapped together at position 12 cM and six markers together at position 14 cM and eight markers together at position 15 cM on the current genetic map of 1110 × 810) (Supplementary Table S2). On the wheat physical map, the identified genomic region with the two QTLs corresponded to a region between 0.3 and 15 Mb (Figure 3). In this region, there were seven genes with known functions on the current assembly of wheat chromosome 3B hosted at URGI-Jbrowse database on traes3bPseudomoleculeV1 (Figure 3). These seven genes interestingly were all related to photosynthesis and heat stress and are therefore, considered as potential candidate genes identified from the QHst.cph-3B.1 and QHst.cph-3B.2 regions (Table 3). This includes frk2 (fructokinase 2) and bglu26 (beta-glucosidase 26) that are involved in carbohydrate metabolism. Further, there were two genes, ndhB2 [chloroplastic NAD(P)H-quinone oxidoreductase subunit 2B] and psaC (photosystem I iron-sulfur center), both having a direct role in the photosynthetic light reaction (Table 3). Besides, there were a gene (BUD31/G10 related) with a conserved site, and two genes encoding chloroplastic 3-isopropylmalate dehydrogenase 2) having known function in metal binding.

In the 1110 × 1313 population the identified QTL (QHst.cph-3B.3) for day 3 F_v/F_m was located downstream of both QHst.cph-3B.1 and QHst.cph-3B.2 identified in the 1110 × 810 population on same chromosome, as indicated by the position of common markers mapped between the two mapping populations (compare bold markers in Figures 3, 4). This QTL was highly significant with a LOD score of 8.9 explaining 35% of the phenotypic variation for F_v/F_m in the population. The heat tolerant parent 1313 contributed the positive allele for increasing F_v/F_m in the population with an additive effect of 0.003 and dominance effect of 0.001 (Table 2). QHst.cph-3B.3 was positioned at 37 cM flanked by markers 1175840p (PAV) and 1127409s (SNP) with one-LOD drop off interval between 35 and 39 cM on the 1110 × 1313 linkage map. This interval consisted of 11 markers including the binned ones. On the wheat physical map, this QTL region corresponded to a physical location approximately between 26 and 30 Mb, a region that holds a single potential candidate genes that have a known function and interestingly, it is psb28, class 1 that encodes a PSII reaction center protein (Figure 4). This protein is known to have a direct involvement in the oxygen evolving complex, biogenesis, assembly, stabilization, and repair of PSII complex (Table 3).

In the 1110 × 1039 population, the identified QTL (QHst.cph-1D) showed a LOD score of 5 and explained about 13% of the phenotypic variation for day 3 F_v/F_m in the population. Similar to the other two populations, the source of heat tolerance originates from the tolerant parent 1039, contributing the positive additive effect of 0.002 and dominance effect of 0.001 (Table 2). The QTL peak is positioned at 49 cM on chromosome 1D with one-LOD drop off interval at 49–51 cM. The QTL is flanked by two PAV markers 985618p and 1698203p mapped at 48.9 and 52.2 cM, respectively (Figure 5). This interval including the binned markers holds four markers mapped in the 1110 × 1039 linkage map (Supplementary Table S2). The genetic location of QHst.cph-1D corresponded to an approximate physical location somewhere between 13 and 32 Mb on the assembly of wheat chromosome 1D hosted at Ensembl Plants (release 29, December 2015). This region holds four potential candidate genes considering only the genes having known functions related to photosynthesis and heat stress (Figure 5). These genes include heme peroxidase known to be involved in the response to oxidative stress; αgalactosidase that is involved in carbohydrate metabolism; psbK encoding a PSII reaction center protein having a direct involvement in the assembly, stability, and repair of PSII complex; a DNAJ hsp—a heat shock protein (Table 3).