Unravelling consensus genomic regions conferring leaf rust resistance in wheat via meta-QTL analysis

Leaf rust, caused by the fungus Puccinia triticina Erikss (Pt), is a destructive disease affecting wheat and a threat to food security. Developing resistant varieties represents a useful method of disease control, and thus, understanding the genetic basis for leaf rust resistance is required. To this end, a comprehensive bibliographic search for leaf rust resistance quantitative trait loci (QTLs) was performed, and 393 QTLs were collected from 50 QTL mapping studies. Afterwards, a consensus map with a total length of 4567 cM consisting of different types of markers (SSR, DArT, Chip-based SNP markers and SNP markers from GBS) was used for QTL projection, and meta-QTL analysis was performed on 320 QTLs. A total of 75 genetic map positions (gmQTLs) were discovered and refined to 15 high confidence mQTLs (hcmQTLs). The candidate genes discovered within the hcmQTL interval were then checked for differential expression using data from three transcriptome studies, resulting in 92 differentially expressed genes (DEGs). The expression of these genes in various leaf tissues during wheat development was explored. This study provides insight into leaf rust resistance in wheat and thereby provides an avenue for developing resistant varieties by incorporating the most important hcmQTLs.

have been performed in wheat; however, the detected QTLs often do not overlap either in part or in whole due to a different combination of parental lines or studies in multiple environments (Rong et al. 2007). Therefore, there is a need to detect the most promising consensus QTL found among studies using different parents that is stable across environments.
Furthermore, recognizing robust and reliable QTLs and refining their intervals can be achieved by meta-QTL (MQTL) analysis without using expensive resources (Goffinet and Gerber, 2000).
The free software BioMercator (Arcade et al. 2004), used for MQTL analysis, allows the compilation of a vast number of genetic maps from various sources and can project QTL to a consensus or reference map (Veyrieras et al. 2007). Thus, MQTL analysis could identify the consensus QTL associated with the trait in multiple environments and genetic backgrounds (Goffinet and Gerber, 2000). QTLs for similar traits may be combined synergistically into MQTL traits by a single MQTL analysis. The reported MQTL method can be used for MAS (Yu et al. 2014;Maccaferri et al. 2015).
Several studies on MQTL-QTL analysis for disease resistance in wheat have been carried out successfully, including those for tan spot resistance (Liu et al. 2020b), Fusarium head blight resistance (Liu et al. 2009;Löffler et al. 2009;Venske et al. 2019;Zheng et al. 2020), leaf rust resistance (Soriano and Royo, 2015), and stem rust resistance (Yu et al. 2014). The first meta-QTL analysis on leaf rust focused only on identifying consensus genomic regions for the trait. In addition, the authors did not use SNP-based genetic maps due to the lack of QTL studies in leaf rust using high-throughput SNP arrays several years ago. However, in this study, we aimed to delve deeper into the genetic architecture underlying leaf rust disease by discovering putative candidate genes from the newly released wheat genome sequence (IWGSC, 2018), incorporating transcriptomic studies, and studying the function of these in genes in different tissues.

Bibliographic search for wheat leaf rust resistant QTLs
Using Google Scholar (https://scholar.google.com/) and the Web of Science (http:// www.webofknowledge.com/), an exhaustive search for publications containing QTLs conferring leaf rust resistance in wheat was performed. For each study, the information collected during the

QTL projection on the consensus map
To project the largest number of QTLs, the high-density consensus map developed by Venske et al. (2019) was used. This consensus map incorporated three marker types: SNPs, DArTs, and SSR markers. The SNPs were sourced from chip-based markers and genotype by sequencing (GBS) (Cavanagh et al. 2013;Saintenac et al. 2013;Wang et al. 2014). SSR markers, including functional markers, were provided from three genetic maps (Wheat, Consensus SSR 2004, Wheat Composite 2004, and Wheat Synthetic x Opata) obtained from the GrainGenes database (https://wheat.pw.usda.gov/GG3/node/976). Wheat consensus map version 4.0, which contains more than 100000 DArTseq markers and nearly 4000 DArT markers developed from over 100 genetic maps, was downloaded from https:/www.diversityarrays.com/technology-andresources/genetic-maps. The initial QTLs were projected following the approach described in Chardon et al. (2004) (Darvasi and Soller, 1997;Guo et al. 2006), where N is the number of genotypes in the mapping population, and PVE is the phenotypic variance explained by the QTL.

Meta-QTL analysis and validation using GWAS
The meta-QTL analysis was conducted using the software BioMercator (Arcade et al. 2004;Sosnowski et al. 2012) by incorporating two approaches for the analysis. The first approach proposed by Goffinet and Gerber (2000) is used when the QTL count for a chromosome is below 10. The second approach proposed by Veyrieras et al. (2007) is used when the QTL count for a chromosome is above 10. For the first approach, the lowest AIC value was selected as the best fit model. However, the best model was selected from the AIC, AICc, AIC3, Bayesian information criterion, and average evidence weight models from the second approach. Therefore, the model with the lowest criteria in at least three of the models was selected and regarded as the best fit model. The MQTL naming convention used by Zheng et al. (2020) was adopted for this study.
The MQTLs identified from the meta-QTL analysis were named based on their genetic map positions (gmQTLs). Afterwards, the sequences of the flanking markers of each gmQTL was BLASTed against the Chinese spring reference genome at https://wheat.pw.usda.gov/blast/ (IWGSC, 2018), and their corresponding physical positions were identified, thus resulting in a sequence-based mQTL (smQTL). In addition, the smQTLs found in this study were validated using recent association studies with different genetic backgrounds and environments aimed at identifying loci/QTLs related to leaf rust resistance in wheat.

Establishment of hcmQTLs and candidate gene mining
To further refine the smQTLs, those with at least five overlapping QTLs having a physical distance lower than 20 Mb and a genetic distance lower than 10 cM were selected and called high confidence mQTLs (hcmQTLs). The annotated reliable genes (HighConfidenceGenesv1.1) within the interval of each hcmQTL were obtained, and their functional annotations were examined at https://wheat-urgi.versailles.inra.fr/Seq-Repository/Annotations. For the DEGs discovered within the refined hcmQTL, Gene Ontology (GO) analysis was performed using the GENEDENOVO cloud platform (https://www.omicshare.com/tools/). dough: flag leaf blade (senescence); and ripening: flag leaf blade (senescence). Transcripts per million (TPM) values were used to assess the candidate genes' level of expression within the hcmQTLs displayed on the heat map using log 2 (TPM+1).

QTL compilation and projection on the consensus map
A comprehensive search for QTLs conferring resistance to leaf rust resulted in 50 articles published from 1999 to 2020 (Online Resource 1). A total of 393 QTLs widely distributed across the genome were collected (Online Resource 2).
Of the 393 QTLs found, only 320 QTLs had flanking markers and were thus used for MQTL analysis, leaving 73 QTLs with no flanking markers on the consensus map (Online Resource 2).
The QTLs were then projected onto the consensus map constructed by Venske et al. (2019). The highest number of QTLs (264) was projected on the B-genome, while the D-genome harboured the lowest number of QTLs (59) (Figure 1). For the A-genome, chromosome 2A had the highest number of projected QTLs, while chromosomes 4A and 7A equally had the lowest (8).
Chromosome 1B had the highest number of projected QTLs (41) for the B-genome, while chromosomes 4B and 5B both had the lowest (14). The overall number of QTLs projected on the D-genome was relatively low compared to other genomes, with the highest number of QTLs (22) projected on chromosome 2D, while 6D did not harbour any projected QTLs. When the trait was considered, a large proportion of the projected QTLs were for disease severity (43%), followed by AUDPC (24%) (Figure 1).

Unravelling consensus regions via MQTL analysis
For the MQTL analysis, the Veyrieras approach was used to analyse all the linkage groups except for chromosomes 1A, 4A, 5A, 7A, 3D, 4D, and 5D because they had fewer than 10 QTLs; thus, the approach of Goffinet and Gerber (2000) was used for their analysis. Overall, 75 gmQTLs were discovered and distributed across all the chromosomes ( Figure 1 and Table 1).
For the reported gmQTLs, the CI ranged from 0.03 cM to 25.23 cM, with an average of 5.36 cM.
For the B genome, a total of 33 gmQTLs were found, representing the genome with the highest number of gmQTLs. Chromosome 2B reported the highest number of gmQTLs (8), followed by 6 gmQTLs on 1B and 6B, and chromosomes 3B and 4B had the lowest number with 3 gmQTLs.
For the D genome, a total of 18 gmQTLs were found, with 1D, 2D, and 7D harbouring the most gmQTLs (4), while chromosomes 3D, 4D, and 5D harboured the fewest gmQTLs (2) . Interestingly, none of the smQTLs whose physical intervals overlapped were included in the genetic map.

Association study validation and colocalization with leaf rust genes
The loci resistant to leaf rust across different genetic backgrounds and environments identified in recent association studies (Online Resource 3) were used to validate the MQTLs discovered in this study. A total of 51 MTAs identified were colocalized with 29 MQTLs; thus, some MQTLs integrated more than one MTA (Online Resource 4, Figure 2). Furthermore, eight leaf rust genes colocalized with some smQTLs found in this study (Online Resource 5, Figure 2), such as smQTL1D.2 colocalizing with Lr60 and Lr42 and smQTL7B.3 colocalizing with the Lr68 and Lr14a genes. Additionally, smQTL1B.4 and smQTL2B.5 both colocalized with two leaf rust genes, Lr46 and Lr13, respectively, conferring adult plant resistance.

Candidate gene mining of established high confidence mQTLs (hcmQTLs)
To further improve the quality of the smQTLs discovered, they were further refined to regions termed high confidence mQTLs (hcmQTLs Owing to the lack of transcriptomic data repositories for leaf rust in wheat, we decided to use transcriptomic expression data for three fungal diseases in wheat, stripe rust, powdery mildew, and Septoria tritici blotch, to explore the DEGs within these hcmQTLs. The ERP013983  of 92 genes were found to be differentially expressed across the three expression datasets used. The 92 DEGs within the hcmQTL interval were analysed for Gene Ontology (GO) enrichment (Table 2). The most significantly enriched GO terms associated with biological processes were for metabolic (44 genes) and cellular processes (31 genes) (Online Resource 10, Figure 4). The most significantly enriched GO terms associated with molecular function were for catalytic activities (41 genes) and binding (39 genes). In terms of cellular components, the genes were enriched mainly in the cell membrane and its components.

Tissue-specific expression profile of DEGs within hcmQTLs
To understand the expression of the DEGs within the hcmQTL intervals across different tissues at different wheat development stages, we used DEGs expressed in all the transcriptomics datasets SRP041017, ERP013983, and ERP009837. Row clustering was applied, and as a result, the 92 DEGs fell into two classes (I-II) based on their expression patterns ( Figure 5). Genes in class I showed moderate to high expression in the flag leaf blade and fifth leaf blade at the anthesis stage when compared to other stages of growth. Moreover, for class II, more genes were highly expressed in the first leaf sheath at the tillering stage. The DEGs in class II accounted for 1 1 shown by genes in class I. For class I, at the seedling stage, the genes TraesCS1D02G003700 (hcmQTL1D.2) and TraesCS7B02G421100 (hcmQTL7B.1) showed moderate expression in the stem axis, while at the adult stage, TraesCS6A02G072500 (hcmQTL6A.2) was highly expressed in the fifth leaf blade at the anthesis stage, and TraesCS1D02G017700 (hcmQTL1D.2) showed high expression in the flag leaf at the dough stage. Furthermore, at the seedling stage, more class II genes were moderately expressed in the radicle and roots, with only TraesCS7D02G211000 (hcmQTL7D.1) showing expression in the shoot apical meristem. At the adult stage, the gene TraesCS7D02G220300 (hcmQTL7D.1) was highly expressed in the fifth leaf blade at the anthesis stage.

Establishment of gmQTL and smQTL regions
To gain deeper insight into the control of leaf rust resistance in wheat, a meta-QTL analysis was performed based on the numerous QTLs conferring leaf rust resistance identified in the literature from various independent studies. The first step to identifying consensus regions via meta-QTL analysis is the projection of the original QTLs onto a consensus or reference map.
A feature of the consensus map and QTL database was that the B genome reported the highest marker saturation, and thus, the highest number of QTLs was mapped to this genome, which is in agreement with previous studies characterizing genetic diversity and unravelling complex traits for disease resistance in bread wheat (Li et al. 2015;Soriano and Royo, 2015;Wang et al. 2014).
The D genome presented a lower number of QTLs, as previously found in other meta-QTL analyses for disease resistance in wheat (Soriano and Royo, 2015;Venske et al. 2019;Liu et al. 2020b;Zheng et al. 2020). Furthermore, no QTLs were found on chromosome 6D, as discovered in previous meta-QTL analysis studies on leaf rust and Fusarium head blight diseases in wheat (Soriano and Royo, 2015;Zheng et al. 2020 Soriano and Royo (2015). Consequently, the number of genomic regions (gmQTLs) discovered in this study was higher than those reported in Soriano and Royo (2015), at 75 and 48, respectively. For the gmQTLs discovered in this study, the confidence interval ranged from 0.03 cM to 25.23 cM, which was significantly reduced by more than 50% from the confidence interval of the original QTLs, ranging from 1.14 cM to 173.11 cM. In addition, in the present study, the physical position of the gmQTLs was reported due to the release of the wheat genome sequence (IWGSC, 2018), improving the mapping resolution of the genome regions (smQTLs) and helping the identification of candidate genes.
These analyses enhance the results provided by studies published prior to the release of the genome sequence. In this study, we discovered 7 smQTLs incorporating at least five original QTLs and having a confidence interval of less than 10 Mb, making them the most promising for candidate gene identification.

Colocalization of smQTLs with leaf rust resistance genes and traits
To strengthen the location of smQTLs discovered in this study, a search for colocalization of leaf rust resistance genes and smQTLs was performed. More than 61 leaf rust genes have been mapped and documented in wheat (Kim et al. 2020), and 4 of them have been cloned (Hafeez et al. 2021). A total of six leaf rust genes (Lr13, Lr14a, Lr46, Lr68, Lr63, Lr60, Lr42, and Lr41) were found to colocalize with smQTLs. Interestingly, the colocalization of Lr13, Lr14a, and Lr46 with smQTL2B.5, smQTL7B.3, and smQTL1B.4 on chromosomes 2B, 7B, and 1B, respectively, in this study was in agreement with the results obtained by Soriano and Royo (2015). As reported by these authors, Lr68 was found to have a tight association with MQTL33 (colocalized with Lr14a) on chromosome 7B; however, in this study, smQTL7B.3 colocalized with both leaf rust genes (Lr14a and Lr68), thus confirming the usefulness of using highly saturated consensus maps for meta-QTL analysis. The gene Lr14a, known to confer adult plant resistance, is thought to have evolved from emmer wheat cv. Yaroslav (McFadden, 1930) and is associated with the stem rust and powdery mildew resistance genes Sr17 and Pm5. Lr68, on the other hand, confers seedling resistance to the majority of P. triticina isolates with low to medium infection types and is linked to small but noticeable leaf tip necrosis (Herrera-Foessel et al. 2012).
Consequently, the smQTL7B.3 region not only confers seedling and adult plant resistance to leaf rust but also constitutes a region of multiple disease resistance. Additionally, smQTL1D.2 colocalized with two leaf rust resistance genes (Lr60 and Lr42). Hiebert et al. (2008) found that Lr60 is 13.5 cm distal to Lr21, which would position Lr60 and Lr42 approximately 40 cm apart (Huang et al. 2003;Somers et al. 2004). The association between Lr60 and Lr42 has not been confirmed, but in this study, we discovered that both genes were located in the same smQTL, with a confidence interval of 8.8 Mb. This supports possible linkage between the two genes; however, a genetic linkage test needs to be carried out to corroborate this claim. Furthermore, Lr60 is known to confer seedling resistance, while Lr42 confers adult plant resistance. Thus, the smQTL1D.2 region has the potential to provide more qualitative resistance against leaf rust.
Additionally, smQTL1B.4 colocalized with Lr46, a gene known to increase the latent period and reduce the frequency of infection and uredinial size in a similar manner to Lr34 (Drijepondt and Pretorius, 1989;William et al. 2003). There is also a tight linkage between Lr46 and a stripe rust gene (Yr29), which is similar to the linkage between Lr34 and Yr18 (McIntosh, 1992;Singh 1992). Consequently, the smQTL1D.2 region confers resistance to both leaf and stripe rust in wheat, thus making this region a hotspot for selecting multiple disease resistance in wheat.
Most of the smQTLs discovered in this study clustered QTLs conferring two or more resistance traits. This phenomenon was also discovered by Kolmer et al. (2018) using RILs, suggesting that they experienced higher disease severity levels and leaf rust responses. In another study, Ren et al. (2012) also discovered that maximum disease severity had a significant association with the area under the disease progress curve (AUDPC) across diverse environments, and this finding was in agreement with previous studies (Wang et al. 2005b;Liang et al. 2006;Lan et al. 2009).
Consequently, this result indicates the possibility of replacing AUDPC with MDS. A possible explanation for this could be that when two or more traits are mapped to the same region, they are most likely under the same genetic control, as suggested by Lu et al. (2017). Furthermore, effects arising from tight linkage and pleiotropism could also be a possible explanation.

Candidate genes within hcmQTLs and their role in leaf rust resistance
The search for candidate genes was extended to hcmQTLs within 20 Mb, thus yielding 15 hcmQTLs. hcmQTLs also have a small confidence interval compared to smQTLs, thus making them more reliable and useful for QTL selection in breeding programmes. The hcmQTLs 1 4 harboured a total of 2240 genes, after which 92 DEGs were narrowed down. Two main types of disease resistance are used in breeding programmes: seedling resistance and adult plant resistance. Thus, the analysis of the expression of the candidate genes across different tissues and developmental stages can inform us of their potential role in seedling or adult plant resistance.
Five out of the 92 genes expressed across the three transcriptomic data sets, TraesCS7D02G212800, TraesCS6A02G073300, TraesCS2B02G104200, TraesCS1D02G003700, and TraesCS2D02G021300, showed moderate expression in the first leaf sheath at the seedling stage. TraesCS7D02G212800 and TraesCS2B02G104200 encode a receptor-like kinase (RLK) and protein kinase family protein, respectively, and both proteins play a crucial role in contributing to disease resistance in wheat. Plant protein kinases, as well as receptor-like kinases, govern the detection and activation of diverse developmental and physiological signals, particularly those involved in defence and symbiosis (Rentel et al. 2004;AbuQamar et al. 2008;Fu et al. 2009;Garcia et al. 2012). Prior studies found that various RLK genes coding wheat leaf rust kinases (WLRKs) were conserved in wheat, with the most studied member of the WLRK family being LRK10, which is genetically linked to the Lr10 locus (Feuillet et al. 1997(Feuillet et al. , 1998(Feuillet et al. , 2001. meta-QTL analysis can be used to identify regions that confer resistance to more than one disease, and the marker information can be used for MAS (Ali et al. 2013). In this study, the hcmQTLs1B.4 region was identified to confer resistance to leaf and stripe rusts, thus making it a potential region to exploit for multiple disease resistance in wheat. Furthermore, breeding for durable resistance is desired in major breeding programmes. Durable resistance remains effective against a pathogen for a significant number of years (Johnson, 1981;Johnson, 1984). The combination of seedling resistance and adult plant resistance has been proven to confer prolonged resistance over several years (Kolmer and Oelke, 2006). In addition, various studies have ascribed durable leaf resistance to adult plant resistance rather than to seedling resistance (Figlan et al. 2020). Therefore, hcmQTL1D.2 discovered in this study can be harnessed to confer durable resistance in wheat, as it incorporates genes conferring both seedling and adult plant resistance. Another useful approach that could be harnessed in breeding for leaf rust resistance in wheat is gene pyramiding. Gene pyramiding involves incorporating multiple desired genes into a single variety. Gene pyramiding is broadly acknowledged by breeders, plant pathologists and farmers to improve disease resistance in wheat (Chen and Kang, 2017). A major requirement for gene pyramiding is to identify various QTLs or genes conferring resistance and then incorporate them into a high-yielding cultivar (Singh 1992). In several instances, this technique has been utilized in crops. For instance, long-term resistance was conferred when diverse genes were pyramided with leaf rust genes (Kolmer, 1996;Bhawar et al. 2011;Aboukhaddour et al. 2020;Babu et al. 2020). Additionally, in barley, MAS combined with gene pyramiding has been used to introgress resistance loci against stripe rust into numerous lines (Toojinda et al. 1998(Toojinda et al. , 2000Castro et al. 2003a, b;Richardson et al. 2006). To this end, the hcmQTLs discovered in this study can be utilized and exploited for gene pyramiding via MAS to bolster the resistance of wheat against leaf rust.

Concluding remarks
One of the most effective methods for analysing the wealth of QTL information available from various studies is MQTL analysis. In this study, we delineated the genetic architecture of leaf rust in wheat via MQTL analysis and by integrating genomics studies. In comparison with initial QTL reports, meta-analysis allowed us to reduce the MQTL CI, thereby facilitating the search for candidate resistance genes in the databases available. The result was the discovery of 15 hcmQTLs, with each having a potential role in MAS. This result will be useful for developing resistance to leaf rust through the introgression of desirable hcmQTLs that could confer a high level of resistance during cultivar development. Last, this study can also help better define the various mechanisms associated with leaf rust resistance in wheat.