NLR Genes Related Transcript Sets in Potato Cultivars Bearing Genetic Material of Wild Mexican Solanum Species

: The long history of potato breeding includes the numerous introgressions of resistance genes from many wild species of South and Central America as well as from cultivated species into the breeding genepool. Most R genes belong to the NLR family with nucleotide-binding site–leucine-rich repeat. The aim of this research concerns an evaluation of NLR genes expression in transcriptomes of three potato cultivars (Evraziya, Siverskij, Sudarynya), which combine genetic material from wild and cultivated potato species, and each bears intragenic markers of RB/Rpi-blb1/Rpi-sto1 genes conferring broad-range resistance to late blight. The transcriptomes of the cultivars were compared before and 24 h after the Phytophthora infestans inoculation. The induction of RB/Rpi-blb1/Rpi-sto1 transcript after 24 h of inoculation was detected in the resistant cultivars Siverskij and Sudarynya but not in susceptible cv. Evraziya. This demonstrates the importance of transcriptomic assay for understanding the results of marker-assisted selection and phenotyping. Interestingly, assembling the transcriptomes de novo and analysis with NLR-parser tool revealed signiﬁcant fractions of novel NLR genes with no homology to the reference genome from 103 (cv. Siverskij) to 160 ( S. stoloniferum, 30514/15). Comparison of novel NLRs demonstrated a relatively small intersection between the genotypes that coincided with their complex pedigrees with several interspeciﬁc hybridization events. These novel NLRs may facilitate the discovery of new efﬁcient R genes.


Introduction
Marker-assisted selection (MAS) is efficiently exploited in selecting desirable genotypes with resistance genes (R genes) and stacking several R genes to facilitate the breeding of new cultivars with durable resistance [1][2][3]. MAS can increase the efficiency of conventional breeding, and it is applicable for traits conferred by both major qualitative genes and quantitative trait loci (QTLs). Occasionally, MAS can be less efficient (or predictable) because of the insufficient linkage between a marker and a gene (QTL), the complexity of the genetic control of a resistance trait, or the instability of QTL effects. In addition, in many cases the genetic factors providing resistance against a specific pathogen remains unidentified.
Systemic comparison of 1301 A. thaliana natural accession showed the presence of variants of defense-related genes in nearly all populations investigated [37].
The long history of potato breeding includes numerous introgression of resistance genes from wild and cultivated species into the breeding genepool. The number of introgression events may vary and reach up to one hundred [38]. An investigation of NBS domains in the genomes of 91 potato cultivars revealed 587 segments mapped on the reference genome of the diploid (doubled monoploid) S. tuberosum Group Phureja clone (clonal cultivar) DM 1-3 516 R44 [39]. As expected, a number of R genes detected in modern cultivars was absent in the reference genome. However, an evaluation of potato genotypes with genome sequencing approaches will reveal both functional and non-functional NLR genes because of their rapid evolution and high recombination rate in the gene clusters [38].
The aim of this research concerns an evaluation of NLR genes expression in the transcriptomes of three potato cultivars (Evraziya, Siverskij, Sudarynya), which combine genetic material from wild and cultivated potato species: Solanum stoloniferum Schltdl., S. demissum Lindl., S. acaule Bitter, and S. tuberosum L. ssp. andigena Hawkes [40] and may bear different subsets of resistance genes. The selected cultivars possessed several allele-specific markers of well-known broad-spectrum late blight resistance genes (RB/Rpi-blb1/Rpi-sto1). According to our knowledge, the RB/Rpi-blb1/Rpi-sto1 homologue in potato cultivars developed by conventional breeding through interspecific crosses was first described in our work [40,41]. The RB/Rpi-blb1 gene was first identified in wild Mexican potato species S. bulbocastanum [42][43][44]. The functional homolog of the RB/Rpi-blb1 gene was detected in wild Mexican potato species S. stoloniferum-Rpi-sto1 [45], as well as in accessions of S. papita, S. cardiophyllum, S. brachistotrichum [46,47]. With the exception of transgenic varieties of potato obtained by cisgenesis, the RB/Rpi-blb1/Rpi-sto1 is extremely rare in the genomes of potato cultivars [48]. To reveal NLR genes' variability, we compared the de novo-assembled transcriptomes of cultivars Evraziya, Siverskij, Sudarynya before and 24 h after inoculation with P. infestans.

Plant Materials
Three potato cultivars chosen for this study, Evraziya, Siverskij, Sudarynya, were bred in Leningrad Research Institute for Agriculture "Belogorka" (LenNIISKh 'Belogorka') located in the Northwestern zone of the Russian Federation where prolonged periods of humid and cool summer weather are favorable for late blight. Cultivar Sudarynya showed late blight resistance under field conditions of 2010-2017, cv. Evraziya demonstrated moderate level of resistance or susceptibility depending on the season, whereas cv. Siverskij was released quite recently [40]. These cvs. originated from the hybrid clone 8889/3, which included genetic material from S. stoloniferum, S. demissum, S. tuberosum ssp. andigena (the detailed genetic pedigrees of these cultivars were published earlier and are available in [40] (Supplement 2)). These three cultivars were selected for the present study based on the results of marker-assisted selection-five allele-specific markers of broad-spectrum late blight resistance genes RB/Rpi-blb1/Rpi-sto1 were detected in all three cultivars (for details, see [40,41]). The sequences of the Rpi-sto1 and BLB1 F/R-amplicons were identical in analyzed genotypes and showed 99 and 100% similarity to the corresponding fragments of the Rpi-blb1 and Rpi-sto1 genes from the GenBank database, respectively [41].
Samples of cvs. Evraziya, Siverskij, and Sudarynya were provided by the cultivar author Dr. Z. Evdokimova from LenNIISKh 'Belogorka'; cv. Bintje, and clone 30514/15 of wild species Solanum stoloniferum Schltdl. were obtained from VIR collection. Cultivar Bintje and S. stoloniferum, clone 30514/15, were used as susceptible and resistant controls, respectively.  [49]. Whole plants were inoculated by suspension with a concentration of approximately 50,000 sporangia/mL. Before inoculation, the suspension was incubated at 12 • C for 3 h to release zoospores. Plants were inoculated on both the underside and the topside of the leaves using a hand sprayer. Inoculated plants were covered with plastic boxes for two days to create a humid chamber (>95% relative humidity) at 23 • C during the day and 15 • C at night with a 14 h photoperiod.
Each cultivar was represented by 22-24 plants in double replication-six plants (3 at 0 hpi and 3 at 24 hpi) per cultivar/per replication were used for RNA extraction, and the remaining plants were cultivated further for the infection assay.
Plant resistance was evaluated on the 7th, 10th, and 13th days after inoculation. Each plant was evaluated on a 9-point scale, in which a score of 9 corresponded to resistance (absence of disease symptoms) and a score of 1 corresponded to susceptibility (90-100% area of leaf was covered with blight-infected lesions) [50]. Plants with a score between 9 and 6 were considered resistant A score of 5 corresponds to moderate resistance; scores of 4 or below indicate susceptibility. In all experiments, control uninfected plants of all potato genotypes had no symptoms of disease.

RNA-seq
Three biological replicates were used for each cultivar analyzed, and the leaf samples were taken before and 24 h after inoculation with P. infestans. Total RNA was extracted with an RNeasy Plant Mini Kit (Qiagen, Hilden, Germany).
The quality of the RNA samples was evaluated using a Bioanalyzer 2100 (Agilent, Santa Clara, CA, USA). All samples had RIN 7.8 or higher. RNAseq library preparations were carried out with 1.5 µg of total RNA fraction using TruSeq ® Stranded mRNA Library Prep Kit (Illumina, San Diego, CA, USA) according to the manufacturer's instructions for barcoded libraries with small modifications (4 min RNA fragmentation time and 12 PCR cycles were used). Final libraries quantification was performed with a Bioanalyzer 2100 and a DNA High Sensitivity Kit (Agilent). After normalization, barcoded libraries were pooled and sequenced on a NextSeq 550 sequencer 2 × 150 bp using NextSeq ® 500/550 High Output v2.5 Kit 300 cycles (Illumina).
Raw sequencing data were deposited in NCBI Sequence Read Archive (PRJNA765923).

Bioinformatic Analysis
FeatureCount [55] was used to estimate the libraries quality. These values were further used to perform multidimensional scaling analysis using plotMDS function in EdgeR v. 3.26.8 [56].

Transcriptome Assembly De Novo
We implemented the combined transcriptome assembly strategy suitable for accurate transcriptome reconstruction [57,58]. Transcripts were assembled independently by de novo-and genome-guided Trinity v. 2.11.0 modes [59] and with the aid of Trans-ABySS v. 2.0.1 [60] for each genotype from their specific paired read libraries. Read alignments onto potato genome obtained with the aid of Dart were further used for Trinity genome-guided transcriptome assembly.
We used Trans-ABySS assembler with four different k-mer length values (67, 87, 107, and 127) for each of the paired-end read libraries independently yielding 24 de novo transcriptome assemblies per genotype. These transcriptome sequences were merged for each genotype with the transabyss-merge tool with parameter k-range set varying from 67 to 127.
Finally, Trinity de novo-and genome-guided assemblies and Trans-ABySS assembly for each of the five potato genotypes were concatenated into a single file. To remove redundancy for genotype-specific transcriptome we used tool uclust from the usearch package v. 10.0.240 [61] with identity threshold equal to 100% to remove the duplicated transcrips.
Nucleotide sequences from non-redundant transcriptome assemblies were aligned to the S. tuberosum reference genome SolTub_3.0 by Gmap v.20 March 2018 [63] at 90% identity and 90% transcript coverage thresholds (only the best alignment for each transcript was taken). The transcript matches a gene if it overlaps with at least half of gene exons.
To find orthologous amino acid sequences for protein-coding transcripts, S. tuberosum reference genome [39] and S. lycopersicum reference genome [64] were taken from Ensembl Plants database v. 49 [53], and OrthoFinder v. 2.5.1 [65] was used as outgroup. The phylogenetic tree of the genotypes investigated was reconstructed with concatenated multiple sequence alignments of amino acid sequences from single-copy orthologous groups. The alignment was made with Mafft 7.480 (L-INS-i, max. 2 iterations) [66], the tree was reconstructed using IQ-TREE v. 1.6.12 (parameters -st AA m TESTNEW -bb 1000 -alrt 1000; JTT+F+I, substitution model was chosen automatically) [67].
Quantification of transcripts abundances was used with the aid of kallisto v. 0.46.1 [68]. Transcripts with low abundancies determined by filterByExpr function of EdgeR package v. 3.26.8 [56] were excluded from analysis. The expression threshold value for each mRNA was determined using its distribution in all transcripts and all libraries [69]. Differentially expressed genes (DEGs) were identified by EdgeR using generalized linear model (glmQLFit function) and F-test of quasi-likelihood (glmQLFTest function) at the 'FDR < 0.01' threshold.
To identify NBS-LRR domains in the amino acid sequences, the NLR-parser tool was used [70].

Compilation of the Marker Genes Associated with Wide-Range Resistance to P. infestans
Several genes conferring late blight resistance (broad-spectrum and race-specific) were selected as known markers for transcriptome-Rpi-blb1/Rpi-sto1, R3a, R3b, R8 (Table S1). A search for homologs of transcripts sequences with no similarity to reference genome was made with NCBI amino acid sequences database and ublast tool from the usearch package v. 10.0.240 [61]. Best-hit homologs were identified at 95% (or higher) sequence identity and coverage at the protein level.

Phylogenetic Analysis of the Novel NBS-LRR Sequences
A number of transcripts encoding amino acid sequences with NBS-LRR domain signatures did not show significant similarity to the potato reference genome. To characterize these putative NBS-LRR proteins, we performed phylogenetic tree reconstruction for novel and 438 known NBS-LRR potato amino acid sequences taken from Jupe et al. [13]. Five protein sequences from Table S1 were added to the dataset as well. Hmmscan from HMMER package v.3.1b2 [71] was used to identify Pfam domain signatures [72] in these sequences and extract fragments of the NB-ARC domains (Pfam id PF00931.23). Mean length of the NB-ARC domains is 209 residues; therefore, sequences with NB-ARC domain length less than 160 amino acids were removed. NB-ARC domain sequences were aligned with Mafft v. 7.480 (L-INS-i, iterative refinement method, max. 1000 iterations) [66]. Phylogenetic tree was reconstructed by maximum likelihood method using IQ-TREE v. 1.6.12 [67] (parameters -st AA m TESTNEW -bb 1000 -alrt 1000; JTT+F+R7 substitution model was chosen automatically). SeqKit v0.16.1 [73] and HmmPy (https://github.com/EnzoAndree/HmmPy; accessed on 7 April 2021) tools were used for sequence data processing at various stages of this analysis.

Assessment of Resistance to P. infestans
The results of evaluation of potato cultivars resistance to P. infestans isolate VZR18 are presented in Table 1. On the 7th day after inoculation (DAI), plants of the susceptible cultivar Bintje in two independent replications were severely affected, and on the 13th day, the almost died (score 1), indicating sufficient infection load and favorable conditions for disease development, allowing a reliable assessment of resistance. Cultivars Sudarynya and Siverskij showed resistance to isolate VZR18. S. stoloniferum clone 30514/15 and had no symptoms of the disease. On the 13th day, the cultivar Evraziya was affected at the level of susceptible cv. Bintje (Table 1).
The multidimensional scaling analysis demonstrates that libraries of the same genotype and conditions are clustered on the MDS plot ( Figure 1). The transcriptomes corresponding to control (0 hpi) and P. infestans infection (24 hpi) were clustered separately. As expected, the libraries from S. stoloniferum clone 30514/15 differ from cultivated potato genotypes. The changes in transcriptomes of cultivated potato genotypes revealed two clusters: smaller changes in response to P. infestans were characteristic for cv. Evraziya and cv. Bintjie, and more pronounced changes were detected for cv. Siverskij and cv. Sudarynya. This difference could be explained by a genotype-specific response to pathogen invasion.

Libraries Comparison and Quality Evaluation
The libraries comprise approximately 2.02 billion paired reads (298 billion bases; statistical metrics presented in more details in Table S2). After filtering by quality, 1.81 billion short reads (90% of raw libraries) comprising 269 billion bases remained. The number of bases with q > 30 in filtered reads varies from 93.4% (Bintje, 0 h, library 1) to 94.9% (Sudarynya, 24 hpi, library 3). The average number of paired reads per library is approximately 60.4 million, the average lengths of reads is 148 bases, and the average number of bases with Q > 30 is 94.
The multidimensional scaling analysis demonstrates that libraries of the same genotype and conditions are clustered on the MDS plot ( Figure 1). The transcriptomes corresponding to control (0 hpi) and P. infestans infection (24 hpi) were clustered separately. As expected, the libraries from S. stoloniferum clone 30514/15 differ from cultivated potato genotypes. The changes in transcriptomes of cultivated potato genotypes revealed two clusters: smaller changes in response to P. infestans were characteristic for cv. Evraziya and cv. Bintjie, and more pronounced changes were detected for cv. Siverskij and cv. Sudarynya. This difference could be explained by a genotype-specific response to pathogen invasion.

De novo Transcriptome Assembly and Reconstruction of the Orthologous Groups
Transcriptomes were assembled de novo (  (Table S4).
A search of orthologous groups in the reconstructed transcriptomes and in the reference transcriptomes of S. tuberosum and S. lycopersicum identified 53,775 orthogroups. A total of 172 orthogroups contain exactly one sequence in each of the seven transcriptomes (single-copy orthogroups). Based on the alignments of protein sequences in these orthogroups, the phylogenetic tree was reconstructed with the aid of maximum likelihood method ( Figure 2). nent 1 (X axis) and 2 (Y axis). Libraries from five genotypes are marked by distinct symbol shapes as shown on the legend in the top right part of figure.

De novo Transcriptome Assembly and Reconstruction of the Orthologous groups
Transcriptomes were assembled de novo (Table S3). The assembly parameter N50 varies from 2337 to 2444. Transcript lengths median values vary from 1691 to 1742 between the five reconstructed transcriptomes. After redundancy filtering, 32% of transcripts were removed on average.
A search of orthologous groups in the reconstructed transcriptomes and in the reference transcriptomes of S. tuberosum and S. lycopersicum identified 53,775 orthogroups. A total of 172 orthogroups contain exactly one sequence in each of the seven transcriptomes (single-copy orthogroups). Based on the alignments of protein sequences in these orthogroups, the phylogenetic tree was reconstructed with the aid of maximum likelihood method ( Figure 2). As expected, three cultivars of common origin, Evraziya, Siverskij, and Sudarynya, form a well-separated subcluster. It is known that the pedigree of cv. Bintje (bred in the Netherlands in 1904) does not include crosses with wild potato species. S. stoloniferum 30514/15 represent a sister group for the cultivated potato genotypes and the reference genome. As expected, three cultivars of common origin, Evraziya, Siverskij, and Sudarynya, form a well-separated subcluster. It is known that the pedigree of cv. Bintje (bred in the Netherlands in 1904) does not include crosses with wild potato species. S. stoloniferum 30514/15 represent a sister group for the cultivated potato genotypes and the reference genome.

Identification and Analysis of the NBS-LRR Proteins
The analysis of amino acid sequences encoded by the de novo-reconstructed transcripts with the aid of NLR-parser tool [70] predicted the proteins containing complete or partial NBS-LRR domains ( Table 2).
The number of amino acid sequences with complete or partial NBS-LRR domain signatures is similar between four potato genotypes and varies from 960 (cv. Sudarynya) to 1188 (cv. Siverskij). The number of sequences with complete NBS-LRR domains is two to three times lower and varies from 334 (cv. Sudarynya) to 504 (S. stoloniferum 30514/15). Most of transcripts encoding complete NBS-LRR domains (>80% for each genotype) have Agronomy 2021, 11, 2426 9 of 19 an expression level above 1 TPM. Detailed information on predicted NBS-LRR proteins is available in File S1.

Novel Transcripts Encoding Proteins with NBS-LRR Domains
De novo-assembled transcripts encoding proteins with complete NBS-LRR domains ( Table 2) were analyzed for similarity to the reference S. tuberosum genome. It was found that a fraction of predicted NLR proteins has no counterparts in the genome ( Table 3). The origin of these novel transcripts may relate to the introgression of resistance genes from wild and cultivated potato species as well as the intense recombination within the NBS-LRR genes clusters [74,75]. Transcriptomes of cvs. Siverskij, Evraziya, Sudarynya have similar numbers of the novel NBS-LRR transcripts, including the CNL/TNL classes, and differ from the cv. Bintje or S. stoloniferum 30514/15 transcriptomes ( Table 3). Comparison of transcriptomes taken before (0 hpi) and after P. infestans inoculation (24 hpi) revealed subsets of transcripts with either increased or decreased levels in response to the pathogen invasion (Table 3; data for novel transcripts are presented in Files S2 (description) and S3 (sequences)). As expected, S. stoloniferum 30514/15 was characterized with much higher number of P. infestans-induced novel NLR transcripts.
The results shown in Table 4 demonstrate that the transcriptomes of genotypes Siverskij, Sudarynya, and Euraziya encode proteins nearly identical to the R3b gene from S. demissum (identity > 99%). The transcriptomes of cv. Siverskij and cv. Sudarynya bear transcripts highly similar to Rpi-blb1 from S. bulbocastanum (or Rpi-sto1 from S. stoloniferum), identity >99% at amino acid level. Interestingly, cv. Euraziya transcriptome does not contain mRNA for RB/Rpi-blb1/Rpi-sto1 either at 0 or at 24 hpi despite the genome of this cultivar containing corresponding intragenic markers [40,41]  To characterize the novel NLR-encoding transcripts, the phylogenetic tree was reconstructed for their NB-ARC domains. The initial dataset contained 1052 sequences, including 608 novel NBS-LRR genes found in this research (description is available in File S3), 438 sequences that were taken from [13], and a few sequences for wide-range resistance genes against P. infestans (description is available in Table S1). After filtering (extraction of partial NB-ARC domains, removal of identical sequences), 812 non-redundant NB-ARC domains were used for tree reconstruction (File S4). The tree was rooted at CNL-R node. The part of tree with phylogenetic relationships in CNL6, 7, and 8 classes is shown in Figure 3. The full phylogenetic tree for NB-ARC domains from known and novel sequences and marker genes is presented in Newick, PDF, and SVG formats in File S4.
The distribution of NB-ARC domains of novel NLR proteins between the ten classes from [13] is non-uniform ( Figure 4). The CNL-8 and TNL classes contain many novel sequences, while CNL-2 has only a few, and the CNL-R class does not contain novel sequences at all. The CNL-8 class includes the R3a protein from S. tuberosum and R3b-like protein from S. stoloniferum 30514/15. The CNL-6 class includes Rpi-sto1 of S. stoloniferum and Rpi-blb1 of S. bulbocastanum with identical NB-ARC domains. The R8 gene from S. demissum does not fall into these classes but into the large cluster sister to the CNL-1 class.
The novel NLR proteins of different potato cultivars under investigation (Table 3) were compared with the aid of the CD-hit program (proteins with identity >99% were combined).
One may see ( Figure 5) that the intersections are rather small; that is, most of the novel NLR-like proteins are genotype specific. As expected, related genotypes (cv. Sudarynya, cv. Siversky, cv. Evraziya) have a larger number of common NLR proteins than cv. Bintje or S. stoloniferum 30514/15.   The novel NLR proteins of different potato cultivars under investigation (Table 3) were compared with the aid of the CD-hit program (proteins with identity >99% were combined).
One may see ( Figure 5) that the intersections are rather small; that is, most of the novel NLR-like proteins are genotype specific. As expected, related genotypes (cv. Sudarynya, cv. Siversky, cv. Evraziya) have a larger number of common NLR proteins than cv. Bintje or S. stoloniferum 30514/15.

Discussion
NLR genes play a vital role in plant defense against various pathogens. Rapid coevolution of pathogens and host plants results in a huge natural diversity of NBS-LRR receptors intensively screened for new valuable R genes for crop improvement. In addition, some incompatible NLR gene combinations might decrease plant fitness (e.g.,  The novel NLR proteins of different potato cultivars under investigation (Table 3) were compared with the aid of the CD-hit program (proteins with identity >99% were combined).
One may see ( Figure 5) that the intersections are rather small; that is, most of the novel NLR-like proteins are genotype specific. As expected, related genotypes (cv. Sudarynya, cv. Siversky, cv. Evraziya) have a larger number of common NLR proteins than cv. Bintje or S. stoloniferum 30514/15.

Discussion
NLR genes play a vital role in plant defense against various pathogens. Rapid coevolution of pathogens and host plants results in a huge natural diversity of NBS-LRR receptors intensively screened for new valuable R genes for crop improvement. In addition, some incompatible NLR gene combinations might decrease plant fitness (e.g.,

Discussion
NLR genes play a vital role in plant defense against various pathogens. Rapid coevolution of pathogens and host plants results in a huge natural diversity of NBS-LRR receptors intensively screened for new valuable R genes for crop improvement. In addition, some incompatible NLR gene combinations might decrease plant fitness (e.g., [26,27,29]). Plant breeding has been conducted for many decades, and cultivars may bear rather different subsets of NLR genes. In addition, NLR loci frequently containing gene clusters are characterized with a high recombination rate; these evolution hot spots may facilitate the selection of new efficient resistance genes in response to newly arising strains of pathogens. Here, we investigated the transcriptomes of potato cultivars combining genetic material from wild and cultivated potato species for the presence of NLR-related transcripts.
The first 24 h period after inoculation is considered the most critical for the interaction between potato and P. infestans (e.g., [76]); thus, the samples were analyzed either before or 24 h after inoculation with P. infestans. The transcriptomes were assembled de novo to take into account the pool of transcripts with no homology to reference genome, and NBS-LRR proteins were predicted computationally with the aid of an NLR-parser tool [70].
Three related cultivars, Evraziya, Siverskij, and Sudarynya, with the interspecific hybrid with wild Mexican species S. stoloniferum and S. demissum in their pedigrees, and two controls (susceptible cv. Bintje and resistant clone 30514/15 of the wild species S. stoloniferum) were taken for this analysis. Cultivars Sudarynya and Siverskij were evaluated as resistant to late blight in field trials in the 2016-2017 season (score 7-8) and to P. infestans isolate VZR18 with virulence formula 1.2.3.4.6.7.10.11 in laboratory assay (Table 1). Cv. Evraziya was evaluated as susceptible to isolate VZR18. The transcriptomes of four potato cultivars (Evraziya, Siverskij, Sudarynya, Bintjee) with different response to late blight were compared for the presence of NLR genes-related mRNAs.
The comparison of cultivars demonstrated the difference between transcriptomes taken before and after P. infestans inoculation. Interestingly, the S. tuberosum transcriptomes at 24 hpi were divided into two clusters: those susceptible to late blight cvs. Evraziya and Bintje versus resistant cvs. Sudarynya and Siverskij, which may reflect their differential response to P. infestans (Figure 1). The reconstructed phylogenetic tree for orthologous groups of genes also revealed a separate cluster containing a close group from cv. Sudarynia and cv. Siverskij together with cv. Evraziya, whereas cv. Bintje was closer to the reference genome ( Figure 2). This coincided with the fact that the pedigree of cv. Bintje (bred in the Netherlands in 1904) does not include crosses with wild potato species.
Actually, these results demonstrated the importance of transcriptomic evaluations for resolving the contradictory data between the marker-assisted selection and phenotyping. In some cases, the cultivar phenotypic resistance may differ from the DNA-marker-based prediction because of inappropriate expression pattern of the targeted R gene. For example, all the cultivars (Siverskij, Evraziya, Sudarynya) originating from hybrid clone 8889/3 (containing genetic material from S. stoloniferum) have the diagnostic intragenic markers of the RB/Rpi-blb1/Rpi-sto1 gene in their genomes (for detail, see [40]). It may be expected that this NLR (if functional) will contribute to late blight resistance in all three genotypes. However, despite the presence of the DNA-markers, cv. Evraziya transcriptome does not contain RB/Rpi-blb1/Rpi-sto1 mRNA, and this cultivar was susceptible to P. infestans in both field and laboratory evaluations. In our opinion, MAS results could be more precise if accompanied with transcriptomic assay or at least targeted RT-qPCR.
Despite the resistance of tested cultivars to the P. infestans isolate VZR18 and their field resistance to late blight in 2016-2017, it was later shown that in 2020 all these cultivars (Siverskij, Evraziya, Sudarynya) were severely affected by P. infestans in the field trials under conditions of high disease severity; they also appeared to be susceptible to P. infestans isolate MP1841 with virulence formula 1.2.3.4.5.6.7.8.9.10.11 (score 1-3, data not shown). Thus, the resistance of these cultivars turned out to be race-specific, which might be explained by insufficient expression level of RB/Rpi-blb1/Rpi-sto1. The last data indicate that for reaching the durable late blight resistance the stacks from several favorable genes and their simultaneous expression is important.
Indeed, the introgression of the RB/Rpi-blb1/Rpi-sto1 gene per se does not necessarily provide a potato cultivar with high-level resistance against P. infestans. It has been previously reported that transgenic potato cultivars bearing the RB/Rpi-blb1 gene from S. stoloniferum in its native form (8.5 kb long segment with promoter, introns, etc.) were resistant against P. infestans. However, the resistance level of transgenic lines varied greatly (susceptible, moderate susceptible, moderate resistant, resistant) and correlated with transgene copy numbers and expression rates [77]. Similar results were reported by another group [78].
Another outcome of this study concerns the considerable abundancy of NLR transcript absent in the reference genome. Interestingly, the transcriptomes of the potato cultivars tested contained significant numbers of novel NLR transcripts. The reconstructed phylogenetic tree (File S4, Figure 3) showed the distribution of novel transcripts according to the classification of NBS-LRR genes based on their domain structures [13]. These potential novel NLR proteins are distributed non-uniformly but are present in almost all groups except CNL-R. In fact, their structure is close to well-known functional NLRs, and they may be involved in regulation of various biological processes. Indeed, the functionality of novel NLRs needs further experimental verification.
The origin of novel NLR transcripts may relate to introgression of resistance genes from other species as well as to intense recombination within the NBS-LRR genes clusters [74,75]. Different subsets of NLR alleles were introgressed into the genepool of modern cultivars from many wild species of Central and South America as well as from cultivated species. A comparison of proteins encoded by novel NLR transcripts demonstrated a relatively small intersection between the genotypes ( Figure 5). This observation coincided well with their complex pedigrees, which contain hybridizations with introgressive lines originating from hybrid clones with S. stoloniferum, S. demissum, and S, acaule, S. tuberosum ssp. tuberosum [40]. The novel cultivar-specific NLR transcripts may also result from intense recombination inside the gene loci. It is likely that NLR gene loci represent some kind of playground for the rapid evolution of resistance genes to provide a reservoir for selection under the severe pressure of newly arising pathogen strains.
Investigation of genotype-specific NLR genes at the level of de novo-assembled transcriptome may considerably supplement genomic data because a cluster structure of NLR loci makes prediction at the genome level complicated. Taking novel NLRs into account may considerably extend the interpretation of major resistance genes mapping experiments.
One of the aims of this study was to compile a set of new candidate R genes (File S2) for further investigation. Association of a candidate R gene with the cultivar specific resistance against some pathogen provides opportunity for further study and inclusion in potato breeding programs or other biotechnological approaches for plant improvement (e.g., [3,[79][80][81][82][83][84][85]). Further investigation of candidate R genes may be conducted with the aid of reverse genetics (suppression or increased expression of a targeted gene under testing combined with phenotype evaluation). This way may be of use to reveal novel R genes in selected resistant clones from the genetic collections. In this work, the late blight resistance of analyzed potato cultivars was race specific. However, S. stoloniferum clone 30514/15 is highly resistant to late blight, and its transcriptome contained 41 novel mRNA encoding NLR proteins. Further detailed investigation of these mRNAs could reveal new NLRs that are potentially useful for the further improvement of S. tuberosum commercial cultivars. P. infestans is rapidly evolving [84] and may gain the resistance to fungicides [85], and some new strains are rather dangerous [86], which makes molecular breeding quite pertinent [87]. The transcriptomic approach has proved itself as a useful tool to reveal new potential R genes for further introgression.

Conclusions
Three potato cultivars (Evraziya, Siverskij, Sudarynya) combining genetic material from wild and cultivated potato species and bearing the same intragenic markers of the RB/Rpi-blb1/Rpi-sto1 gene were investigated. It was found that: 1. The transcriptomic assay may resolve the contradictive results of marker-assistedselection and phenotypic data. Here, both resistant and susceptible cultivars have the same intragenic markers of the RB/Rpi-blb1/Rpi-sto1; however, the induction of RB/Rpi-blb1/Rpi-sto1 transcript was found only in cultivars resistant to P. infestans isolate VZR18.
2. The significant fractions of transcriptomes analyzed (cvs. Evraziya, Siverskij, Sudarynya, Bintjee, S. stoloniferum clone 30514/15) represent novel NLR genes with no homology to the reference genome. Taking these NLRs into account may extend the limitations of the reference genome and increase the efficiency of novel R genes mapping.
Further functional analysis of novel NLR genes may provide a valuable resource for breeding. The computational data analysis was conducted using the resources of the ICG "Bioinformatics" Joint Computational Center. Data Availability Statement: Raw sequencing data were deposited in the NCBI Sequence Read Archive (PRJNA765923).