To obtain comprehensive NLR from spinach, we identified NLR across 19 representative Spinacia assemblies using NLR-Annotator v. 2.1b (Steuernagel et al., 2020) with default parameters. The NLRs in spinach could be further divided into six subfamilies: CC-NBARC-LRR, CC-NBARC, NBARC-LRR, NBARC, TIR-NBARC, TIR-NBARC-LRR. To assess the distribution of NLRs, we visualized them using the telomere-to-telomere genome (Sp_YY_v2) as an example.

Based on the result from NLR-Annotator, we also obtained each NLR and its flanking sequences (2 kb) for each genome, generating pan-NLRome sequences.

To further determine NLRs and their corresponding protein-coding genes, we performed automatic structural gene annotations for 19 Spinacia assemblies using Helixer v0.3.6 (Holst et al., 2023), as previous studies have demonstrated its exceptional accuracy in annotating NLR genes (Belinchon-Moreno et al., 2025). We employed an in-house script to identify NLR and its overlapped protein-coding genes, which were then subjected to further analysis.

We constructed a pan-NLRome using the NLR protein-coding genes from 19 Spinacia assemblies. First, we clustered these genes using OrthoFinder (Emms and Kelly, 2015) with the default parameters, and then divided these gene families into three categories, core, dispensable, and private based on their frequency. We defined the NLR families as core, dispensable, and private if they were present in all 19 accessions, 3–18 accessions, and 1–2 accessions, respectively.

For comparison of the K_A/K_S value between the core and dispensable NLR families, we first selected the longest length NLR protein-coding gene from each accession within both the core and dispensable NLR families as the representative gene. Subsequently, TBtools (Chen et al., 2020) was employed to calculate the K_A/K_S value for homologous gene pairs between each pair of genomes within both the core and dispensable NLR families.

We utilized 116 spinach accessions (She et al., 2024a), comprising 30 resistant and 86 susceptible to Peronospora effuse race 9 (Pfs9), to identify causal variants associated with resistance to Pfs9 based on our optimized k-mer-based GWAS. First, we filtered low-quality short reads from 116 accessions using fastp v0.23.4 (Chen et al., 2018) with the parameter “-q 20”. Then, the clean reads were aligned to pan-NLRome sequences using BWA 0.7.17-r1188 (Li, 2013) with default parameters. After removing duplicated reads, we obtained the paired-end (PE) reads that were located at pan-NLRome sequences (termed NLR-reads) using SAMtools (Li et al., 2009) with the parameters “view –F 12 –q 30”.

We created the k-mers table as described in Jaegle et al. (2025). Specifically, we firstly extracted k-mers (31 bp) from NLR-reads for each accession using KMC (Kokot et al., 2017). Then, we combined and filtered lists of k-mers using the script list_kmers_found_in_multiple_samples with parameters “--mac 5 –p 0.2”. Last, we obtained the k-mers table, containing the presence/absence pattern of each k-mers in the 116 spinach accessions. Combining the phenotype and k-mer table, we run k-mer-based GWAS with a permutation-based threshold for 5% family-wise error rate.

To retrieve the coordinates of the significant k-mers, we first extracted the PE reads containing the significant k-mers using fetch_reads_with_kmers (https://github.com/voichek/fetch_reads_with_kmers). Subsequently, we quantified the number of PE reads containing significant k-mers on pan-NLRome sequences to identify the target NLR or determined the coordinates of the significant k-mers by aligning the corresponding PE reads against the target genome using BWA 0.7.17-r1188 (Li, 2013) with default parameters.

We aligned the clean reads against the Sp_YY_v2 genome using HISAT2 v2.1.0 (Kim et al., 2015), followed by calculating the read count using featureCounts v2.0.1 with the parameters ‘-T 10 -p -t exon -g ID’ (Yang et al., 2014). Gene expression was normalized using the transcripts per million (TPM) method with an in-house script. The expression patterns are shown using R v4.1.1.

The 116 spinach accessions mentioned above were used for SNP-based GWAS analysis based on the Sp_YY_v2 assembly. The clean reads were aligned to Sp_YY_v2 assembly using BWA 0.7.17-r1188 (Li, 2013) with default parameters. The SNPs were identified using GATK v4.3 (McKenna et al., 2010), and then filtered using the GATK with “QD < 2.0 | | FS > 60.0 | | MQ < 40.0| | SOR > 3.0 | | MQRankSum < −12.5 | | ReadPosRankSum < −8.0”, and indels with “QD < 2.0 | | FS > 200.0 | | SOR > 10.0 | | MQRankSum < −12.5 | | ReadPosRankSum < −8.0; (2)”, and VCFtools v0.1.16 (Danecek et al., 2011) with the parameters “--max-missing 0.85 --mac 4 --minQ 30 --maf 0.05 --min-alleles 2 --max-alleles 2”. A total of 1,087,947 high-quality SNPs were obtained. The GWAS was performed using EMMAX software (Kang et al., 2010).

To determine the dynamics in different Spinacia genomes, we predicted NLR genes in 19 Spinacia genomes using NLR-Annotator (Steuernagel et al., 2020). Fourteen spinach genomes were from our previous studies (She et al., 2023, 2024a, 2024b); the remaining five assemblies (Monoe-Viroflay, SpoV3, Cornell-No.9, SOL_r2, and Sp75) were used in others (Xu et al., 2017; Cai et al., 2021; Hirakawa et al., 2021; Hulse-Kemp et al., 2021; Ma et al., 2022). These accessions included 14 cultivated spinach species, S. oleracea, two of its closest wild relatives, S. turkestanica (0.8 million years ago (Mya)), and three more distant relatives, S. tetrandra (~6.3 Mya) (Figure 1A; Supplementary Table S1).

In total, 2,549 NLR genes were annotated in the 19 spinach genomes, with a range of 128 to 143 genes per accession (Supplementary Table S2). The 12 cultivated spinach plants and their five wild relatives shared a similar number of NLR genes (Wilcoxon test, p = 0.89) (Supplementary Figure S1). Furthermore, these NLR genes were categorized into six subfamilies, and the most frequent NLR class was CC-NBARC-LRR (60.58%), followed by NBARC-LRR (23.90%) and NBARC (9.96%) (Figure 1A, Supplementary Table S2). TIR-NBARC and TIR-NBARC-LRR are rarely detected in Spinacia species, particularly in S. tetrandra species (Supplementary Table S2), suggesting that the two subfamilies probably evolved after splitting from S. tetrandra and S. turkestanica/S. oleracea. NLRs exhibited an uneven distribution along chromosomes, invariably clustering together and predominantly located near subtelomere regions (Figure 1B). Applying the definition of NLR clusters as genes within 200 kb of each other in the genome (Holub, 2001), 47.73%–57.86% of NLRs in each accession were located in such clusters (Supplementary Table S3, Supplementary Figure S2), consistent with a previous report on Arabidopsis thaliana (Van de Weyer et al., 2019).

To understand the variation in NLR content, we inferred 2846 NLR protein-coding genes that overlapped with NLR loci in the corresponding assembly. 2813 (99%) of these genes were annotated and associated with disease resistance using the Non-redundant protein Sequence (NR), Swiss-Prot, Pfam, TrEMBL, and A. thaliana databases (Supplementary Table S4). Then, we identified 186 pan-NLRome by clustering NLR protein-coding genes (see Methods) in the 19 Spinacia accessions (Supplementary Table S5). The number of NLR families exhibits a positive correlation with genome number, with their abundance stabilizing at a genome number of 12 (Figure 2A), indicating that the 19 Spinacia accessions are diverged and that the pan-NLRome closely captures the NLRs of spinach. Furthermore, only one and none additional NLR family was found when the 12^th and 17^th accession was added, respectively (Supplementary Figure S3).

Based on the frequency of NLR gene families in the 19 Spinacia genomes, we further divided the pan-NLRome into core, dispensable, and private families. In total, we obtained 72 (38.7%) core NLR families present in 19 accessions, 95 (51.1%) dispensable NLR families present in 3–18 accessions, and 19 (10.2%) NLR families present in 1–2 accessions (Figure 2B). The proportion of pan-NLRome gene families across 19 genomes revealed that dispensable NLRs are more likely to be present in S. oleracea/S. turkestanica compared to S. tetrandra, whereas private NLRs showed the opposite trend (Figure 2C), suggesting high diversity between S. oleracea/S. turkestanica and S. tetrandra. Moreover, the 19 spinach accessions shared an average of 58.99% of core NLR genes, exceeding the proportion found in dispensable NLR genes (39.93%) (Figure 2D, Supplementary Table S6). These dispensable NLR genes were significantly prevalent in S. oleracea/S. turkestanica (Wilcoxon test, p < 0.01); conversely, the private NLR genes (1.48%) were predominantly present in the S. tetrandra (Wilcoxon test, p < 0.05) (Supplementary Figure S4), further confirming that S. tetrandra represents a more distant spinach wild relative.

Among the core, dispensable, and private NLR genes, the most frequent NLR class was CC-NBARC-LRR (52.38%–70.24%), followed by the NBARC-LRR (22.29%–24.16%) and NBARC (3.25–9.52%) (Figure 2E). The proportion of CC-NBARC in private NLR genes (14.29%) was higher than that in core (2.82%) and dispensable (2.29%) NLR genes (Figure 2E; Supplementary Table S7). Moreover, the core NLR genes exhibited longer CDS length (Wilcoxon test, p < 2e-5) and lower K_A/K_S values (Wilcoxon test, p < 0.0097) compared to dispensable NLR genes, suggesting that core NLR genes possess more conserved functions, consistent with previous findings in soybean (Liu et al., 2020) and broomcorn millet (Chen et al., 2023).

K-mer-based GWAS is an efficient approach to identify causal variants associated with phenotypes, particularly disease resistance in plants (Voichek and Weigel, 2020; Jaegle et al., 2025). To enhance the efficiency of disease-resistance gene discovery, we optimized this approach by integrating the pan-NLRome, which primarily consists of three steps: identification of NLR and flanking sequences (Figure 3A), identification of NLR reads and construction of k-mer table (Figure 3B), and conduction of k-mer-based GWAS (Figure 3C). First, we obtained the nucleotide-binding site leucine-rich repeat (NLR) and its flanking sequences (2 kb) for each of the 19 Spinacia genomes, generating a pan-NLRome sequences. Secondly, we identified NLR reads by aligning short reads of spinach accessions to pan-NLRome sequences, followed by extracting k-mers (31 bp) and constructing a k-mers table. Last, we performed k-mer-based GWAS based on phenotype and k-mer table as described by Voichek and Weigel, 2020.

To demonstrate the usefulness of our approach, we utilized 116 spinach accessions (Supplementary Table S8) (She et al., 2024a), comprising 30 resistant and 86 susceptible to Peronospora effuse race 9 (Pfs9), to identify downy mildew resistance gene RPF1 loci that resist Pfs9 (Correll et al., 2011; She et al., 2018). Among these accessions, a total of 111.02 Gb of short reads were obtained from the NLR and flanking sequences, generating a k-mer table comprising 10,001 k-mers. In our k-mer-based GWAS, we detected 50 significant k-mers associated with RPF1, with the Te17S24XX_Chr1_nlr42 gene enriching the highest number of k-mers and exhibiting high significance. (Figure 4A). Notably, our previous study identified Te17S24XX_Chr1_nlr42 (formerly designated as Spo12903) as a candidate gene for RPF1 using BC₁ and F₂ populations (She et al., 2018).

Otherwise, to assess the effectiveness of k-mer-based GWAS in a single genome, we mapped the PE reads containing the significant k-mers against the Sp_YY_v2 assembly. The vast majority of these reads were located within the 781 kb–1101 kb interval on chromosome 1 (Figures 4B, C, Supplementary Table S9), consistent with previous observations (She et al., 2018; Bhattarai et al., 2020; Cai et al., 2021). Within the candidate region, we identified four NLR genes, three of which carry the canonical CC-NBARC-LRR architecture, and one gene, YYT2T_Chr1_005275.1, a homolog of Te17S24XX_Chr1_nlr42 mentioned above, lacks the N-terminal coiled-coil domain (NBARC-LRR only) (Figure 4D). The expression analysis indicates that all four genes exhibited elevated expression levels at 7 days post-inoculation, with only the YYT2T_Chr1_005275.1 gene maintaining sustained expression at both 14 and 24 days (Figure 4E). Therefore, we confirmed that Spo12903 is probably RPF1.

Furthermore, we performed the SNP-based GWAS using the Sp_YY_v2 reference genome, which successfully detected significant SNPs encompassing the YYT2T_Chr1_005275.1 locus (Supplementary Table S10). However, the SNP-GWAS also identified numerous other significant signals across the genome, making it difficult to unequivocally prioritize this specific gene (Supplementary Figure S5). Taken together, compared to a single genome, pan-NLRome enables faster and more accurate anchoring of candidate genes.