Selection of Candidate Genes Conferring Blast Resistance and Heat Tolerance in Rice through Integration of Meta-QTLs and RNA-Seq

Due to global warming, high temperature is a significant environmental stress for rice production. Rice (Oryza sativa L.), one of the most crucial cereal crops, is also seriously devastated by Magnaporthe oryzae. Therefore, it is essential to breed new rice cultivars with blast and heat tolerance. Although progress had been made in QTL mapping and RNA-seq analysis in rice in response to blast and heat stresses, there are few reports on simultaneously mining blast-resistant and heat-tolerant genes. In this study, we separately conducted meta-analysis of 839 blast-resistant and 308 heat-tolerant QTLs in rice. Consequently, 7054 genes were identified in 67 blast-resistant meta-QTLs with an average interval of 1.00 Mb. Likewise, 6425 genes were obtained in 40 heat-tolerant meta-QTLs with an average interval of 1.49 Mb. Additionally, using differentially expressed genes (DEGs) in the previous research and GO enrichment analysis, 55 DEGs were co-located on the common regions of 16 blast-resistant and 14 heat-tolerant meta-QTLs. Among, OsChib3H-c, OsJAMyb, Pi-k, OsWAK1, OsMT2b, OsTPS3, OsHI-LOX, OsACLA-2 and OsGS2 were the significant candidate genes to be further investigated. These results could provide the gene resources for rice breeding with excellent resistance to these 2 stresses, and help to understand how plants response to the combination stresses of blast fungus and high temperature.


Introduction
Rice is a significant staple food worldwide that provides more than 20% of the daily caloric needs for at least 50% of the global population [1]. To meet the demand of a growing global population, rice yields need a yearly increase of 0.6 to 0.9% [2]. However, rice frequently suffers from various biotic and abiotic stresses in nature [3,4]. As a consequence of global warming, the combination of pathogens and high temperatures (HT) frequently exists in the cultivation of cereal crops [5][6][7][8]. As the "cancer" of rice, M. oryzae is widely distributed, and terribly destructive under favorable conditions, which can cause severe yield losses in rice [9][10][11]. It is estimated that yearly, the rice blast fungus can destroy enough rice to feed 60 million people [12]. Likewise, HT also poses a serious threat to rice growth and development that harshly affects rice yields and quality [13,14]. Currently, introgression of the resistant and tolerant genes has been proven to improve resistance and tolerance in existing rice cultivars [15][16][17][18]. Therefore, it is essential to mine genes conferring blast resistance and heat tolerance (BR-HT).
Previous studies have reported that heat induction could enhance plant resistance to pathogen stresses. On the one hand, HT induction can inhibit the growth of pathogens. For example, cucumber seedlings increased resistance to Cladosporium cucumerinum after HT preheating [19]. Susceptible barley varieties enhanced resistance to powdery mildew after preheating [20]. Likewise, melon also strongly resisted Botrytis cinerea after HT (LOD), phenotypic variance (R 2 ), flanking markers and physical interval of target QTLs. In addition, the QTLs without CIs were calculated according to the published calculation formula [58,59], as follows: CI = 530/(N × R 2 ) (1) where N represents the size of mapping population, and R 2 is for the phenotypic variance of target QTLs. Equation (1) is applicable to BILs, DH and F2 populations, and Equation (2) is applicable to RILs.

Meta-QTL Analysis
The physical positions of flanking markers along with target QTLs were determined by Gramene Marker Search (https://archive.gramene.org/db/markers/marker_view, accessed on 9 January 2022). Secondly, we prepared two separate input files (map file and qtl file) in txt format for each study. Following that, a consensus map was constructed, and all the maps with the markers and original QTLs were iteratively projected on a reference map (rice physical map) by Biomercatorv4.2 (https://sourcesup.renater.fr/projects/ biomercator/, accessed on 9 January 2022) [60,61]. Meta-QTLs of individual chromosomes were determined based on AIC [52]. Meta-QTLs projected by at least two initial QTLs were selected as results with high reliability, which were left for subsequent analysis [48,50,55]. Then, we utilized Mapchart software (https://www.wur.nl/en/show/Mapchart.htm, accessed on 9 January 2022) [62] to output the vector map of the consensus map. According to the physical positions and sequences of markers flanked at meta-QTL regions, the physical distances of meta-QTL intervals and adjacent genes were determined by NCBI BLAST alignment.

Difference Analysis of RNA-Seq Data
Onaga et al. (2017) [5] provided transcriptome data on CO and LT at seedling stage, which were inoculated with M. oryzae isolate, TAN211.16, after preheating 7 days under 35 • C (this treatment was namely 35 • C + M. oryzae). The two cultivars both showed stronger resistance to rice blast under the above treatment in marked contrast to 28 • C + M. oryzae. There were 6454 DEGs from CO and 5666 DEGs from LT under 35 • C + M. oryzae in supplement files.

Integration Analysis of Meta-QTLs and RNA-Seq
A Venn diagram was used to obtain the common genes of meta-QTL interval genes and DEGs from RNA-seq. Via the singular enrichment analysis (SEA) of AgriGo (http: //bioinfo.cau.edu.cn/agriGO/analysis.php, accessed on 9 January 2022), Go terms of target genes were obtained, and then, based on FDR (false discovery rate) <0.05, GO enrichment analysis were conducted. The top GO items were visualized by ggplot2 and GOplot R software packages [63].

Compilation and Characterization of QTL Studies Regarding Blast Resistance in Rice
We updated the collection of blast-resistant QTLs in rice in our previous research (783 blast-resistant QTLs from 43 publications) [64]. A total of 839 blast-resistant QTLs in rice were collected from 51 publications in this study (Table S1). These studies used different parent lines, population size, marker type and mapping method. The overwhelming majority of those studies employed resistant and susceptible cultivars as parent lines. Furthermore, the type of mapping population included RIL, DH, BIL, F2 and F3. The size of the assayed population ranged from 63 to 587. In addition, the mapping methods included interval mapping (IM), composite interval mapping (CIM), single marker analysis (SMA), inclusive composite interval mapping (ICIM), analysis of variance (ANOVA), multiple interval mapping (MIM), simplified composite interval mapping (sCIM), mixed linear model (MLM), and generalized linear mode (GLM). The number of identified QTLs ranged from 1 to 83 in those studies.

Meta-Analysis Results for Blast-Resistant QTLs in Rice
Through meta-analysis of 839 blast-resistant QTLs, the consensus map contained 1706 markers, with an average distance of 0.22 Mb, and a total of 71 blast-resistant meta-QTLs were identified on 12 chromosomes (Figure 1). Of those, 67 blast-resistant meta-QTLs were projected by at least two original QTLs, which would be subsequently analyzed ( Table 1). The maximum number of original QTLs projected to individual meta-QTLs was up to 24. The interval distance of these meta-QTLs ranged from 0.04 to 3.12 Mb, with an average value of 1.00 Mb. Additionally, 42 blast-resistant meta-QTLs had an interval distance of ≤1.0 Mb, and 23 meta-QTLs had an interval distance of ≤0.5 Mb. Furthermore, the number of interval genes in 67 blast-resistant meta-QTLs ranged from 9 to 357. A total of 7054 interval genes were obtained from 67 blast-resistant meta-QTLs (Table S3). In addition, 47 cloned blast-resistant genes were found in 23 meta-QTL regions, among which Metaq11-6 contained 9 blast-resistant genes. different parent lines, population size, marker type and mapping method. The overwhelming majority of those studies employed resistant and susceptible cultivars as parent lines. Furthermore, the type of mapping population included RIL, DH, BIL, F2 and F3. The size of the assayed population ranged from 63 to 587. In addition, the mapping methods included interval mapping (IM), composite interval mapping (CIM), single marker analysis (SMA), inclusive composite interval mapping (ICIM), analysis of variance (ANOVA), multiple interval mapping (MIM), simplified composite interval mapping (sCIM), mixed linear model (MLM), and generalized linear mode (GLM). The number of identified QTLs ranged from 1 to 83 in those studies.

Meta-Analysis Results for Blast-Resistant QTLs in Rice
Through meta-analysis of 839 blast-resistant QTLs, the consensus map contained 1706 markers, with an average distance of 0.22 Mb, and a total of 71 blast-resistant meta-QTLs were identified on 12 chromosomes (Figure 1). Of those, 67 blast-resistant meta-QTLs were projected by at least two original QTLs, which would be subsequently analyzed ( Table 1). The maximum number of original QTLs projected to individual meta-QTLs was up to 24. The interval distance of these meta-QTLs ranged from 0.04 to 3.12 Mb, with an average value of 1.00 Mb. Additionally, 42 blast-resistant meta-QTLs had an interval distance of ≤ 1.0 Mb, and 23 meta-QTLs had an interval distance of ≤ 0.5 Mb. Furthermore, the number of interval genes in 67 blast-resistant meta-QTLs ranged from 9 to 357. A total of 7054 interval genes were obtained from 67 blast-resistant meta-QTLs (Table  S3). In addition, 47 cloned blast-resistant genes were found in 23 meta-QTL regions, among which Metaq11-6 contained 9 blast-resistant genes.

Compilation and Characterization of QTL Studies Regarding Heat Tolerance in Rice
We collected 308 heat-tolerant QTLs in rice from 32 studies published from 2002 to 2021 (Table S2). The parent lines, population size, marker type and mapping method differed in these studies. They employed heat-tolerant and heat-sensitive cultivars as parent population. The size of assayed population ranged from 37 to 1027. Moreover, the mapping methods included IM, CIM, SMA, ICIM, etc. in the previous works. The identified heat-resistant QTLs in those studies ranged from 1-53.

Meta-Analysis Results for Heat-Tolerant QTLs in Rice
With the meta-analysis of 308 heat-tolerant QTLs in rice, 43 heat-tolerant meta-QTLs were detected on 12 chromosomes which contained 1385 markers, with an average distance of 0.27 Mb, (Figure 2). The number of initial QTLs projected to meta-QTLs was 2-7. Among these, 40 heat-tolerant meta-QTLs were derived from at least two initial QTLs ( Table 2). The interval distances of these meta-QTLs ranged from 0.03 to 7.17 Mb, with an average value of 1.49 Mb. In addition, ten heat-tolerant meta-QTLs had an interval distance of ≤0.5 Mb, and 20 meta-QTLs had an interval distance of ≤1.0 Mb. Furthermore, the number of interval genes in 40 heat-tolerant meta-QTLs ranged from 6 to 887. A total of 6425 genes were obtained in these meta-QTLs (Table S4). In addition, 21 heat-tolerant genes were found in 10 meta-QTL regions, among which Metah11-1 contained eight heat-tolerant genes.

Integration Analysis Results for Meta-QTLs and RNA-Seq
By integrating interval genes from meta-QTLs regarding the two traits, we identified 1058 common genes (Figure 3a and Table S5). There were 6454 DEGs (  [5]. By integrating those DEGs, we identified 118 common DEGs (Table S8), which were located on blast-resistant and heat-tolerant meta-QTL regions (Figure 3b). In addition, 118 common DEGs were involved in 14 terms of biological processes (blue column), 4 terms of molecular function (red column), and 10 terms of cellular component (green column) (Figure 3c). Furthermore, assigned to GO enrichment analysis with FDR value < 0.05, the top 5 GO term enrichments were involved in 5 pathways of molecular function, containing 55 genes (Figure 3d and Table S9). These 55 genes were co-located on the common regions of 16 blast-resistant meta-QTLs and 14 heat-tolerant meta-QTLs in 9 chromosomes (Table S10)

Discussion
In recent years, both HT and rice blast have damaged the growth and development of rice, causing serious losses in rice production. Therefore, breeding new rice varieties with blast resistance and heat tolerance has become an increasingly urgent research task. Over the past few decades, advances in molecular genetics have led to the identification and utilization of QTLs related to yield, abiotic and biotic stress resistance [120][121][122][123][124]. However, the genetic inconsistency of these QTLs hinders their application in MAS breeding [60,124]. Meta-analysis can overcome the limitations of a single study and identify "consistent" QTLs from previous research in different genetic backgrounds and environments. It has been applied in rice [49,60], maize [53] and soybean [51,54]. In addition, compared with genetic maps with lower molecular marker density, physical maps can cover almost all markers from previous genetic maps to the greatest extent. Courtois et al. (2009) [60] used a rice physical map as reference map to identify drought-tolerant meta-QTLs through meta-analysis. Although the genetic map is not consistent with the physical map, the low recombination rate of rice mainly affects the region around the centromere, and

Discussion
In recent years, both HT and rice blast have damaged the growth and development of rice, causing serious losses in rice production. Therefore, breeding new rice varieties with blast resistance and heat tolerance has become an increasingly urgent research task. Over the past few decades, advances in molecular genetics have led to the identification and utilization of QTLs related to yield, abiotic and biotic stress resistance [120][121][122][123][124]. However, the genetic inconsistency of these QTLs hinders their application in MAS breeding [60,124]. Meta-analysis can overcome the limitations of a single study and identify "consistent" QTLs from previous research in different genetic backgrounds and environments. It has been applied in rice [49,60], maize [53] and soybean [51,54]. In addition, compared with genetic maps with lower molecular marker density, physical maps can cover almost all markers from previous genetic maps to the greatest extent. Courtois et al. (2009) [60] used a rice physical map as reference map to identify drought-tolerant meta-QTLs through meta-analysis. Although the genetic map is not consistent with the physical map, the low recombination rate of rice mainly affects the region around the centromere, and the inactive region of rice recombination is limited to a small interval [125]. Therefore, the physical map of rice was used as the reference map in this study. Except for some QTLs with lost data, most QTLs can be anchored to the physical map to minimize the loss of data and increase the number of QTLs available for meta-analysis, which would help improve the accuracy of meta-analysis. In addition, Ballini et al. (2008) [124] conducted a meta-analysis of 347 blastresistance QTLs of rice in 18 papers, and the results demonstrated that the average interval distance of meta-QTLs was 3.3 Mb, and the average number of original QTLs mapped to individual meta-QTLs was around 1.9. In this study, we conducted a meta-analysis of blast-resistance 839 QTLs from 51 studies in the literature published from 1994 to 2021. The results indicated that the average interval distance of meta-QTLs was 1.00 Mb, and individual meta-QTLs were averagely obtained from 4.8 original QTLs. Likewise, compared with Raza et al. (2020) [126], the average interval distance of heat-tolerance meta-QTLs was smaller in this study. Besides, 47 blast-resistant and 21 heat-tolerant genes published were found in meta-QTL regions (Tables 1 and 2), such as Pi27(t) [65], Pi-Da(t) [66], Pi14 [67], Pi16(t) [68], Pid1(t) [68], OsHsfA7 [104], OsHsfC1a [104], OsWRKY11 [105], OsGR2 [106], OsRb1 [107] and OsHSP58.7 [112]. Therefore, this study not only narrowed the confidence interval, but also improved the credibility of meta-QTLs related to these two traits, which would lay a foundation for further mining of crucial resistant genes.
Previous studies have conducted integrating analysis of RNA-seq and meta-QTLs to screen out key candidate genes for target traits such as cold tolerance in rice [56], seed storage components of soybean [51] and veraison time in grapevine [57]. Although confidence intervals of the above meta-QTLs were further narrowed, the meta-QTLs still contained a large number of genes in this study. Onaga et al. (2017) [5] found that rice varieties CO and LT displayed complete resistance to blast with 35 • C preheating treatment in marked contrast to 28 • C, demonstrating that HT preheating could enhance rice defense against blast fungus. In our study, the RNA-seq data was downloaded from Onaga et al. (2017) and was analyzed with blast-resistant and heat-tolerant meta-QTLs. Totally, we identified 118 common genes, which not only were positioned on meta-QTL regions related to these two traits, but also responded to rice blast and heat stresses. Furthermore, based on GO term enrichment, 55 significant candidate genes were selected (Table S10). Among them, 34 genes were up-regulated, while 21 genes were down-regulated under 35 • C and blast stresses.
To identify blast-resistant and heat-tolerant genes, we preferred to choose genes that were published to be involved in blast resistance or heat tolerance in the previous studies. Consequently, 24 significant BR-HT candidate genes were selected in Table 3. To withstand various stresses, plants have evolved a battery of complicated immune system to protect themselves from various stresses. Especially, PTI (PAMPs triggered immunity) and ETI (effector triggered immunity) play vital roles in defensing against pathogens. Among defense mechanisms, defense response genes (such as chitinases), and resistance (R)-genes are most significant. Five genes, including C10923, OsChib3H-c, C10150, Os11g0701500, and Os11g0702100, belong to glycosyl hydrolase. Chitinases are just one kind of glycosyl hydrolase. OsChib3H-c was co-located on Metab11-6 and Metah11-4 and up-regulated with 3.25-3.36 fold change in CO and LT under 35 • C + M. oryzae. OsChib3H-c had been identified as a novel chitinase gene, which could enhance resistance to sheath blight pathogen in rice [127]. Interestingly, plant hormones and abiotic stresses also regulated the expression and activity of chitinases [128][129][130][131]. OsChib3H-c was also up-regulated to respond to heat [132], drought [133] and jasmonic acid [131]. Thus, we inferred OsChib3H-c could confer tolerance to both blast fungus and heat stress. Likewise, it was found that C10150 [132,134], C10923 [132,133], Os11g0702100 [132], and Os11g0701500 [134,135] were up-regulated to respond to pathogen and abiotic stresses. According to the structure characterizes of proteins encoded by R genes [136], nine genes, including OsJAMyb, OsWRKY59, OsMYBAS1, Pi-k, WAK124, OsiWAK1, WAK3, OsRLCK198, and SDK6, belong to resistant gene analogues. OsJAMyb is an R2R3 MYB transcription factor, which was located in the same region of Metab11-6 and Metah11-4 and was up-regulated with 5.39-4.98 fold change in CO and LT under 35 • C + M. oryzae. It had been reported that OsJAMyb overexpressed in Suyu variety to enhance rice defense against blast, suggesting that OsJAMyb was involved in resistance to rice blast [98]. Pi-k [95], a resistant gene, was also up-regulated with 2.28-2.65 fold change in CO and LT under 35 • C + M. oryzae. Besides, WAK kinases (wall associated kinases) are cell wall-associated receptor kinases and had been found to be involved in pathogen resistance and abiotic stress tolerance of rice [99]. OsWAK1 (Os11g0690066) was co-located in Metab11-6 and Metah11-4 regions but not differentially expressed in rice under 35 • C + M. oryzae, but WAK124, OsiWAK1, and WAK3 were all up-regulated in CO and LT under 35 • C + M. oryzae. Previous studies had found that rice resistant genes could be induced to express by preheating [7,137]. For instance, TaRPM1 [138] under HT stress was up-regulated 6 folds higher than that under normal temperature, which actively regulated the resistance of wheat to wheat stripe rust. CaWRKY40 had a positive regulatory effect on the single stress of HT and Fusarium wilt [29]. Further studies showed that the interaction between CaWRKY40 and CaWRKY6 in pepper regulated the resistance to bacterial blight under HT [139]. Under hormone induction, TaWRKY70 [137], TaWRKY49 [140] and TaWRKY62 [140] participated in the regulation of wheat resistance to the dual stress of HT and stripe rust. Thus, OsJAMyb, Pi-k, OsWAK124, OsiWAK1, and WAK3 could be up-regulated by preheating to enhance CO and LT defense against M. oryzae.
On other hand, when plants are suffering from abiotic and biotic stresses, a change in redox state controlled by oxidordeuctases is a common outcome, due to the production and accumulation of reactive oxygen species (ROS) [141,142]. Eight genes were related to oxidordeuctases, including OsMT2b, OsTPS3, OsNDB3, OsHI-LOX, OsLOX8, ACLA3, OsACLA-2 and OsGS2. It depends on the level of ROS to determine whether it will be a defensive or destructive molecule [142,143]. Signal transduction pathways can regulate the level of ROS to protect plants from adverse effects of ROS [142]. In the present study, the metallothionein encoded by OsMT2b was co-located on Metab5-1 and Metah5-1 regions and that was down-regulated in CO and LT under 35 • C + M. oryzae. Previous study illustrated that OsMT2b was down-regulated by the small GTPase OsRac1 in rice to scavenge ROS to increase resistance to bacterial blight and blast fungus [71]. OsHI-LOX is a chloroplastlocalized type lipoxygenase gene in rice. A previous study had found that OsHI-LOX participated in insect-induced JA synthesis and enhanced resistance to BPH (brown planthopper) by scavenging BPH-induced H 2 O 2 [144]. Likewise, OsTPS3 (caryophyllene synthase) [145], OsACLA-2 (ATP-citrate lyases) [146] and OsGS2 (glutathione synthetase) [147] had been reported to negatively regulate cell death and disease resistance in rice. Thus, OsMT2b, OsHI-LOX, OsTPS3, OsACLA-2 and OsGS2 could balance the ROS level to defend against blast and HT stresses in CO and LT under 35 • C + M. oryzae.
Another one, OsUBC6, encoding one of the ubiquitin-conjugating enzymes, E2, was also up-regulated with 2.44-2.62 fold change in CO and LT under 35 • C + M. oryzae. Previous research showed that OsHTAS, encoding a ubiquitin ligase, interacted with components of the ubiquitin/26S proteasome system to enhance heat tolerance through modulation of hydrogen peroxide-induced stomatal closure [110]. Therefore, OsUBC6 was inferred to play a role in responding to blast and HT stresses through interaction with components of the ubiquitin/26S proteasome system.
Overall, integration analysis of meta-QTLs and RNA-seq provided new insight for further screening of candidate genes conferring both blast resistance and heat tolerance in rice. Of course, the genes identified in meta-QTLs should be further investigated.

Conclusions
In this study, 67 blast-resistant and 40 heat-tolerant meta-QTLs in rice were obtained from 839 and 308 QTLs, respectively. Combined with RNA-seq and GO enrichment analysis, 24 significant genes were mined, which would be the gene resources for functional verification and rice breeding with double tolerance.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/genes13020224/s1, Table S1: Information of new publications on blast-resistant QTL mapping in rice, Table S2: Significant information of 32 previous studies related to heat-tolerance QTL mapping, Table S3: Interval genes in 67 blast-resistant meta-QTLs, Table S4. Interval genes in 40 heattolerant meta-QTLs, Table S5: Commen genes from blast-resistant meta-QTLs and heat-tolerant meta-QTLs, Table S6: Differently expressed genes (Gene ID converted from MSU to RAP type) in CO under 35 • C and blast stresses, Table S7: Differently expressed genes (Gene ID converted from MSU to RAP type) from LT under 35 • C and blast stresses, Table S8: Commen genes from meta-QTLs and DEGs in rice, Table S9: Genes in the top 5 enrichenment GO terms, Table S10