2. What is the NLR birth-and-death process?

The Red Queen hypothesis posits that rapidly evolving pathogens that make use of multiple, redundant and overlapping mechanisms to colonize their hosts lead to continuous cycles of reciprocal adaptation, without either the host or the pathogen prevailing in the long term. As the host evolves resistance, the pathogen jettisons genes for molecules recognized by the host, deploying instead alternative molecules for successful infection. In response, the host evolves new means that can suppress the pathogen until the pathogen again changes its genetic makeup. The defining feature is negative frequency-dependent selection (FDS) generating allele frequencies that fluctuate over time. In the Red Queen scenario (or trench warfare), genetic diversity is maintained in both host and pathogen (Box 2). In contrast, under the arms race situation, both host and pathogens rely on continuous innovation rather than recycling of existing variation. While the timescales over which alleles persist in a population differ under the two scenarios (Box 2; Brown & Tellier, Reference Brown and Tellier2011; Tellier et al., Reference Tellier, Moreno-Gámez and Stephan2014), common to both situations is that new genetic diversity needs to be generated. This variability can be generated by various mechanisms (Figure 1): 1) Mutation at the nucleotide level (point mutation and small insertion–deletions). 2) Recombination reshuffling variability between homologous copies of NLRs, with unequal recombination reshuffling variation between recently diverged copies in the same gene cluster, and gene conversion doing the same, but between more distant paralogous copies. 3) Gene duplication or deletion occurring by DNA break and repair, sometimes associated with the activity of transposable elements (TEs), generating presence/absence variation (PAV) at single genes or duplicated copies. 4) Unequal recombination between sets of paralogous copies further enhancing CNV and driving the expansion of the gene family. Variation is removed primarily by genetic drift but also by costs such as autoimmunity, which can be the consequence of NLR gene truncations, NLR overexpression or inappropriate interactions between NLR and other immune genes (Freh et al., Reference Freh, Gao, Petersen and Panstruga2021).

Three types of selection occur at a given NLR locus (Box 2). Positive selection, due to negative indirect FDS (Tellier & Brown, Reference Tellier and Brown2007), generates selective sweeps when an NLR variant provides an advantage, for example, by activating effective defence against a common pathogen. FDS is indirect, as the fitness of host alleles depends on the parasite allele frequencies and vice versa, and negative, as rare alleles have an advantage against common alleles. This yields over several cycles of selection the so-called arms race dynamics (Bergelson et al., Reference Bergelson, Kreitman, Stahl and Tian2001; Woolhouse et al., Reference Woolhouse, Webster, Domingo, Charlesworth and Levin2002). Negative selection (or purifying selection) removes deleterious variants, for example, if mutations at a given NLR locus incur a cost of resistance (Freh et al., Reference Freh, Gao, Petersen and Panstruga2021; Tian et al., Reference Tian, Araki, Stahl, Bergelson and Kreitman2002) or activate auto-immunity (Freh et al., Reference Freh, Gao, Petersen and Panstruga2021). Balancing selection, due to stable coevolutionary cycles (fluctuating selection) driven by negative indirect FDS in large populations, can generate the so-called trench warfare dynamics at the NLR locus (Gladieux et al., Reference Gladieux, van Oosterhout, Fairhead, Jouet, Ortiz, Ravel, Shrestha, Frouin, He, Zhu, Morel, Huang, Kroj and Jones2024; Stahl et al., Reference Stahl, Dwyer, Mauricio, Kreitman and Bergelson1999). While costs of resistance and virulence are necessary for balancing selection to occur, these are not sufficient and additional characteristics of natural systems (negative direct FDS; Tellier & Brown, Reference Tellier and Brown2007; Box 2) are required, such as spatial and temporal heterogeneous selection (review in Brown and Tellier, Reference Brown and Tellier2011), diffuse coevolution (Karasov et al., Reference Karasov, Kniskern, Gao, DeYoung, Ding, Dubiella, Lastra, Nallu, Roux, Innes, Barrett, Hudson and Bergelson2014) or epidemiological feedbacks (Boots et al., Reference Boots, White, Best and Bowers2014). Negative direct FDS occurs when the fitness of an allele depends negatively on its own frequency (Tellier & Brown, Reference Tellier and Brown2007), so that balancing selection can occur via the maintenance of long-lived polymorphism in host and pathogen populations (against the effect of genetic drift; Tellier et al., Reference Tellier, Moreno-Gámez and Stephan2014).

These evolutionary processes and mechanisms drive the birth-and-death process of NLR genes, thereby shaping variation in the composition of the NLR gene family between individuals, populations and species. When the birth of an NLR occurs via gene duplication, the fate of each paralogous copy depends on these neutral and selective processes (e.g., Innan & Kondrashov, Reference Innan and Kondrashov2010). Balancing selection at one NLR locus for several resistance alleles can favour gene duplication and the subsequent subfunctionalization of each copy to recognize different pathogen effectors or effector alleles (Van der Hoorn et al., Reference Van der Hoorn, De Wit and Joosten2002). Neo-subfunctionalization can also occur through selection at one copy for recognition of new pathogen variants (arising by mutation, recombination or gene conversion). In both cases, diversifying selection occurs and is observed across the different NLR paralogous gene copies (Innan & Kondrashov, Reference Innan and Kondrashov2010; Otto et al., Reference Otto, Zheng and Wiehe2022). Pseudogenization (by mutation, recombination or gene conversion) may inactivate one NLR copy, which becomes free to evolve a new function or may be lost. Ultimately, an NLR copy undergoes gene death either by deletion or by non-functional pseudogenization following genetic drift or selection against deleterious copies/alleles (Michelmore & Meyers, Reference Michelmore and Meyers1998). Copies/alleles become deleterious when they exhibit high fitness costs that are not compensated by providing ad hoc resistance to a pathogen (Stahl et al., Reference Stahl, Dwyer, Mauricio, Kreitman and Bergelson1999; Tellier & Brown, Reference Tellier and Brown2007), or when they inadvertently trigger autoimmune responses (Freh et al., Reference Freh, Gao, Petersen and Panstruga2021; Wan et al., Reference Wan, Kim, Castel, Charoennit and Chae2021). In this framework, the number of NLR gene copies in the genome and their polymorphism (at the CNV levels, a.k.a. the pan-NLRome) depends thus on the balance between the rates of gene birth and death. For the plant host to cope with the fast adaptation of their pathogens, it is classically expected that higher mutation, recombination and gene duplication rates may be advantageous in generating more variable NLR genes, which in turn would promote higher levels of CNV reflecting the NLR gene turnover rate.

3. How to identify structural NLR variation in the genome?

The first task to study NLR diversity and dynamics is to identify NLR variation. While whole-genome short read sequencing or NLR Ren-seq coupled with short or long reads (e.g., SMRT-RenSeq; Wu et al., Reference Wu, Jones, Lin, Kawasaki and Yoshioka2026) provides some insights into NLR CNVs (e.g., Gladieux et al., Reference Gladieux, van Oosterhout, Fairhead, Jouet, Ortiz, Ravel, Shrestha, Frouin, He, Zhu, Morel, Huang, Kroj and Jones2024; Li et al., Reference Li, Liu, Zhang, Chen and Shao2025; Seong et al., Reference Seong, Seo, Witek, Li and Staskawicz2020; Silva-Arias et al., Reference Silva-Arias, Gagnon, Hembrom, Fastner, Khan, Stam and Tellier2025; Stam et al., Reference Stam, Nosenko, Hörger, Stephan, Seidel, Kuhn, Haberer and Tellier2019; Van de Weyer et al., Reference Van de Weyer, Monteiro, Furzer, Nishimura, Cevik, Witek, Jones, Dangl, Weigel and Bemm2019), such analyses suffer from mapping uncertainty and reference biases. The clustering of NLR genes in the genome and the quantification of the family members and the nucleotide diversity can be assessed based on their NB-ARC domain (Prigozhin & Krasileva, Reference Prigozhin and Krasileva2021; Silva-Arias et al., Reference Silva-Arias, Gagnon, Hembrom, Fastner, Khan, Stam and Tellier2025); but the resolution of CNV remains error-prone (Silva-Arias et al., Reference Silva-Arias, Gagnon, Hembrom, Fastner, Khan, Stam and Tellier2025). Long-read assemblies are much better at disentangling paralogous sequences, with the caveat that near-identical copies, especially in heterozygous and polyploid genomes, continue to be challenging (Cheng et al., Reference Cheng, Concepcion, Feng, Zhang and Li2021; Rautiainen et al., Reference Rautiainen, Nurk, Walenz, Logsdon, Porubsky, Rhie, Eichler, Phillippy and Koren2023; Rhie et al., Reference Rhie, McCarthy, Fedrigo, Damas, Formenti, Koren, Uliano-Silva, Chow, Fungtammasan, Kim, Lee, Ko, Chaisson, Gedman, Cantin, Thibaud-Nissen, Haggerty, Bista, Smith and Jarvis2021). The use of new bioinformatic methods for pangenomes including massive pair-wise comparisons (Igolkina et al., Reference Igolkina, Vorbrugg, Rabanal, Liu, Ashkenazy, Kornienko, Fitz, Collenberg, Kubica, Morales, Jaegle, Wrightsman, Voloshin, Bezlepsky, Llaca, Nizhynska, Reichardt, Bezrukov, Lanz and Nordborg2025) and genome graphs (Garrison et al., Reference Garrison, Guarracino, Heumos, Villani, Bao, Tattini, Hagmann, Vorbrugg, Marco-Sola, Kubica, Ashbrook, Thorell, Rusholme-Pilcher, Liti, Rudbeck, Golicz, Nahnsen, Yang, Mwaniki and Prins2024) has become instrumental in accounting for genome rearrangements occurring in the NLR gene family (see the recent review by Bao & Weigel, Reference Bao and Weigel2025). However, while genome graphs are useful representations, they describe only the outcome of evolution and do not allow directly for the inference of the causal mechanisms and evolutionary forces underpinning a given graph. We suggest that an important aim of future research must be to integrate genome graphs with population genomics theory – which should also improve parameter choices that have to be made when building genome graphs. Furthermore, the reconstruction of gene ancestry within gene clusters and more generally the question of orthology in gene families need to be developed.

An example of how long-read assemblies improve our understanding of the genomic landscape of NLR genes comes from a comparison of 17 A. thaliana accessions (Teasdale et al., Reference Teasdale, Murray, Collenberg, Contreras-Garrido, Schlegel, van Ess, Jüttner, Lanz, Deusch, Fitz, Mencia, van Velthoven, Drost, Weigel and Shirsekar2025). This study introduced the concept of an NLR neighbourhood, defined as a region that contains at least one NLR in at least one of the accessions. It revealed a higher-than-expected density of TEs, which may indicate a role of TEs in gene duplication and/or deletion. So far, this is, however, merely a correlation: are NLR loci more prone to attracting TE insertions or allowing for TE mobilization, or is TE fixation driven by hitch-hiking (see below linked selection) due to evolutionary (selective) changes at NLR genes? Furthermore, NLR gene families appear varied in their diversity (Teasdale et al., Reference Teasdale, Murray, Collenberg, Contreras-Garrido, Schlegel, van Ess, Jüttner, Lanz, Deusch, Fitz, Mencia, van Velthoven, Drost, Weigel and Shirsekar2025), and there do not appear to be single metrics that fully capture NLR diversity. For example, the RPP13 gene is highly diverse in sequence but not in isoforms or copy number, whereas the DAR4 gene exhibits the opposite pattern, while the RPP4 gene consistently shows extreme values for most metrics. This suggests that certain NLR genes and families have perhaps primary roles in dealing with slowly changing pathogen threats, while others have higher potential for adaptability in the face of fast evolving pathogens (Bakker et al., Reference Bakker, Toomajian, Kreitman and Bergelson2006; Laflamme et al., Reference Laflamme, Dillon, Martel, Almeida, Desveaux and Guttman2020; Prigozhin & Krasileva, Reference Prigozhin and Krasileva2021; Sutherland et al., Reference Sutherland, Prigozhin, Grey Monroe and Krasileva2024 in A. thaliana; Stam et al., Reference Stam, Nosenko, Hörger, Stephan, Seidel, Kuhn, Haberer and Tellier2019 in Solanum chilense). A central aim of the present study is to provide a roadmap for revealing and quantifying this heterogeneity in NLR gene evolutionary trajectories and the underlying mechanisms.

4. Is the high genetic diversity of NLR genes due to a higher mutation rate?

It is usually difficult to estimate the mutation rate from polymorphism data. As a proxy for it, the observed substitution rate is commonly used. Based on phylogenetic reconstruction of A. thaliana NLR families from 62 long-read assemblies of NLR loci, Sutherland et al. (Reference Sutherland, Prigozhin, Grey Monroe and Krasileva2024) define two categories of NLR genes: 1) highly variable NLR genes (hvNLRs) encoding proteins with unusually high amino acid diversity and 2) low variability NLR genes (non-hvNLRs), which presumably experience stronger purifying selection. Although the exact definition of high and low is arbitrary and based on a statistic sensitive to sample size, the observed categories seem still to capture different and relevant genomic signatures here. The authors find that hvNLRs exhibit higher expression levels, lower gene body cytosine methylation and closer proximity to TEs than non-hvNLRs. As expected, hvNLRs have elevated synonymous and nonsynonymous nucleotide diversity (i.e., $\pi_\mathrm{S} $ and $\pi_\mathrm{N} $ , respectively, see Table 2), higher Tajima’s D values (see Table 2), and are preferentially found in chromatin states associated with an increased probability of per-site mutation (based on the estimates in Monroe et al., Reference Monroe, Srikant, Carbonell-Bejerano, Becker, Lensink, Exposito-Alonso, Klein, Hildebrandt, Neumann, Kliebenstein, Weng, Imbert, Ågren, Rutter, Fenster and Weigel2022). While purifying (negative) selection to maintain a specific function is likely weaker at hvNLRs than non-hvNLRs (along with the associated background selection, see below), these analyses do not tell us whether balancing selection or positive selection drives hvNLR diversity, as relaxed purifying selection (and weak linked selection) or neutral evolution would generate similar footprints of nucleotide diversity. An analysis of the paralogous assignments and clustering of hvNLRs is required (for which genome graphs could help, see above) to assess the duplication (paralogous) status of the hvNLRs, as duplicated genes should more easily accumulate mutations due to (neo-)subfunctionalization (Innan & Kondrashov, Reference Innan and Kondrashov2010). In addition, as more complete genomes become available, it will become possible to ask whether hvNLR loci feature individual alleles (haplotypes) that undergo purifying selection, that is, whether diversity (and the speed of evolution) at these loci is discrete or continuously distributed. We also suggest that quantitative analyses of the selection strength can now be performed both for entire NLR subfamilies (e.g., TIR-NLRs, CC-NLRs and CCr-NLRs) or categories (hvNLRS vs. non-hvNLRs, single genes or clustered genes, or canonical NLRs and truncated NLRs) via the estimation of the distribution of fitness effects (DFE; Keightley & Eyre-Walker, Reference Keightley and Eyre-Walker2007) and the percentage of adaptive mutations (method in Sendrowski and Bataillon, Reference Sendrowski and Bataillon2024 and applications in Grandaubert et al., Reference Grandaubert, Dutheil and Stukenbrock2019). While DNA mutations are generally considered to be non-targeted with respect to their genomic location and their functional consequence, there are rare cases in which targeted rearrangements of the DNA sequence occur, which enhance the probability for adaptive changes (review in Hanlon et al., Reference Hanlon, Cagan and den Akker2025). Nonetheless, based on published results (Prigozhin & Krasileva, Reference Prigozhin and Krasileva2021; Sutherland et al., Reference Sutherland, Prigozhin, Grey Monroe and Krasileva2024; Van de Weyer et al., Reference Van de Weyer, Monteiro, Furzer, Nishimura, Cevik, Witek, Jones, Dangl, Weigel and Bemm2019), we cannot yet conclude whether or not the latter applies to NLR genes. We thus caution against the general claim that NLR genes exhibit higher per-site mutation rates than the rest of the genome.

Table 2

Summary of population genomics analyses to interrogate the evolution of NLR genes and clusters

Table 2 Long description

The table consists of nine rows detailing various genomic analyses.

* Genetic diversity for SNP s or structural variants: Requires genomic data for several individuals. Measures include pairwise nucleotide diversity (pi) and Watterson Theta (Theta w). It interprets diversity as a function of effective population size (N sub e) and recombination.

* Recombination: Measures Linkage Disequilibrium (L D) as a correlation coefficient (r super 2) between loci. It depends on N sub e and per-site recombination rates.

* Neutrality tests: Requires S N P calls with outgroup species. Measures include Tajima’s D, Fu and Li’s D, and Fay and Wu’s H to identify deviations in the site frequency spectrum (S F S).

* Gene conversion: Requires paralog genomic data. Uses joint S F S for both paralogs to infer gene conversion and assess selection.

* Neutrality tests with outgroup: Requires coding region S N P calls. Measures include pi N, pi S, the ratio pi N / pi S, McDonald-Kreitman (M K) test, Hudson-Kreitman-Aguade (H K A) test, and omega. It assesses selective or relaxed constraints.

* Selection across phylogeny: Requires genomic data for several species. Measures K A (d N), K S (d S), and the ratio K A / K S (d N / d S) to indicate long-term deviations from neutral evolution.

* Genetic differentiation between populations: Requires S N P or structural variant calls. Measures include F S T, Jost’s D, d X Y, and V S T to identify local adaptation or balancing selection.

* Inference of neutral demographic and spatial processes: Based on nucleotide diversity, L D, and population differentiation to predict thresholds for outlier loci.

* Genome scans for positive selection: Uses nucleotide diversity and L D per genomic region to detect genes under positive selection based on demographic inference.

Note: To keep the list of references manageable, we refrain from citing all publications and method papers referenced in the table. We direct readers to textbooks (Charlesworth & Charlesworth, Reference Charlesworth and Charlesworth2010; Hartl & Clark, Reference Hartl and Clark2023; Wiuf & Hein, Reference Wiuf and Hein1997) and the reference manual of state-of-the-art population genomics methods (Dutheil, Reference Dutheil2020). Recommendations to perform population genomics analyses rigorously are found in Bourgeois and Warren (Reference Bourgeois and Warren2021) and Johri et al. (Reference Johri, Aquadro, Beaumont, Charlesworth, Excoffier, Eyre-Walker, Keightley, Lynch, McVean, Payseur, Pfeifer, Stephan and Jensen2021) and in ad hoc chapters of Dutheil (Reference Dutheil2020). Software recommendations can be found in Bourgeois and Warren (Reference Bourgeois and Warren2021) and we only indicate here recent additions: ngsLD (Fox et al., Reference Fox, Wright, Fumagalli and Vieira2019), eSMC2 (Sellinger et al., Reference Sellinger, Patrick, Awad and Tellier2021) and SM $\!\! \beta $ C (Korfmann et al., Reference Korfmann, Sellinger, Freund, Fumagalli and Tellier2024).

5. Are meiotic recombination and gene conversion specifically elevated at NLR loci?

Following the tenets of the Red Queen hypothesis (Barton & Charlesworth, Reference Barton and Charlesworth1998; Van Valen, Reference Van Valen1973), sexual reproduction, and especially meiotic recombination (introducing chromosome arm exchanges and gene conversions), is expected to be an important evolutionary mechanism for generating new genotypes in hosts (and pathogens). The shuffling of existing nucleotide diversity generates novel multi-locus haplotypes with every generation. Moreover, if recombination occurs within clusters or genes, it can even create new NLR cluster or allelic variants. The situation for NLRs remains, however, unclear, as recent insights challenge the idea that increased meiotic recombination is underlying the high diversity of NLR genes (or gene clusters) and thereby drive accelerated evolution of NLR genes. Studies in A. thaliana have shown that the frequency of meiotic crossovers in NLR regions does not deviate greatly from genome-wide average (Choi et al., Reference Choi, Reinhard, Serra, Ziolkowski, Underwood, Zhao, Hardcastle, Yelina, Griffin, Jackson, Mézard, McVean, Copenhaver and Henderson2016). Although meiotic crossovers occur in NLR gene clusters, occasionally even forming recombination hotspots, the overall crossover rates within these clusters typically align with genome-wide averages (Choi et al., Reference Choi, Reinhard, Serra, Ziolkowski, Underwood, Zhao, Hardcastle, Yelina, Griffin, Jackson, Mézard, McVean, Copenhaver and Henderson2016). Notably, structural rearrangements, which are hallmarks of complex NLR gene clusters (Jiao et al., Reference Jiao, Patel, Klasen, Liu, Pecinkova, Ferrand, Gy, Camilleri, Effgen, Koornneef, Pecinka, Loudet and Schneeberger2021), inhibit recombination, resulting in extensive recombination suppression across these regions. The lack of elevated crossover frequency, coupled with the likely suppression of crossovers due to massive structural variation, suggests that crossover recombination alone may not be the main driver behind the rapid diversification of NLR genes in A. thaliana. Alternative mechanisms involving gene duplication and CNV as prerequisites, such as (meiotic or somatic) gene conversion (Hörger et al., Reference Hörger, Ilyas, Stephan, Tellier, Renier, der Hoorn and Rose2012 in wild tomato) between NLR paralogs or unequal recombination (e.g., as seen at MHC loci in animals; Otto et al., Reference Otto, Zheng and Wiehe2022), may be additional drivers of NLR gene evolution. Yet, their occurrence and frequency have to be quantified with high-quality genome assemblies, preferably using a collection of genomes from both more closely and more distantly related accessions, since excessive structural variation will reduce opportunities for unequal crossovers. We note that population genetics models of gene conversion (Innan, Reference Innan2003; Thornton, Reference Thornton2007) and unequal recombination (Otto et al., Reference Otto, Zheng and Wiehe2022) already exist and can be used to generate predicted polymorphism patterns (SNPs) across paralogs, which in turn can be used to infer the rate of gene conversion (Hörger et al., Reference Hörger, Ilyas, Stephan, Tellier, Renier, der Hoorn and Rose2012; Table 2). In a recent review, Metzger et al. (Reference Metzger, Abu-Awad, Sellinger and Tellier2026) summarized available sets of genome-wide statistics and inference methods (Sequential Markovian Coalescent, Approximate Bayesian Computation and Deep Learning) that can be applied to long-read pangenomes to infer past demographic history and past changes in life-history traits. We suggest training and applying these methods to infer rates of gene conversion (Thornton, Reference Thornton2007) and unequal recombination (Otto et al., Reference Otto, Zheng and Wiehe2022), and other neutral processes shaping NLR polymorphism data. While available data suggest that crossover recombination is not likely to be the sole driver behind the rapid diversification of NLR genes, such new inference methods applied to pangenome polymorphism data would allow us to assess the importance of meiotic or somatic gene conversion and unequal recombination in species with different life-history traits and genome characteristics.

6. How many different types of selection act at NLR genes?

Since the discovery of plant resistance loci encoded by NLR genes, numerous studies have utilized polymorphism data to evaluate the occurrence of positive selection (arms race) versus balancing selection (trench warfare dynamics) both at functionally defined loci as well as NLRs genome wide. To keep the list of references to a reasonable number, we mention here only a few systems and studies. Early efforts based on PCR amplification and Sanger sequencing have been conducted in A. thaliana (Bakker et al., Reference Bakker, Toomajian, Kreitman and Bergelson2006; Bergelson et al., Reference Bergelson, Kreitman, Stahl and Tian2001; Stahl et al., Reference Stahl, Dwyer, Mauricio, Kreitman and Bergelson1999), wild tomato species (Hörger et al., Reference Hörger, Ilyas, Stephan, Tellier, Renier, der Hoorn and Rose2012; Rose et al., Reference Rose, Michelmore and Langley2007) or common bean (De Meaux et al., Reference De Meaux, Cattan-Toupance, Lavigne, Langin and Neema2003). These initial studies found clear footprints of selection only at a few loci. With the advent of new sequencing technologies and improved methods for carrying out selection scans (see Table 2), this picture has not significantly changed (see Stam et al., Reference Stam, Nosenko, Hörger, Stephan, Seidel, Kuhn, Haberer and Tellier2019; Wei et al., Reference Wei, Silva-Arias and Tellier2023; Wu et al., Reference Wu, Han, Chen, Chen, Zou, Li, Yong-Chao and Guo2017). As predicted by theory (Brown & Tellier, Reference Brown and Tellier2011; Tellier et al., Reference Tellier, Moreno-Gámez and Stephan2014), arms races may not be that common in natural systems, which is very different from agricultural crop systems. Looking at this from the other side, it was found that pathogen population of wild plants might indeed be stable over very long time scales, with diversity among them driven as much by pathogen–pathogen antagonism as by plant–pathogen interactions (Backman et al., Reference Backman, Latorre, Symeonidi, Muszyński, Bleak, Eads, Martinez-Koury, Som, Hawks, Gloss, Belnap, Manuel, Deutschbauer, Bergelson, Azadi, Burbano and Karasov2024; Shalev et al., Reference Shalev, Karasov, Lundberg, Ashkenazy, Ayutthaya and Weigel2022). However, typical polymorphism footprints of balancing selection are more subtle to uncover (than selective sweeps), as they depend on the local recombination rate and the time during which selection acted to maintain variation (Charlesworth & Jensen, Reference Charlesworth and Jensen2021; Stahl et al., Reference Stahl, Dwyer, Mauricio, Kreitman and Bergelson1999; Tellier et al., Reference Tellier, Moreno-Gámez and Stephan2014). By a detailed population genetic analysis of the NLR genes of modern rice varieties and rice landraces, Gladieux et al. (Reference Gladieux, van Oosterhout, Fairhead, Jouet, Ortiz, Ravel, Shrestha, Frouin, He, Zhu, Morel, Huang, Kroj and Jones2024) demonstrate the occurrence of footprints of balancing selection at some NLR loci. However, in general, balancing selection is more easily uncovered if it underlies variation that is maintained across species after speciation. On the flip side, such signals can be difficult to detect, because recombination over a long time will reduce the size of the segments generating the signals; in extreme cases, only a single nucleotide polymorphism might be left, which can be difficult to distinguish from technical artefacts and recurrent mutations. A special situation was described for an outcrossing species, Capsella grandiflora, and a sister species, C. rubella, that had only very recently diverged from it, such that shared blocks of polymorphism were still large (Bachmann et al., Reference Bachmann, Tedder, Laenen, Fracassetti, Désamoré, Lafon-Placette, Steige, Callot, Marande, Neuffer, Bergès, Köhler, Castric and Slotte2019; Gos et al., Reference Gos, Slotte and Wright2012; Koenig et al., Reference Koenig, Hagmann, Li, Bemm, Slotte, Neuffer, Wright and Weigel2019). Trans-species SNPs are very unevenly distributed along the genome and strongly enriched in immune-related loci, strongly suggesting balancing selection at NLR and other immune genes. Many of the trans-species SNPs could be reproduced in a comparison with a third, also selfing species, C. orientalis, that had apparently undergone a very recent, very strong genetic bottleneck (Koenig et al., Reference Koenig, Hagmann, Li, Bemm, Slotte, Neuffer, Wright and Weigel2019). While intriguing, a note of caution comes from a study that showed how the heterogeneous distribution of genomic regions with high nucleotide diversity (and ancestral polymorphism) in the genome of selfing species could also be attributed to the process of transition from outcrossing to selfing (Metzger et al., Reference Metzger, Abu-Awad, Sellinger and Tellier2026; Strütt et al., Reference Strütt, Sellinger, Glémin, Tellier and Laurent2023). This means that such footprints of trans-species polymorphism may also be generated by other processes, such as a transition in reproductive mode, without the action of balancing selection. The original analysis (Koenig et al., Reference Koenig, Hagmann, Li, Bemm, Slotte, Neuffer, Wright and Weigel2019) was based on short-read data, and reconstructing entire haplotypes with long reads should greatly help to further resolve the mutational history and structural variation at the loci in question (Metzger et al., Reference Metzger, Abu-Awad, Sellinger and Tellier2026) and confirm what type of selection has acted on these. Following theoretical predictions (Brown & Tellier, Reference Brown and Tellier2011), scans for selection tend to confirm the predominance of footprints of positive selection (arms race) in crops (and their pathogens) and of balancing selection (trench warfare) in natural systems (but see Gladieux et al., Reference Gladieux, van Oosterhout, Fairhead, Jouet, Ortiz, Ravel, Shrestha, Frouin, He, Zhu, Morel, Huang, Kroj and Jones2024). However, these scans are usually not accurate for gene families (with several paralogs) due to bias of SNP calls for short-read data. Improved scans for selection based on new Deep Learning or Approximate Bayesian Computation methods (Metzger et al., Reference Metzger, Abu-Awad, Sellinger and Tellier2026) and including information from genome graphs (Bao & Weigel, Reference Bao and Weigel2025) would be expected to help confirm or refute these trends. Furthermore, recent advances in bioinformatic prediction of protein–protein interactions using the Deep Learning method AlphaFold (Jumper et al., Reference Jumper, Evans, Pritzel, Green, Figurnov, Ronneberger, Tunyasuvunakool, Bates, Žídek, Potapenko, Bridgland, Meyer, Kohl, Ballard, Cowie, Romera-Paredes, Nikolov, Jain and Adler2021) are now used to predict potential interaction strength between plant and pathogen molecules (Fick et al., Reference Fick, Fick, Swart and van den Berg2025; Homma et al., Reference Homma, Lyu, Renier and der Hoorn2024; Lyu & Van der Hoorn, Reference Lyu and Van der Hoorn2026). We suggest connecting these predictions with evolutionary analyses to enhance the accuracy of the inferred selective pressures (or lack thereof). This can be realized, for example, by integrating scores of protein–protein interaction strength (from AlphaFold) for different host and pathogen alleles as statistics into selection scan methods.

7. What is the effect of linked selection at and around NLR genes?

A consequence of the heterogeneity in recombination rates (Choi et al., Reference Choi, Reinhard, Serra, Ziolkowski, Underwood, Zhao, Hardcastle, Yelina, Griffin, Jackson, Mézard, McVean, Copenhaver and Henderson2016; Jiao et al., Reference Jiao, Patel, Klasen, Liu, Pecinkova, Ferrand, Gy, Camilleri, Effgen, Koornneef, Pecinka, Loudet and Schneeberger2021) and in the strength of positive, negative and balancing selection along the genome is that the genome is replete with regions that experience linked selection due to hitch-hiking and background selection, which in turn affects local nucleotide diversity and population genetics statistics beyond the focal loci that are under selection (Charlesworth & Jensen, Reference Charlesworth and Jensen2021; Table 2). The effect of linked selection is more pronounced in species with long-range linkage disequilibrium, such as selfing species (A. thaliana, Arabis ssp. and Capsella rubella), compared to outcrossers (Solanum chilense and Capsella grandiflora). Linked selection has several often overlooked consequences for the study of NLR genes, gene clusters and duplicates.

First, linked selection leads to per-site recombination rate and nucleotide diversity becoming correlated (demonstrated in Begun and Aquadro, Reference Begun and Aquadro1993 and in A. thaliana by Lian et al., Reference Lian, Solier, Walkemeier, Durand, Huettel, Schneeberger and Mercier2022). This correlation stems from meiotic recombination reducing the effect of linked selection (decreased linkage disequilibrium) around the selected sites (Begun & Aquadro, Reference Begun and Aquadro1993; Charlesworth & Jensen, Reference Charlesworth and Jensen2021). The effect of linked selection can be, for example, uncovered by analysing measures of synonymous to non-synonymous polymorphisms ( $ \pi $ _N/ $ \pi $ _S, Table 2) while accounting for local genomic recombination rates (Campos et al., Reference Campos, Halligan, Haddrill and Charlesworth2014). When recombination is low, higher linkage disequilibrium and stronger linked selection generate so-called genetic draft around selected loci by which a genomic region exhibits footprints of neutral evolution (Table 2), albeit with reduced nucleotide diversity (Achaz & Schertzer, Reference Achaz and Schertzer2023). Furthermore, reduced recombination enhances the Hill–Robertson effect (Hill & Robertson, Reference Hill and Robertson2007) by which selection at the focal locus influences the efficacy of selection at neighbouring loci/genes. In other words, a locus under positive selection or negative selection can influence adjacent loci that are under the opposite type of selection. The polymorphism pattern in the adjacent loci would thus not produce the expected footprints of selection (Table 2). Deciphering the evolutionary history of NLR genes, which can be found in tandem arrangements or in gene clusters, may prove to be particularly difficult in species with long-range LD. The correlation in nucleotide diversity or other neutrality tests (Table 2) between tandem or duplicated NLRs should be interpreted with caution, as it may not reflect co-adaptation of genes but rather correlation due to linked selection. The same caution applies to study polymorphism within NLR domains. These may exhibit different evolutionary trajectories (and selective pressures), and the resulting polymorphism patterns depend on the intra-NLR recombination rate. When analysing NLR gene evolution over phylogenetic timescales, recombination and linked selection also play a role, affecting the accuracy of phylogenetic reconstruction of the history of tandem NLR genes or members of NLR gene clusters. We note that the sensitivity of phylogenetic methods to infer selection via the d_N/d_S (Table 2) ratio does depend on recombination rate (Kryazhimskiy & Plotkin, Reference Kryazhimskiy and Plotkin2008). To our knowledge, it has not yet been fully assessed how sensitive to recombination (and potential biases) the tools are that are used to infer episodic measures of selection along a phylogeny (as used in Sutherland et al., Reference Sutherland, Prigozhin, Grey Monroe and Krasileva2024). Accounting for recombination rate (and linked selection) is thus a prerequisite for studying the possible co-adaptation of NLR genes in tandem or in clusters. For example, the strength, prevalence and consequences of linked selection and co-adaptation between NLR genes are likely to differ between A. thaliana, a selfing species, and the outcrossing wild tomato S. chilense (see discussion in Silva-Arias et al., Reference Silva-Arias, Gagnon, Hembrom, Fastner, Khan, Stam and Tellier2025).

8. How many NLR genes are pseudogenized?

A key feature of the birth-and-death process is the intermediate step of pseudogenization of NLR genes. The assessment of pseudogenization faces, however, several issues. First, there are many examples of both recently truncated NLR alleles or of bona fide NLR genes that lack many, or even most of the canonical NLR domains, such as TIR-only genes, that are clearly functional (Barragan & Weigel, Reference Barragan and Weigel2021). In addition, while the progressive accumulation of mutations might lead to eventual complete removal of a gene from the genome, pseudogenized NLR genes could act as a reservoir of diversity (for subsequent neofunctionalization, e.g., see Garewal et al., Reference Garewal, Pathania, Bhatia and Singh2022). It has been investigated whether pseudogenized NLR genes are likely to evolve neutrally by comparing the ratio of non-synonymous to synonymous changes in genes and pseudogenes in A. thaliana and cultivated rice (Zou et al., Reference Zou, Lehti-Shiu, Thibaud-Nissen, Tanmay Prakash, Buell and Shiu2009). While not all NLR genes could be resolved in this study due to the lack of long-read sequencing data, the authors observed 1) that many genes annotated as pseudogenes are actually subject to purifying selection, 2) that some apparent pseudogenes are expressed (albeit at a low level), and 3) that there is a positive correlation between the size of gene families with different functional domains and the number of pseudogenes. We suggest that by combining the methods described above and in Table 2, including AlphaFold predictions of host NLR–pathogen effector interactions (Homma et al., Reference Homma, Lyu, Renier and der Hoorn2024; Lyu & Van der Hoorn, Reference Lyu and Van der Hoorn2026), we can dissect the various evolutionary forces driving the evolution of NLR genes that may or may not be true pseudogenes and infer the rate of neofunctionalization.

9. What is the relative importance of neutral versus selective forces?

The efficacy of any type of selection depends on the effective population size, that is, the amount of genetic drift. This is encapsulated in the drift barrier concept (Lynch, Reference Lynch2007), which states that the evolution of mutation rates is constrained by the opposing effects of natural selection and genetic drift. High mutation rates are deleterious to the genome, and low mutation rates decrease a population’s adaptability (Hanlon et al., Reference Hanlon, Cagan and den Akker2025; Lynch, Reference Lynch2007). Therefore, it can be expected that an optimal mutation rate is reached for the genome of a given species, and possibly for different types of genes, if genetic drift is not too strong. In A. thaliana, the mutation rate may be only marginally higher for NLR genes than for other genes in the genome, if at all (see above), but this rate remains unknown for other plant species. We speculate that the small local population size and high selfing rate in A. thaliana generate relatively strong genetic drift (and LD), possibly limiting the ability of fine-tuned mutation rates at NLR genes to evolve. We also note that the efficacy of selection for positive mutations and against deleterious mutations (purifying selection) is likely to be weaker in A. thaliana compared to outcrossing plant species. This could mean that mildly deleterious (costly) NLR gene copies may remain in the genome for extended periods, appearing thus as segregating PAVs/CNVs in A. thaliana, whereas they would be more quickly lost from the entire population in outcrossers with stronger purifying selection. In addition, advantageous NLR genes that confer pathogen resistance may not always become rapidly fixed (and remain polymorphic as CNVs) for a long period if local selection is weak compared to drift. In addition, certain combinations of alleles at unlinked NLRs can be deleterious (Bomblies & Weigel, Reference Bomblies and Weigel2007), which could increase the strength of purifying selection against NLRs that have otherwise no positive function in the organism. In species that mostly self, this effect is expected to be negligible, but it might be more relevant in obligate outcrossers.

A second process to consider is meiotic recombination. As for mutation, there are advantages and disadvantages to recombination (Barton & Charlesworth, Reference Barton and Charlesworth1998), and the optimum rate of homologous recombination at NLR genes must be determined by the strength of genetic drift (as it constrains the efficacy of selection). However, contrary to the Red Queen hypothesis, meiotic recombination in A. thaliana does not appear to be higher at NLR genes than in the rest of the genome (Choi et al., Reference Choi, Reinhard, Serra, Ziolkowski, Underwood, Zhao, Hardcastle, Yelina, Griffin, Jackson, Mézard, McVean, Copenhaver and Henderson2016), and thus other mechanisms such as gene duplication, deletion and unequal recombination may be necessary to generate the observed NLR gene clusters and large NLR gene families. NLR gene duplicates need to diverge by mutation to undergo subfunctionalization or neosubfunctionalization (see above), a process that is slowed down by gene conversion and by the low effective population size and thus mutational input in A. thaliana. These various mechanisms (Figure 1) collectively provide a general explanation for the large number of CNVs found at NLR genes in A. thaliana (Lee & Chae, Reference Lee and Chae2020; Teasdale et al., Reference Teasdale, Murray, Collenberg, Contreras-Garrido, Schlegel, van Ess, Jüttner, Lanz, Deusch, Fitz, Mencia, van Velthoven, Drost, Weigel and Shirsekar2025; Van de Weyer et al., Reference Van de Weyer, Monteiro, Furzer, Nishimura, Cevik, Witek, Jones, Dangl, Weigel and Bemm2019) versus outcrossing species (such as S. chilense; Silva-Arias et al., Reference Silva-Arias, Gagnon, Hembrom, Fastner, Khan, Stam and Tellier2025).

10. How do NLR genes evolve in the context of gene clusters and networks?

Ultimately, we aim to integrate the previous insights from population genomics with recent molecular insights on the function of NLR genes (Figure 1). As NLR genes function as part of gene networks, one can predict heterogeneity in selective constraints across genes. For example, more conserved NLR genes (such as the genes coding for the RNL class; Saile & El Kasmi, Reference Saile and El Kasmi2023; or the NRC helper clade; Sakai et al., Reference Sakai, Contreras, Martinez-Anaya, Lüdke, Kamoun, Chih-Hang and Adachi2024; Stam et al., Reference Stam, Nosenko, Hörger, Stephan, Seidel, Kuhn, Haberer and Tellier2019; Wu et al., Reference Wu, Abd-El-Haliem, Bozkurt, Belhaj, Terauchi, Vossen and Kamoun2017) are likely to be more pleiotropic and have more connections to downstream signalling (and be under stronger purifying selection), whereas sensor NLR genes involved in effector recognition may be more peripheral in such networks (and under weaker purifying selection). NLR gene network structures are likely organized around host genes coding for key guardee proteins that are targeted by pathogen effectors for their function in plant physiology, development or defence response. The type of genes involved in defence response, such as those coding for guardee (an effector target, often having an important role in the host, but not a resistance protein), guard (resistance protein monitoring this guardee) and decoy (resistance protein, which chiefly recognizes pathogen effector and distracts these from the guardee, hence the name decoy) could thus be understood in the context of gene networks, multilocus epistatic interactions and trade-offs for effector binding and downstream signalling (Weiner et al., Reference Weiner, Märkle, Laderman, Demirjian and Bergelson2025). This epistasis between NLR genes defines the strength and type of selection occurring at these loci (positive, negative and balancing) as well as the inherent evolutionary constraints. We speculate that neofunctionalization of gene copies of more connected NLR genes in the networks may explain the evolution of new functions, like guard and decoy coding genes, as gene duplication can free one copy for further specialization and refinement of the effector detection function (Adachi et al., Reference Adachi, Derevnina and Kamoun2019; Dodds & Rathjen, Reference Dodds and Rathjen2010; Van der Hoorn et al., Reference Van der Hoorn, De Wit and Joosten2002). This process is arguably most visible in the evolution of integrated domains (IDs), where duplicated guardee domains are integrated directly into the NLR architecture (Li et al., Reference Li, Wang, Jia, Ishikawa, Kosami, Ueba, Tsujimoto, Yamanaka, Yabumoto, Miki, Sasaki, Fukao, Fujiwara, Kaneko-Kawano, Tan, Kojima, Wing, Sebastian, Nishimura and Kawano2024; Yang et al., Reference Yang, Kim, Cevik, Jacob, Wan, Furzer and Dangl2022). Duplicated gene copies of peripheral genes of the NLR networks may also (neo-)subfunctionalize if selection pressure is strong enough. Duplicated pairs of phylogenetically distant paralog NLRs, which subfunctionalize, could be found as genes coding for sensor and helper. In other words, the topology of these functional networks would determine the physical arrangement of NLRs in the genome, for example, co-functioning pairs (sensor/helper) would be maintained in tight physical clusters (head-to-head orientation), which seems to be supported by the detected enrichment of head-to-head pairs of rice NLRs under balancing selection (Gladieux et al., Reference Gladieux, van Oosterhout, Fairhead, Jouet, Ortiz, Ravel, Shrestha, Frouin, He, Zhu, Morel, Huang, Kroj and Jones2024). Due to the relatively high linkage disequilibrium and linked selection in A. thaliana, the sensor and helper genes are expected to exhibit similar population genetics statistics (Table 2; Van de Weyer et al., Reference Van de Weyer, Monteiro, Furzer, Nishimura, Cevik, Witek, Jones, Dangl, Weigel and Bemm2019). In line with this idea, there is evidence for the evolution of NLR gene clusters and duplicates starting from a duplicated pair of sensor and helper NRCs in the Asteroid family (Sakai et al., Reference Sakai, Contreras, Martinez-Anaya, Lüdke, Kamoun, Chih-Hang and Adachi2024).

Moving beyond single genes and isolated defence pathway responses, it becomes evident that NLRs are parts of complex and diverse gene regulatory networks (Wu et al., Reference Wu, Abd-El-Haliem, Bozkurt, Belhaj, Terauchi, Vossen and Kamoun2017; Sakai et al., Reference Sakai, Contreras, Martinez-Anaya, Lüdke, Kamoun, Chih-Hang and Adachi2024) that integrate pathogen recognition and downstream defence response in a more quantitative way than classically acknowledged (Weiner et al., Reference Weiner, Märkle, Laderman, Demirjian and Bergelson2025). The large body of theory on the evolution of gene regulatory networks (e.g., Pouzet & Le Rouzic, Reference Pouzet and Le Rouzic2025; Salathé & Soyer, Reference Salathé and Soyer2008; Schrom et al., Reference Schrom, Prada and Graham2018; Wagner, Reference Wagner1996) allows us to spell out the selective and mutational constraints determining robustness, inducibility and evolvability of the plant immune system gene network. The immune gene network is robust and inducible by mounting an efficient defence response to various pathogens (despite parasite interference), and concomitantly by avoiding triggering defences in the absence of pathogens (avoiding auto-immunity). Furthermore, the immune gene network requires some degree of evolvability to incorporate new defence traits, recognize new pathogens and fend off parasite interference, enabling plants to keep up in the coevolutionary Red Queen race with pathogens (Salathé & Soyer, Reference Salathé and Soyer2008; Schrom et al., Reference Schrom, Prada and Graham2018). Gene regulatory networks evolve through regulatory or coding point mutations, gene duplications and deletions, generating possible redundancy (Salathé & Soyer, Reference Salathé and Soyer2008). Large-scale effects on the network robustness and inducibility can occur by gene duplications and deletions (Salathé & Soyer, Reference Salathé and Soyer2008; Schrom et al., Reference Schrom, Prada and Graham2018) and mutations in transcription factors (affecting a substantial amount of target genes), while smaller localized effects are expected to occur following cis-acting point mutations (Pouzet & Le Rouzic, Reference Pouzet and Le Rouzic2025). The network structure appears to have a greater influence on the robustness and output response than the type of mutation or gene duplication or loss (Pouzet & Le Rouzic, Reference Pouzet and Le Rouzic2025). Further, larger adaptive steps (with larger positive selection coefficients) may occur first at more connected hub genes in the network, followed by fine-tuning of the response via mutations at peripheral, less connected genes (Pouzet & Le Rouzic, Reference Pouzet and Le Rouzic2025; Wei et al., Reference Wei, Sharifova, Zhao, Sinha, Nakayama, Tellier and Silva-Arias2024). The latter observations contradict in part the alternative expectations that more connected hub genes are very conserved, and adaptation occurs only by changes at the peripheral, less connected genes (e.g., the omnigenic model; Boyle et al., Reference Boyle, Li and Pritchard2017). We suggest two hypotheses that can be tested using evolutionary analyses of NLR gene families and of the plant immune gene network. 1) The change of hub (well-connected and pleiotropic) NLR genes (such as NRC helper genes) in networks would allow large steps of adaptation, for example, recognition of new pathogens with new host target gene (guardee), but at the price of possible loss of signalling competence for upstream sensors and/or of control of the resistance pathway (yielding autoimmunity). 2) The change of peripheral (less connected) NLR genes would result in recognizing additional pathogen effectors without destabilizing the networks, but these changes would be limited to altering recognition specificity (e.g., via LRR variation or ID acquisition). NLR gene duplication, mutation, recombination and CNV would then be under the selective constraint defined by the structure and robustness of the NLR defence networks and the quantitative nature of NLR recognition and activation (Figure 1; Weiner et al., Reference Weiner, Märkle, Laderman, Demirjian and Bergelson2025). We also suggest extending AlphaFold predictions to interactions between different host proteins (NLRs, guardee, guard, decoy, sensor, helper and downstream defence signalling) to generate a priori expectations on the structure of the gene network and the selective constraints at the level of NLR (or at the level of their protein domains). These predictions can thereafter be tested using polymorphism from pangenome data (Figure 1).