Ferjaoui et al. BMC Genomics (2022) 23:372 https://doi.org/10.1186/s12864-022-08560-2 RESEARCH Open Access Deciphering resistance to Zymoseptoria tritici in the Tunisian durum wheat landrace accession ‘Agili39’ Sahbi Ferjaoui1,2†, Lamia Aouini3,4,5†, Rim B. Slimane1,6, Karim Ammar7, Suzanne Dreisigacker7, Henk J. Schouten8, Suraj Sapkota9,10, Bochra A. Bahri1,9, Sarrah Ben M’Barek11, Richard G. F. Visser8, Gert H. J. Kema3,12 and Sonia Hamza1* Abstract Background: Septoria tritici blotch (STB), caused by Zymoseptoria tritici (Z. tritici), is an important biotic threat to durum wheat in the entire Mediterranean Basin. Although most durum wheat cultivars are susceptible to Z. tritici, research in STB resistance in durum wheat has been limited. Results: In our study, we have identified resistance to a wide array of Z. tritici isolates in the Tunisian durum wheat landrace accession ‘Agili39’. Subsequently, a recombinant inbred population was developed and tested under green- house conditions at the seedling stage with eight Z. tritici isolates and for five years under field conditions with three Z. tritici isolates. Mapping of quantitative trait loci (QTL) resulted in the identification of two major QTL on chromo- some 2B designated as Qstb2B_1 and Qstb2B_2. The Qstb2B_1 QTL was mapped at the seedling and the adult plant stage (highest LOD 33.9, explained variance 61.6%), conferring an effective resistance against five Z. tritici isolates. The Qstb2B_2 conferred adult plant resistance (highest LOD 32.9, explained variance 42%) and has been effective at the field trials against two Z. tritici isolates. The physical positions of the flanking markers linked to Qstb2B_1 and Qstb2B_2 indicate that these two QTL are 5 Mb apart. In addition, we identified two minor QTL on chromosomes 1A (Qstb1A) and chromosome 7A (Qstb7A) (highest LODs 4.6 and 4.0, and explained variances of 16% and 9%, respectively) that were specific to three and one Z. tritici isolates, respectively. All identified QTL were derived from the landrace acces- sion Agili39 that represents a valuable source for STB resistance in durum wheat. Conclusion: This study demonstrates that Z. tritici resistance in the ‘Agili39’ landrace accession is controlled by two minor and two major QTL acting in an additive mode. We also provide evidence that the broad efficacy of the resist- ance to STB in ‘Agili 39’ is due to a natural pyramiding of these QTL. A sustainable use of this Z. tritici resistance source and a positive selection of the linked markers to the identified QTL will greatly support effective breeding for Z. tritici resistance in durum wheat. Keywords: Durum wheat, Landraces, Zymoseptoria tritici, Mycosphaerella graminicola quantitative trait loci (QTL), Gene efficacy, QTL epistasis, Gene pyramiding †Sahbi Ferjaoui and Lamia Aouini contributed equally. Background *Correspondence: hamza.sonia@inat.agrinet.tn Wheat has been, for centuries, the prime food and feed 1 Laboratory of Bioaggressors and Integrated Protection in Agriculture (BPIA), crop especially in the Mediterranean basin [1]. This sta- National Institute of Agronomy of Tunisia (INAT), 43 Avenue Charles Nicolle, ple crop supplies 20% of the human calorie intake, and is 1082 El Mahrajène, Tunis, Tunisia thereby a major component for global food security [2, 3]. Full list of author information is available at the end of the article © The Author(s) 2022. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/. The Creative Commons Public Domain Dedication waiver (http:// creat iveco mmons. org/ publi cdoma in/ zero/1. 0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data. Ferjaoui et al. BMC Genomics (2022) 23:372 Page 2 of 20 The genus Triticum L. comprises several wheat species of re-emerging and upcoming threats into a manage- with various ploidy levels, but global wheat production is able problem [21, 22], such as the stem rust caused by almost entirely based on bread wheat, T. aestivum L. em. the Ug99 strain [23–25]. Before modern plant breed- Thell. (2n = 6x = 42, sub-genomes AABBDD), and durum ing, improved crops frequently resulted from farmers’ wheat, T. turgidum L. var. durum (2n = 4x = 28, sub- selections of outperforming genotypes in terms of yield genomes AABB), also known as pasta wheat [4]. Durum stability. Often, such so-called landraces contained a vari- wheat accounts for about 8% to the global wheat produc- ety of closely related lines that quenched biotic threats. tion, and its cultivation is concentrated in latitudes rang- During the onset of breeding, these landraces were often ing from 55°N to 40°S [5, 6], corresponding mostly to the the starting material for targeted efforts to improve for Mediterranean Basin, the North American Great Plains, instance disease resistance [26–28]. India and the former USSR [6]. Durum wheat is also Several studies have revealed that durum wheat lan- produced in sub-Saharan Africa (SSA) where Ethiopia is draces are a valuable source of resistance alleles against a leading country and considered as one of the biggest fungal pathogens [8, 12, 29–31]. STB resistance sources durum wheat producers with approximately 0.6 million on durum wheat were identified in many countries such ha [7], and a center of diversity for tetraploid wheat [8]. as Tunisia [12], Ethiopia [8], Iran [32] and Spain [33]. Northern Africa has been also the cradle of wheat pro- Until now, up to 22 septoria tritici blotch (Stb) resistance duction for centuries and was the breadbasket for the genes have been identified and mapped [34, 35]. How- Roman Empire [9, 10] with locations such as Dougga in ever, due to the apparent dichotomy in natural Z. tritici Tunisia, as exquisite trading zones for wheat and other populations for either bread wheat or durum wheat [36– commodities until the late 500’s AD [11]. In Tunisia, 38], the presence of these mapped Stb genes in durum durum wheat occupies 725 Mha approximates represent- wheat cannot be determined using well characterized Z. ing 49% of the total annual cereal area [12], with an aver- tritici strains originating from bread wheat. Thus far, the age yield estimated at 1.7 tons per hectare between the substantial research progress is mainly based on the Z. cropping seasons 2014/2015 and 2019/2020 [13]. tritici – bread wheat pathosystem [34, 39]. Therefore, and Alike bread wheat, durum wheat production is sig- albeit the growing efforts to dissect the durum wheat – Z. nificantly affected by abiotic stress conditions—mostly tritici interactions [12, 33, 40, 41], resistance breeding to drought—and by the emergence of more aggressive path- Z. tritici in durum wheat has slowly progressed over the ogens [14]. Throughout Maghreb region, the foliar blight last 25 years compared to bread wheat [34]. This affects septoria tritici blotch (STB), caused by the hemibiotroph many small growers who produce this wheat as a staple Zymoseptoria tritici (Desm.) Quaedvlieg & Crous (for- crop in an area that is severely struck by septoria tritici merly Mycosphaerella graminicola (Fuckel) J. Schröt. in blotch. Cohn), is among the major threats [15]. Estimated yield In this study, we have embarked on increasing the losses amounted up to 385 kg.ha−1 in 2008–2009, which understanding of the Z. tritici—durum wheat patho- is more than 30% in most regions [16]. Recent research system. Here, we have first screened eight durum wheat increased the general understanding of the Z. tritici epi- Tunisian landrace accessions for Z. tritici resistance. Sub- demiology in the Maghreb. Hamada [17] reported the sequently, we developed a mapping population between occurrence of the teleomorph of the fungus in Tunisia, the resistant landrace accession ‘Agili39’ and the suscep- despite the arid conditions in the region, and Neddaf tible modern cv. Khiar to identify the genetic basis of et al. [18] determined an equal distribution of both mat- resistance to Z. tritici and to map the underlying genes ing types in Algeria, indicating regular sexual reproduc- under greenhouse and field conditions. tion, which likely contributes to the vast genetic diversity in this region. Similar results have recently been reported Results in Tunisia on durum wheat [19]. The use of fungicides Phenotyping of RILs, landrace accessions and modern has been more slowly adopted by durum wheat growers cultivars at the seedling stage as compared to bread wheat producers in Europe, and The 20 Z. tritici isolates grew successfully under labo- the first occurrence of strobilurin resistance have been ratory conditions enabling appropriate inoculum pro- reported in Tunisia and Algeria [18, 20]. duction and phenotyping. None of the tested durum One of the best management strategies for all plant landraces and cultivars was resistant to the entire suite diseases is the generation of new disease resistant germ- of Z. tritici isolates (Table 1), but the landrace accessions plasm through plant breeding. The huge genetic diversity showed a broader efficacy compared to the cvs. Khiar and in wheat and its ancestors has provided new varieties for Karim, resulting in a significant ‘line x isolate’ interac- almost a century [4]. Releasing new resistant germplasm tion, indicating specific gene action (Additional Table 1). has proven its efficacy and has turned the potential havoc Interestingly, necrosis values were high and ranged F erjaoui et al. BMC Genomics (2022) 23:372 Page 3 of 20 between 72 and 97% (data not shown). The parents of the least population mean necrosis coverage (34.5%) was the developed RILs, ‘Agili39’ and cv. Khiar showed highly observed for RILs inoculated with Tun1 isolate. Pycnidia significantly different pycnidia values of 6% and 36%, population mean scores were also variable and were respectively (Table  1), and henceforward, we selected a high for lines inoculated with Tun6 isolate (23.9%), but set of eight Z. tritici isolates that discriminated between relatively low for lines inoculated with IPO95052 iso- ‘Agili39’ and cv. Khiar for subsequent phenotyping of the late (9.4%). For all tested isolates, necrosis scores ranged developed F6 RILs population (Table 2). between 0 and 100%, however a maximum of 90% of The seedling screening of the RILs with the selected pycnidia score was registered for lines inoculated with eight Z. tritici isolates resulted in non-symmetric fre- IPO95052 isolate (Table 4). All seedling phenotypic traits quency distributions skewed towards the resistance phe- were repeatable with the highest repeatability registered notype for all tested isolates (Fig.  1 panel A, Additional for isolates IPO92003 and Tun6 for necrosis and pyc- Fig. 1). Subsequent analyses of variance of the split-plot nidia scores, respectively and with an equal repeatability design seedling experiment revealed that the ‘RIL’ term of 0.96 for both seedling phenotypic traits and isolates was highly significant for necrosis and pycnidia AUDPC (Table 4). scores at p = 0.0001 (Table  3). This result indicates that the observed variation in the data is accounted for the variable genetic make-up of the tested lines. Phenotyping of RILs at the adult plant stage Differentiation between the isolates was observed for During all field trials, STB developed well after the necrosis and pycnidia scores in the RIL population tested inoculations, but only pycnidia coverage was assessed. at the seedling stage (Table  4). The highest population The susceptible parent cv. Khiar showed high dis- mean necrosis coverage (60.60%) was registered for RILs ease severities throughout the trials, whereas ‘Agili39’ inoculated with IPO92003 isolate (Table  4). However, remained free of disease (0 pycnidia) (Fig. 1 panel B). Table 1 Percentage of pycnidia on the sprayed-inoculated primary leaves of durum wheat landraces and cultivars with 20 Zymoseptoria tritici isolates at 21 days-post inoculations. Coloured cells indicate least significant differences (LSDs; P = 0.05) with resistant in green (not significantly different from 0% Pycnidia), intermediate in yellow (significantly different from 0% Pycnidia and 100% Pycnidia) and susceptible in red (not significantly different from 100% Pycnidia) Ferjaoui et al. BMC Genomics (2022) 23:372 Page 4 of 20 Table 2 Origin of 20 Zymoseptoria tritici isolates that were isolated from durum wheat in the Mediterranean Basin and that were used for phenotyping in the seedling and adult plant stage Experiment Region Isolate ID Country Location Year 1 2 3 Middle-East IPO91004 Syria Lattakia 1991 + + IPO95002 Syria Lattakia 1995 + IPO95003 Syria Lattakia 1995 + North Africa IPO91009 Tunisia Bejá 1991 + + IIB-123 Tunisia Bejá 2005 + + + Tun1 Tunisia Qued bagrat - + + + Tun6 Tunisia Sidi Nsir - + + + IPO91018 Morocco Jenica Shaim 1991 + + IPO95052 Algeria Berrahal 1995 + + Europe IPO92003 Portugal - 1992 + + IPO13001 Italy Emilia Romagna 2013 + IPO13003 Italy Emilia Romagna 2013 + IPO13006 Italy Emilia Romagna 2013 + IPO13007 Italy Emilia Romagna 2013 + IPO13008 Italy Emilia Romagna 2013 + IPO13018 Italy Sicily 2013 + IPO13019 Italy Sicily 2013 + IPO13023 Italy Sicily 2013 + IPO13024 Italy Sicily 2013 + IPO13056 Italy Tuscany 2013 + Experiment 1 = Pre-screening of the “Agili39” landrace and the cv. Khiar with 20 Zymoseptoria tritici isolates under controlled conditions Experiment 2 = Screening of the F6 “Agili39”/Khiar recombinant inbred lines with 8 Zymoseptoria tritici isolates under controlled conditions Experiment 3 = Screening of the F6 “Agili39”/Khiar recombinant Inbred population with 3 Zymoseptoria tritici isolates under field conditions The three-way analysis of variance revealed that the coverage was relatively low at the F10 generation inocu- ‘genotype’, ‘isolate’ and their interaction ‘genotype lated with the IIB-123 isolate (6.5% of pycnidia), with a x isolate’ terms are highly significant at p = 0.0001 maximum pycnidia coverage of 19.4% (Table 4). (Table  5). This result indicates a ‘genotype x isolate’ Adult pycnidia coverage heritability (H2) was high for specificity at the adult plant. The term ‘year’ was sig- the Tun1 and the Tun6 isolates, with a higher heritability nificant at p = 0.05. However, the term ‘block’, the two- for the adult pycnidia coverage caused by the Tun6 iso- way interaction terms ‘genotype x block’ and ‘genotype late (H2 = 0.98) compared to the heritability of the pyc- x year’, and the three-way interaction term ‘genotype x nidia coverage caused by the Tun1 isolate (H2 = 0.88) year x isolate’ were not significant, indicating no varia- (Table 4). Field data generated by the inoculation of the tion in the micro-environment and the homogeneity of F10 RIL with the IIB123 isolate were also repeatable the field inoculation (Table 5). (0.97) (Table 4). Overall, adult plants F9 showed the highest average of pycnidia coverage compared to the F7 inoculated with Correlations between the seedling and the adult plant the Tun1 isolate (21.5% of pycnidia for the F9 compared assays to 15.4% for the F7), and to the F6 and the F8 generations Low to high correlations were observed between the inoculated with the Tun6 isolate (29.9% of pycnidia for different traits (Fig. 2, Additional Table 2). Phenotypic the F9 compared to 19.6% and 19.0% for the F6 and the scores obtained by the IPO92003 isolate were the least F8, respectively) (Table  4, Fig.  1 panel B). The field dis- correlated with all phenotypic scores obtained by the ease severity on the F9 generation inoculated with Tun1 tested isolates at the seedling and the adult plant stages and Tun6 isolates was also variable with a higher sever- (Fig.  2, Additional Table  2). Nonetheless, necrosis and ity for Tun6 isolate (29.9%) compared to the Tun1 isolate pycnidia AUDPC scores generated by the IPO92003 (21.5%), which indicates once more a specific ‘genotype isolate on the tested RILs at the seedling stage were x isolate’ interaction. Interestingly, average pycnidia highly correlated (r = 0.6). Z. tritici isolates IPO95052, F erjaoui et al. BMC Genomics (2022) 23:372 Page 5 of 20 Fig. 1 Frequency distributions of the disease severity assessed as pycnidia percentage in seedlings and adult plants of the F6-F10 recombinant inbred lines of Agili39/Khiar with three Zymoseptoria tritici under field and controlled conditions. ‘A’ and ‘K’ are referring to ‘Agili39’ and cv. Khiar parents, respectively Table 3 Analysis of variance of necrosis and pycnidia AUDPC (N-AUDPC and P-AUDPC) seedling data from the ‘Agili39’/Khiar recombinant inbred lines (RILs) tested following a split-plot design Source of variation dfa Mean of Square F value Pr(> F) N-AUDPC P-AUDPC N-AUDPC P-AUDPC N-AUDPC P-AUDPC Replicate 2 5,303,251 834,974 2.2345 2.8532 0.14966 0.09689 Isolate 7 2,092,146 167,005 0.8815 0.5707 0.54841 0.7668 Error a (Main plot) 12 2,373,376 292,643 RIL 160 217,971 125,635 12.3639 13.2979 < 2.00E-16*** < 2.00E-16*** Isolate x RIL 1099 16,572 8909 0.94 0.943 0.87956 0.86706 Error b (Sub-plot) 2183 17,630 9448 Signifiance codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘’ 1 a Degrees of freedom Ferjaoui et al. BMC Genomics (2022) 23:372 Page 6 of 20 Table 4 Mean, range (minimum and maximum) and repeatability / heritability (H2) of Necrosis and Pycnidia values of the durum parents ‘Agili39’ and cv. Khiar and their derived recombinant inbred lines at the seedling and the adult plant stages Isolate Parent Mean (± Stda) RILs Population mean (± Stda) Range (Min–Max) Repeatability / H2b Physiological stage Trait Necrosis % Pycnidia % Trait Necrosis % Pycnidia % Necrosis % Pycnidia % Necrosis Pycnidia Tun1 Agili39 0.0 ± 0.5 0.0 ± 0.8 Adult-F7 15.40 ± 1.7 0.00- 51.60 0.88 b Khiar 70.0 ± 5.0 46.7 ± 3.5 Adult-F9 21.50 ± 2.2 0.00—82.40 Seedling 34.50 ± 2.2 12.80 ± 3.2 0.00—90.00 0.00—70.00 0.85 0.82 Tun6 Agili39 6.7 ± 4.2 0.0 ± 1.8 Adult-F6 19.60 ± 2.1 0.0—55.3 0.98 b Adult-F8 19.00 ± 1.9 0.00- 60.30 Khiar 90.0 ± 0.5 0.0 ± 2.5 Adult-F9 29.90 ± 3.8 0.00—87.30 Seedling 35.30 ± 3.0 23.90 ± 3.2 0.00—100.00 0.00—93.30 0.93 0.96 IIB123 Agili39 1.7 ± 0.2 0.0 ± 3.5 Adullt-F10 6.50 ± 0.94 0.00—19.40 0.97 Khiar 90.0 ± 5.0 70.0 ± 2.6 Seedling 42.10 ± 2.9 14.00 ± 2.2 0.00—100.00 0.00 -100.00 0.92 0.91 IPO91004 Agili39 23.3 ± 0.6 0.0 ± 2.3 Seedling 40.50 ± 2.7 17.00 ± 2.5 0.00—90.00 0.00—86.70 0.91 0.93 Khiar 80.0 ± 1.5 80.0 ± 5.3 IPO91009 Agili39 1.7 ± 0.5 0.0 ± 3.0 Seedling 30.50 ± 2.4 9.70 ± 1.9 0.00—93.30 0.00—80.00 0.88 0.85 Khiar 66.7 ± 3.5 56.7 ± 1.5 IPO91018 Agili39 11.7 ± 5.3 0.0 ± 0.2 Seedling 41.80 ± 2.2 10.80 ± 1.8 0.00—90.00 0.00—76.70 0.92 0.89 Khiar 80.0 ± 3.2 70.0 ± 3.5 IPO92003 Agili39 66.7 ± 5.3 13.3 ± 0.2 Seedling 60.60 ± 1.5 18.40 ± 1.1 15.00—100.00 0.00—63.30 0.96 0.87 Khiar 50.0 ± 3.4 5.0 ± 1.8 IPO95052 Agili39 5.3 ± 2.6 0.0 ± 1.5 Seedling 36.40 ± 3.1 9.40 ± 2.1 0.00—100.00 0.00—90.00 0.92 0.89 Khiar 70.0 ± 1.6 80.0 ± 2.5 a Standard deviation b Broad sense heritability (Fig.  2, Additional Table  2). The highest correlation Table 5 Analysis of variance for adult plant disease severity coefficient (r) of 0.9 was registered between pycnidia scores in the F6-F9 ‘Agili39’/Khiar recombinant inbred lines (RILs) AUDPC scores generated on lines inoculated by the that were inoculated with Zymoseptoria tritici isolates Tun1 and IPO91004 and the II-B123 Z. tritici isolates at the seed- Tun6 ling stage. In contrast, necrosis and pycnidia scores Source of Variance Dfa Mean of Square F value Pr(> F) were rather moderately correlated between IPO91018, IIB-123, IPO95052 and IP91004 isolates with a maxi- Genotype 174 4283 9.439 < 2e-16 *** mum r of 0.5 (Additional Table 2). Field data generated Isolate 1 20,847 45.94 2.18e-09 *** by Tun6 isolate were highly correlated across the F6 – Year 3 1721 3.793 0.013581 * F9 generation (r = 0.8) (Additional Table  2). Interest- Block 4 300 0.661 0.620963 ingly, positive correlations were also observed between Genotype x Block 51 262 0.578 0.980775 adult and seedling disease scores induced by Tun6 iso- Genotype x Isolate 167 955 2.105 0.000165 *** late (Additional Table 2). Similarly, positive correlations Genotype x Year 465 162 0.356 1 were registered between field scores generated by Tun6, Genotype x Year x 653 365 0.89 0.80147 Isolate and II-B123 isolates, with a correlation coefficient Residuals 77 454 r = 0.4 (Fig. 2, Additional Table 2). Signifiance codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘’ 1 a Degrees of freedom The ‘Agili39’/ Khiar linkage genetic map and the identification of QTL for Z. tritici resistance The ‘Agili39’/Khiar genetic linkage map consisted of IPO91009, Tun6, IIB-123 and IPO91004 tested at the 1959 SNP markers assigned to 30 linkage groups repre- seedling stage were correlated between each other sentative of the 14 durum wheat chromosomes (Addi- for necrosis and pycnidia AUDPC scores (0.3 < r > 0.5) tional Table 3). The total length of the genetic map was F erjaoui et al. BMC Genomics (2022) 23:372 Page 7 of 20 Fig. 2 Heatmap correlation between the Zymoseptoria tritici isolates tested under field and controlled conditions on the ‘Agili39’/Khiar mapping population performed using the “Corr” function and visualized using the “corrplot” package in the R environment [42, 43]. Different traits are named by the isolate followed by the Necrosis (N) or pycnidia (P) development at the seedling stage (S). Isolates tested at the adult plant stage are followed by the term field, and the RILs advanced generation (F6 – F10). Colours gradient is an indication of high (red colour) and low (blue colour) correlation between traits 4220  cM. In average, linkage groups consisted of 80.83 IPO91009, IPO95052 and IPO91004, but not for isolates SNP markers spanning an average length of 140.7 cM, in Tun1, IPO91018 and IPO92003—and the adult plant which adjacent loci were in average 2.3 cM distant (Addi- stage, where it provided resistance to Z. tritici isolate tional Table 3). Linkage group sizes ranged from 28.7 cM Tun6, but not to Tun1 and IIB-123 isolates. Qstb2B_1 and 338.0  cM, corresponding to linkage groups 12 and was mapped between the two SNP markers Tag_1056626 7, respectively, both representative of chromosome 1B and Tag_111757 (Additional Table  4) with a confidence (Additional Table 3). interval ranging between 69.5 and 75.5 cM and explained For the QTL analysis, the permutation test used to up to 61.6% of the phenotypic variance at the adult plant define the significant threshold LOD resulted into a LOD stage for Tun6 isolate, and up to 38.0% for IPO91004 value of 3.5, hence only QTL with a LOD ≥ 3.5 were con- isolate at the seedling stage (Table 6). LOD values of the sidered, which excluded all detected QTL with isolate Qstb2B_1 QTL were variable depending on the isolate IPO92003. In total, we identified four significant QTL on (Table  6). The highest LOD for the Qstb2B_1 was reg- three chromosomes 1A, 7A and 2B. None of these QTL istered for Tun6 at the F9 generation (33.9) at the adult was mapped with every tested Z. tritici isolate, which plant stage, and for IPO91004 isolate for the necrosis underscores specificity of the interaction between Z. score (17.3) at the seedling stage (Table 6). tritici and durum wheat. Two QTL were identified on A second QTL, designated as Qstb2B_2, was mapped chromosome 2B (Table 6; Fig. 3) designated as Qstb2B_1 on chromosome 2B. This QTL was solely mapped at the and Qstb2B_2. Qstb2B_1 was effective in both the seed- adult plant stage for isolates Tun1 and Tun6, but not for ling stage—particularly against isolates Tun6, IIB123, the IIB123 isolate. This QTL is flanked by the two SNP Ferjaoui et al. BMC Genomics (2022) 23:372 Page 8 of 20 Table 6 Detected quantitative trait loci (QTL) for Zymoseptoria tritici necrosis and pycnidia resistance in the seedling and the adult plant stages in the ‘Agili39’/Khiar recombinant inbred lines Chromosome QTL ID Isolate Trait – LOD Peak Flanking markers Confidence PVEa (%) Additif effect Physiological position Interval stage (cM) (cM) 1A Qstb1A Tun1 Necrosis-Seedling 4.6 63 Tag_1694925— 62.5—63 16.0 36.6 Tag_1127081 Tun6 Necrosis-Seedling 4.1 63 Tag_1694925— 62.5—63 9.1 50.9 Tag_1127081 Pycnidia-Seedling 4.4 63 Tag_1694925— 62.5—63 9.9 54.1 Tag_1127081 IPO91009 Necrosis-Seedling 4.1 63 Tag_1694925— 62.5—63 9.8 38.0 Tag_1127081 7As Qstb7A IIB123 Necrosis-Seedling 4.0 53 Tag_2277193— 52.5—60.5 9.0 46.2 Tag_100009953 2B Qstb2B_1 Tun6 Necrosis-Seedling 11.6 74 Tag_1056626 72.5—75.5 29.3 91.4 -Tag_111757 Pycnidia-Seedling 11.9 74 Tag_1056626— 72.5—76.5 29.6 93.7 Tag_1117572 Pycnidia-Adult (F7:8) 23.4 75 Tag_1056626— 74.5—76.5 52.9 1.5 Tag_1117572 Pycnidia—Adult 33.9 74 Tag_105662— 72.5—75.5 61.6 24.4 (F8:9) Tag_1117572 IIB123 Necrosis-Seedling 12.9 74 Tag_1056626 72.5—76.5 33.7 87.6 -Tag_1117572 Pycnidia-Seedling 5.4 74 Tag_1056626— 70.5—76.5 17.5 38.2 Tag_1117572 IPO91004 Necrsosis-Seedling 17.3 75 Tag_1056626— 72.5—76.5 38.0 96.9 Tag_1117572 Pycnidia-Seedling 7.2 74 Tag_1056626— 69.5—76.5 21.2 52.5 Tag_1117572 IPO91009 Necrosis-Seedling 9.5 74 Tag_1056626— 73.5—76.5 24.6 60.1 Tag_1117572 IPO95052 Necrsosis-Seedling 7.2 75 Tag_1056626— 71.5—76.5 21.1 77.2 Tag_1117572 Pycnidia-Seedling 3.2 74 Tag_1056626— 72.5—76.5 10.0 36.0 Tag_1117572 Qstb2B_2 Tun1 Pycnidia-Adult (F6:7) 4.8 117 Tag_100031118— 107.5—123 14.9 0.6 Tag_3027184 Pycnidia-Adult (F8:9) 5.5 113 Tag_100031118— 106.5—118.5 17.1 9.1 Tag_3027184 Tun6 Pycnidia-Adult (F5:6) 32.8 116 Tag_100031118— 114.5—118.5 54.3 1.9 Tag_3027184 a Phenotypic Variance Explained markers Tag_100031118 and Tag_3027184 (Additional showed specificity for Tun1, Tun6 and IPO91009 isolates Table  4) at a confidence interval of 106.5 and 123  cM for Qstb1A, and for IIB-123 isolate for Qstb7A (Table 6, (Table 6). The highest LOD and explained variance were Fig. 3). These QTL were uniquely detected at the seedling registered for the Tun6 isolate at the F9 generation with plant stage. Qstb1A was flanked by the two SNP markers 32.8 and 54.3%, respectively (Table  6).Comparing the Tag_1694925 and Tag_1127081 (62.5  cM—63  cM), and physical positions of the flanking markers linked to the explained up to 16.0% of the phenotypic variance with Qstb2B_1 and Qstb2B_2 QTL indicates that these two a maximum LOD of 4.6 for Tun1 isolate (Table 6, Addi- QTL are 5 Mb apart (Additional Table 4). tional Table  4). Qstb7A, mapped on chromosome 7A at Two additional QTL with lower LODs and explained 53  cM, was solely detected with IIB-123 isolate at the variances were mapped on chromosomes 1A and 7As seedling plant stage for necrosis scores (Table  6). This designated as Qstb1A and Qstb7A, respectively and QTL is flanked by the SNP markers Tag_2277193 and F erjaoui et al. BMC Genomics (2022) 23:372 Page 9 of 20 Fig. 3 ‘Agili39’/Khiar linkage groups associated with Zymospetoria tritici resistance. Panel a represents resistance QTL mapped at chromosomes 1A and 7A. Panel b represents the wide spectrum of resistance to Zymoseptoria tritici at the seedling and the adult plant stages detected in the ‘Agili39’/ Khiar population on chromosome 2B. Marker intervals linked to these QTLs are underlined and written in red. The centiMorgan (cM) distances between marker loci and the marker loci IDs are on the left and right sides of the linkage map, respectively. Different colours are associated to Zymoseptoria tritici isolates. Empty and dashed square boxes refer to seedling and adult plant stages, respectively. Trait names are marked on top of each box Tag_100009953 (Additional Table 4), with a LOD of 4.0 QTL detected for resistance to Z. tritici in the ‘Agili39’/ and explained 9% of the phenotypic variance (Table 6). Khiar population were derived from the ‘Agili39’ parent Overall, Qstb2B_1 and Qstb2B_2 had the widest effi- (Table 6). cacy compared to Qstb1A and Qstb7A (Fig.  3, Table  6). We also noticed that only Qstb1A and Qstb2B_1 QTL Identification of epistatic QTL were effective to both necrosis and pycnidia, whereas The QTL interaction analysis has disclosed an epistatic Qstb7A was only effective to necrosis (Table 6). All four QTL mapped on chromosome 5B, designated as Qstb5B Ferjaoui et al. BMC Genomics (2022) 23:372 Page 10 of 20 (Table  7). Epistatic interactions were observed between In our study, all data indicated and confirmed sig- the Qstb5B QTL and QTL Qstb1A, Qstb7A, Qstb2B_1 nificant ‘isolate’ and ‘line x isolate’ interactions as deter- and Qstb2B_2. Qstb5B was flanked between the two SNP mined in earlier studies [37, 38, 46, 47]; and recently markers Tag_1125523 and Tag_4909926 and was mapped proven in the bread wheat – Z. tritici pathosystem where at variable genetic position depending on the pairwise both Stb6 and AvrStb6 genes were cloned [48–50]. More- interaction ranging between 45 and 65  cM (Table  7). over, and comparably to other cereal diseases, namely Interestingly, the epistatic interactions were revealed to rust [51–54], we determined QTL that are detected with several isolates and at both stages (Table 7). Epistatic for either seedling (Qstb1A and Qstb7A) or adult plant LODs and explained variances were variable and ranged stage (Qstb2B_2), as well as a QTL that was detected at between 6.8—29.1, and 0.9%—54.1%, respectively. The both stages (Qstb2B_1). Our findings confirm that spe- highest pairwise LOD was detected between Qstb2B_2 cific plant physiological stage resistances are commonly and Qstb5B for the IPO95052 seedling pycnidia AUDPC observed in the Z. tritici – wheat pathosystem [34]. score (LOD = 29). However, the Tun6 adult pycnidia score Specific plant physiological stage resistances were also at F6 generation revealed the minimum epistatic LOD confirmed for other fungal diseases such as the pow- (6.8) among all detected pairwise interactions. Moreover, dery mildew and the leaf rust diseases [55, 56]. In fact, pairwise QTL interactions were mostly revealed by pyc- some Z. tritici resistance genes are uniquely effective at nidia AUDPC scores, except for the Qstb2B_1/Qstb5B the seedling stage, such as the Stb7 gene mapped in the interaction that was detected for necrosis and pycnidia spring wheat cultivar ST6 [57], or at the adult plant stage, AUDPC scores generated on the inoculated RILs by the such as the Stb17 gene [58, 59]. In contrast, other resist- Tun6 isolate at the seedling stage (Table 7). ance genes have proven to be effective at both seedling Finally, we selected the pairwise epistasis Qstb2B_1 and adult plant stages alike the Stb4 and Stb5 qualitative / Qstb5B and Qstb2B_2 / Qstb5B showing the highest genes [60, 61]. epistatic LODs for a two-way interaction test to exam- Subsequently, we compared the identified QTL to ine pycnidia AUDPC means of the variant allele combi- formerly identified Z. tritici genes using the reported nations linked to QTL Qstb2B_1, Qstb2B_2 and Qstb5B literature. This comparison has revealed that a puta- (Fig. 4). The first allele refers to the SNP markers linked tive QTL for resistance to Z. tritici was mapped on to Qstb2B_1 or Qstb2B_2, however the second allele chromosome 1A at 68  cM at the adult plant stage by refers to the Qstb5B SNP marker (Fig. 4). Kidane et  al. [8] through a genome-wide association Pycnidia AUDPC means were reduced when pairing study conducted on an Ethiopian durum wheat lan- the resistant alleles (RR) linked to Qstb2B_1 /Qstb5B drace population. Two other QTL mapped on chro- QTL and linked to Qstb2B_2/Qstb5B with pycnidia mosome 1A were also revealed by Goudemand et  al. AUDPC scores of 7.4 and 0.8 for lines inoculated with the [62] in the bread wheat Apache/Balance population, IPO91004 and IPO95052 Z. tritici isolates at the seedling and by Risser et  al. [63] named as QStb.lsa_fb-1A in stage, respectively (Fig.  4, panel A and B). Lines carry- the bread wheat bi-parental mapping population Flo- ing the susceptible allele linked to the Qstb5B QTL (RS) rett/Biscay. These QTL were mapped at the adult showed an increase in pycnidia AUDPC scores compared plant stage between 56 and 69 cM and could thus co- to lines carrying both resistant allele (RR) (Fig.  4, panel localize with the QTL mapped in the ‘Agili39’/Khiar A and B). Lines carrying solely the resistant allele linked population. However, and in contrast to the above- to Qstb5B (SR) were rather susceptible showing high mentioned studies, the Qstb1A QTL mapped in the pycnidia AUDPC scores for both pairwise interactions ‘Agili39’/Khiar population was solely detected at the Qstb2B_1/Qstb5B and Qstb2B_2/Qstb5B. This observa- seedling plant stage. tion could be explained by the minor effect of the Qstb5B The Qstb7A QTL particularly conferred reduced QTL in controlling pycnidia development, compared to necrosis values to Z. tritici isolate IIB123 and co-local- the major effect of the Qstb2B_1 and Qstb2B_2 QTL in izes with the major Stb3 gene that was mapped in the reducing pycnidia development. bread wheat cultivar Israel 493 [64, 65]. Two QTL for STB resistance were mapped on chro- Discussion mosome 2B in the ‘Agili39’/khiar population. Other Zymoseptoria tritici is a major threat to European and studies have also revealed genomic regions on chro- Mediterranean bread and durum wheat production [44]. mosome 2B associated with the STB resistance [8, 62, Despite the increasing efforts to elucidate the genetic 66–69]. basis of tetraploid wheat resistance to STB [8, 40, 45], The Qstb2B_1 QTL identified in the ‘Agili39’/Khiar more studies are required for an effective breeding strat- population likely co-localized with the known major egy in durum wheat for STB resistance. gene Stb9 that was mapped in the French bread wheat F erjaoui et al. BMC Genomics (2022) 23:372 Page 11 of 20 Table 7 Epistatic analysis between quantitative traits loci in the Agili39/Khiar bi-parental mapping population using the inclusive composite interval mapping of digenic epistatic QTL (ICIM-EPI) method in QTL the IciMapping software Isolate-Trait- Chromosome1 QTL ID Position Flanking markers Chromosome QTL ID Position Flanking markers LODe PVE Add1 g Add2 h AddbyAddi Physiological Stage a 1 (cM) b (Position 1) 2 c 2 (cM) d (Position 2) (%)f IPO95052-Pycnidia- 1A Qstb1A 60 Tag_1123459—Tag_1694925 5B Qstb5B 55 Tag_1125523— 19.1 0.9 74.7 70.5 76.6 Seedling Tag_4909926 IIB123-Pycnidia- 7As Qstb7A 55 Tag_2277193Tag_100009953 5B Qstb5B 55 Tag_1125523— 12.9 2.0 56.4 49.8 62.9 Seedling Tag_4909926 IIB123-Pycnidia-Adult Qstb7A 75 Tag_4411932—Tag_4005129 Qstb5B 65 Tag_1125523— 9.0 1.9 -0.1 0.0 0.8 (F9:10) Tag_4909926 IPO91009-Pycnidia- Qstb7A 55 Tag_2277193—Tag_100009953 Qstb5B 55 Tag_1125523— 19.5 1.1 59.4 -59.0 -56.8 Seedling Tag_4909926 IPO95052-Pycnidia- Qstb7A 80 Tag_4411932—Tag_4005129 Qstb5B 50 Tag_1125523— 19.8 0.9 74.3 73.1 75.0 Seedling Tag_4909926 Tun6-Necrosis- Seed- 2B Qstb2B_1 75 Tag_1056626—Tag_1117572 5B Qstb5B 60 Tag_1125523— 7.3 31.0 94.4 -63.2 -76.5 ling Tag_4909926 Tun6-Pycnidia- Qstb2B_1 75 Tag_1056626—Tag_1117572 Qstb5B 55 Tag_1125523— 15.8 10.6 104.2 -83.8 -69.0 Seedling Tag_4909926 IPO91004-Pycnidia- Qstb2B_1 70 Tag_100012058—Tag_1056626 Qstb5B 45 Tag_1125523— 17.1 2.9 66.8 58.4 60.2 Seedling Tag_4909926 Tun6-Pycnidia-Adult Qstb2B_2 115 Tag_100031118—Tag_3027184 Qstb5B 55 Tag_1125523— 6.8 54.6 1.9 -0.4 -0.4 (F5:6) Tag_4909926 IIB123-Pycnidia- Qstb2B_2 115 Tag_100031118—Tag_3027184 Qstb5B 55 Tag_1125523— 15.2 1.8 26.8 54.4 57.3 Seedling Tag_4909926 IPO91009-Pycnidia- Qstb2B_2 120 Tag_100031118 -Tag_3027184 Qstb5B 55 Tag_1125523— 23.0 1.1 56.5 -57.7 -60.0 Seedling Tag_4909926 IPO95052-Pycnidia- Qstb2B_2 100 Tag_1077600—Tag_100008169 Qstb5B 55 Tag_1125523— 29.1 1.0 49.6 78.7 79.7 Seedling Tag_4909926 a Chromosome at the first scanning position b Scanning position in cM of the first flanking marker pair c Chromosome at the second scanning position d Scanning position in cM of the second flanking marker pair e LOD score caused by epistatic effects f Phenotypic variation explained by epistatic effects g Estimated additive effect of position 1 h Estimated additive effect of position 2 i Additive by additive epistatic effect at the two scanning positions. A positive value indicates that the effect of the parents’ effect is larger than the recombinant effect, and a negative value means that the recombinant effect is larger than the parents’ effect Ferjaoui et al. BMC Genomics (2022) 23:372 Page 12 of 20 Fig. 4 Epistatic interactions between QTL Qstb2B_1 and Qstb5B (Panel A), and between Qstb2B_2 and Qstb5B (Panel B). Box plots illustrate the significant effects of allele variants on pycnidia AUDPC on seedlings inoculated with the IPO91004 and the IPO95052 isolates. “R” and “S” denote the resistant and the susceptible alleles at each locus, respectively. The first allele refers to the SNP markers linked to Qstb2B_1 or Qstb2B_2, however the second allele refers to the Qstb5B SNP marker. Number of genotypes and mean values of Pycnidia AUDPC are indicated under each allele classes. The black horizontal lines in the middle of the boxes are the median values for the Pycnidia AUDPC value in the respective allele classes. The vertical size of the boxes represents the inter-quantile range. The upper and lower whiskers represent the minimum and maximum values of data cv. Courtot [66]; with the qSTB.2 QTL mapped in the Thus far, in durum wheat, only partial resistance to Z. Ethiopian durum wheat landrace population [8]; with tritici was reported [16, 71]. Here, we derived Qstb2B_1 the QStb.ihar-2B.2 QTL mapped in the Liwilla/Begra from ‘Agili39’ that provides resistance to five Z. tritici bread wheat doubled-haploid population [68]; and with isolates at the seedling stage, and to two isolates at the the QStb.lfl-2B.1 mapped in the eight-founder MAGIC adult plant stage. Qstb2B_1 explains up to 61.6% of the population of winter wheat [69]. However, the Qstb2B_1 observed phenotypic variance and was characterized by QTL is different from the QTL mapped at the long arm a high heritability (0.98) with a dual action at the seed- of chromosome 2B in the Nimbus/Stigg bread wheat ling and the adult plant stages. ‘Agili39’ is also the origin mapping population [67]. of Qstb2B_2, QTL providing a major adult resistance The Qstb2B_2 QTL likely co-localized with the QTL explaining up to 54.3% of the observed phenotypic vari- identified on chromosome 2B in the mapping popula- ance. Our findings confirm that Tunisian durum lan- tions Apache/Balance and FD3/Robigus [62] associated draces harbor highly effective Z. tritici resistance QTL. with both necrosis and pycnidia resistance in the adult In fact, the initial screening of the Tunisian landrace plant stage. However, due to the unavailability of marker accessions showed a remarkable genetic diversity for STB sequences, we cannot conclude that Qstb2B_2 derived resistance, as claimed also by Ouaja et  al. [12] proving from ‘Agili39’ is the same locus that was mapped in the that Tunisian durum wheat landraces encompass diverse aforementioned bread wheat mapping populations. and valuable sources of resistance to Z. tritici. Eight lan- Thus, the identified QTL in ‘Agili39’ co-localized with drace accessions (Agili37; Agili38; Agili39, Sbei99; Der- previously mapped QTL for STB resistance in bread and bessi 12, Mahmoudi101, JK85 and Azizi27) were highly durum wheat populations, hence we cannot claim a new resistant and one landrace showed an intermediate Stb gene in the ‘Agili39’ landrace accession. However, we response (‘Agili41’). The different ‘Agili’ landrace acces- clearly have identified QTL conferring resistance to a sions reacted differently to the deployed Z. tritici isolates, wide range of Z. tritici isolates under artificial inocula- suggesting a different genetic background, which is in tion conditions in seedlings and adult plants, known as accord with Ferjaoui et al. [72] study hypothesizing that field resistance [70]. the tested ‘Agili’ accessions most likely carry different Stb F erjaoui et al. BMC Genomics (2022) 23:372 Page 13 of 20 genes combinations. Hence, and alike other durum wheat interaction between the two major QTL Qstb2B_1 and landrace populations [8], the Tunisian durum wheat lan- Qstb2B_2 mapped on chromosome 2B. draces uncover untapped allelic diversity that is of a great value to support effective breeding strategies to enhance Conclusion STB resistance in durum wheat. Our study deciphered ancient broad-based resistance to Our data demonstrated that the broad efficacy of Z. tritici in a durum wheat landrace accession. A positive the observed STB resistance in ‘Agili39’ is due to sev- selection for the QTL linked markers may result in new eral stacked QTL, both at seedling as well as adult plant high yielding durum wheat cultivars with wide resistance stages, which was also commonly observed in inherit- to Z. tritici reminiscent of the durable resistance to STB ance studies in bread wheat [34, 58]. Pyramiding genes in landraces. Given the overall high susceptibility to STB for disease resistance has been an effective strategy in in modern durum wheat cultivars, our data shed new preventing boom-and-bust cycles, and is now amenable light on disease resistance breeding in durum wheat. through marker assisted breeding as a strategy to main- tain disease resistance durability, such as for wheat stem Methods rust where various resistance gene combinations have Plant materials and study layout well controlled the disease since the mid-1950s and more Eleven durum wheat accessions (Table 8) and a bi-paren- recently to the devastating Ug99 race [56, 73–75]. A con- tal recombinant inbred (RIL) population derived from a crete illustration for Z. tritici is the effective resistance to single seed descent cross between the Tunisian landrace a wide range of isolates in the bread wheat germplasm accession ‘Agili39’ and the high yielding commercial cv. ‘KK4500’ and ‘TE11’ which is conferred by stacking sev- Khiar, were screened for resistance to septoria tritici eral known Stb genes [76–78] and also in other germ- blotch. The RIL population was generated by crossing the plasm several QTL have contributed to broad efficacy of resistant parent ‘Agili39’ to the susceptible Khiar follow- resistance [59]. Our data also confirm that stacking QTL ing a single seed descent approach (SSD). The F1 plants in durum wheat results in broad efficacy of STB resist- were selfed to generate the F2 seeds. One head row of ance. This study has identified genotypes harboring each F2 plants were then randomly selected and sown in diverse resistance loci entailing dual actions at the differ- one row to produce the F2 derived F3 plants (F3). This ent physiological stages constituting thus potential effec- procedure was followed for all subsequent generations up tive sources for Z. tritici resistance and will thus support to the F10 plants. Therecombinant population (Agili39’/ sustainable breeding approach for Z. tritici resistance in Khiar) is maintained and available upon request at the durum wheat. National Institute of Agronomy -Tunisia (INAT) and at Finally, we explored QTL epistasis and identified four the CIMMYT gene banks. significant pairwise interactions of the identified QTL For STB resistance screening, we performed three with an epistatic QTL mapped on chromosome 5B, experiments (Table 2). The first experiment was repeated designated as Qstb5B. Hence, the epistasis analysis has three times and comprised the screening of the 11 Tuni- revealed other QTL that affects the expression of Z. trit- sian landraces at the seedling stage (Z13.3/21) with a ici resistance in the ‘Agili39’/Khiar population. In fact, panel of 20 Z. tritici isolates (Table  2). The first experi- epistatic interactions between QTL are an important ment was performed to understand overall resistance factor that affects the phenotypic expression of genes patterns to STB and to select potential isolates for further and genetic variation in populations [79–81]. Similarly, screening of the RILs derived from the ‘Agili39’/Khiar to many other studies [82–84], our data demonstrated cross, which consisted of the second experiment (experi- interaction between QTL having main effect (Qstb2B_1 ment 2), performed thrice at the seedling stage. Finally, in and Qstb2B_2) that are involved in epistasis with the the third experiment we tested the ‘Agili39’/Khiar popu- Qstb5B QTL affecting the same trait. The epistasis lation under field conditions in Oued-Bejá, located in analysis showed an additive-by-additive effect between North-Western Tunisia, over a period of five years, 2011 the Qstb2B_1/ Qstb5B QTL and the Qstb2B_2/Qstb5B – 2014 and 2016 with three different Z. tritici isolates. QTL, with a major effect of QTL mapped on chromo- some 2B (Qstb2B_1 and Qstb2B_2) over the Qstb5B QTL. Interestingly, the epistasis analysis showed that Screening landraces and RILs population at the seedling the Qstb2B_2 QTL, proven to control pycnidia develop- stage for resistance to Zymoseptoria tritici ment at the adult plant stage by QTL analysis, has also an Experimental design, Plants management and growth effect in controlling pycnidia development at the seedling conditions stage when interacting with the epistatic Qstb5B QTL. A first experiment that consisted of the pre-screening Nonetheless, the epistasis analysis did not indicate an of 11 Tunisian landrace accessions with 20 Z. tritici Ferjaoui et al. BMC Genomics (2022) 23:372 Page 14 of 20 Table 8 Nine Tunisian durum wheat landraces and two cultivars Zymoseptoria tritici isolates and inoculation procedures that were investigated for resistance to Zymoseptoria tritici Pre-cultures of each isolate were prepared in an auto- Name Habitus Source Empiric evaluation of septoria claved 100  ml Erlenmeyer flask containing 50  ml yeast tritici blotch under field glucose (YG) liquid medium (30  g glucose, 10  g yeast conditions per litre distilled water). The flasks were inoculated with Agili 37 landrace INATb Resistant frozen isolate samples that were selected from different Agili 38 landrace INAT Resistant durum wheat growing countries and maintained at -80 °C Agili 39 (P)a landrace INAT Resistant in a Z. tritici isolate collection at WUR. The inoculated Agili 41 landrace INAT Resistant flasks were subsequently placed in an incubated rotary Azizi 27 landrace INAT Resistant shaker (Innova 4430, New Brunswick Scientific, USA) Derbessi 12 landrace INAT Resistant set at 125 rpm and 15 °C for 5–7 days. These pre-cultures Jneh Khotifa 85 landrace INAT Resistant were then used to inoculate two 1L Erlenmeyer flasks Mahmoudi 101 landrace INAT Resistant containing 500 ml YG media per isolate that were incu- Sbei 99 landrace INAT Resistant bated under the aforementioned conditions to provide sufficient inoculum for the seedling inoculation assays at Khiar (P)a cultivar INRATc Susceptible growth stage (GS) 11 [85]. Spores were subsequently col- Karim cultivar INRAT Susceptible lected after overnight settling in static cultures, concen- a Parents of the recombinant inbred population trated by decanting the supernatant medium, adjusted to b National Institute of Agronomy-Tunis, Tunisia 107 c spores ml-1 in a total volume of 40 ml for a set of 18 National Institute of Agronomical Research-Tuni plastic pots or 24 Jiffy® pots and supplemented with two drops (µl/ml) of Tween 20 surfactant (MERCK®, Not- isolates was conducted at the seedling stage. The pre- tingham, UK). screening assay followed a complete block design in Prior to the inoculation, plant development was three replicates and included the susceptible parent cv. allowed for 10 days in a greenhouse adjusted at a temper- Khiar, the susceptiblecv. Karim and the resistant par- ature of 18/16 °C (day/night rhythm) and relative humid- ent ‘Agili39’ as checks (Table 2 and 8). This experiment ity (RH) of 70%. Inoculations were conducted by spraying enabled the selection of eight Z. tritici isolates that dis- the inoculum over the seedlings that were placed in an criminated between the ‘Agili39’ and the cv. Khiar par- inoculation cabinet on a rotary table, adjusted at 15 rpm, ents and were subsequently used to screen the ‘Agili39’/ which is equipped with interchangeable atomizers and Khiar derived recombinant inbred lines (RILs) at the a water cleaning device to avoid cross- contamination. seedling stage. Infected plants were incubated in transparent plastic bags For the seedling assay of the RIL population (experi- for 48  h under 100% RH. Post-inoculation conditions ment 2), we followed a split plot design with isolates as were set at a temperature of 22/ ± 2 °C and RH of ≥ 95%. whole plots. Each whole plot consists of three neighbour- Light intensity (son-T Agro 400 W lamps) and day length ing trays of fifty-four pots. The tested isolates were ran- (16/8  h light/dark) were similar during pre- and post- domly allocated to the whole plots. Individual pots were inoculation conditions. Ten days after inoculation, seed- the experimental units, and they were randomly arranged lings were trimmed for the second and subsequent leaves for each isolate/replicate combination on separate par- to enable sufficient light on the inoculated primary leaves allel greenhouse tables. ‘Replicate’, ‘isolate’, ‘line’ and the for appropriate disease development. Fertilizer (Sporu- ‘isolate x line’ interaction were the fixed terms of our split mix PG®, Rotterdam, Netherlands; 0.5  g.l1) was applied plot model. However, the random term of the model con- to support plant growth. sists of the ‘replicate x whole plot’ interaction. In all three replicates, eleven checks several checks were included Disease severity scoring in the seedling assay with both parents ‘Agili39’ and cv. Khiar (Table 8). In the seedling assays, disease severities were evalu- Five seeds per pot per accessions were grown in VQB ated at 15, 18 and 21 days post-inoculation (dpi). These 7 × 7x8 cm plastic pots (TEKU®, Lohne, Germany), multiple observations enabled Area Under the Disease whereas 157 F6 RILs of the ‘Agili39’/Khiar mapping pop- Progress Curve (AUDPC) calculations for quantitative ulation were planted in round peat pots (Jiffy, Moerdijk, analyses of temporal differences in disease progress. We Netherlands), also five seeds per pot, using a special mix- estimated the quantitative presence of necrosis and pyc- ture for growing seeds (Substraat Zaai) provided by the nidia on the inoculated seedling leaves in percentages greenhouse facility Unifarm of Wageningen University [36, 37, 86]. AUDPC calculations for seedling scores fol- and Research (WUR), The Netherlands. lowed the trapezoidal method, which approximates the F erjaoui et al. BMC Genomics (2022) 23:372 Page 15 of 20 time variable and calculates the average disease intensity Genotyping and construction of linkage map between each pair of adjacent time points [87]. Genomic DNA was extracted from fresh leaves using a modified CTAB (cetyltrimethylammonium bromide) Screening RILs population at the adult plant stage method and quantified using NanoDrop 8000 spec- for resistance to Zymoseptoria tritici trophotometer V 2.1.0. Whole-genome profiling was Experimental design and inoculation procedures performed using DArT-Seq™ technology by Diversity Adult plant assays of the ‘Agili39’/Khiar population were Arrays Technology Pty Ltd, Australia, as described by conducted at Oued-Bejá, located in North-Western Kilian et  al. [89]. In brief, the DArT-Seq™ technology Tunisia (36° 46′ 27.516’’ N, 10° 3′ 36.432’’ E), for five years, was optimized by selecting the most appropriate com- 2011 – 2014 and 2016, with Tun1 – Tun6 and IIB-123 Z. plexity reduction method for wheat (PstI-MseI restric- tritici isolates, respectively. This region belongs to the tion enzymes). DNA fragments digested with restriction sub-humid bioclimatic zone of Tunisia with an average enzymes were ligated with PstI adaptors and unique rainfall ranging from 500 to 850  mm and a daily mean barcodes, and then amplified following PCR. Ampli- temperature between 10–28 °C [72]. cons were pooled and sequenced in a 96-multiplex on a For the field trials, we adopted an augmented rand- HiSeq2000 (Illumina, USA) resulting in a total of 5,891 omized complete block design. Five blocks with 1.5  m raw DArTSeq SNP markers. width and spaced 1.5 m were linearly drilled with 30 to The SNP markers were first filtered according to their 35 RILs, parents and four modern durum cultivars per polymorphism between the parents ‘Agili39’ and ‘Khiar’. block. Each line was sown as one spike per row of 1.5 m A total of 3,459 filtered SNPs was subsequently tested length and spaced 25  cm. We randomized all RILs, the for the segregation distortion that was determined by parents and four additional checks modern durum wheat the calculation of the Chi-Square (X 2). Hence, SNPs cvs. Karim, Nasr, Maali and Salim, important in Tunisian with a major distortion of p < 0.001 were removed. breeding programs and showing different levels of sus- Moreover, identical SNPs, which should be mapped to ceptibility in each block. For the 2016 field trial, we used the same position on the linkage group, were identified a complete random block design with three replicates (loci similarity = 1). Hence, only one marker of ‘similar with both parents ‘Agili39’ and ‘Khiar’ as checks. loci’ was retained on the linkage map to reduce the cal- Field inoculations were conducted with three isolates culation burden. The final set of filtered SNPs was then (Tun1, Tun6 and IIB123) across the F5:6-F9:10 RILs gen- used to generate a genetic linkage map using JoinMap erations. We used Z. tritici isolate Tun6 for three years to ® 4 software [90]. The regression mapping algorithm screen the F5:6 (N = 164), the F7:8 (N = 158) and the F8:9 was used to construct the genetic map. Linkage between (N = 157) RILs in 2011, 2013 and 2014, respectively and markers, recombination rate (Θ), and map distances isolate Tun1 to screen the F6:7 (N = 158) in 2012 and the was calculated using the Kosambi mapping function in F8:9 (N = 157) in 2014. In 2016, we screened the F9:10 JoinMap [91]. Markers were grouped using a minimum (N = 155) with Z. tritici isolate IIB123. In all field trials and independence LOD (logarithm of the odds) score of across all generations, plants were inoculated twice. The 10.0 and linkage groups were established at a minimum first inoculation occurred at the tillering stage (GS 21–26) LOD score of 3.0. Markers were linearly aligned in each in order to initiate the disease infection, and the second linkage group, converting the recombination rates into inoculation was applied after 25–30  days when all RILs the Kosambi’s map distance in centimorgans. A sequen- reached approximately the elongation stage (GS37) [85]. tial map builds up strategy was followed to determine The second inoculation was applied to ensure an increase the best marker position [92]. The best fitting posi- in the disease pressure. Field inoculations were conducted tion of markers was examined based on the goodness- with a spore suspension adjusted to 107 spores/ml of the of-fit test (chi-square) for the resulting map. The final corresponding Z. tritici isolate using a CO2-pressurized map included bin markers that excluded similar SNP knapsack sprayer with a 1 m hand-held boom till run-off. markers. Linkage groups were subsequently assigned to chro- Disease severity scoring in the field tests mosomes by aligning the SNP marker sequences of the For the field evaluations, we scored pycnidial classes at linkage groups to the reference genome of the Triticum 28  days post the second inoculation (GS 51) [85] that turgidum subsp. durum Svevo.v1 and using the BLASTn were later transformed to pycnidia percentages (0 = no function in the publicly available sequence database pycnidia, 1 = 12%, 2 = 25%, 3 = 50%, 4 = 75% and 5 = 87%) ‘Ensemble Plants’(https:// plants. ensem bl. org/ Triti cum_ [60, 88]. turgi dum/ Tools/ Blast? db= core) [93]. Ferjaoui et al. BMC Genomics (2022) 23:372 Page 16 of 20 Statistical analyses and r is the number of replicates which is three per iso- Seedling data statistical analyses late for all seedling assays. Pycnidia AUDPC scores (P-AUDPC) of the pre-screen- ing experiment (experiment 1) were analysed for their variance using the ‘aov’ function in the R environment Field data statistical analyses [42] to test for the effect of ‘line’, the effect of ‘isolate’ and Field pycnidia scores generated by the inoculation of the effect of any ‘line x isolate’ interactions. The following adult plants with Tun1 and Tun6 Z. tritici isolates were linear model was fitted to the observed P-AUDPC scores: tested for their variance to check for the effect of ‘gen- otype’, ‘block’, ‘year’, ‘isolate’ and for any two-way and Yijk = µ+ αi + βj + (αβ)ij + ǫijk three-way interactions between the independent variable where Yijk is the P-AUDPC score in the Kth replicate ‘genotype’ and the independent factors ‘block’, ‘year’ and with isolate i and line j. μ is the overall mean. αi is the ‘isolate’. A linear model was fitted to the observed pyc- isolate i main effect. βj is the line j main effect. (αβ)ij is nidia at the adult plant stage as follow: the interaction effect between isolate i and line j and ϵijk Yijkt =  + i + j + k + t + ()ij + ()kj + ()tj + ()ijt + ijkt is the unexplained error. For the split-plot seedling experiment performed on where Yijkt is the pycnidia coverage score of RIL j inocu- the RIL population (experiment 2), an analysis of vari- lated with isolate i in the Kth block at year t. μ is the overall ance for the necrosis and the pycnidia AUDPC scores mean. αi is the main effect of the isolate i. βj is the main (N-AUDPC and P-AUDPC) was performed using the sp. effect of the RIL j. γk is the main effect of the block k. ηt is plot function available in the Agricolae package in the R the main effect of the year t. (αβ)ij, (γβ)kj, (ηβ)tj are the two- environment [42, 94].Isolates were the whole-plot fac- way interaction effects of the RIL j with the isolate i, the RIL tors and RILs were the sub-plot factors. We fitted the j with the block k and the RIL j with the year t, respectively. observed N-AUDPC and P-AUDPC scores to the follow- (αβη)ijt is the three-way interaction effect of RIL j with iso- ing linear model: late i and with year t and finally, εijkt is the residual. For the Tun1 and the Tun6 isolates tested on multiple Y IJK = µ+ αi + ηk(i)+ βj + (αβ)ij + ǫijk , consecutive years, a broad sense heritability (H2) was estimated as follow: where Yijk is the N-AUDPC or P-AUDPC score in the Kth replicate of a plot with isolate i  and RIL j. μ is the σ2G overall mean. αi is the fixed effect of isolate i. ηk(i) is the H2 = σ 2 E σ2 2G σ G ε + E + r whole-plot error. βj is the fixed effect of RIL j. (αβ)ij is the interaction effect between isolate i and RIL j and ϵijk is where σ2G is the variance component due to genotypes. the sub-plot error. σ2GE is the variance component due to the interaction Significant differences between isolates were deter- between the genotype and the isolate. E is the number of mined using the least significant difference (LSD) of years which is 2 and 3 for Tun1 and Tun6 isolates, respec- N-AUDPC and P-AUDPC scores and using the Agricolae tively. σ2ε is the variance component due to the unex- package in the R environment. ‘RIL x isolate’ grouping of plained error and r is the number of blocks which is five necrosis and pycnidia AUDPC scores was defined based for all adult assays. on the Bonferroni test at p < 0.05. Homogeneity of the The least-square means (Lsmeans) of pycnidia cover- seedling replicates was checked and homogeneous data ages were derived from the individual year trials using across replications were subsequently averaged and used the SAS software [96]. Transformed field data were sub- for the seedling QTL analysis [95]. sequently considered for the field QTL analysis. Repeatability was estimated for the necrosis and pyc- nidia AUDPC generated at the seedling stage as follow: Correlation between seedling and adult plant assays / ( ) and frequency distribution σ2G σ E 2G σ2G Repeatability = + E A heatmap correlation matrix was calculated for the gen- σ2ε + r erated phenotypes at the seedling and the adult plant where σ2G is the variance component due to genotypes, stages using the “Corr” function and visualized using the σ2GE is the variance component due to the interaction “corrplot” package in the R environment [43]. The “Corr” between the genotype and the isolate, E is the number function was set to calculate the Pearson correlation. of isolates which is eight in the seedling experiment, σ2ε Positive and negative correlations were displayed in the is the variance component due to the unexplained error heatmap matrix by a color gradient from red, indicat- ing a positive correlation, to blue, indicating a negative F erjaoui et al. BMC Genomics (2022) 23:372 Page 17 of 20 correlation between the phenotypic traits. Frequency dis- Table 3. Linkage groups, correspondent durum wheat chromosome, the tribution figures were generated using the “hist” function average length indicated in centimorgans (cM), number of SNPs and the in the R environment. average inter-loci distance (cM) in the ’Agili39’/khiar genetic linkage map. Table 4. Sequences, genetic and physical positions of the flanking mark- ers linked to the detected QTL on the ’Agili39’/Khiar mapping population. QTL analysis procedure An inclusive composite interval mapping of additive Acknowledgements (ICIM-ADD) functionality in QTL IciMapping v4.1 [97] We would like to thank Els Verstappen and all members of Bio-interaction and Plant health (Wageningen University and Research) for excellent technical was used. Additive QTL were detected using a walk speed support, Nicolas Letreux and Bertus Van der Laan (Unifarm, Wageningen of 1.0  cM and the probability used in stepwise regres- University) for greenhouse assistance and maintenance. We, also would like to sion for additive QTLs was 0.001. The logarithm of odds thank all members of the Genetics and Plant Breeding group at the National Institute of Agronomy, Tunisia , and staff members at the Regional Field Crops (LOD) value of 3.0 was chosen to declare significant QTL, and Research Center—Beja (Tunisia) for field assistance, and CIMMYT staff and the LOD value was calculated from 1000 permuta- members (Mexico) for genotyping support. We, finally, highly acknowledge tions with type I error of 0.01. The phenotypic variance Paul Keizer and Pieter Vereijken (Biometris, Wageningen University and Research) for providing statistical support. explained (PVE) and additive effect of individual QTL at the LOD peaks were also obtained. Identified QTL were Authors’ contributions plotted against their corresponding linkage groups using GHJK and SH designed the study. SF produced the ‘Agili39’/Khiar RIL population,provided the plant materials, and generated the field phenotypic the MapChart © software version 2.3 [98]. Subsequently data. LA generated the seedling phenotypic data. KA and SD genotyped the and for the identified QTL contributing to resistance, we ‘Agili39’/Khiar mapping population. HJS, LA and RBS generated the linkage aligned the linked SNP markers to the reference genome genetic map. LA performed the statistical analysis. BB and SS performed the QTL and the epistasis analyses. LA, SH and SF wrote the manuscript. All of the Triticum turgidum subsp. durum Svevo.v1 and authors read and reviewed the manuscript. using the BLASTn function (https:// plants. ensem bl. org/ Triti cum_ turgi dum/ Tools/ Blast? db= core) in the publicly Funding The research was partly supported by the Japan International Cooperation available sequence database ‘Ensemble Plants’ [93]. Agency (JICA)–Japan Science and Technology Agency (JST)’s Science and Epistatic interactions between QTL were identified by Technology Research Partnership for Sustainable Development (SATREPS) the inclusive composite interval mapping of digenic epi- under the project entitled “Valorization of Bio-resources in Semi-Arid and Arid Land for Regional Development”. static QTL (ICIM-EPI) method implemented in QTL LA was financially supported by the Monsanto’s Beachell-Borlaug International IciMapping software v4.1 [97]. The LOD threshold was Scholars Program (3,340,030,501) and The ‘Bourse d’alternance’ from the Minis- set at 5.00 to declare significant epistatic QTL and LOD try of Higher Education and Scientific Research—Tunisia. value was calculated from 1000 permutations at the sig- Availability of data and materials nificance level of 0.05. The used plant materials are available at the CIMMYT gene bank upon a rea- The highest epistasis interactions detected between sonable request to the co-author Karim Ammar < k.ammar@cgiar.org > The sequencing raw data files generated during this study are uploaded to the QTL (highest LOD) were subsequently selected for the CIMMYT data repository (https:// hdl. handle. net/ 11529/ 10548 618). The a two-way interaction test to examine pycnidia AUDPC analyzed data are available as additional files to this article. means of the variant allele combinations linked to the epistatic QTL. Alleles linked to the detected QTL were Declarations named by ‘R’ and ‘S’ denoting the resistant and suscepti- Ethics approval and consent to participate ble alleles, respectively. Hence, pycnidia AUDPC means The used plant material comprises either registered Triticum durum cultivars of the four allele combinations ‘RR’, ‘RS’, ‘SR’ and ‘SS’ and landrace accessions or experimental lines developed by the correspond- were examined and Boxplots were generated using the ing author Professor Sonia Hamza at the National Institute of Agronomy (INAT) and by Dr. Sahbi Ferjaoui at the Regional Field Crops Research Center of Beja ‘ggplot2’ package in the R environment [99]. (CRRGC). Experimental research and field studies, including the collection of plant material complies with relevant institutional, national and international Supplementary Information guidelines and legislation. The online version contains supplementary material available at https:// doi. Consent for publication org/ 10. 1186/ s12864- 022- 08560-2. Not applicable. Additional file 1: Fig. S1. Frequency distributions of the disease severity Competing interests assessed as percentage pycnidia in seedlings of the F6 recombinant The authors declare that they have no competing interests. inbred lines of the ‘Agili39’/Khiar population. ‘A’ and ‘K’ are referring to the ‘Agili39’and cv. Khiar parents, respectively. Author details 1 Laboratory of Bioaggressors and Integrated Protection in Agriculture (BPIA), Additional file 2: Table 1. Analysis of variance of pycnidia percent of National Institute of Agronomy of Tunisia (INAT), 43 Avenue Charles Nicolle, nine Tunisian durum landraces and two modern varieties inoculated 1082 El Mahrajène, Tunis, Tunisia. 2 Present Address Field Crops Labora- with a diverse range of twenty durum derived Zymoseptoria tritici isolates. tory, Regional Field Crops Research Center of Beja (CRRGC), P.O. Box 9000, Table 2. Pearson correlation between the different tested isolates on the Beja, Tunisia. 3 Bio-Interaction and Plant Health, Wageningen University ’Agili39’/khiar population RILs at the seedling and the Adult plant stages. and Research, PO Box 16, 6700AA Wageningen, The Netherlands. 4 The Gradu- ate School ‘Experimental Plant Sciences’ (EPS), Wageningen Campus, 6708 Ferjaoui et al. BMC Genomics (2022) 23:372 Page 18 of 20 PB Wageningen, The Netherlands. 5 Present Address Center for Desert Agri- 18. Meamiche Neddaf H, Aouini L, Bouznad Z, Kema GH. Equal distribu- culture, Biological and Environmental Science and Engineering Division, King tion of mating type alleles and the presence of strobilurin resist- Abdullah University of Science and Technology, Thuwal, Saudi Arabia. 6 Present ance in algerian zymoseptoria tritici field populations. Plant Dis. address Higher Institute of Agronomy of Chott Meriam (ISA-CM), 4042 Sousse, 2017;101(4):544–9. Tunisia. 7 International Maize and Wheat Improvement Center (CIMMYT), 19. Hassine M, Siah A, Hellin P, Cadalen T, Halama P, Hilbert J-L, Hamada W, Apdo. Postal 6-641, 06600 Mexico, D.F., Mexico. 8 Plant Breeding, Wageningen Baraket M, Yahyaoui A, Legrève A. Sexual reproduction of Zymosepto- University and Research, P.O. Box 386, 6700 AJ Wageningen, The Netherlands. ria tritici on durum wheat in Tunisia revealed by presence of airborne 9 Institute of Plant Breeding, Genetics and Genomics, Department of Plant inoculum, fruiting bodies and high levels of genetic diversity. Fungal Pathology and Institute of Plant Breeding, University of Georgia, Griffin, GA Biol. 2019;123(10):763–72. 30223, USA. 10 Present Address United States Department of Agriculture USDA, 20. Boukef S, McDonald BA, Yahyaoui A, Rezgui S, Brunner PC. Frequency Crop Genetics and Breeding Research Unit, Tifton, GA, USA. 11 CRP-Wheat of mutations associated with fungicide resistance and population Septoria Phenotyping Platform (CIMMYT-IRESA), Regional Field Crops Research structure of Mycosphaerella graminicola in Tunisia. Eur J Plant Pathol. Center of Beja (CRRGC), BP 350, 9000 Beja, Tunisia. 12 Laboratory of Phytopa- 2012;132(1):111–22. thology, Wageningen University and Research, PO box 16, 6700AA Wagenin- 21. Singh RP, Herrera-Foessel S, Huerta-Espino J, Singh S, Bhavani S, Lan C, gen, The Netherlands. Basnet BR. Progress towards genetics and breeding for minor genes based resistance to Ug99 and other rusts in CIMMYT high-yielding Received: 20 September 2021 Accepted: 14 April 2022 spring wheat. J Integr Agric. 2014;13(2):255–61. 22. Singh R, Huerta-Espino J, Bhavani S, Herrera-Foessel S, Singh D, Singh P, Velu G, Mason R, Jin Y, Njau P. Race non-specific resistance to rust diseases in CIMMYT spring wheats. Euphytica. 2011;179(1):175–86. 23. Saintenac C, Zhang W, Salcedo A, Rouse MN, Trick HN, Akhunov E, References Dubcovsky J. Identification of wheat gene Sr35 that confers resistance 1. Arzani A, Ashraf M. Cultivated Ancient Wheats (Triticum spp): A Potential to Ug99 stem rust race group. Science. 2013;341(6147):783–6. Source of Health-Beneficial Food Products. Compr Rev Food Sci Food 24. Chen S, Rouse MN, Zhang W, Jin Y, Akhunov E, Wei Y, Dubcovsky J. Fine Safety. 2017;16(3):477–88. mapping and characterization of Sr21, a temperature-sensitive diploid 2. FAOSTAT. 2017. http:// www. fao. org/ faost at/ en/# compa re wheat resistance gene effective against the Puccinia graminis f. sp. 3. Shiferaw B, Smale M, Braun H-J, Duveiller E, Reynolds M, Muricho G. Crops tritici Ug99 race group. Theor Appl Genet. 2015;128(4):645–56. that feed the world 10. Past successes and future challenges to the role 25. Bajgain P, Rouse M, Bulli P, Bhavani S, Gordon T, Wanyera R, Njau P, played by wheat in global food security. Food Sec. 2013;5(3):291–317. Legesse W, Anderson J, Pumphrey M. Erratum to: Association mapping 4. Charmet G. Wheat domestication: lessons for the future. CR Biol. of North American spring wheat breeding germplasm reveals loci 2011;334(3):212–20. conferring resistance to Ug99 and other African stem rust races. BMC 5. Palamarchuk A. Selection strategies for traits relevant for winter and Plant Biol. 2016;16(1):1–9. facultative durum wheat. Durum wheat breeding: current approaches 26. Kingsbury N. Hybrid: the history and science of plant breeding: Univer- and future strategies. New York: Food Products Press; 2005. p. 599–644. sity of Chicago Press; 2009. p. 512. 6. Royo C, Nazco R, Villegas D. The climate of the zone of origin of Mediter- 27. Mondal S, Rutkoski JE, Velu G, Singh PK, Crespo-Herrera LA, Guzmán ranean durum wheat (Triticum durum Desf ) landraces affects their C, Bhavani S, Lan C, He X, Singh RP. Harnessing diversity in wheat agronomic performance. Genet Resour Crop Evol. 2014;61(7):1345–58. to enhance grain yield, climate resilience, disease and insect pest 7. Negisho K, Shibru S, Pillen K, Ordon F, Wehner G. Genetic diversity of resistance and nutrition through conventional and modern breeding Ethiopian durum wheat landraces. PLoS ONE. 2021;16(2):e0247016. approaches. Front Plant Sci. 2016;7:991. https:// doi. org/ 10. 3389/ fpls. 8. Kidane YG, Hailemariam BN, Mengistu DK, Fadda C, Pè ME, Dell’Acqua M. 2016. 00991. Genome-Wide Association Study of Septoria tritici Blotch Resistance in 28. Lopes MS, El-Basyoni I, Baenziger PS, Singh S, Royo C, Ozbek K, Aktas Ethiopian Durum Wheat Landraces. Front Plant Sci. 2017;8:1586. https:// H, Ozer E, Ozdemir F, Manickavelu A: Exploiting genetic diversity from doi. org/ 10. 3389/ fpls. 2017. 01586. landraces in wheat breeding for adaptation to climate change. J Exp 9. Fabricant F. Rome’s glory is now Tunisia’s. New York Times 1998. https:// Bot. 2015;66(12):3477–86. https:// doi. org/ 10. 1093/ jxb/ erv122. www. nytim es. com. 29. Aoun M, Kolmer JA, Rouse MN, Elias EM, Breiland M, Bulbula WD, Chao 10. Oliveira HR, Campana MG, Jones H, Hunt HV, Leigh F, Redhouse DI, Lister S, Acevedo M. Mapping of novel leaf rust and stem rust resistance DL, Jones MK. Tetraploid wheat landraces in the Mediterranean basin: genes in the Portuguese durum wheat landrace PI 192051. G3: Genes, taxonomy, evolution and genetic diversity. PLoS ONE. 2012;7(5):e37063. Genomes, Genet. 2019;9(8):2535–47. 11. Davis DK: Resurrecting the granary of Rome: environmental history and 30. Alemu SK, Huluka AB, Tesfaye K, Haileselassie T, Uauy C. Genome-wide French colonial expansion in North Africa: Ohio University Press; 2007. association mapping identifies yellow rust resistance loci in Ethiopian 12. Ouaja M, Aouini L, Bahri B, Ferjaoui S, Medini M, Marcel TC, Hamza durum wheat germplasm. PLoS ONE. 2021;16(5):e0243675. S. Identification of valuable sources of resistance to Zymoseptoria 31. Royo C, Ammar K, Soriano JM, Villegas D. Agronomic, physiological and tritici in the Tunisian durum wheat landraces. Eur J Plant Pathol. genetic changes associated with evolution, migration and modern 2020;156(2):647–61. breeding in durum wheat. Front Plant Sci. 2021;12:1318. 13. Kamel S, Cherif M: Tan spot of wheat in Northern Tunisia: distribution, 32. Ghaneie A, Mehrabi R, Safaie N, Abrinbana M, Saidi A, Aghaee M. prevalence, incidence and severity. Cereal Research Communications Genetic variation for resistance to septoria tritici blotch in Iranian 2021. tetraploid wheat landraces. Eur J Plant Pathol. 2012;132(2):191–202. 14. McDonald BA, Mundt CC. How knowledge of pathogen population 33. Porras R, Pérez-de-Luque A, Sillero J, Miguel-Rojas C. Behavior of Span- biology informs management of Septoria tritici blotch. Phytopathol- ish durum wheat genotypes against Zymoseptoria tritici: resistance ogy. 2016;106(9):948–55. and susceptibility. Span J Agric Res. 2021;19(3): e1002. 15. Abdedayem W. M’barek S, Souissi A, Laribi M, Araar C, Kouki H, Fakhfakh 34. Brown JK, Chartrain L, Lasserre-Zuber P, Saintenac C. Genetics of M, Yahyaoui A: Septoria tritici blotch disease progression and physi- resistance to Zymoseptoria tritici and applications to wheat breeding. ological traits variation in durum wheat variety mixtures. J New Sci. Fungal Genetic Biololgy. 2015;79:33–41. 2021;80:4664–74. 35. Yang N, McDonald MC, Solomon PS, Milgate AW. Genetic mapping of 16. Berraies S, Ammar K, Salah Gharbi M, Yahyaoui A, Rezgui S. Quantitative Stb19, a new resistance gene to Zymoseptoria tritici in wheat. Theor inheritance of resistance to Septoria tritici blotch in durum wheat in Appl Genet. 2018;131(12):2765–73. Tunisia. Chilean J Agri Res. 2014;74(1):35–40. 36. Kema GH, van Silfhout CH. Genetic variation for virulence and resist- 17. Hamada W. First isolation of the Mycosphaerella graminicola tele- ance in the wheat-Mycosphaerella graminicola pathosystem III. omorph stage causing Septoria leaf blotch on wheat in Tunisia. New Comparative seedling and adult plant experiments. Phytopathology. Dis Rep. 2014;29:18–18. 1997;87(3):266–72. F erjaoui et al. BMC Genomics (2022) 23:372 Page 19 of 20 37. Kema G, Annone JG, Sayoud R, Van Silfhout CH, Van Ginkel M, De 58. Ghaffary SMT, Faris JD, Friesen TL, Visser RG, van der Lee TA, Rob- Bree J. Genetic variation for virulence and resistance in the wheat- ert O, Kema GH. New broad-spectrum resistance to septoria tritici Mycosphaerella graminicola pathosystem I. Interactions between path- blotch derived from synthetic hexaploid wheat. Theor Appl Genet. ogen isolates and host cultivars. Phytopathology. 1996;86(2):200–12. 2012;124(1):125–42. 38. Kema G, Annone J, Sayoud R, Van Silfhout C, Van Ginkel M, De 59. Ghaffary SMT, Robert O, Laurent V, Lonnet P, Margalé E, van der Lee TA, Bree J. Genetic variation for virulence and resistance in the wheat- Visser RG, Kema GH. Genetic analysis of resistance to septoria tritici blotch Mycosphaerella graminicola pathosystem. I: Interactions between path- in the French winter wheat cultivars Balance and Apache. Theor Appl ogen isolates and host cultivars. Phytopathology. 1996;86(2):200–12. Genet. 2011;123(5):741–54. 39. Tabib Ghaffary SM. Efficacy and mapping of resistance to Mycosphaerella 60. Adhikari TB, Cavaletto JR, Dubcovsky J, Gieco JO, Schlatter AR, Goodwin graminicola in wheat. PhD dissertation Wageningen University; 2011. SB. Molecular mapping of the Stb4 gene for resistance to septoria tritici https:// edepot. wur. nl/ 169465. blotch in wheat. Phytopathology. 2004;94(11):1198–206. 40. Ballini E, Tavaud M, Ducasse A, Sanchez D, Paux E, Kitt J, Charmet G, 61. Arraiano LS, Worland AJ, Ellerbrook C, Brown JKM. Chromosomal loca- Audigeos D, Roumet P, David J. Genome wide association mapping for tion of a gene for resistance to septoria tritici blotch (Mycosphaerella resistance to multiple fungal pathogens in a panel issued from a broad graminicola)in the hexaploid wheat ’Synthetic 6x’. Theor Appl Genet. composite cross-population of tetraploid wheat Triticum turgidum. 2001;103(5):758–64. Euphytica. 2020;216:1–17. 62. Goudemand E, Laurent V, Duchalais L, Tabib Ghaffary S, Kema G, Lonnet 41. El Haddad N, Kabbaj H, Zaïm M, El Hassouni K, Tidiane Sall A, Azouz M, P, Margalé E, Robert O: Association mapping and meta-analysis: two Ortiz R, Baum M, Amri A, Gamba F. Crop wild relatives in durum wheat complementary approaches for the detection of reliable Septoria tritici breeding: Drift or thrift? Crop Sci. 2021;61(1):37–54. blotch quantitative resistance in bread wheat (Triticum aestivum L.). Mol 42. Team RC: R. A language and environment for statistical computing. 2013. Breeding 2013, 32(3). 43. Wei T, Simko V, Levy M, Xie Y, Jin Y, Zemla J. Package ‘corrplot.’ Statistician. 63. Risser P, Ebmeyer E, Korzun V, Hartl L, Miedaner T. Quantitative trait loci 2017;56:316–24. for adult-plant resistance to Mycosphaerella graminicola in two winter 44. Fones H, Gurr S. The impact of Septoria tritici Blotch disease on wheat: An wheat populations. Phytopathology. 2011;101(10):1209–16. EU perspective. Fungal Genet Biol. 2015;79:3–7. 64. Goodwin SB, Cavaletto JR, Hale IL, Thompson I, Xu SS, Adhikari TB, Dub- 45. Aouini L. Durum wheat and septoria tritici blotch: genes and prospects covsky J. A New Map Location of Gene for Resistance to Septoria Tritici for breeding. PhD dissertation Wageningen University 2018. https:// doi. Blotch in Wheat. Crop Sci. 2015;55(1):35–43. org/ 10. 18174/ 443719. 65. Adhikari TB, Wallwork H, Goodwin SB. Microsatellite Markers Linked to 46. Arraiano L, Brading P, Brown J. A detached seedling leaf technique to the and Genes for Resistance to Septoria Tritici Blotch in Wheat. Crop Sci. study resistance to Mycosphaerella graminicola (anamorph Septoria 2004;44(4):1403–11. tritici) in wheat. Plant Pathol. 2001;50(3):339–46. 66. Chartrain L, Sourdille P, Bernard M, Brown J. Identification and location 47. Berraies S, Gharbi MS, Belzile F, Yahyaoui A, Hajlaoui MR, Trifi M, Jean of Stb9, a gene for resistance to septoria tritici blotch in wheat cultivars M, Rezgui S. High genetic diversity of Mycospaherella graminicola Courtot and Tonic. Plant Pathol. 2009;58(3):547–55. (Zymoseptoria tritici) from a single wheat field in Tunisia as revealed by 67. Odilbekov F, Armoniené R, Koc A, Svensson J, Chawade A. GWAS-Assisted SSR markers. Afr J Biotechnol. 2013;12(12):1344–9. Genomic Prediction to Predict Resistance to Septoria Tritici Blotch in 48. Kema GHJ, Mirzadi Gohari A, Aouini L, Gibriel HAY, Ware SB, van den Nordic Winter Wheat at Seedling Stage. Front Genet. 2019;10:1224. Bosch F, Manning-Smith R, Alonso-Chavez V, Helps J, Ben M’Barek S, 68. Radecka-Janusik M, Czembor PC. Genetic mapping of quantitative trait et al. Stress and sexual reproduction affect the dynamics of the wheat loci (QTL) for resistance to Septoria tritici blotch in a winter wheat cultivar pathogen effector AvrStb6 and strobilurin resistance. Nat Genet. Liwilla. Euphytica. 2014;200(1):109–25. 2018;50(3):375–80. 69. Stadlmeier M, Hartl L, Mohler V. Usefulness of a multiparent advanced 49. Zhong Z, Marcel TC, Hartmann FE, Ma X, Plissonneau C, Zala M, Ducasse generation intercross population with a greatly reduced mating design A, Confais J, Compain J, Lapalu N. A small secreted protein in Zymosepto- for genetic studies in winter wheat. Front Plant Sci. 1825;2018:9. ria tritici is responsible for avirulence on wheat cultivars carrying the Stb6 70. Arraiano L, Balaam N, Fenwick P, Chapman C, Feuerhelm D, Howell P, resistance gene. New Phytol. 2017;214(2):619–31. Smith S, Widdowson J, Brown J. Contributions of disease resistance and 50. Saintenac C, Lee W-S, Cambon F, Rudd JJ, King RC, Marande W, Powers escape to the control of Septoria tritici blotch of wheat. Plant Pathol. SJ, Bergès H, Phillips AL, Uauy C, et al. Wheat receptor-kinase-like protein 2009;58(5):910–22. Stb6 controls gene-for-gene resistance to fungal pathogen Zymoseptoria 71. Tuberosa R: Mapping QTLs for Partial Resistance to Zymoseptoria tritici in tritici. Nat Genet. 2018;50(3):368–74. Durum Wheat. In: Plant and Animal Genome XXII Conference: 2014: Plant 51. Lin Y, Gnanesh BN, Chong J, Chen G, Beattie AD, Mitchell Fetch JW, and Animal Genome; 2014. Kutcher HR, Eckstein PE, Menzies JG, Jackson EW. A major quantitative 72. Ferjaoui S, M’Barek S, Bahri B, Slimane R, Hamza S. Identification of trait locus conferring adult plant partial resistance to crown rust in oat. resistance sources to septoria tritici blotch in old Tunisian durum wheat BMC Plant Biol. 2014;14(1):250. germplasm appliad for the analysis of the Zymoseptoria tritici-durum 52. Hou L, Chen X, Wang M, See DR, Chao S, Bulli P, Jing J. Mapping a Large wheat interaction. J Plant Pathol. 2015;97(3):471–81. Number of QTL for Durable Resistance to Stripe Rust in Winter Wheat 73. Singh RP, Hodson DP, Huerta-Espino J, Jin Y, Bhavani S, Njau P, Herrera- Druchamp Using SSR and SNP Markers. PLoS ONE. 2015;10(5):e0126794. Foessel S, Singh PK, Singh S, Govindan V. The emergence of Ug99 races 53. Liu Y, Qie Y, Li X, Wang M, Chen X. Genome-Wide Mapping of Quantita- of the stem rust fungus is a threat to world wheat production. Annu Rev tive Trait Loci Conferring All-Stage and High-Temperature Adult-Plant Phytopathol. 2011;49:465–81. Resistance to Stripe Rust in Spring Wheat Landrace PI 181410. Int J Mol 74. Mundt CC. Durable resistance: a key to sustainable management of Sci. 2020;21(2):478. pathogens and pests. Infect Genet Evol. 2014;27:446–55. 54. Tahir S, Zia I, Dilshad I, Fayyaz M, Noureen N, Farrakh S. Identification 75. Gautam T, Dhillon GS, Saripalli G, Singh VP, Prasad P, Kaur S, Chhuneja P, of stripe rust resistant genes and their validation in seedling and adult Sharma P, Balyan H, Gupta P. Marker-assisted pyramiding of genes/QTL plant glass house tests. Genet Resour Crop Evol. 2020;67(4):1025–36. for grain quality and rust resistance in wheat (Triticum aestivum L.). Mol 55. Han G, Liu S, Wang J, Jin Y, Zhou Y, Luo Q, Liu H, Zhao H, An D. Identifica- Breeding. 2020;40:1–14. tion of an elite wheat-rye T1RS· 1BL translocation line conferring high 76. Chartrain L, Brading P, Brown J. Presence of the Stb6 gene for resistance resistance to powdery mildew and stripe rust. Plant Dis. 2020. https:// doi. to Septoria tritici blotch (Mycosphaerella graminicola) in cultivars used in org/ 10. 1094/ PDIS- 02- 20- 0323- RE. wheat-breeding programmes worldwide. Plant Pathol. 2005;54(2):134–43. 56. Prasad P, Savadi S, Bhardwaj SC, Gupta PK. The progress of leaf rust 77. Chartrain L, Berry S, Brown J. Resistance of wheat line Kavkaz-K4500 L. research in wheat. Fungal Biol. 2020;124(6):537–50. https:// doi. org/ 10. 6. A. 4 to Septoria tritici blotch controlled by isolate-specific resistance 1016/j. funbio. 2020. 02. 013. genes. Phytopathology. 2005;95(6):664–71. 57. McCartney C, Brûlé-Babel A, Lamari L, Somers D. Chromosomal location 78. Chartrain L, Brading P, Widdowson J, Brown J. Partial resistance to Septo- of a race-specific resistance gene to Mycosphaerella graminicola in the ria tritici blotch (Mycosphaerella graminicola) in wheat cultivars Arina and spring wheat ST6. Theor Appl Genet. 2003;107(7):1181–6. Riband. Phytopathology. 2004;94(5):497–504. Ferjaoui et al. BMC Genomics (2022) 23:372 Page 20 of 20 79. Li ZK, Jiang XL, Peng T, Shi CL, Han SX, Tian B, Zhu ZL, Tian JC. Map- ping quantitative trait loci with additive effects and additive x additive epistatic interactions for biomass yield, grain yield, and straw yield using a doubled haploid population of wheat (Triticum aestivum L.). Genet Mol Res. 2014;13(1):1412-24. https:// doi. org/ 10. 4238/ 2014. 80. Farokhzadeh S, Fakheri BA, Nezhad NM, Tahmasebi S, Mirsoleimani A, McI- ntyre CL. Genetic control of some plant growth characteristics of bread wheat (Triticum aestivum L.) under aluminum stress. Genes Genom. 2020;42(3):245–61. 81. Li Z, Pinson S, Stansel J, Park W. Identification of quantitative trait loci (QTLs) for heading date and plant height in cultivated rice (Oryza sativa L.). Theor Appl Genet. 1995;91(2):374–81. 82. Lecomte L, Duffé P, Buret M, Servin B, Causse M. Marker-assisted intro- gression of five QTLs controlling fruit quality traits into three tomato lines revealed interactions between QTLs and genetic backgrounds. Theor Appl Genet. 2004;109(3):658–68. 83. Wang D, Zhu J, Li Z, Paterson A. Mapping QTLs with epistatic effects and QTL× environment interactions by mixed linear model approaches. Theor Appl Genet. 1999;99(7–8):1255–64. 84. Lin H, Yamamoto T, Sasaki T, Yano M. Characterization and detection of epistatic interactions of 3 QTLs, Hd1, Hd2, and Hd3, controlling heading date in rice using nearly isogenic lines. Theor Appl Genet. 2000;101(7):1021–8. 85. Zadoks JC, Chang TT, Konzak CF. A decimal code for the growth stages of cereals. Weed Res. 1974;14(6):415–21. 86. Kema G, Sayoud R, Annone J, Van Silfhout C. Genetic variation for virulence and resistance in the wheat-Mycosphaerella graminicola patho- system. II: analysis of interactions between pathogen isolates and host cultivars. Phytopathology. 1996;86(2):213–20. 87. Madden LV, Hughes G, Bosch F. The study of plant disease epidemics: American Phytopathological Society (APS Press); 2007. https:// doi. org/ 10. 1094/ 97808 90545 058. 88. Eyal Z, Brown MB. A quantitative method for estimating density of Sep- taria tritici pycnidia on wheat leaves. Phytopathology. 1976;66:11–4. 89. Kilian A, Wenzl P, Huttner E, Carling J, Xia L, Blois H, Caig V, Heller-Uszynska K, Jaccoud D, Hopper C, Aschenbrenner-Kilian M, Evers M, Peng K, Cayla C, Hok P, Uszynski G. Diversity arrays technology: a generic genome pro- filing technology on open platforms. Methods Mol Biol. 2012;888:67–89. https:// doi. org/ 10. 1007/ 978-1- 61779- 870-2_5. 90. Van Ooijen J. JoinMap 4. Software for the calculation of genetic linkage maps in experimental populations Kyazma BV, Wageningen, Netherlands: Scietific research; 2006. p. 33. 91. Kosambi DD. The estimation of map distances from recombination val- ues. Ann Eugenic. 1943;12(1):172–5. https:// doi. org/ 10. 1111/j. 1469- 1809. 1943. tb023 21.x. 92. Stam P. Construction of integrated genetic linkage maps by means of a new computer package: Join Map. Plant J. 1993;3(5):739–44. 93. Maccaferri M, Harris NS, Twardziok SO, Pasam RK, Gundlach H, Span- nagl M, Ormanbekova D, Lux T, Prade VM, Milner SG, et al. Durum wheat genome highlights past domestication signatures and future improve- ment targets. Nat Genet. 2019;51(5):885–95. 94. De Mendiburu F, Simon R: Agricolae-Ten years of an open source statisti- cal tool for experiments in breeding, agriculture and biology. In.: PeerJ PrePrints; 2015. 95. Chu C-G, Chao S, Friesen T, Faris J, Zhong S, Xu S. Identification of novel tan spot resistance QTLs using an SSR-based linkage map of tetraploid wheat. Mol Breeding. 2010;25(2):327–38. 96. Institute S: Base SAS 9.4 procedures guide: SAS Institute; 2015. Ready to submit your research ? Choose BMC and benefit from: 97. Meng L, Li H, Zhang L, Wang J. QTL IciMapping: integrated software for • fast, convenient online submission genetic linkage map construction and quantitative trait locus mapping in biparental populations. Crop J. 2015;3(3):269–83. • thorough peer review by experienced rese archers in your field 98. Voorrips R. MapChart: software for the graphical presentation of linkage • rapid publication on acceptance maps and QTLs. J Hered. 2002;93(1):77–8. • support for research data, including large and complex data types 99. Wickham H. Elegant graphics for data analysis. Media. 2009;35(211):10.1007. • gold Open Access which fosters wider collaboration and increased citations • maximum visibility for your research: over 100M website views per year Publisher’s Note At BMC, research is always in progress. Springer Nature remains neutral with regard to jurisdictional claims in pub- lished maps and institutional affiliations. Learn more biomedcentral.com/submissions