Lathyrus cicera L. (chickling pea) is a dual-purpose cool season legume and a source of protein for animals (Hanbury et al., 2000) and humans (Peña-Chocarro and Peña, 1999) nutrition. This species-genus belongs to the tribe Fabeae (syn. Vicieae) along with Vicia, Lens, Pisum, and Vavilovia (reviewed in Smýkal et al., 2011). Lathyrus sativus (grass pea) is the most closely related relative to L. cicera and these two Lathyrus species have a close phylogenetic relationship with pea (Pisum sativum), so close that there are suggestions that the genus Pisum should be included in the genus Lathyrus (Schaefer et al., 2012).

Lathyrus cicera owns important agronomic traits as resistance to biotic and abiotic stresses (Vaz Patto and Rubiales, 2014; Hammer et al., 2019; Lambein et al., 2019). It is therefore an attractive choice for sustainable feed and food production, mainly in more marginal environments, and because of pathogen sharing, it could act as a promising alternative source of resistance to related species, such as grass pea and pea (Vaz Patto et al., 2006; Vaz Patto and Rubiales, 2014; Hammer et al., 2019).

Fungal diseases are major constraints for yield stability in legumes (Rubiales et al., 2015; Martins et al., 2020). Rusts are among the most important diseases recorded in grain and forage legumes (Sillero et al., 2006; Rubiales et al., 2011). Several rust species can infect legumes, most of them belonging to the Uromyces genus (Rubiales et al., 2011). A good example is the wide host range biotrophic Uromyces pisi that infects species of Lathyrus, Pisum, Lens, and Vicia genera (Farr, D.F., and Rossman, A.Y. Fungal Databases) (Barilli et al., 2012). Chemical control of rust is possible (Emeran et al., 2011), but the use of host plant resistance is the most economical and ecological desired means of control (Rubiales et al., 2011). Due to the reduced selective pressure imposed on the pathogen, the use of plant partial resistance is a potentially more durable approach than complete resistance (McDonald and Linde, 2002; Niks and Rubiales, 2002). Indeed, rusts are among the pathogens with the highest risk of breaking down the effectiveness of major resistance genes (R-genes) due to their effective air dispersal and the coexistence of sexual and asexual reproduction cycles (McDonald and Linde, 2002).

Partial resistance not associated with host cell necrosis (hypersensitive reaction) is common in major grain legumes against rusts (Sillero et al., 2006) and was already identified in Lathyrus spp. (Vaz Patto and Rubiales, 2009; Vaz Patto et al., 2009) or in Pisum spp. (Barilli et al., 2009a,b,c, 2012) against U. pisi and U. vicia-fabae. In the particular case of Pisum fulvum, a wild relative of pea, one to three quantitative trait loci (QTLs) controlling partial resistance to U. pisi have been identified (Barilli et al., 2010, 2018). In contrast, little is known about the genetic control of L. cicera resistance to U. pisi.

In an Iberian collection of L. cicera accessions, microscopic and macroscopic variable levels of resistance were identified against U. pisi. Resistant accessions partially restricted the formation of haustoria, resulting in a high percentage of early aborted fungal colonies, a decreased number of haustoria per colony, and a reduced intercellular growth of infection hyphae compared to susceptible accessions (Vaz Patto et al., 2009). A segregating recombinant inbred line (RIL) population was developed from the cross of the most contrasting L. cicera accessions of this Iberian collection. The RILs were later on used for the development and refinement of the first L. cicera linkage map with transcriptome based SNPs and e-SSR markers (retrieved from the RIL parental lines leaf RNAseq response to U. pisi infection) as well as genotype-by-sequencing based markers (Santos et al., 2018, 2020). The mentioned L. cicera RIL parental lines transcriptomic study also highlighted upregulated genes in response to rust infection involved in hormone metabolism, cell wall degradation, secondary metabolism, ROS production, signalling and regulation of transcription of defence (Santos et al., 2018). In spite of these recent efforts, the causative variants controlling partial resistance to U. pisi remain elusive, leading to a poor understanding of the partial genetic resistance and molecular mechanisms of Lathyrus spp. against rust disease.

Classical approaches to unveil quantitative resistance genetic basis in non-model species include QTL mapping in segregating populations and syntenic analyses using model/related species to search for orthologous genes within the QTL detected regions. Whereas much progress has been made in plant QTL mapping controlling phenotypic natural variation, this approach has been hampered by the complex interrelation of genetic variants and expression regulators (Hansen et al., 2008; Wallace et al., 2013; Albert and Kruglyak, 2015). In species with no sequenced or no fully assembled genome, as the Lathyrus spp. (Emmrich et al., 2020), alternative approaches may add value over classical QTL mapping. Genetical-genomics is considered a very powerful tool to improve our knowledge on the genetic architecture of complex traits, including disease response (Lima et al., 2018; Fauteux et al., 2019). In this approach, transcript expression levels are treated as quantitative phenotypes in a segregating population and the genomic variants that influence expression levels of each transcript are identified by conventional QTL analysis (Li and Burmeister, 2005). The found genomic regions controlling gene expression are referred to as expression-QTLs (eQTLs). Previous studies have reported that distant or trans-eQTLs may explain a higher proportion of expression variance than local (at the same locus as the structural gene) or cis-eQTLs (Liu et al., 2017; Carrasco-Valenzuela et al., 2019; Fauteux et al., 2019; Miculan et al., 2021). Hotspots of trans-eQTL may act as key regulators of phenotypes, whereas cis-eQTLs display local gene expression regulation, with co-regulated gene clusters (Wang et al., 2018; Miculan et al., 2021). Several studies have been using a hybrid approach (pQTL and eQTL analyses) to better understand the gene networks underlying traits of interest in plants (Lima et al., 2018; Carrasco-Valenzuela et al., 2019; Fauteux et al., 2019; Miculan et al., 2021).

The main aim of this work was to elucidate the genetic basis and putative molecular strategies of chickling pea partial resistance against rust disease. Using L. cicera RIL population phenotypic response to U. pisi infection (mainly at the seedling stage, under controlled growth chamber, and complemented by an exploratory assay at the adult plant stage, under semi-controlled field conditions), in the expression analysis of genes related to rust resistance and genomic data, we performed a combined pQTL/eQTL linkage mapping analysis. This will allow prioritising candidate genes for validation and future use in precision breeding, and advancing our understanding on the molecular mechanisms underlying partial resistance to U. pisi in L. cicera.

Most of the rust resistance reactions described so far in cool season legumes are incomplete and so potentially durable, but in most cases, the genetic basis of these resistances is still largely unknown, hampering their use in precision breeding (Rubiales et al., 2015). In this study, L. cicera RIL population segregating for rust resistance was used to elucidate the genetic basis and putative molecular strategies of chickling pea partial resistance against rust disease. By integrating the phenotypic response to U. pisi infection (at two different developmental plant stages and growing conditions), with the expression analysis of genes related to rust resistance and genomic data, we performed a combined pQTL/eQTL analysis. This allowed us to prioritise candidate genes that after validation, may be relevant for resistance precision breeding and advance our understanding of the molecular mechanisms underlying partial resistance to U. pisi in L. cicera.

All analysed L. cicera RIL individuals showed a compatible reaction with U. pisi, at the seedling stage, under a controlled growth chamber, as well as at the adult plant stage under semi-controlled field conditions in one-season exploratory analysis, characterised by well-formed pustules with no associated chlorosis or necrosis (IT = 4). As previously described in L. cicera (Vaz Patto et al., 2009), the most resistant genotypes presently identified, showed a low DS despite this compatible infection type, confirming their partial resistant nature (low DS, high IT) (McDonald and Linde, 2002; Niks and Rubiales, 2002). In the well-studied wheat-rust pathosystem, partial resistance commonly has a polygenic nature (controlled by adult plant resistance or APR genes), being expressed only in adult plants (except under very specific conditions) by a reduced and slow pathogen growth, without hypersensitive response. In contrast, pathogen race- or strain-specific major resistance genes (R-genes), generally conferring complete resistance, mostly function from seedling to adult growth stages (Dakouri et al., 2013; Ellis et al., 2014; Zegeye et al., 2014). In this study, DS evaluated in seedlings and adult plants were weakly correlated, indicating that also in L. cicera as in cereals, a different genetic basis, with multiple “minor effect” genes (explaining 7.1–19% observed variance) depending on the plant developmental stage may be involved in resistance. Indeed, it is widely acknowledged that partial resistance to rust is better identified in polycyclic infections, and even further, clearer on the adult plant stage than on seedlings (Sillero et al., 2006; Barilli et al., 2009a). This also seems to be the case in L. cicera. By considering even just one-season exploratory experiment at the adult plant stage, under semi-controlled field conditions, we allowed polycyclic infection vs. the monocyclic infection to occurr in the well-replicated growth chamber experiments at the seedling stage. Therefore, also in L. cicera there might be valuable small “adult plant” factors, not seen in the accurate monocyclic infections’ evaluations in seedlings. These findings are of utmost importance as currently, breeders put a higher emphasis on the discovery, characterisation, and complementary use of genes for partial, more durable resistance than on using only major effect R-genes, due to their potential lack of durability (Ellis et al., 2014).

From all the highlighted pQTL candidate genes, only the PPR (candidate for the UpDSIV_chamber pQTL), DGK, and mitochondria carrier genes (UpDSIIa_field and UpDSIIc_field pQTLs, respectively) were also identified differentially expressed between the RIL parental lines in the previous RNAseq L. cicera-U. pisi transcriptomic study (Fold Change of 0.94, 1.48, and 2.56, respectively) (Santos et al., 2018). For the remaining pQTL candidate genes, it was not possible to identify any DEG in the mentioned transcriptomic study, possibly due to the short nucleotide sequence available from DartSeq markers (65 bp), hampering a precise alignment between marker and transcriptomic sequences. Thus, the discussion of putative function in the variation of DS L. cicera response against U. pisi will focus on PPR, DGK, and mitochondria carrier genes, as the most promising candidate resistance genes considering the available information.

Pentatricopeptide proteins have been identified as playing important roles in organellar RNA metabolism, organ development, and in abiotic and biotic stresses (Lurin et al., 2004; Koussevitzky et al., 2007; Laluk et al., 2011; Xing et al., 2018). The function of PPR proteins has been reported in plant response to necrotrophic fungi and pathogenic bacteria (Laluk et al., 2011; Park et al., 2014). More recently, genes from the PPR gene family were identified underlying QTLs for partial resistance to the biotrophic Erysiphe pisi and E. trifolii powdery mildew pathogens in L. cicera (Santos et al., 2020) and in L. sativus (Martins et al., 2022⁵). All these studies reporting PPR as a gene involved in Lathyrus spp. response against biotrophic pathogens, instigate further analysis to the function of genes encoding for members of the PPR gene family on chickling pea response to rust fungi.

Concerning DGK genes, some studies have revealed their involvement in the modulation of plant growth and adaptation to both biotic and abiotic stresses. DGKs are the main moderators of lipid signalling in plants, and this enzymatic activity is increased upon pathogen infection or elicitor treatment in different species (Foka et al., 2020). In pea, for instance, the inhibition of DGK activity promoted an elicitor-mediated accumulation of the phytoalexin pisatin, inducing phenylalanine ammonia-lyase expression (Toyoda et al., 2000). In the present study, UpDSIIa_field pQTL (for which the DGK gene was proposed) showed a positive additive effect indicating that alleles for increased DS came from susceptible L. cicera parental line, and so, resistance may be caused by inhibition of DGK activity, as suggested for pea. However, in the previous L. cicera transcriptomic study (Santos et al., 2018), the partial resistant parental line showed higher expression of this DGK gene than the susceptible line (Fold Change of 1.48). Therefore, the involvement of DGK genes in the L. cicera response against U. pisi remains unclear.

The mitochondrial carrier protein family is over-represented among the stress-responsive genes, suggesting that stress induces altered needs for metabolite transport across the mitochondrial inner membrane (Van Aken et al., 2009). Mitochondrial carriers are highly expressed in stress conditions, such as application of cadmium or auxin, exposure to cold, and induction of cell death. Abscisic acid application, on the other hand, decreased the expression of some mitochondrial carriers (Millar and Heazlewood, 2003). The potential role of L. cicera mitochondrial carriers against U. pisi infection here indicated is in accordance with the previous transcriptomic study (Santos et al., 2018).

To better understand the complex interrelation of genetic variants and expression regulators, the identification of eQTLs for some genes related to L. cicera resistance against U. pisi was integrated in this study. From the initial nine genes selected from the transcriptomic data obtained for L. cicera RIL parental lines (BGE023542 and BGE008277) inoculated with U. pisi (Santos et al., 2018), Gluc (Glycosyl hydrolase family 17) was the only analysed gene which RILs’ expression levels correlated with RILs’ DS scorings in the present study. Glycosyl hydrolases family 17 serves diverse roles in plant defence and development, since it comprises degrading enzymes of 1,3-β-glucan polysaccharides found in the cell wall matrix of plants and fungi (Thomas et al., 2000; Gaudioso-Pedraza and Benitez-Alfonso, 2014). The significant negative correlation (P < 0.01) observed between DS at adult plant stage under semi-controlled field conditions and Gluc relative expression support that this gene expression may increase L. cicera resistance against U. pisi. Indeed, the partial resistant BGE023542 parental line showed higher transcript abundance than BGE008277 for this gene, suggesting that its expression may increase resistance to the pathogen. Positive and negative additive effects were found in the eQTLs for Gluc, indicating that both parental lines may harbour alleles for L. cicera resistance to rust.

Two eQTLs were detected for Gluc RIL expression variation. Strong eQTLs are typically cis-regulated (Võsa et al., 2021), and this was also observed in the present study as GlucIVa eQTL (the only detected cis-eQTL) showed the highest LOD score (9.29) and the highest percentage of explained expression variation (32.2%). Nevertheless, the influence of distal regulations is not trivial and generally numerous plant genes are controlled by distant-factors (Swanson-Wagner et al., 2009; Wang et al., 2010; Hammond et al., 2011; Cubillos et al., 2012). Indeed, a distant-eQTL (trans-eQTL) was detected in the present study associated with the expression variation of Gluc but explaining a smaller percentage of variation.

The candidate genes for the GlucIVb trans-eQTL were proposed to be a SCO1 protein and a rust resistance kinase Lr10. SCO1 protein has been described as playing a role in cellular copper homeostasis and mitochondrial redox signalling (Attallah et al., 2011). Regarding the Lr10 leaf rust resistance gene, in wheat encodes a coiled coil–nucleotide-binding site–leucine-rich repeat (CC–NBS–LRR) (Loutre et al., 2009). The majority of disease resistance genes (R-genes) isolated from plants, conferring resistance to pathogens, encode proteins containing an NBS–LRR domain, used for pathogen perception triggering host response (Jones and Dangl, 2006). To the best of our knowledge, no direct relation between glycosyl hydrolases and SCO1 and Lr10 genes were described so far. Therefore, further studies are crucial to unveil if and how the rust resistance kinase Lr10 and SCO1 may regulate Glycosyl hydrolases family 17 members in response to rust infection in L. cicera.

Noteworthy, the function and position of genes underlying the identified pQTLs and eQTLs should be confirmed when the scaffolding of the assembly of the L. sativus genome (Emmrich et al., 2020) to the pseudochromosome level and gene annotation becomes completed.

Disease resistance genes are commonly clustered in genetic regions, conferring resistance to different pathogens and/or to different races of the same pathogen (Michelmore and Meyers, 1998; Loridon et al., 2005). This has been previously reported in Pisum spp. for fungal and oomycete pathogens (Barilli et al., 2018; Jha et al., 2021; Wu et al., 2021). In L. cicera, a similar situation seems to occur to a certain extent. One of the recently identified QTLs for partial resistance to E. pisi (EpDSIV in Santos et al., 2020) locates within a close distance to pQTLs and eQTLs here identified for partial resistance to U. pisi: about 13.8 cM to the UpDSIV_chamber, 22 cM to the UpDSIVa_field and GlucIVa, and 17.5 cM to the GlucIVb eQTL. Therefore, a QTL hotspot for partial disease resistance to different biotrophic fungi (U. pisi and E. pisi) may be suggested in L. cicera LGIV (considering UpDSIV_chamber, UpDSIVa_field, GlucIVa, and GlucIVb and EpDSIV; Santos et al., 2020).

In the present study, only the Gluc and Pi49 genes showed significant differences in gene expression by RT-qPCR between the two L. cicera parental lines. A motive that might be influencing our results is that different biological material was used on the RNAseq (a pool of equally mixed RNA individually extracted from 24 biological replicates per genotype) vs. the RT-qPCR (3 biological replicates individual RNA samples per genotype). Some residual heterozygosity potentially present in the RIL parental lines might have been exposed during different processing and contribute to the different fold changes (and consequently significant differences) obtained between gene expression of the RIL parental lines using RNA-seq and RT-qPCR techniques. The limited number of DEG from the RIL parental lines RNAseq data screened in the RIL individuals hampered the scope of the integration approach for prioritising candidate genes for chickling pea rust resistance. Nevertheless, the information obtained from the eQTL linkage mapping allowed a better understanding of the complex interrelation of genetic variants and expression regulators of the DEG Gluc.

Different pQTLs for U. pisi partial resistance were identified in L. cicera RIL, suggesting that a different genetic control may be involved in different stages of plant development. The integration of these linkage mapping results and the previously obtained RIL parental lines transcriptomic data (Santos et al., 2018), helped to prioritise QTLs candidate genes (if DEG between the parental lines) as a cross-validation approach. This integration has materialised also by the selection of à priori genes related to resistance to rust from the RIL parental lines transcriptomic study (Santos et al., 2018) to “validate” at the RIL level by RT-qPCR. Taking all the above into consideration, the L. cicera response to U. pisi candidate gene prioritised list for future validation and use in precision breeding is constituted by PPR, DGK, Mitochondrial carrier, and Gluc genes. Furthermore, the presence of a putative hotspot of resistance-related genes in the L. cicera LGIV is suggested. Candidate genes underlying the identified pQTLs/eQTLs hotspot will be useful for a better understanding of the complex interrelation of genetic variants and regulators of expression related to the partial resistance of L. cicera against rust.

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

CS, DCM, and MJG-B contributed to phenotyping and methodology. CS contributed to gene expression analysis, QTL and eQTL analysis, software, and data curation. DR contributed to RIL population development. MCVP and CS contributed to conceptualisation, formal analysis, and writing—original draft preparation. CS, DCM, DR, and MCVP contributed to writing, reviewing, and editing. MCVP and DR contributed to funding acquisition, supervision, and resources. MCVP contributed to project administration. All authors have read and agreed to the published version of the manuscript.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.