Meta-QTLs, ortho-MQTLs, and candidate genes for thermotolerance in wheat (Triticum aestivum L.)

Meta-QTL analysis for thermotolerance in wheat was conducted to identify robust meta-QTLs (MQTLs). In this study, 441 QTLs related to 31 heat-responsive traits were projected on the consensus map with 50,310 markers. This exercise resulted in the identification of 85 MQTLs with confidence interval (CI) ranging from 0.11 to 34.9 cM with an average of 5.6 cM. This amounted to a 2.96-fold reduction relative to the mean CI (16.5 cM) of the QTLs used. Seventy-seven (77) of these MQTLs were also compared and verified with the results of recent genome-wide association studies (GWAS). The 85 MQTLs included seven MQTLs that are particularly useful for breeding purposes (we called them breeders’ MQTLs). Seven ortho-MQTLs between wheat and rice genomes were also identified using synteny and collinearity. The MQTLs were used for the identification of 1,704 candidate genes (CGs). In silico expression analysis of these CGs permitted identification of 182 differentially expressed genes (DEGs), which included 36 high confidence CGs with known functions previously reported to be important for thermotolerance. These high confidence CGs encoded proteins belonging to the following families: protein kinase, WD40 repeat, glycosyltransferase, ribosomal protein, SNARE associated Golgi protein, GDSL lipase/esterase, SANT/Myb domain, K homology domain, etc. Thus, the present study resulted in the identification of MQTLs (including breeders’ MQTLs), ortho-MQTLs, and underlying CGs, which could prove useful not only for molecular breeding for the development of thermotolerant wheat cultivars but also for future studies focused on understanding the molecular basis of thermotolerance.


Introduction
Wheat (Triticum aestivum L.) is grown in the tropical and sub-tropical regions of the world, and is plagued with various abiotic stresses . Heat, drought, floods, salinity, cold, and toxic chemicals are all prominent abiotic stresses, which adversely affect crop productivity (Lesk et al. 2016). At the global level, heat and drought are recognized as the two most widely occurring abiotic stresses reducing wheat yield significantly . According to a global climate change model, the average ambient temperature is expected to rise by 6 °C by the end of the twenty-first century (De Costa 2011). In India, the temperature during the wheat growing season increased at a rate of 0.04 °C per decade during the twentieth century   (Chauhan et al. 2014). A recent estimate suggested that wheat yield can sometimes be reduced up to 60% in the field due to extreme heat stress (Zampieri et al. 2017). According to several studies on the effect of rising temperature on wheat grain yield, it is now known that every 1 °C rise in average temperature reduces global wheat yield by 6% (Asseng et al. 2011). This rise of 1 °C during the grain filling period can cause a yield loss of 0.25 tons ha −1 (~ 8%, assuming a global average yield of 3,200 kg ha −1 ) in winter wheat (Kristensen et al. 2011). In another study, a rise of 2 °C was estimated to reduce global wheat yields by 11% (Zaveri and Lobell 2019). Also, based on simulation studies involving agronomic trials conducted during 1985 to 2011, a 1 °C rise in mean temperature during reproductive stages was shown to cause 21% reduction in wheat yield (Barkley et al. 2014).
Studies have also been conducted on the effect of heat stress on physiological traits associated with yield losses. Based on these studies, reduction in wheat yield due to heat stress has been shown to be associated with the following traits: poor seed germination, inactivation of photosynthetic enzymes, reduced photosynthetic capability, slower nutrient transport, lower chlorophyll content, premature leaf senescence, poor pollen viability/germination, shorter grain filling period, fewer and lighter grains and poor starch/protein content resulting into poor bread-making quality (Akter and Rafiqul Islam 2017;Bheemanahalli et al. 2019;Jacott and Boden 2020;Jagdish 2020;Djanaguiraman et al. 2020;Poudel and Poudel 2020).
In view of the above, a breeding strategy that enables plants to become resilient to high-temperature stress appears to be an effective and sustainable solution. Therefore, the development of heat-tolerant, high-yielding wheat varieties, and/or pre-breeding materials for future breeding programs aimed at breeding for heat tolerance is crucial for food and nutritional security. Therefore, a complete understanding of the genetics of thermotolerance, which depends on a number of interacting QTLs/genes, should certainly help in the development of heat-tolerant cultivars.
Hundreds of QTLs are now known in wheat {including both tetraploid (4x) and hexaploid (6x)}, which control various agronomic and quality traits under heat stress (WheatQTLdb by Singh et al. 2021; http:// www. wheat qtldb. net/). However, limited efforts have been made for the discovery of meta-QTLs (MQTLs) for heat tolerance in wheat (Acuña-Galindo et al. 2015;Liu et al. 2020). BioMercator (Arcade et al. 2004;Sosnowski et al. 2012) and MetaQTL (Veyrieras et al. 2007) are the two most versatile software packages available for meta-QTL analysis. Algorithms furnished in BioMercator v2.0 estimate the most appropriate MQTLs using the available information. Meta-QTL analysis may be followed by the selection of promising MQTLs (breeders' QTLs) for marker-assisted selection (MAS) (Goffinet and Gerber 2000;Yu et al. 2014;Maccaferri et al. 2015) and may also be used for identification of candidate genes (CGs) for future studies.
Two studies, which have already been conducted on MQTL analysis for tolerance against heat stress in wheat (Acuña-Galindo et al. 2015;Liu et al. 2020), gave only limited information. In one study, two MQTLs and in the other study, 36 MQTLs for heat tolerance alone were reported (Acuña-Galindo et al. 2015;Liu et al. 2020). Acuña-Galindo et al. (2015) reported 43 MQTLs, which provided tolerance against both heat and drought stresses. The recent study by Liu et al. (2020) reported 59 MQTLs using QTLs associated with grain yield and its major components only, thus leaving out several other important traits such as stay green habit, chlorophyll content, and chlorophyll fluorescence. Since a large number of studies involving both interval mapping and GWAS for heat stress tolerance have been reported in wheat during the past decade (see WheatQTLdb by Singh et al. 2021), further MQTL analysis for heat-responsive traits is warranted.
The present study on MQTL and ortho-MQTL analysis was undertaken using the information on available QTLs for a number of traits related to heat tolerance (or thermotolerance) in wheat. MQTLs in rice along with information on synteny and collinearity between wheat and rice genomes (Ahn et al. 1993;Kumar et al. 2009) were also utilized for identification of ortho-MQTLs. The MQTLs identified in the present study were also verified using available information on GWAS and known genes for thermotolerance in wheat. Important genes including high confidence CGs relevant for thermotolerance in wheat were also identified using the information on MQTLs and the reference sequence of the hexaploid wheat genome (Zhu et al. 2021).

Materials and methods
Bibliographic search and collection of data on QTLs A comprehensive search using Google Scholar (https:// schol ar. google. com/) and PubMed (https:// pubmed. ncbi. nlm. nih. gov/) was made for reports on QTLs associated with thermotolerance in wheat {including both 4 × and 6x (spring and winter) wheats}; this information was supplemented by the information available in WheatQTLdb (see Singh et al. 2021; http:// www. wheat qtldb. net/). The details of parental crosses, populations type/size, traits and markers associated with QTL are given in Online Resource 1. The following information for each QTL was collected: (i) traits/parameters relevant to thermotolerance, (ii) most closely associated markers, (iii) chromosomal location of associated markers and QTLs, (iv) logarithm of the odds (LOD) value for each QTL, and (v) percent phenotypic variation explained (PVE%) or R 2 value for individual QTLs (Online Resource 2). All QTLs were treated as single QTLs, even if they were discovered for more than one trait or detected in more than one environment. Some studies reported QTLs separately under normal (control) and heat-stress conditions; however, we considered only QTLs that were identified under heat stress conditions.

Construction of the consensus map
A saturated consensus genetic map (described as Wheat_Consensus_map_2021) was generated using the following integrated genetic maps: (i) the "Wheat Consensus SSR 2004" with 1258 marker loci (Somers et al. 2004); (ii) an integrated map for durum wheat with 3618 markers (Marone et al. 2013); (iii) the SNP-based consensus map of tetraploid wheat with 30,144 markers (Maccaferri et al. 2015); (iv) SNP-based genetic map constructed with 46,977 loci from the Illumina iSelect 90K SNP Array (Wang et al. 2014); individual genetic maps for different mapping populations reported in previous studies on QTL analysis (Online Resources 1 and 2). The number of common and unique markers among the four individual maps was calculated and illustrated using custom Venn diagrams webtool (https:// bioin forma tics. psb. ugent. be/ webto ols/ Venn/). The R package LPMerge (Endelman and Plomion, 2014) was used for the construction of the Wheat_Consensus_map_2021. This map was used as a reference map for QTL projection and identification of MQTLs.

Projection of QTLs on the consensus map
The projection of QTLs onto the Wheat_Consen-sus_map_2021 was performed by BioMercator v4.2 (https:// urgi. versa illes. inra. fr/ Tools/ BioMe rcator-V4; Arcade et al. 2004) using the QTLProj command. QTLs were projected on the consensus map following the homothetic approach proposed by Chardon et al. (2004). When confidence interval (CI) was not reported in the original studies, the genetic distance between flanking markers was taken as the CI. When the QTL peak position was available, CI (95%) was calculated for each locus using different population-specific equations (Darvasi and Soller 1997;Guo et al. 2006;Liu et al. 2009), which included the following: (i) for backcross and F 2 populations: CI = 530/NR 2 ; (ii) for doubled haploid (DH) lines: CI = 287/NR 2 , and (iii) for recombinant inbred lines (RILs), CI was calculated as 163/NR 2 , where N is the population size and R 2 is the PVE% by the QTL.
Meta-analysis of the QTLs QTL meta-analysis was performed for each wheat chromosome, separately, using algorithms implemented in BioMercator v4.2 (Goffinet and Gerber 2000;Veyrieras et al. 2007). Depending on the number of projected QTLs per chromosome, the following two distinct approaches were applied: (i) when the number of QTLs per chromosome was 10 or less, the methods proposed by Goffinet and Gerber (2000) was applied, and (ii) when the number of QTLs per chromosome was > 10, the method proposed by Veyrieras et al. (2007) was utilized. In the first approach, the model with the lowest Akaike Information Criterion (AIC) value was used to select the best QTL model for determining the number of MQTLs on each consensus chromosome; this approach was followed in several earlier meta-QTL studies (Lu et al. 2018;Venske et al. 2019;Liu et al. 2020). In the second approach, the following five parameters were used to ascertain the correct number of MQTLs (Veyrieras et al. 2007): (i) Akaike Information Criterion (AIC), (ii) AIC correction (AICc), (iii) AIC3 candidate model, (iv) Bayesian Information Criterion (BIC), and (v) Approximate Weight of Evidence Criterion (AWE). The model with the lowest value in at least three of the five parameters was chosen as the best fit ). In our meta-analysis, the QTLs identified from the best-fit model were considered as the MQTL. The number of MQTLs on each wheat chromosome was determined based on the genetic positions and 95% CI. Detailed information on statistical procedures and algorithms implemented in the BioMercator are available elsewhere (Arcade et al. 2004;Veyrieras et al. 2007;Sosnowski et al. 2012).
The physical positions for individual MQTLs were obtained using the nucleotide sequences of the markers flanking the MQTLs; these nucleotide sequences were obtained from published literature or retrieved from GrainGenes (https:// wheat. pw. usda. gov/ GG3/) and/or CerealDB (https:// www. cerea lsdb. uk. net/ cerea lgeno mics/ Cerea lsDB/ index NEW.php) databases. The marker sequences were used for BLAST against the hexaploid wheat reference genome (IWGSC RefSeq v1.0). In some cases, the physical positions of SNP markers were directly obtained from the JBrowse wheat genome browser (https:// wheat-urgi. versa illes. inra. fr/ Tools/ Jbrow se). The genomic region between the flanking markers was treated as the physical position of the corresponding MQTL.

Comparison of MQTLs with known GWAS-MTAs and genes
MTAs reported in five earlier genome-wide association studies (GWAS; published during 2019-2020) on heat stress-responsive traits were available; these were utilized for a comparison with MQTLs detected in the present study. The details of the size of association mapping panel, traits-parameters, SNP markers used, and MTAs identified in these studies are given in Online Resource 3. The physical positions of the MTAs obtained from different databases (such as CerealDB and JBrowse wheat genome browser) were compared with the physical coordinates of the MQTL regions on the same chromosome; an individual MTA occurring within a specific MQTL region was considered co-located. In this manner, MQTLs detected in this study were verified by the MTAs reported in at least one of the five GWAS reports.
In addition to the above, the previously known heat stress-responsive gene sequences were subjected to BLAST analysis against the IWGSC RefSeq v1.0 and the physical coordinates of known genes were extracted, and compared with the physical positions of the MQTL regions. The MQTLs matching with sequences of heat-responsive genes were treated as coincident.

Identification of genes and their expression analysis
The genes underlying individual MQTL regions were identified in the 1-Mb interval on either side (upstream and downstream) of the peak position of the MQTL (total 2 Mb region; Jan et al. 2021); for this purpose, the BioMart tool of EnsemblPlants (https:// plants. ensem bl. org/ bioma rt/ martv iew) was used. In silico expression analysis of genes was carried out using transcriptome data available at wheat expression browser (expression and visualization platform) powered by expVIP (http:// www. wheatexpre ssion. com; Ramírez-González et al. 2018). For this purpose, the relevant expression dataset including expression data related to heat stress was utilized (Liu et al. 2015). This RNA-seq dataset consisted of differential expression of genes in a heat and drought tolerant wheat cv. TAM107, grown under heat-stressed condition, with leaf samples collected separately at 1 h and 6 h after heat stress treatment (40 °C); seedlings grown in normal growth condition (22 °C, wellwatered) were taken as control (Liu et al. 2015). The data on expression was available as log2 transformed TPM (transcripts per million) values. Only genes showing fold change (FC ≥ 2 or FC ≤ − 2 TPM values) relative to control were considered as differentially expressed. Gene ontology (GO) analysis was conducted using the BioMart tool and the heat maps of the expression patterns were generated using the Morpheus program (https:// softw are. broad insti tute. org/ morph eus/).

Identification of ortho-MQTLs
Seven wheat MQTLs were used for the identification of ortho-or syntenic MQTLs in rice. Following steps were involved: (i) genes detected in the regions of the above seven MQTLs were subjected to BLAST analysis against rice genome database at EnsemblPlants to identify rice orthologues. (ii) The rice orthologues were extracted with their physical positions from the database. (iii) Physical positions of genes underlying wheat MQTLs were compared with the orthologous genes underlying rice MQTLs for heat tolerance (Raza et al. 2020), and (iv) MQTLs harboring similar genes located on known syntenic genomic positions of wheat and rice were considered as ortho-MQTLs.

Results and discussion
The QTL data used for conducting meta-QTL (MQTL) analysis in the present study were largely similar to those used in two previous studies (Acuña-Galindo et al. 2015;Liu et al. 2020). This also facilitated a comparison across three studies and interpretation of results. A detailed comparison with respect to the data used and the results obtained are presented in Online Resource 4. For identification of MQTLs involved in thermotolerance, integration of GWAS (for MQTL verification), heat-responsive genes, genes underlying MQTLs and their expression, and homology comparison with rice (ortho-MQTL) offer an interesting contrast between the present study and two earlier studies. The chromosomes in the consensus map were more densely populated with markers and the QTLs used in the present study. In our analysis, Wheat_Consensus_map_2021 was prepared with 50,310 markers and 441 QTLs were projected for 154 heat-responsive traits belonging to 31 categories (Online Resource 5), which was higher relative to those in previous two studies ( Acuña-Galindo et al. 2015;Liu et al. 2020). Further, out of 85 MQTLs detected in our study, the genetic position of 18 MQTLs matched with the MQTLs reported in above two studies (Acuña-Galindo et al. 2015;Liu et al. 2020), whereas the remaining 67 MQTLs were perhaps novel which were not reported earlier.
Interestingly, two of the 18 MQTLs were exclusively associated with thermotolerance, 14 MQTLs were associated either only with drought or heat/drought tolerance. However, the remaining two MQTLs were neither associated with heat tolerance nor with drought tolerance. Overall, the results suggest that ~ 20% (16/85) of MQTLs involved in thermotolerance as well as drought were common among the three studies including our own study and therefore provide a genetic basis for future studies on abiotic stress tolerance in wheat. Moreover, above 16 MQTL regions associated with heat and drought toleranceresponsive traits can be targeted for future studies involving fine mapping, cloning, functional analysis, and development of markers for use in markerassisted selection.
QTLs associated with thermotolerance Information on QTLs from 31 interval mapping reports (published during 2008-2020) for heatresponsive traits was collected (Online Resource 1). Only 28 of these studies (25 for hexaploid and three for tetraploid wheats) had complete information (marker positions and PVE values) needed for conducting meta-QTL analysis; 826 individual QTLs for 154 relevant traits were available from these 28 studies (Online Resources 1 and 2). The populations involved in these 28 studies included 13 recombinant inbred line (RIL) populations, 6 doubled haploid (DH) populations, 5 F 2 /backcross populations, and one family-based population (rarely same population was used in more than one study). The population size in different interval mapping studies varied from 64 (Mason et al. 2010) to 420 (Sharma et al. 2017) and the number of QTLs reported in individual studies ranged from 2 (Beecher et al. 2012) to 118 QTLs (Shirdelmoghanloo et al. 2016) (Fig. 1a).
The number of QTLs on individual chromosomes ranged from a minimum of 13 QTLs on chromosome 1D to a maximum of 90 on chromosome 3B (Online Resource 2). Distribution of QTLs on individual homoeologous groups also differed, with homoeologous group 2 having a maximum of 197 QTLs, and group 6 having a minimum of 66 QTLs. This distribution of QTLs on the seven homoeologous groups may be attributed to higher gene density of stressrelated expressed sequence tags (ESTs) on group 2 chromosomes and the lowest gene density on group 6 chromosomes as reported in an earlier study in wheat (Ramalingam et al. 2006). The number of QTLs for different traits ranged from 2 (spike length; SL) to 164 (bread-making quality; BMQ), which comprise 25 individual traits-parameters, with an average of ~ 26 QTLs per category of traits (Fig. 1b). The distribution of QTLs in three sub-genomes also differed with sub-genome B carrying the maximum of 350 QTLs, followed by sub-genomes A with 282 and sub-genome D with 194 QTLs (Fig. 1c). This distribution of QTLs on homoeologous groups and subgenomes is also in agreement with previous MQTL studies in wheat for thermotolerance (Acuña-Galindo et al. 2015;Liu et al. 2020) and drought tolerance . The 826 QTLs were further categorized based on LOD scores and PVE%. The LOD score of the individual QTLs ranged from 2 to 60 with an average of 6.22 (Fig. 1d). Nearly half (46.7%) of all QTLs had their LOD scores ranging from 2 to 4. The PVE% for individual QTLs ranged from 1 to 60.1% with a mean PVE of 12%. Approximately 75% (619/826) QTLs had up to 15% PVE, whereas about Traits abbreviations include GY grain yield, GW grain weight, GN grain number, DTH days to heading, DTA days to anthesis, DTM days to maturity, GFD grain filling duration, HI harvest index, BM biomass, PT productive tiller, KD kernel diameter, SL spike length, SD spike density, SP spikes per plant, PH plant height, BMQ bread-making quality traits, SG stay green, CF chlorophyll fluorescence, CC chlorophyll content, CT canopy temperature, NDVI normalized difference vegetation index, HSI heat susceptibility index, TT thermotolerance, FLD flag leaf dimension, WX waxiness, TMD thylakoid membrane damage, KH kernel hardness, PL length of the ear peduncle, PMD plasma membrane damage, WSC water soluble carbohydrates, and PC proline content 5% (41/826) QTLs exhibited > 30% PVE (Fig. 1e), suggesting the availability of both major and minor QTLs for thermotolerance.

Consensus map of wheat (Wheat_Consensus_ map_2021)
The marker types on the Wheat_Consensus_ map_2021 included RFLPs, AFLPs, SSRs, DArTs, SNPs, and some well-known gene loci, such as Vrn, Ppd, Rht, and Glu. The map contained 50,310 marker loci with a total genetic length of 12,280.7 cM, giving a density of 4.1 marker loci/cM ( Fig. S1 and S2, Online Resource 6). The length of the linkage maps for the individual chromosomes ranged from 280.1 cM (4D) to 894.1 cM (4A), and the number of markers in individual linkage groups ranged from 392 on 4D to 4,121 on chromosome 1B. The marker density for individual chromosomes ranged from 1.4 markers/cM for 4D to 11.1 markers/cM for 2A. Sub-genome A covered 4,102.7 cM with 19,518 markers (4.76 markers/cM), sub-genome B covered 4,579.9 cM with 23,058 markers (5.03 markers/cM), and sub-genome D spanned 3,598.1 cM size with 7,734 markers (2.15 markers/cM) (

QTL projection and MQTL analysis
The distribution of the number of available QTL, the QTLs utilized for projection, and the MQTLs identified during the present study among different wheats (durum wheats, spring 6 × wheats, and winter 6 × wheats) are summarized in Table 1. From the 826 individual QTLs, 441 (53.4%) QTLs for different heat-responsive traits were successfully projected onto the above Consensus map. These 441 QTLs included 27 QTLs derived from tetraploid wheat (durum) and 414 from hexaploid wheat (involving 244 from spring wheats, 65 from winter wheats, 105 from spring × winter crosses). The remaining 385 (46.6%) QTLs could not be projected due to the following reasons: (i) use of marker types that were not used in consensus map, (ii) lack of commonly associated markers between consensus and initial maps, and (iii) very long CI.
As a result of meta-analysis, as many as 85 MQTL regions were identified from 411 QTL with at least one MQTL on each of the 21 wheat chromosomes ( Table 2 In addition to reduction in the number of QTLs, MQTL analysis also leads to reduction in the size of CI by accumulating QTLs information from different populations of diverse genetic origins. In the present study, the CI of the MQTLs ranged from 0.11 to 34.9 cM with a mean of 5.6 cM. The average length of CI of MQTLs was reduced 2.95-fold relative to the CIs in QTLs, which is quite similar to those reported in previous meta-QTL studies in wheat (Acuña-Galindo et al. 2015;Liu et al. 2020;   (1), spring wheat (7), winter wheat (2) (2), winter wheat (4) Yang et al. 2021). There were significant differences in average CIs of MQTLs that mapped on different chromosomes; CI of MQTLs on 6D was reduced 14.6 times and that on 5A was reduced 8.9 times, followed by 8.8 times for 7A and 8.1 times for 7B (Fig. S3).  Yang et al. (2021) showing that the CI was proportional to the number of QTLs on which the MQTL was based. These results and those of the present study suggest that an increase in the number of QTLs used for meta-QTL analysis could significantly reduce the size of the CI of MQTL relative to those in the initial QTLs. As expected, a solitary MQTL based exclusively on tetraploid (durum) wheat and those based on both tetraploid and hexaploid wheats belong to only subgenomes A and B (none belonging to sub-genome D), while those belonging exclusively to hexaploid wheat (spring and/or winter wheats) carried MQTLs belonging to all the three sub-genomes. It is interesting to note that out of 85 MQTLs identified in the present study, 25 MQTLs contained original QTLs exclusively from spring wheat and two MQTLs were such which contained original QTLs belonging exclusively to winter wheat; a solitary MQTL was however based exclusively on the original QTL from durum wheat (Table 2). Thus, our findings clearly dissected specific genomic regions that contribute to heat tolerance in spring and winter wheats, separately. The distribution of MQTLs on 21 individual wheat chromosomes ranged from 1 on 6D to 7 on 7A (Fig. 2, Online Resource 7) and this number did not depend on the number of QTLs on individual chromosomes utilized for MQTL analysis. For instance, chromosomes 2A, 2D, and 3B each had many more QTLs but relatively fewer MQTLs. Several MQTLs had a clustered distribution in the genome, with some having overlapping regions. One of these clusters, consisting of 4 MQTLs, was identified on chromosome 1A physically located between 2.54 and 12.09 Mb (1A.1-1A.4). All MQTLs (except 4B.2) on chromosome 4B were physically clustered in 5.19-255.38 Mb region. Similarly, five MQTLs (5B.1-5B.5) were found  7). Such a clustered pattern of distributions of MQTLs in wheat genome was also reported in three recent studies Saini et al. 2021;Yang et al. 2021). The number of individual QTLs covered within an individual MQTL ranged from 2 to 22 (Table 2) and a large proportion of MQTLs (65%, 55/85) was each based on > 3 QTLs derived from different mapping populations. The following six MQTLs comprised ≥ 12 QTLs: MQTL2B.2 (18), MQTL2B.3 (12), MQTL2D.2 (22), MQTL3B.2 (14), MQTL5A.2 (12), and MQTL7B.3 (13). This is in sharp contrast to earlier studies, where a single MQTL associated with thermotolerance was based on no more than 10 QTLs (Acuña-Galindo et al. 2015;Liu et al. 2020). The MQTLs that harbored original QTLs derived from more than one mapping population and consistently appeared in several environments may be the most reliable, robust, and stable for wheat breeding.
The number of traits controlled by an individual MQTL ranged from a solitary trait controlled by each of many MQTLs to as many as 14 traits controlled by a solitary MQTL (e.g., MQTL2D.1); 76 of the 85 MQTLs were multi-trait MQTLs (Table 2), each controlling more than one trait indicating either the occurrence of pleiotropic genes or tight linkage among genes for different traits. It could also be due to the use of related traits as also shown in the case of MQTL analysis in wheat for grain yield and related traits (Saini et al. 2021) and root-related traits in wheat (Soriano and Alvaro, 2019).
The LOD score of individual MQTLs ranged from 2.7 to 29.8 with an average of 5.5. Similarly, the PVE ranged from 2.5 to 26.5% with a mean of 11.5%. The MQTLs were subjected to further selection for identification of some promising MQTLs that we liked to call "breeders' MQTLs," each fulfilling the following conditions: (i) CI < 2 cM, (ii) average PVE > 10%, (iii) average LOD > 4, and (iv) involvement of at least 7 QTLs within the MQTL. This exercise resulted in the selection of seven promising MQTLs (2A.4, 3B.2, 5B.3, 6D.1, 7A.1, 7B.1, and 7B.3) that contained original QTLs belonging to both spring and winter wheats. These seven promising MQTL regions overlapped genomic regions carrying MTAs reported in earlier GWA studies (see later). Although three more MQTLs (viz., 1B.1, 6A.2, and 6B.2) each had high PVE (> 20%), these were not included in the list of breeders' MQTLs, due to their association with large CI or due to the fact that these were each based on fewer QTLs. The so-called breeders' MQTLs should be subjected to further scrutiny (including validation) for their utility in MAS and may also be used as a resource for fine mapping and cloning.

Verification of MQTLs using MTAs from GWAS
The physical coordinates of MQTLs were also compared with physical positions of 1,034 MTAs for the traits associated with thermotolerance reported earlier in GWA studies; this resulted in the identification of MQTLs matching with as many as 704 MTAs (involving 139 MTAs, each associated with multiple heat-responsive traits within an individual or across several GWA studies) (Fig. 2, Online Resource 8). As many as 63 MQTLs were verified with MTAs obtained from GWAS involving spring wheats and 67 MQTLs were identified for winter wheat populations. However, 53 MQTLs were verified using GWAS involving both spring and winter wheats. The MQTL6B.5 matched with 110 MTAs and MQTL6B.1 matched with 103 MTAs. There were other MQTLs (e.g., 2B.2, 2B.3, and 5A.3) each involving 10 or more QTLs and matched with more than 60 MTAs (Online Resource 8). Overall, ~ 90% of the total MQTLs matched GWAS results. These findings are in sharp contrast to recent studies by Aduragbemi and Soriano (2021) and Yang et al. (2021), where only 38.66% and 62.68% MQTLs, respectively, overlapped MTAs obtained from GWA studies. This might be due to the availability of a comparatively large number of MTAs for comparison in the present study. The MQTLs, which matched GWAS-based MTAs, provide a basis for accurately mining of high confidence candidate genes associated with thermotolerance in wheat.
MQTLs co-located with known genes for thermotolerance A comparison of MQTLs with known genes for thermotolerance in wheat may also facilitate the efforts being made towards molecular dissection of genomic regions involved in thermotolerance. Therefore, association of MQTLs with known thermotolerance genes was also examined. Nine such MQTLs were associated with seven genes known for thermotolerance in wheat (Table 3). These known genes associated with thermotolerance included three different classes: (i) Following four HSP/HSF genes were found in the following seven MQTL regions: TaHSP17   Only the flanking markers are shown for better visualization. MQTLs in blue boxes (right side) were co-located with marker-trait associations (MTAs) identified using genome-wide association studies thermotolerance (Kumar et al. 2012;Campbell et al. 2001;Vishwakarma et al. 2018;Xue et al. 2015).
(ii) The gene TaGF14-JM22 (encoding 14-3-3 protein-like) corresponding to MQTL4A.3 is known to be involved in starch biosynthesis. This gene has a negative association with thermotolerance so that its reduced expression enhances thermotolerance (Guo et al. 2018). (iii) TaFER-5B and TaPEPKR2 genes, respectively, encoding for ferritin and phosphoenolpyruvate carboxylase kinase-related kinase overlap MQTL5B.1. TaFER-5B is known to provide thermotolerance through scavenging reactive oxygen species (Zang et al. 2017) and TaPEPKR2 has been characterized for thermotolerance and drought tolerance in wheat and Arabidopsis (Zang et al. 2018). These data suggest that the nine MQTLs (Table 3), which correspond to seven known genes for thermotolerance, are promising candidates for the improvement of this trait and could be useful in wheat breeding for climate resilience. Additional efforts are needed to identify functional variants of these genes in the targeted MQTL regions.

Identification of candidate genes, in silico expression, and GO analysis
Mining of genes in genomic regions of all the MQTLs (except MQTL1D.1, 4B.2, 5A.3) allowed identification of 1,704 non-redundant CGs (Online Resource 9), which included the following: (i) 165 CGs including 55 CGs with each of the three MQTLs (MQTL1A.4, 2A.1, 7D.4) and a solitary CG each with two MQTL (MQTL4B.5 and 4D.4). (ii) 260 genes, for which no description for functions was available, and another 78 CGs, which were either uncharacterized or with no known encoded proteins (Online Resource 9). (iii) A large number of CGs, each encoding a variety of proteins products in more than one MQTL region showing redundancy of genes in the wheat genome. These included 102 CGs for protein kinases, 55 CGs for F-box-like domains, 43 CGs for cytochrome P450, 38 CGs for proteins with leucine-rich repeats, 37 CGs for zinc finger RING/ FYVE/PHD-type, and 27 CGs for pentatricopeptide repeats (Fig. S4).
Functional relevance of the CGs/gene families listed above in category (iii) with thermotolerance has also been documented previously (Park et al. 2014;Agarwal and Khurana 2018;Chen et al. 2018;Pandian et al. 2020;Kumar et al. 2021;Guérin et al. 2021). The following are some examples: (i) a mitogen-activated protein kinase (MAPK) was recently shown to be associated with oxidative stress tolerance and carbon assimilation in wheat under terminal heat stress ); (ii) F-box-like proteins are known to play a major role in plant responses to diverse developmental and stress conditions in wheat (Guérin et al. 2021); (iii) cytochrome P450s are the largest enzyme family involved in O 2 and/or NADPHdependent hydroxylation reactions. The expressions of cytochrome P450 genes have been shown to play a prominent role in the crosstalk between abiotic and biotic stress responses (Pandian et al. 2020); (iv) leucine-rich repeats are known to act as receptors of external signals, e.g., abiotic stresses such as heat, salt, and osmotic stresses (Park et al. 2014).
The above CGs were also subjected to an in silico expression analysis under 1 h and 6 h of heat stress using RNA-seq data available in a public database. As many as 182 CGs, belonging to 71 MQTLs, exhibited differential expression; these differentially expressed CGs (DECGs) were categorized into upregulating (FC ≥ 2) and down-regulating (FC ≤ − 2) DECGs. From these DECGs, 92 were differentially expressed at 1 h of heat stress (16 up-regulated and 76 down-regulated), 116 were differentially expressed at 6 h of heat stress (55 up-regulated and 61 downregulated), and only 31 were differentially expressed at both conditions (10 up-regulated, 15 down-regulated, 6 up-regulated at one time-point and downregulated at the other). A representative heat map of some important DECGs is shown in Fig. 3. The genes with uncharacterized proteins need to be subjected to further investigations for determining their roles in thermotolerance.
The above 182 DECGs (belonging to 71 MQTLs) were also analyzed for gene ontology (GO) enrichment (Online Resource 9). The GO terms associated with biological processes belonged to metabolic (12 CGs) and cellular processes (6 CGs) (Fig. S5). Similarly, GO terms associated with molecular function were for heterocyclic compound binding (25 CGs) and catalytic activity (19 CGs). In relation with cellular components, the CGs were enriched mainly in the cell membrane and its components such as intracellular anatomical structure and cell periphery.
Among the above 1,704 genes, thirty-six (36) high confidence CGs were also detected (Table 4), which either exhibited differential expressions under heat stress or were known to be important for traits related to thermotolerance. The role of some of the important high confidence CGs in providing thermotolerance can be summarized as follows: (i) CGs encoding UDP-glucosyltransferase redirect the metabolic flux from lignin biosynthesis to flavonoid biosynthesis under abiotic stress (heat, drought, and salt stress) and the accumulation of flavonoid glycosides providing protection against abiotic stress leading to improvement in grain size (Dong et al. 2020). (ii) Over-expression of a CG encoding a zinc finger transcription factor conferred improved thermotolerance and grain yield, as demonstrated in Arabidopsis thaliana (Agarwal and Khurana 2018), (iii) CGs encoding soluble N-ethylmaleimide-sensitive factor attachment protein receptors (SNAREs) are also known to be involved in abiotic and biotic stresses, as reported in several plant species (Kwon et al. 2020), (iv) over-expression of MYB genes provided improved thermotolerance and drought tolerance; this has been attributed to enhanced levels of cellular abscisic acid in wheat , (v) a KH domain-containing putative RNA-binding protein was reported to be crucial for the regulation of heat stress-responsive CG and thermotolerance in A. thaliana (Guan et al. 2013), and (vi) a CG encoding C2 domain protein was reported to play a key role in the improvement of A. thaliana thermotolerance, mainly through the enhancement of antioxidant enzyme activity, leading to enhanced expression of genes for response to heat stress response (Ye et al. 2020). While these potential CGs are promising targets for future breeding and genetic studies, it is not apparent whether they are the actual regulators behind the MQTLs, since many other genes occur within the CIs of the MQTLs. Further studies are warranted to validate the role of these CGs.

Ortho-MQTLs for wheat and rice
Chromosomal regions carrying MQTLs that are conserved across more than one cereal were also identified. For this purpose, MQTLs for different traits relevant to thermotolerance in the present study in wheat were compared with MQTLs reported in rice (Raza et al. 2020). Ortho-MQTLs could be identified for only seven wheat MQTLs positioned on chromosomes 1B, 2B, 2D, 3B, 7B, and 7D that correspond to rice-MQTLs (rMQTLs) located on rice chromosomes 1, 4, 5, 6, 7, and 8. This confirms the occurrence of wheat MQTLs that are orthologous to rMQTLs. Detailed information related to ortho-MQTLs and underlying genes along with physical positions and functional description is presented in Table 5 and Online Resource 10. Among ortho-MQTLs, two ortho-MQTLs (ortho-MQTL2B.3 and ortho-MQTL2D.1) were considered most promising for revealing gene collinearity within a syntenic region between wheat and rice. For ortho-MQTL2B.3, four underlying CGs on chromosome 2B were found to be syntenous and collinear with similar CGs belonging to rMQTL4.4 on rice chromosome 4 ( Fig. 4a and Table 5). The wheat ortho-MQTL2D.1 on chromosome 2D was syntenous with rMQTL7.1 on rice chromosome 7 ( Fig. 4b and Table 5). The orthologous genes located at the two most promising syntenic regions (viz., ortho-MQTL2B.3 and 2D.1) were subjected to further analysis. Ortho-MQTL2B.3 contained at least two heat-responsive genes encoding for protein kinase (Sinha et al. 2011) and amino acid/polyamine transporter (Shen et al. 2016). Ortho-MQTL2D.1 included at least 5 heat-responsive genes encoding for the following proteins/products: NAC domain-containing protein (Ren et al. 2021), nucleoside phosphatase (Escobar et al. 2001), PLD-regulated protein1-like (Mishkind et al. 2009), RNA recognition motif domain (Lee et al. 2017), and helicase (Gong et al. 2005). The crucial roles of some of these genes associated with ortho-MQTLs in conferring thermotolerance can be summarized as follows: (i) mitogen-activated protein kinase cascades have been shown to be involved in signalling pathways activated by different abiotic stresses including heat stress in different plant species (Sinha et al. 2011); (ii) polyamine transporter is required for the uptake of extracellular polyamines and is known to play a key role in stabilizing the mRNAs of several important heat stress-responsive genes under heat stress in Arabidopsis (Shen et al. 2016); (iii) CRISPR/Cas9 based genome editing and transgenesis demonstrated a transcriptional regulatory network with two NAC genes (ONAC127 and ONAC129) coordinating multiple pathways to modulate grain development and heat stress responses at reproductive stages in rice (Ren et al. 2021); (iv) a mitochondrial nucleoside diphosphate kinase (mtNDPK) is known to be involved in heat stress response in pea, possibly as a modulator of a protein (86 kilo-dalton protein) synthesized upon exposure to heat (Escobar et al. 2001); (v) phospholipase D and phosphatidylinositol 4,5-bisphosphate are accumulated at the plasma membrane and nucleus, thereby regulate the ion channels and the cytoskeleton under heat stress in Arabidopsis and rice (Mishkind et al. 2009); and (vi) the RNA recognition motif of Arabidopsis splicing factor affects the alternative splicing   (Lee et al. 2017). Among the remaining five ortho-MQTLs, only one ortholog matched between wheat and rice MQTLs (Table 5).
Ortho-MQTLs in wheat in particular and cereals in general have been examined sparingly; at least three such studies are available (Khahani et al. 2020(Khahani et al. , 2021Saini et al. 2021). The discovery of ortho-MQTLs in closely related species increases their significance and verifies their stability, as well as the confidence of associated genes (Khahani et al. 2020, Table 5 Identification of ortho-MQTLs between wheat and rice genomes a Genomic position of meta-QTL (MQTL) peak region (total 2 Mb region), b Raza et al. (2020), MB mega base

Fig. 4
Representing syntenic regions of two ortho-MQTLs between wheat and rice. The chromosome number, physical position, and common genes between the wheat and rice are indicated. More details are given in Online Resource 10 2021). The orthologous genes detected in the present study may be functionally validated and gene-based functional markers can be developed for their use in cereal breeding programs (Quraishi et al. 2009).

Conclusions
The present study involved identification of MQTLs using information about the known QTLs for thermotolerance. Improved meta-analysis strategy was used leading to the identification of 85 most stable and robust MQTLs; 77 of these MQTLs were confirmed using available results of recent GWA studies. Seven promising MQTLs (breeders' MQTLs) were recommended for use in MAB for improvement of thermotolerance in wheat. Ortho-MQTLs for wheat and rice identified in the present study allowed identification of conserved genomic regions associated with heat-responsive traits in rice. As many as 1,704 CGs underlying 82 MQTLs were also identified which included 182 CGs with differential expression under heat stress. In addition, 36 high confidence CGs underlying 24 MQTLs with known functions were identified as being reliable for thermotolerance. These CGs may be validated using one or more of the available approaches (including transgenesis, gene editing, and gene knockout strategies). The results may also be used for studies involving proteomics and metabolomics. On validation, these CGs can be used in breeding programs for improving thermotolerance and may also be used as a resource for designing research programs for understanding the mechanism for thermotolerance in wheat.