\documentclass[10pt]{article}

\usepackage{fullpage}
\usepackage{setspace}
\usepackage{parskip}
\usepackage{titlesec}
\usepackage[section]{placeins}
\usepackage{xcolor}
\usepackage{breakcites}
\usepackage{lineno}
\usepackage{hyphenat}





\PassOptionsToPackage{hyphens}{url}
\usepackage[colorlinks = true,
            linkcolor = blue,
            urlcolor  = blue,
            citecolor = blue,
            anchorcolor = blue]{hyperref}
\usepackage{etoolbox}
\makeatletter
\patchcmd\@combinedblfloats{\box\@outputbox}{\unvbox\@outputbox}{}{%
  \errmessage{\noexpand\@combinedblfloats could not be patched}%
}%
\makeatother


\usepackage{natbib}




\renewenvironment{abstract}
  {{\bfseries\noindent{\abstractname}\par\nobreak}\footnotesize}
  {\bigskip}

\titlespacing{\section}{0pt}{*3}{*1}
\titlespacing{\subsection}{0pt}{*2}{*0.5}
\titlespacing{\subsubsection}{0pt}{*1.5}{0pt}


\usepackage{authblk}


\usepackage{graphicx}
\usepackage[space]{grffile}
\usepackage{latexsym}
\usepackage{textcomp}
\usepackage{longtable}
\usepackage{tabulary}
\usepackage{booktabs,array,multirow}
\usepackage{amsfonts,amsmath,amssymb}
\providecommand\citet{\cite}
\providecommand\citep{\cite}
\providecommand\citealt{\cite}
% You can conditionalize code for latexml or normal latex using this.
\newif\iflatexml\latexmlfalse
\providecommand{\tightlist}{\setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}}%

\AtBeginDocument{\DeclareGraphicsExtensions{.pdf,.PDF,.eps,.EPS,.png,.PNG,.tif,.TIF,.jpg,.JPG,.jpeg,.JPEG}}

\usepackage[utf8]{inputenc}
\usepackage[ngerman,greek,english]{babel}



\usepackage{float}








\begin{document}

\title{Phenotypic resistance diversity underpinned by a diverse repertoire of
candidate NLR loci and genotype-specific expression patterns}



\author[1]{Pezhman Safdari}%
\author[1]{Layla Höckerstedt}%
\author[1]{Mikael Brosche}%
\author[2]{Jarkko Salojärvi}%
\author[1]{Anna-Liisa Laine}%
\affil[1]{University of Helsinki}%
\affil[2]{Nanyang Technological University}%


\vspace{-1em}



  
  \date{May 08, 2024}


\begingroup
\let\center\flushleft
\let\endcenter\endflushleft
\maketitle
\endgroup





\selectlanguage{english}
\begin{abstract}
High levels of phenotypic variation in resistance appears to be nearly
ubiquitous across natural host populations. To date, the molecular
processes associated with this variation in nature are poorly known,
although theory predicts resistance to evolve at specific loci in
response to pathogen -imposed selection. Nucleotide-binding leucine-rich
repeat (NLR) genes play an important role in pathogen recognition,
downstream defense responses and defense signaling. Exploring the depth
of NLR natural variation has the potential to increase our understanding
of how NLR diversity is generated and maintained, and how to manage
disease resistance. Here, we sequenced the transcriptomes of five
different Plantago lanceolata genotypes in response to inoculation by
the same strain of obligate fungal pathogen Podosphaera plantaginis. A
de novo transcriptome assembly of RNA-sequencing data yielded 24332 gene
models with N50 value of 1329 base pairs and gene space completeness of
66.5\%, suggesting a high-quality assembly. The gene expression data
showed highly varying responses, with each plant genotype having a
unique gene expression profile in response to the pathogen, regardless
of the resulting resistance phenotype. We discovered a diverse NLR
repertoire in P. lanceolata which is consistent with high phenotypic
resistance diversity and high genetic diversity of the species. We find
evidence of selection generating diversity at some of the NLR loci.
Jointly, our results demonstrate that phenotypic resistance diversity
results from a crosstalk between different defense mechanisms. In
conclusion, characterizing the architecture of resistance in natural
host populations may shed unprecedented light on the potential of
evolution to generate variation.%
\end{abstract}%



\sloppy


\section*{Key words:}

{\label{key-words}}

Phenotypic variation, pathogen imposed selection, phenotypic resistance
diversity, gene expression profile

\section*{INTRODUCTION}

{\label{introduction}}

Pathogens and herbivores play an important role in maintaining
biodiversity in natural populations (Bever, Mangan, \& Alexander, 2015;
Kursar et al., 2009). The threats imposed by pathogens on humans and on
managed, food production systems have motivated research that aims to
predict where pathogens will occur and how risks of infection evolve
(Gilligan, 2002; Koff, 1992; Woolhouse, Taylor, \& Haydon, 2001).
Pathogens can only occur where they have susceptible hosts, and hence,
resistance diversity is the key determinant of disease dynamics. Thus,
our ability to understand how diversity in resistance is generated and
maintained underlies our ability to predict and prevent disease
emergence and epidemics. In agriculture increasing the diversity of
crops - even from a monoculture to a mixture of two cultivars - has been
shown to reduce disease levels significantly (Mundt, 2002b; Zhu et al.,
2000). Natural host populations typically support diversity in
resistance phenotypes (Laine, Burdon, Dodds, \& Thrall, 2011; Salvaudon,
Giraud, \& Shykoff, 2008), and limited data available to date show that
increasing resistance diversity decreases disease risk also in the wild
(Jousimo et al., 2014a).

Hosts and pathogens are assumed to coevolve through Red Queen dynamics,
where the pathogen overcomes host's defenses and the host in turn
responds with new counter-defenses (Hamilton, 1980; Jaenike, 1978).
Theory predicts such reciprocal coevolutionary selection to be a
powerful mechanism for maintaining diversity in both host and parasite
populations, as the selection rate for resistance depends on the
frequency of parasite alleles, and vice versa, in a negative indirect
frequency-dependent manner (Bergelson, Kreitman, Stahl, \& Tian, 2001;
Leonard, 1977). There are numerous examples of pathogens overcoming host
resistance mechanisms, both from agriculture and from the wild (Mundt,
2002a, 2014). While evidence of resistance evolving under pathogen
attack in the wild is scarce (Laine, 2006), there is ample support for
coevolution from local adaptation studies where parasite/host fitness is
measured in sympatry vs. allopatry (Greischar \& Koskella, 2007;
Hoeksema \& Forde, 2008). To date, a handful of ground-breaking studies
have demonstrated that fluctuations in resistance and infectivity in
natural systems match the predictions of coevolutionary selection
(Decaestecker et al., 2007; G\selectlanguage{ngerman}ómez \& Buckling, 2011; Thrall et al.,
2012).

The interaction between plants and their pathogens is mediated by
complex defense mechanisms, incorporated in several layers of defense.
Thick and waxy cell walls are an example of a mechanical defense barrier
against pathogen invasion (Miedes, Vanholme, Boerjan, \& Molina, 2014).
Upon the arrival of a pathogen, the pathogen-associated molecular
patterns (PAMPs) trigger the so-called PAMP-triggered immunity (PTI)
response, which can stop the pathogen infection even before it begins.
If the pathogen succeeds in overcoming these first two defense layers,
the next layer is effector triggered immunity (ETI), involving the
recognition of pathogen virulent factors (effector proteins), either
directly or indirectly by recognizing the modifications that they make
to the plant proteins (Jones \& Dangl, 2006). Many of the proteins
involved in intracellular pathogen recognition belong to
nucleotide-binding--leucine-rich repeat (NLR) receptors (Monteiro \&
Nishimura, 2018). After pathogen recognition a multitude of different
signaling pathways, including production of reactive oxygen species,
elevated Ca\textsuperscript{2+} and MAP kinases lead to activation of
plant defenses. These defenses include the induction of stress hormones
salicylic acid, jasmonic acid and ethylene, as well as extensive
transcriptional re-programming ultimately resulting in the production of
defensive compounds, such as toxic secondary metabolites and
antimicrobial chemicals and enzymes. Against some pathogens, plants also
activate the hyper-sensitive response, programmed cell death, in which
the plant rapidly kills the cells surrounding the infection to prevent
the spread to nearby tissues (Coll, Epple, \& Dangl, 2011; Egorov \&
Odintsova, 2012).

Within these defense responses, NLR proteins play a key role. They are
involved in recognition of the pathogen's effector proteins both
directly and indirectly, as well as in triggering the plant immune
response (Meunier \& Broz, 2017). NLRs have also been shown to be
involved in signaling and transcript regulation (Chisholm, Coaker, Day,
\& Staskawicz, 2006; Jones \& Dangl, 2006). Moreover, NLR genes play an
important role in local adaptation and habitat expansion of plants (Stam
et al., 2017; Thrall et al., 2012). The antagonistic interaction between
plant NLR and pathogen effector proteins is considered to have a
profound effect on the evolution of both organisms, shaping their
genomes and gene repertoire (Upson, Zess, Bialas, Wu, \& Kamoun, 2018).
NLRs usually form large tandemly arrayed gene families and hence
questions regarding their origins and evolutionary history have been
under active research in both plants and animals (Andolfo et al., 2019;
Borrelli et al., 2018). The numbers of identified NLRs differ
substantially within and between plant families (Baggs, Dagdas, \&
Krasileva, 2017), for example \emph{Arabidopsis thaliana} (Arabidopsis)
contains between 165 to 251 NLR genes (Shao et al., 2016; Van de Weyer
et al., 2019) and crop species such as wheat, barley, rice, tomato and
potato contain 627, 224, 438,137 and 309 NLRs, respectively (Sarris,
Cevik, Dagdas, Jones, \& Krasileva, 2016). In \emph{A. thaliana} there
is evidence of widespread positive selection in the core NLRs shared
among accessions, especially in the canonical NLR domains (Van de Weyer
et al., 2019), while a pioneering study on wild tomato revealed high NLR
diversity with a small subset of NLR genes driving local adaptation to
pathogens (Stam et al., 2017).

A current key challenge in molecular ecology is to understand the role
of pathogen-imposed selection on generating NLR diversity. Exploring the
breadth and depth of plant NLR natural variation has the potential to
increase our understanding of how NLR diversity is generated and
maintained, and to establish a toolbox of deployable disease resistance
traits (Monteiro \& Nishimura, 2018). In natural plant populations,
neither pathogen epidemiology nor host resistance is under human
management in contrast to agricultural systems where disease is managed
both via resistance breeding and fungicides. Hence, natural populations
have the potential to offer unique insight into the processes generating
NLR diversity. Our study is focused on the interaction
between\emph{Plantago lanceolata} and its fungal pathogen
\emph{Podosphaera plantaginis} . Previous studies have detected
considerable phenotypic variation in \emph{P. lanceolata} resistance
against \emph{P. plantaginis} (Laine, 2004); diversity is shown to
accumulate in the well-connected populations across the landscape
(Hockerstedt, Siren, \& Laine, 2018), and has a direct negative impact
on disease dynamics (Jousimo et al., 2014b). Moreover, there is evidence
of coevolution in this interaction (Laine, 2005, 2006, 2008).

Here, we carried out a controlled experiment where five \emph{P.
lanceolata} individuals were inoculated with the same \emph{P.
plantaginis} strain. The aim of our study was to characterize the
transcriptional responses and the expression pathways activated in
response to the pathogen treatment, as well as in resistant versus
susceptible phenotypes. For this purpose, we carried out the
first\emph{de novo} transcriptome assembly for \emph{P. lanceolata} and
used it to characterize the induced NLR repertoire in \emph{P.
lanceolata}and test whether any particular NLRs show signs of selection.
We found that each plant genotype demonstrated a unique gene expression
profile in response to the pathogen and discovered a diverse NLR
repertoire in\emph{P. lanceolata} which is consistent with the high
phenotypic resistance diversity uncovered in earlier studies.

\section*{MATERIALS AND METHODS}

{\label{materials-and-methods}}

Study system and plant and fungal material

Ribwort plantain, \emph{Plantago lanceola} L., is a perennial monoecious
plant that reproduces both sexually by wind pollination, and clonally by
producing side rosettes (Sagar \& Harper, 1964). ~\emph{Podosphaerea
plantaginis} (Castagne; U. Braun and S. Takamatsu) (\emph{Erysiphales} ,
Ascomycota) is a specialist obligate biotroph infecting \emph{P.
lanceolata.} As all powdery mildews, it requires living host tissue
throughout its life cycle (Bushnell, 2002), and completes its life cycle
as localized lesions on host leaves. Infected plants suffer significant
stress, and infection may increase host mortality (Laine, 2004). The
interaction between \emph{P. lanceolata} and \emph{P. plantaginis} is
strain-specific suggesting gene-for-gene type control (Laine, 2004,
2007; Thompson \& Burdon, 1992). In some cases~the host can mitigate
pathogen´s reproduction; the putative resistance mechanism includes two
steps, where the host plant first recognizes the attacking pathogen and
blocks its growth (Laine, 2004) and the following infection outcome
depends on both host and pathogen genotypes (Laine, 2004, 2007).

In resistant interactions no infection develops, while in susceptible
interactions there is considerable variation in pathogen development
that is affected by both host and pathogen genotype (Laine, 2007). An
inoculation protocol where conidia from small colonies or individual
chains are placed on detached leaves or intact leaves of plants yields a
robust characterization of resistance-susceptibility. A plant genotype
is characterized as resistant when no pathogen growth is detected
following inoculation, or when the test plant shows rapid cell death
around inoculum source. A genotype is characterized as susceptible when
infection is detected following inoculation. From an earlier large
inoculation study consisting of 2944 host genotype--pathogen genotype
combinations (Hockerstedt et al., 2018), we selected three genotypes
(IDs 193\_2, 2818\_3 and 2818\_6, named Res1, Res2 and Res3 here after)
that were resistant against all tested pathogen strains, and two
genotypes (IDs 313\_6, 1553\_5, named Sus1 and Sus2 here after) that
were susceptible against all tested pathogen strains. Selected genotypes
were cloned into six plants each as described in Laine (2004).

\section*{Inoculation experiment}

{\label{inoculation-experiment}}

A schematic overview of the experiment design is shown in Figure 1a.
When plant clones were two months old, three replicates of each of the
five genotypes, in total 15 plants, were exposed to single \emph{P.
plantaginis} strain (Lammi\_3) by brushing spores gently with a fine
paintbrush onto six test leaves and two positive control leaves.
Similarly, three replicates of each genotype, in total 15 plants, were
mock inoculated by brushing leaves without mildew spores. Both powdery
mildew and mock-inoculated leaves were marked with a piece of masking
tape. Inoculated and mock-inoculated plant clones were placed in two
separate growth chambers (Panasonic MLR-352) at 20 ± 2 °C (day) and 16 ±
2 °C (night) with 16:8 light-darkness (L:D) photoperiod, and were
randomly organized to minimize potential variation in microclimatic
conditions. We collected two inoculated or mock-inoculated leaves from
every plant at 24, 48 and 72 hours post inoculation (hpi), snap froze
the leaves in liquid nitrogen, and stored samples in glassine bags in
-80 °C until RNA was extracted. Control leaves in whole plants were
screened until 14 days post inoculation to detect putative powdery
mildew infection to confirm the phenotypic response of the plants,
resistant or susceptible. Viability of spores used in the experiment was
confirmed by inoculating detached leaves of a susceptible genotype.

\section*{RNA extraction}

{\label{rna-extraction}}

For total RNA extraction, 0.2 g of frozen leaf material was ground in
lysing buffer (2\% CTAB, 2\% PVP K-30, 100 mM Tris-HCl pH 8.0, 2 M NaCl,
25 mM EDTA, \selectlanguage{greek}β-\selectlanguage{english}MeOH (200 \selectlanguage{greek}µ\selectlanguage{english}l/10ml) added in prior to use (Chang, Puryear,
\& Cairney, 1993). Thoroughly vortexed solution was then extracted twice
with equal volume of acid phenol-chloroform-isoamyl-OH (ph 4.5). Prior
to precipitation, 160 \selectlanguage{greek}µ\selectlanguage{english}l of 10M LiCl was added and samples were kept on
ice overnight, followed by 30 min centrifugation (10000 rpm) in +4 \selectlanguage{ngerman}°C.
Pellets were dissolved in 500 \selectlanguage{greek}µ\selectlanguage{english}l of 65 \selectlanguage{ngerman}°C SSTE (1M NaCl, 0.5\% SDS, 10
mmTris-HCl pH 8.0, 1mM EDTA) and RNA was extracted twice with
chloroform-Isoamyl alcohol (24\selectlanguage{english}:1). After EtOH precipitation and 70\%
wash, the pellets were dissolved in 40 \selectlanguage{greek}µ\selectlanguage{english}l H\textsubscript{2}O and RNA
quantity and quality were checked using NanoDrop (Thermo Fischer
Scientific). Potential contamination of genomic DNA in the RNA samples
was removed using DNase I (Thermo Fischer Scientific) and samples were
then reverse-transcribed to cDNA using iScript cDNA Synthesis Kit
(Bio-Rad) according to the manufacturer's instructions.

\section*{Selecting the time point for RNA-Seq using Quantitative
real-time PCR
(qPCR)}

{\label{selecting-the-time-point-for-rna-seq-using-quantitative-real-time-pcr-qpcr}}

In order to select the most informative time point for RNA sequencing to
detect significant transcriptomic responses induced by the inoculation
treatment, we initially performed a gene expression analysis of selected
marker genes using qPCR. We selected inoculated and mock-inoculated
clones of two plant genotypes (resistant 193\_2.1 and susceptible
1553\_5.1), and samples from all three time points (24, 48 and 72 hours
post inoculation), resulting in 12 samples. Primers were designed with
Primer3 (Rozen \& Skaletsky, 1999) based on previously \emph{in
situ}sequenced transcriptome of \emph{P. lanceolata} (unpublished data)
and known disease-induced genes in Arabidopsis\emph{.} We tested seven
putative disease-induced genes (Supplementary File 1). Amplification
efficiencies (E) of the primer pairs were determined with five dilutions
(1 : 1, 1 : 4, 1 : 24, 1 : 124, 1 : 624) of template cDNA, where E =
10\textsuperscript{-1}/slope. The qPCR was performed with three
technical replicates, one water control and a plate control sample in a
384-well plate with 10 \selectlanguage{greek}μ\selectlanguage{english}L volume, using C1000 Thermal Cycler (Bio-Rad).
All samples were tested for genomic DNA contamination with -RT controls
prior to qPCR. Each reaction had 1 \selectlanguage{greek}μ\selectlanguage{english}L of the 1:4 diluted cDNA, 5 \selectlanguage{greek}μ\selectlanguage{english}L of
SYBR\selectlanguage{ngerman}® Green containing master mix (iQ\selectlanguage{english} SYBR(r) Green Supermix for qPCR;
Bio-Rad), 3 \selectlanguage{greek}μ\selectlanguage{english}L of nuclease-free water and 0.5 \selectlanguage{greek}μ\selectlanguage{english}L (10 \selectlanguage{greek}µ\selectlanguage{english}m) of each primer.
The cycle conditions were one cycle initiating at 95\selectlanguage{ngerman}°C for 3 min, 40
cycles at 95°C for 10 s, 60°C 30 s, and ending with melting curve
analysis. To normalize the qPCR data several different reference genes
were tested and Elongation factor\_CL4, GADPH\_28221 and Actin\_34737
were validated with geNorm to have stable expression in the samples used
in this study (Supplementary File 1). Calculation of relative expression
(CNRQ) and normalization was done in qBase+ 3.2.

\section*{RNA sequencing (RNA-Seq)}

{\label{rna-sequencing-rna-seq}}

Several studies of powdery mildew pathogen induced genes in host plants
have found the highest number of differentially expressed genes in later
time points (Li et al., 2016; Li, Dong, et al., 2019; Polonio et al.,
2019). The qPCR results revealed that gene expression of the marker
genes showed elevated levels in time point 72 h post inoculation (hpi;
Supplementary Figure 1). Hence, we selected samples from time point 72
hpi for RNA sequencing. The samples were sequenced using Illumina
NextSeq 500 platform in the Institute of Biotechnology of the University
of Helsinki using paired end sequencing. The length of the forward reads
was 78 and reverse reads 74 bases. The sequenced reads were first
trimmed and low quality reads were removed using Trimmomatic software
(Bolger, Lohse, \& Usadel, 2014), resulting in an average library size
of 14.6 million reads.

\section*{Transcriptome assembly}

{\label{transcriptome-assembly}}

We combined all libraries and assembled them \emph{de novo} using
Trinity software (Grabherr et al., 2011), SOAPdenovo-Trans (k-mer sizes
39, 41) (Xie et al., 2014) and Oases (k-mer sizes 39, 43 and 47)
(Schulz, Zerbino, Vingron, \& Birney, 2012). Next, we combined all the
resulting contigs and ran the EvidentialGene pipeline (Gilbert, 2013).
We then combined okayset and okayalt outputs of the pipeline and
clustered them using RapClust (Srivastava, Sarkar, Malik, \& Patro,
2016), and merged the genes within the clusters using Lace (Davidson,
Hawkins, \& Oshlack, 2017). To remove contamination, the resulting
contigs were aligned against NCBI non-redundant protein database
(Pruitt, Tatusova, \& Maglott, 2007) using BLAST, and only the
transcripts that had a best hit in plant kingdom were retained. We also
removed the transcripts mapping to ribosomal genes and having ambiguous
sites (Ns). Minimum read coverage of three was used for all the
assemblies.

\section*{Differential gene expression
analysis}

{\label{differential-gene-expression-analysis}}

All libraries (30 in total) were mapped to the combined transcriptome
assembly using kallisto software (Bray, Pimentel, Melsted, \& Pachter,
2016) with one hundred bootstrap replicates. The averages of bootstrap
replicates of Transcript per Million (TPM) (Li, Ruotti, Stewart,
Thomson, \& Dewey, 2010) values were used as counts for differential
gene expression analysis. We imported the count tables to R by tximport
package (Soneson, Love, \& Robinson, 2015). Principal Component Analysis
of all samples, both for combined and for each genotype separately was
carried out using plotPCA command in DESeq2 and using rgl package
(Adler, Nenadic, \& Zucchini, 2017) for 3D plot and ggplot2 (Wickham H,
2016) for 2D plots. We analyzed the samples at genotype and phenotype
levels (Res, Sus and Res\_vs\_Sus) by DESeq2 (Love, Huber, \& Anders,
2014) to obtain differentially expressed genes.

\section*{Redundancy analysis of the count
data}

{\label{redundancy-analysis-of-the-count-data}}

We used vegan package (Oksanen et al., 2018) for redundancy analysis
(RDA) of the count data to identify genetic variation significantly
associated with the key conditions in our experiment: Inoculation,
Genotype, and Phenotype as well as the interactions of
Genotype-by-Inoculation and Phenotype-by-Inoculation. The statistical
significance of each of the conditions was tested with a permutation
test (permtest) with 10000 permutations. The genes contributing to
significant differences along the RDA axes were identified using cut off
of three standard deviations (corresponding to two-tailed p-value =
0.0027 in Z test). The overlap among gene sets was analyzed using Venn
diagrams in R using venn package (Dusa, 2018).

\section*{Annotation and gene ontology (GO)
analysis}

{\label{annotation-and-gene-ontology-go-analysis}}

To obtain a functional annotation for the transcripts, we used blastp
(Camacho et al., 2009) to search for the best match among Arabidopsis
representative set of proteins (Berardini et al., 2015) downloaded from
TAIR FTP server
(ftp://ftp.arabidopsis.org/home/tair/Proteins/TAIR10\_protein\_lists/).
The retrieved best hit was assumed to be an Arabidopsis ortholog and the
functional annotation and GO category assignment of the ortholog was
transferred to the \emph{P. lanceolata} query transcript. The
Arabidopsis GO annotations were downloaded from the TAIR FTP server
(ftp://ftp.arabidopsis.org/home/tair/Ontologies/Gene\_Ontology/). Gene
ontology enrichment analysis was carried out using piano software
(Varemo, Nielsen, \& Nookaew, 2013) and the log\textsubscript{2} fold
changes and false discovery rate (FDR) adjusted P-values of
differentially expressed genes were imported from DESeq2 results.

The GO hierarchy enrichment tests were carried out using the GOATools
toolkit (Klopfenstein et al., 2018). To focus on signaling responses, we
selected the responses and signaling sub-branches of the Biological
process category. The GO enrichment results were plotted using R with
fold change values obtained from the piano package. The GO enrichment
analyses of RDA loadings was carried out with piano software (Varemo et
al., 2013) using RDA loadings as gene level statistics and GO categories
as gene set collection and mean statistics. Again, the results were
plotted in R.

\section*{Prediction of candidate
NLRs}

{\label{prediction-of-candidate-nlrs}}

To predict the candidate NLR (nucleotide binding leucine rich repeat)
resistance genes in our reference transcriptome, the NLR-Parser pipeline
was applied (Steuernagel, Jupe, Witek, Jones, \& Wulff, 2015). We picked
the highest scoring domain found per reading frame per transcript and
manually screened the ORFs containing that domain. The transcripts were
left out if the ORF was too short or if they were missing the
appropriate start and stop codons and BLAST queries did not identify
hits from NCBI nr database. The NLR transcripts were considered to be
expressed in a genotype level assembly if the count was greater than one
in at least two replicates. The Venn diagram of expressed and
differentially expressed NLR genes across all genotypes was made using
venn package (https://CRAN.R-project.org/package=venn) in R.

Evolutionary analysis of the NLR transcripts was carried out
using\emph{Antirrhinum majus} L. (snapdragon) proteins and coding
sequences (Li, Zhang, et al., 2019) as outgroup, since it is
evolutionarily the most recently diverged plant from \emph{P.
lanceolata} where full genome assembly is available. Multiple sequence
alignment of all of the NLRs in reference assembly was carried out using
MAFFT, and the phylogenetic tree was estimated using FastTree (Price,
Dehal, \& Arkin, 2010). ClusterPicker (Ragonnet-Cronin et al., 2013)
script with 90 percent similarity and genetic distance of 0.03 with gap
option was used to cut the phylogenetic tree into clusters. The longest
sequence in each cluster was used as BLAST queries against snapdragon
proteins and the closest hit was selected to represent the ancestral
state and added to the cluster. The phylogenetic tree was first
visualized using the newick format output tree produced by ClusterPicker
and the final visualization was carried out with ggtree package (Yu et
al., 2017) in R.

\section*{Neutrality test (dN/dS and H
statistic)}

{\label{neutrality-test-dnds-and-h-statistic}}

To identify the loci under positive selection, we performed dN/dS
analysis. For each cluster obtained from ClusterPicker, we used the most
similar homolog in snapdragon (\emph{Antirrhimum majus} L.; Li, Zhang,
et al. (2019)) as ancestral NLR ortholog sequence for the sequences in
the cluster and conducted multiple sequence alignment using MAFFT,
estimated phylogenetic tree with FastTree and reverse-transcribed the
alignment of the sequences. We calculated the dN/dS ratios (ratio
between non-synonymous mutations and synonymous mutations) using PAML
software (Yang, 2007) within each of the clusters and used ratio
\textgreater{}1 as an indicator of putative positive selection. We also
calculated per base dN/dS ratios since, due to functional constraints on
conserved protein domains, it is much more likely that certain regions
in a gene are under selective pressure rather than the whole gene.

We carried out the Fay \& Wu H statistic analysis (Fay \& Wu, 2000) on
NLR transcripts by mapping the Plantago reads to the reference
transcriptome assembly using ANGSD software (Korneliussen, Albrechtsen,
\& Nielsen, 2014). To identify loci under selective pressure, we
calculated the Fay \& Wu H (H) values within a sliding window of nine
nucleotides. We chose H values less than -3.5 to signify selection.
Functional domains within the transcripts were identified with
InterProScan (Jones et al., 2014) search with the online InterPro
service of the EMBL-EBI website (Mitchell et al., 2019). The significant
H value loci were plotted against the transcript domain structure in R.

\section*{RESULTS}

{\label{results}}

\section*{Inoculation experiment and
qPCR}

{\label{inoculation-experiment-and-qpcr}}

We detected powdery mildew spores growing on the inoculated, susceptible
plant clone leaves on day 14 post inoculation. None of the
mock-inoculated or inoculated resistant plant clones showed visible
disease symptoms at that time. Gene expression of the marker genes
varied considerably in the susceptible inoculated plants and showed
elevated levels only at time point 72 h post inoculation (Supplementary
Figure 1). This time point was chosen for further analysis.

\section*{Transcriptome assembly and expression
analysis}

{\label{transcriptome-assembly-and-expression-analysis}}

On average, the 90 transcriptome libraries produced around 14.6 million
reads. All libraries were pooled to construct a \emph{de novo} assembly
with Trinity, yielding 177,948 isotigs. The assembly was of high
quality, since the Busco score of universal single-copy genes (Seppey,
Manni, \& Zdobnov, 2019) was 87\% (including complete and fragmented
genes), but it likely contained many splice variants and allelic
variants, because the proportion of duplicated gene models was high
(46.3\%). Careful clustering and filtering (see M\&M) narrowed the set
down to 24,332 high quality non-redundant transcripts with an average
length of 1,858 bases. The procedure reduced the Busco score to 66.5\%,
but clearly removed the allelic variants, as only 2\% of gene models
remained duplicated. The filtered gene models had mostly low expression
counts and therefore were of low biological significance to the
experiment.

On average 75 \% of the transcriptome data mapped to the \emph{de
novo}assembly (Supplementary Table 1). Principal Component Analysis
(PCA) of TPM normalized gene expression data showed a clear grouping of
the samples by genotype (Figure 1g) along the first three PCs. These
first three PCs explained altogether 53 percent of the total variation,
illustrating that genotype effect is the main contributor to the
variation between samples. The inoculation treatment had a smaller but
marked effect on gene expression, which was clearly demonstrated in
genotype-specific PCA plots where the effect of the inoculation was
associated with PC1 and PC2 axes (Figure 1 b to f). For example, in
resistant R1 and R2 genotypes the variation explained by PC1 was 65 and
52 percent, respectively, and clearly separated the inoculated and
control plants (Figure 1 b and c).

\section*{Differential gene expression
analysis}

{\label{differential-gene-expression-analysis-1}}

The PCA analyses showed marked differences between the gene expression
profiles of the genotypes and their responses to the inoculation.
Similarly, the genes expressed differentially between mildew inoculated
and mock-inoculated plants showed high genotype-specific variation. The
R1 genotype had the highest number of differentially expressed (DE)
transcripts among all genotypes, 3803, from which about 2000 had
absolute log\textsubscript{2} fold change greater than one. On the other
hand, the S2 genotype had the lowest number of DE transcripts, 43, with
only 20 having absolute log\textsubscript{2} fold change greater than
one (Table 1; Supplementary Table 2).

\section*{Redundancy Analysis (RDA)}

{\label{redundancy-analysis-rda}}

To study the relative effect sizes and their statistical significance,
we carried out multivariate regression using Redundancy Analysis (RDA).
The analysis confirmed that the genotype and phenotype (resistant versus
susceptible) effects were highly statistically significant (P = 0.001
and P = 0.004), describing 35 \% and 9 \% of the overall variation
(Table 2). The effect of the inoculation treatment alone was not
significant (P = 0.238), but the combined effects of
genotype-by-experiment and phenotype-by-experiment were significant (P =
0.001 and P = 0.048, respectively), again suggesting genotype-specific
response profiles. Accordingly, the RDA plots displayed clear separation
when using genotype and phenotype as a covariate but not with
inoculation treatment (Figures 2a, b and c). Venn diagrams of the genes
contributing to the separation in the RDA analyses demonstrate that
while 87 genes contribute significantly to separation according to
genotype, only 7 genes contribute directly to the phenotypic variation
and 109 genes to the joint effect of phenotype-by-experiment (Figures 2
d and e).

\section*{Gene ontology (GO) analysis based on differential
expression}

{\label{gene-ontology-go-analysis-based-on-differential-expression}}

To look for common pathways that are activated by inoculation in all
plants, we carried out GO enrichment analysis of the differentially
expressed genes. The different responses by the genotypes were visible
also at the level of activated pathways (Figure 3; Supplementary Figure
2; Supplementary Table 3). This may be due to differences in plant
defense responses or manipulation of the plant defense mechanisms by the
pathogen.

To explore molecular underpinnings between susceptible and resistant
genotypes we searched for differential activation of defense response
pathways by identifying the GOs with decreased average expression levels
in susceptible phenotypes and elevated levels in resistant phenotypes.
In resistant phenotypes, genes encoding photosynthesis related proteins
and NAD(P)H dehydrogenase complex had increased transcript levels. This
could contribute in defense against the pathogen, as it has been shown
that photosynthesis plays an important role in plant defense against
biotic stress (Gohre, 2015). Genes assigned to photosynthesis functions
showed elevated transcript levels in susceptible phenotypes as well but
not to the same extent. Chlorosis is a hallmark sign of powdery mildew
infection and biotrophic fungi are known to reduce photosynthetic rate
and possibly damage chloroplast structure (Perez-Bueno, Pineda, \&
Baron, 2019), thus the upregulation could be either compensation, plant
defense mechanism or induced by pathogen. In addition, uroporphyrinogen
decarboxylase activity (GO:0004853) was upregulated in resistant
phenotypes. Involved in chlorophyll biosynthesis, it also points towards
acting against the chlorosis induced by the pathogen (Mock, Keetman,
Kruse, Rank, \& Grimm, 1998).

Among both susceptible phenotypes, the GO category with most decreased
expression levels was induction of programmed cell death (GO:0012502),
suggesting that as a biotrophic pathogen, \emph{P. plantaginis} may have
disabled the programmed cell death and is keeping the host cells alive.
However, also the resistant phenotypes showed reduced expression levels,
possibly due to successful manipulation by the pathogen, and therefore
the comparison between susceptible versus resistant did not identify
this process as significantly different between phenotypes (P=0.0559).

In addition to the shared responses, the genotypes showed individual
enrichment of various disease resistance pathways. In susceptible
genotype 1 (S1), the processes with most decreased average expression
levels were tripeptide transporter activity (GO:0042937), tripeptide
transport (GO:0042939) and delta12\selectlanguage{english}-fatty acid dehydrogenase activity
(GO:0016720), whereas S2 demonstrated decrease in Oxazole or thiazole
biosynthetic process (GO: 0018131) and Low\selectlanguage{english}-affinity nitrate transport
(GO:0080054 \& GO:0080055). Fatty acids play a direct role in modulating
the plant defense response to pathogens (Kachroo \& Kachroo, 2009), and
thiazole or thiamine has been shown to play a crucial role in activation
of the defense responses, callose/lignin deposition and stomatal closure
(Zhou, Sun, \& Xing, 2013).

Tripeptide transport includes also nitrate transporters. Interestingly,
powdery mildew causative agent \emph{Erysiphe necator} elevates the
expression levels of nitrate transporters in grapevine and Arabidopsis
(Pike et al., 2014), possibly to acquire nutrients from the host. In
addition to decreased levels of the GO categories related to nitrate
transport in both S1 and S2, we identified homolog of Arabidopsis
nitrate transporter (AtNRT1.5) to be upregulated after inoculation in
susceptible vs resistant comparison. In Arabidopsis, the protein is
responsible for nitrate transport from roots to shoots, and in this
context suggests towards manipulation of host nutrient distribution~by
the pathogen. Mur, Simpson, Kumari, Gupta, and Gupta (2017) have argued
that nitrogen and nitrates and their transportation to different tissues
in the plant during the pathogen infection could be the ``silver
bullet'' of the plant defense. Besides nitrate transport, tripeptide
transport activities also play an important role for defense against
biotic and abiotic stress (Karim et al., 2007), suggesting a reason for
the decreased expression of the tripeptide transporters as a whole in
the susceptible phenotypes.

In resistant phenotypes, the glucosyltransferase (GO:0050284)
upregulation in R1 is a sign of early preparation for pathogen response
(Le Roy, Huss, Creach, Hawkins, \& Neutelings, 2016), and in R2
genotype, the activation of NADH dehydrogenase complex assembly
(GO:0010258) has been shown to be involved in defense signaling
(Wallstrom et al., 2014).

\section*{Gene ontology (GO) analysis based on Redundancy Analysis
(RDA)}

{\label{gene-ontology-go-analysis-based-on-redundancy-analysis-rda}}

We identified the GO categories contributing to the differences in RDA
analyses by averaging the RDA loadings of the genes in each of the
categories, and tested whether the values significantly differ from
random scores obtained by permutation. Genes contributing to the
separation between inoculation and control were enriched for ABA and
cytokine signaling, primary metabolism and chloroplast activity. ABA
induces resistance to powdery mildew in barley (Wiese, Kranz, \&
Schubert, 2004), and repression of ABA biosynthesis as well as genes
regulated by ABA, such as cold/dehydration/salinity responsive genes,
are associated with mildew resistance in nonhost plants in general
(Jensen et al., 2008). Cytokinin suppresses programmed cell death and
plays a role in the synthesis and maintenance of chlorophyll (Walters \&
McRoberts, 2006); both processes were enriched among the differentially
expressed genes. Additionally, cytokinin levels regulate cell division
together with auxin. Interestingly, in Arabidopsis, \emph{Golovinomyces
orontii} inoculation induced cell cycle related genes and
endoreduplication, possibly due to increased metabolic demands of the
pathogen (Chandran, Inada, Hather, Kleindt, \& Wildermuth, 2010). On the
other hand, Choi, Choi, Lee, Ryu, and Hwang (2011) have shown that plant
based cytokinins systematically induce plant resistance against
pathogens by cytokinin and salicylic acid signaling.

Genes associated with the differences between the phenotypes showed GO
enrichments for kinase activity, carbohydrate metabolism, plant cell
wall organization, photosystem II and response to cold GO categories,
whereas the genes contributing to the differences between genotypes were
enriched for tryptophan metabolism, plant cell wall and chloroplast
(Figure 4a and Figure 2d/e; Supplementary Figure 3). In Arabidopsis
(Chandran et al., 2010), the expression of cold/drought responsive genes
were decreased together with ABA biosynthesis after inoculation
with\emph{G. orontii} . Together with the observed induction of ABA
during inoculation, this suggests that the phenotypes may differ in how
strongly ABA activates its targets such as cold responsive genes.

The responses to infection start to play a role in
genotype-by-experiment effect where the genes are mainly enriched for
immune and defense responses, illustrating that the genes in these
processes highly differ by the transcription levels between genotypes.
Overall, the enriched GOs show a clear activation of defense responses
in general, and defense responses to fungi in particular. The GO
category with highest positive average of RDA loadings (and therefore,
high contribution to separation) is aldose 1\selectlanguage{english}-epimerase activity
(GO:0004034) which, based on Sheshukova et al. (2017), is activated
because of the mechanical damage inflicted by the pathogen and results
in methanol emission and priming of the non-infected leaves. The
presence of hydrogen peroxide metabolic process and salicylic acid
mediated signaling pathway are both well-established pathogen-induced
defense mechanisms (Hua, 2009; Kuniak \& Urbanek, 2000; Niu \& Liao,
2016; Sheshukova et al., 2017), further demonstrating the activation of
the defense processes due to the pathogen infection. The GO category
with most negative average RDA loadings is RNA splicing, via
endonucleolytic cleavage and ligation (GO:0000394). An emerging field in
plant molecular biology, it is becoming increasingly clear that plants
use alternative RNA splicing extensively as a means to respond to their
environment and defend against pathogens (Shang, Cao, \& Ma, 2017;
Staiger, Korneli, Lummer, \& Navarro, 2013). Within the
signaling-specific GOs (Figure 4b) the genotype by experiment effect
showed the activation of jasmonic acid (JA) and abscisic acid signaling
(as expected, (Yang et al., 2019)), but again, in a genotype-specific
manner. Further inspection of putative orthologs of marker genes for
different hormonal signaling pathways showed activation of auxin
biosynthesis and signaling, as well as differences in the activation of
JA signaling and NLR signaling through EDS1 ortholog (Supplementary
Table 4). In R1, specific induction of the SA biosynthesis and systemic
acquired resistance was also detected.

The most significant contributor to phenotype-by-experiment is
photosystem II activity. The GO category with highest average RDA
loadings for phenotype by experiment effect is oligopeptide
transmembrane transporter activity (GO:0035673). The perception and
transduction of fungal oligopeptides will trigger multiple defense
responses (Hahlbrock et al., 1995; Nurnberger et al., 1994). Multitude
of photosynthetic processes were also enriched; their role in defense
was discussed above. The categories with most negative average loadings
were response to fungus, and cytokinin biosynthetic process
(GO:0009691).

\section*{NLR transcripts}

{\label{nlr-transcripts}}

To look for the variation in the plant defense responses we carried out
an in-depth study of the resistance NLR genes induced in the experiment.
The transcriptome annotation predicted altogether 543 candidate NLR
transcripts, and 51 of them were upregulated and 35 downregulated in
resistant versus susceptible comparison. The clustering of the NLR
transcripts resulted in 58 clusters containing 457 sequences and 54
singletons (Figure 5a). Cluster 69 with 22 transcripts had the highest
number of transcripts. The BLAST query against Arabidopsis protein
database for the longest transcript in this cluster returned a hit to
the RPP13 gene. RPP13 has the highest amount of amino acid diversity in
Arabidopsis and is involved in defense against \emph{Peronospora
parasitica} (Hall et al., 2009; Rose et al., 2004), an oomycete causing
downy mildew in Brassicaceae. A closer look at the expressed NLR
transcripts revealed that although many of the NLR transcripts are
expressed across all genotypes, the differentially expressed NLRs differ
between genotypes with no common NLR transcripts, as illustrated by Venn
diagrams of all expressed and differentially expressed NLR transcripts
(Figures 5b and c).

\section*{Neutrality test (dN/dS and H
statistic)}

{\label{neutrality-test-dnds-and-h-statistic-1}}

To look for NLR clusters under positive selection, possibly in response
to coevolution with the pathogen, we analysed dN/dS, the ratio between
non-synonymous (amino acid changing) to synonymous mutations (Figure 5a,
outmost ring). None of the NLR transcript clusters had an \selectlanguage{greek}ω \selectlanguage{english}value
greater than one, which would indicate positive selection pressure.
However, site-wise analysis of dN/dS revealed that 23 of the clusters
contained a varying number of one to 12 amino acid positions under
positive selection, based on Bayes Empirical Bayes (BEB) analysis
(P\textgreater{}95\%). The cluster with the highest number of loci under
selection returned a BLAST hit to an Arabidopsis NLR protein responsible
for resistance against Lettuce mosaic virus (Revers et al., 2003).

In order to investigate potential selection pressure by a complementary
method, considering the shortcomings of within population dN/dS analysis
(Kryazhimskiy \& Plotkin, 2008), we also calculated Fay \& Wu H
statistics on the whole NLRs using putative snapdragon orthologs (Li,
Zhang, et al., 2019) as the ancestral state. A positive value of H is an
indication of balancing or purifying selection, suggesting the tendency
to preserve the current protein sequence, whereas high negative values
are an indication of diversifying positive selection, selective sweeps,
or drift, for example from population bottlenecks.

We identified 23 NLR transcripts with regions having H statistics less
than -3.5 (Figure 5a; Supplementary Figure 4; Supplementary Table 5).
This set included two genes from the cluster with homology to RPP13, and
the sites under selection in all but one were located in the conserved
NB-ARC domain. Since this domain is the activation domain (Van Ooijen et
al., 2008), the mutations may affect the strength of the activation.
Additionally, we identified selection in the LRR domains of three genes.
Since they are involved in pathogen perception and induced in our
experiment, these domains may have been under positive selection to
detect the \emph{P. plantaginis} pathogen.

\section*{DISCUSSION}

{\label{discussion}}

Given that pathogens are prevalent across all ecosystems, an
individual's reproductive success and survival will depend on its
ability to resist infection. Natural host populations have been shown to
support considerable diversity in resistance (Laine et al., 2011;
Salvaudon et al., 2008), and theory predicts that this variation is
maintained by pathogen-imposed selection. However, empirical support for
the role of selection in generating resistance diversity is to date
still scarce. With recent advances uncovering the molecular
underpinnings of resistance, it is becoming increasingly feasible to
study resistance also in non-model systems.

Here, we established a high-quality \emph{de novo} transcriptome
assembly of \emph{P. lanceolata} to investigate the gene expression and
processes activated in different plant genotypes in response to
inoculation of the same pathogen strain. In our study, all five plant
genotypes showed unique gene expression patterns. This was clearly
demonstrated in the principal component analysis where the genotypes
separate, while the inoculated and mock-inoculated replicates of the
same genotype cluster together. Significant variation in gene expression
patterns among plant genotypes has also been discovered in other studies
(Burghardt et al., 2017; Muller, Kersten, Fladung, \& Schroeder, 2019).
In the redundancy analysis (RDA), inoculation explains only 4 \% of the
total variation, while genotype-by-inoculation interactions contributes
46 \%, suggesting that the genotypes have highly unique responses to
pathogen attack. While such genotype specificity may be expected between
susceptible and resistant genotypes, the split to resistant versus
susceptible phenotypes explains only 9 \% of the variance, with
considerable transcript expression pattern differences between
phenotypes. Overall, the plant genotypes also differ by the number, fold
change and the function of the transcripts which are differentially
expressed in response to the pathogen. Furthermore, even though the gene
expression shows the known induction of JA, SA and ABA signaling
pathways, they show highly varying activation patterns with JA and ABA
significantly contributing to the genotype by experiment differences in
RDA. This suggests that different plant genotypes use different
strategies in response to the same pathogen, which could be an important
mechanism generating phenotypic resistance diversity. One possible
explanation for the diverse responses is the extremely high genetic
variation within the species; overall, the transcriptome had very high
Watterson \selectlanguage{greek}θ\selectlanguage{english}=0.068 and nucleotide diversity \selectlanguage{greek}π\selectlanguage{english}\textbf{=} 0.077, suggesting
effective population sizes in the order of millions. The high genetic
diversity where, on average, eight nucleotides out of 100 differ between
any two individuals, is likely manifested also in the diverse responses.
The experimental take home message of genotype-specific expression
pathways to treatment is that including multiple genotypes in
experiments and avoiding pooling for RNA-Seq is essential to uncover
variation relevant for phenotypic differentiation.

The pathways induced by the pathogen were also visible in the gene
expression data. The common responses included the induction of specific
nitrate transport genes in susceptible phenotypes as well as elevated
expression of photosynthesis-associated genes and related biological
processes taking part in chloroplast in all genotypes. The powdery
mildew fungi have a contracted carbohydrate metabolism, for example they
are not able to degrade pectin, an essential component of plant cell
walls (Liang et al., 2018), whereas the lipid metabolism is intact,
suggesting that their main source of energy is from lipids. Chlorosis is
another hallmark sign of a successful pathogen attack. In our results,
the fact that specific nitrate transporters show elevated expression
levels as well as chloroplast processes in general, suggests elevated
chlorophyll biosynthesis. Together, this suggests that, at least at the
early stage of infection, P\emph{. plantaginis} may target the
chloroplast lipids of Plantago to obtain its nutrients, most likely the
lipids involved in chlorophyll.

NLRs play an important role in pathogen recognition and downstream
defense responses, defense signaling, as well as activation of a hyper
sensitive response (Monteiro \& Nishimura, 2018). In our study, we
discovered 543 NLR isoforms from a transcriptome that is a combination
of the five different \emph{Plantago} genotype NLR repertoires. Majority
of these transcripts were expressed at least at some level in all five
plant genotypes. Presence-absence polymorphism in a subset of NLRs has
been demonstrated across Arabidopsis accessions (Van de Weyer et al.,
2019), and hence it could contribute to the slight differences in the
numbers of NLRs detected in the genotypes. Overall, the putative number
of NLR genes in \emph{P. lanceolata,} based on the identified 543 NLR
isoforms, seems to be quite high compared to the range reported in other
species (Baggs et al., 2017). The high count could result from
preferential retention of NLR genes after Plantaginaceae specific whole
genome duplication event 46-49 million years ago (Li, Zhang, et al.,
2019), since it has been proposed that numerous and highly diversified
defense genes are particularly favored in long-lived plants to avoid
infection by their rapidly evolving pathogens (Plomion et al., 2018).

The predicted NLR transcripts resulted in 58 clusters of varying sizes,
with 22 transcripts in the largest cluster. We also found considerable
variation in the branch lengths among clusters, which could indicate
variation among clusters in their likelihood to accumulate mutations
(Tucker, Ackerman, Eads, Xu, \& Lynch, 2013). Indeed, NLR genes are
among the fastest evolving gene families in plants. They often form
tandemly arrayed gene clusters, and this is believed to be critical for
the fast pace of their structural and functional diversification
(Meyers, 2003; Michelmore \& Meyers, 1998). Frequent homologous
recombination events and errors produced during the process, followed by
diversifying selection, may generate the structural diversity needed to
match high effector evolution rates characteristics of the microbes they
encounter (Jacob, Vernaldi, \& Maekawa, 2013; McDowell \& Simon, 2006).
NLR genes are also under evolutionary pressure resulting from
inappropriate activation of cell death. If the plant cannot control
NLR-activated cell death, it leads to decreased fitness (Phadnis \&
Malik, 2014). In particular, we found an increase in the number of
homologs of Arabidopsis RPP13 gene, which is involved in defense against
downy mildew (\emph{Peronospora parasitica} ) in Arabidopsis, as well as
other defense processes and signaling (Bittner-Eddy, Crute, Holub, \&
Beynon, 2000; Rentel, Leonelli, Dahlbeck, Zhao, \& Staskawicz, 2008).
While none of the clusters had significant dN/dS values, we found
between one and twelve loci under selective pressure in 23 of the
clusters. The cluster with the highest number of loci under selection
has been suggested to be involved in defense against viruses in
Arabidopsis (Revers et al., 2003). The H statistic analysis confirms the
transcripts identified by the dN/dS analysis (10 transcripts), as well
as identifies 13 other NLR transcripts that may have been under putative
selection pressure, two of which show selection in the LRR domain
responsible for detecting the type of pathogen.

We find that the NLR transcripts are differentially expressed in
response to the pathogen treatment, and that this response varies
according to genotype. Transcripts of many NLR genes are known to
accumulate in response to defense induction or related stimuli (Lai \&
Eulgem, 2018). For example, 75 of the 124 studied Arabidopsis NLR genes
were found to exhibit at least two-fold higher transcript levels in
response to one or more of the 15 implemented defence-related treatments
(Mohr et al., 2010). Up-regulation of NLR transcripts after defence
induction has also been observed in other plant species, such
as\emph{Brassica rapa} , soybean and rice (Brechenmacher et al., 2015;
Chen, Pang, Chen, Zhang, \& Piao, 2015; Ribot et al., 2008). While we
were more likely to observe up-regulation of NLR expression levels in
response to the pathogen treatment, this was not consistent across
transcripts and genotypes. This is in line with recent studies on crop
plants testing different genotypes in response to pathogen infection
(Cruz-Miralles, Cabedo-Lopez, Perez-Hedo, Flors, \& Jaques, 2019; Sari
et al., 2018; Sari, Bhadauria, Vandenberg, \& Banniza, 2017). Plants
have evolved mechanisms to stabilize their basal expression levels, and
to reduce the fitness costs of an overexpressed immune response that
could have more deleterious effects on plant fitness than the infection
(Fei, Xia, \& Meyers, 2013). This may explain the down-regulation of
some of the NLR transcripts we observe in both susceptible and resistant
phenotypes.

Our results here are well in line with the extensive phenotypic
variation, and highly strain-specific disease resistance measured
in\emph{P. lanceolata} in earlier studies (Hockerstedt et al., 2018;
Jousimo et al., 2014a). High levels of variation in resistance seems to
be nearly ubiquitous across natural host populations that experience
pathogen-imposed selection without any human interference, in contrast
to agricultural systems (Laine et al., 2011; Salvaudon et al., 2008). We
show that phenotypic resistance may be generated by different
mechanisms. First, we discovered a large repertoire of candidate NLRs
in\emph{P. lanceolata.} We also find evidence of selection generating
diversity in a subset of the identified NLRs. Moreover, we discovered
that the genotypes have unique expression profiles in response to
pathogens, a mechanism which may further contribute to phenotypic
variation. Indeed, this high level of genetic and expression profile
diversity may be the key to successful defense against pathogens in
sessile plants that lack a long lasting immune memory (Hall et al.,
2009; Roux \& Bergelson, 2016). Finding different mechanisms that
contribute to phenotypic resistance is nontrivial, given how effectively
this variation may be utilized to predict and control disease epidemics
(Mundt, 2002a). Moreover, resistance in agricultural crops is highly
prone to breakdown following pathogen adaptation, and for many
commercially important pathogens, the known effective resistance genes
are becoming limited. Wild plant populations are currently identified as
the most promising source of genes required for development of
sustainable agriculture (Fu et al., 2019). In conclusion, characterizing
the architecture of resistance in natural host populations may yield
unprecedented light on the potential of evolution to generate variation,
and it can have broad and long-lasting impacts in our food production
environments.

\section*{ACKNOWLEDGEMENTS}

{\label{acknowledgements}}

We would like to acknowledge Krista Raveala for maintaining the plant
and fungal material needed for the experiment. This work was funded by
grants from the Academy of Finland (296686) and the European Research
Council (Consolidator Grant RESISTANCE 724508) to A-LL, and Academy of
Finland (318288) to JS. The authors also wish to acknowledge CSC -- IT
Center for Science, Finland, for computational resources.

\section*{REFERENCES}

{\label{references}}

\begin{enumerate}
\tightlist
\item
  Adler, D., Nenadic, O., \& Zucchini, W. (2017). RGL: A R-library for
  3D visualization with OpenGL.
\item
  Andolfo, G., Di Donato, A., Chiaiese, P., De Natale, A., Pollio, A.,
  Jones, J. D. G., . . . Ercolano, M. R. (2019). Alien Domains Shaped
  the Modular Structure of Plant NLR Proteins. \emph{Genome Biology and
  Evolution, 11} (12), 3466-3477. doi:10.1093/gbe/evz248
\item
  Baggs, E., Dagdas, G., \& Krasileva, K. V. (2017). NLR diversity,
  helpers and integrated domains: making sense of the NLR
  IDentity.\emph{Current Opinion in Plant Biology, 38} , 59-67.
  doi:10.1016/j.pbi.2017.04.012
\item
  Berardini, T. Z., Reiser, L., Li, D., Mezheritsky, Y., Muller, R.,
  Strait, E., \& Huala, E. (2015). The Arabidopsis information resource:
  Making and mining the ''gold standard'' annotated reference plant
  genome. \emph{Genesis, 53} (8), 474-485. doi:10.1002/dvg.22877
\item
  Bergelson, J., Kreitman, M., Stahl, E. A., \& Tian, D. (2001).
  Evolutionary Dynamics of Plant Genes. \emph{Science, 292} (5525),
  2281-2285. doi:10.1126/science.1061337 \%J Science
\item
  Bever, J. D., Mangan, S. A., \& Alexander, H. M. (2015). Maintenance
  of Plant Species Diversity by Pathogens. \emph{Annual Review of
  Ecology, Evolution, and Systematics, 46} (1), 305-325.
  doi:10.1146/annurev-ecolsys-112414-054306
\item
  Bittner-Eddy, P. D., Crute, I. R., Holub, E. B., \& Beynon, J. L.
  (2000). RPP13 is a simple locus in Arabidopsis thaliana for alleles
  that specify downy mildew resistance to different avirulence
  determinants in Peronospora parasitica. \emph{The Plant Journal, 21}
  (2), 177-188. doi:10.1046/j.1365-313x.2000.00664.x
\item
  Bolger, A. M., Lohse, M., \& Usadel, B. (2014). Trimmomatic: a
  flexible trimmer for Illumina sequence data. \emph{Bioinformatics, 30}
  (15), 2114-2120. doi:10.1093/bioinformatics/btu170
\item
  Borrelli, G. M., Mazzucotelli, E., Marone, D., Crosatti, C.,
  Michelotti, V., Val\selectlanguage{ngerman}è, G., \& Mastrangelo, A. M. (2018). Regulation and
  Evolution of NLR Genes: A Close Interconnection for Plant
  Immunity.\emph{International Journal of Molecular Sciences, 19} (6).
  doi:10.3390/ijms19061662
\item
  Bray, N. L., Pimentel, H., Melsted, P., \& Pachter, L. (2016).
  Near-optimal probabilistic RNA-seq quantification. \emph{Nature
  Biotechnology, 34} (5), 525-527. doi:10.1038/nbt.3519
\item
  Brechenmacher, L., Nguyen, T. H., Zhang, N., Jun, T. H., Xu, D., Mian,
  M. A., \& Stacey, G. (2015). Identification of Soybean Proteins and
  Genes Differentially Regulated in Near Isogenic Lines Differing in
  Resistance to Aphid Infestation. \emph{Journal of Proteome Research,
  14} (10), 4137-4146. doi:10.1021/acs.jproteome.5b00146
\item
  Burghardt, L. T., Guhlin, J., Chun, C. L., Liu, J., Sadowsky, M. J.,
  Stupar, R. M., . . . Tiffin, P. (2017). Transcriptomic basis of genome
  by genome variation in a legume-rhizobia mutualism. \emph{Molecular
  Ecology, 26} (21), 6122-6135. doi:10.1111/mec.14285
\item
  Bushnell, W. R. (2002). The role of powdery mildew research in
  understanding host-parasite interaction: past, present and
  future.\selectlanguage{english}\emph{In The Powdery Mildews , A Comprehensive Treatise
  (Belanger, R. R., Bushnell, W. R., Dik, A. J., and Carver, T. A. W.,
  eds). APS Press, St. Paul, Minnesota, USA.} , 1-12.
\item
  Camacho, C., Coulouris, G., Avagyan, V., Ma, N., Papadopoulos, J.,
  Bealer, K., \& Madden, T. L. (2009). BLAST+: architecture and
  applications. \emph{BMC Bioinformatics, 10} , 421.
  doi:10.1186/1471-2105-10-421
\item
  Chandran, D., Inada, N., Hather, G., Kleindt, C. K., \& Wildermuth, M.
  C. (2010). Laser microdissection of Arabidopsis cells at the powdery
  mildew infection site reveals site-specific processes and
  regulators.\emph{Proceedings of the National Academy of Sciences of
  the United States of America, 107} (1), 460-465.
  doi:10.1073/pnas.0912492107
\item
  Chang, S., Puryear, J., \& Cairney, J. (1993). A simple and efficient
  method for isolating RNA from pine trees. \emph{Plant Molecular
  Biology Reporter, 11} (2), 113-116.
\item
  Chen, J., Pang, W., Chen, B., Zhang, C., \& Piao, Z. (2015).
  Transcriptome Analysis of Brassica rapa Near-Isogenic Lines Carrying
  Clubroot-Resistant and -Susceptible Alleles in Response to
  Plasmodiophora brassicae during Early Infection. \emph{Frontiers in
  Plant Science, 6} , 1183. doi:10.3389/fpls.2015.01183
\item
  Chisholm, S. T., Coaker, G., Day, B., \& Staskawicz, B. J. (2006).
  Host-microbe interactions: shaping the evolution of the plant immune
  response. \emph{Cell, 124} (4), 803-814.
  doi:10.1016/j.cell.2006.02.008
\item
  Choi, J., Choi, D., Lee, S., Ryu, C. M., \& Hwang, I. (2011).
  Cytokinins and plant immunity: old foes or new friends? \emph{Trends
  in Plant Science, 16} (7), 388-394. doi:10.1016/j.tplants.2011.03.003
\item
  Coll, N. S., Epple, P., \& Dangl, J. L. (2011). Programmed cell death
  in the plant immune system. \emph{Cell Death \& Differentiation, 18}
  (8), 1247-1256. doi:10.1038/cdd.2011.37
\item
  Cruz-Miralles, J., Cabedo-Lopez, M., Perez-Hedo, M., Flors, V., \&
  Jaques, J. A. (2019). Zoophytophagous mites can trigger plant-genotype
  specific defensive responses affecting potential prey beyond
  predation: the case of Euseius stipulatus and Tetranychus urticae in
  citrus. \emph{Pest Managment Science, 75} (7), 1962-1970.
  doi:10.1002/ps.5309
\item
  Davidson, N. M., Hawkins, A. D. K., \& Oshlack, A. (2017).
  SuperTranscripts: a data driven reference for analysis and
  visualisation of transcriptomes. \emph{Genome Biology, 18} (1), 148.
  doi:10.1186/s13059-017-1284-1
\item
  Decaestecker, E., Gaba, S., Raeymaekers, J. A., Stoks, R., Van
  Kerckhoven, L., Ebert, D., \& De Meester, L. (2007). Host-parasite
  'Red Queen' dynamics archived in pond sediment. \emph{Nature, 450}
  (7171), 870-873. doi:10.1038/nature06291
\item
  Dusa, A. (2018). venn: Draw Venn Diagrams. \emph{R package, version
  1.7.}
\item
  Egorov, T. A., \& Odintsova, T. I. (2012). Defense peptides of plant
  immunity. \emph{Russian Journal of Bioorganic Chemistry, 38} (1), 1-9.
  doi:10.1134/s1068162012010062
\item
  Fay, J. C., \& Wu, C. I. (2000). Hitchhiking under positive Darwinian
  selection. \emph{Genetics, 155} (3), 1405-1413.
\item
  Fei, Q., Xia, R., \& Meyers, B. C. (2013). Phased, Secondary, Small
  Interfering RNAs in Posttranscriptional Regulatory Networks. \emph{The
  Plant Cell, 25} (7), 2400-2415. doi:10.1105/tpc.113.114652
\item
  Fu, Y. B., Peterson, G. W., Horbach, C., Konkin, D. J., Beiles, A., \&
  Nevo, E. (2019). Elevated mutation and selection in wild emmer wheat
  in response to 28 years of global warming. \emph{Proceedings of the
  National Academy of Sciences (PNAS) USA, 116} (40), 20002-20008.
  doi:10.1073/pnas.1909564116
\item
  Gilbert, D. (2013). Gene-omes built from mRNA seq not genome
  DNA.\emph{7th annual arthropod genomics symposium. Notre Dame} .
\item
  Gilligan, C. A. (2002). An epidemiological framework for disease
  management. In \emph{Advances in Botanical Research} (Vol. 38, pp.
  1-64): Academic Press.
\item
  Gohre, V. (2015). Immune responses: Photosynthetic
  defence.\emph{Nature Plants, 1} , 15079. doi:10.1038/nplants.2015.79
\item
  Gomez, P., \& Buckling, A. (2011). Bacteria-Phage Antagonistic
  Coevolution in Soil. \emph{Science, 332} (6025), 106-109.
  doi:10.1126/science.1198767 \%J Science
\item
  Grabherr, M. G., Haas, B. J., Yassour, M., Levin, J. Z., Thompson, D.
  A., Amit, I., . . . Regev, A. (2011). Full-length transcriptome
  assembly from RNA-Seq data without a reference genome. \emph{Nature
  Biotechnology, 29} (7), 644-652. doi:10.1038/nbt.1883
\item
  Greischar, M. A., \& Koskella, B. (2007). A synthesis of experimental
  work on parasite local adaptation. \emph{Ecology Letters, 10} (5),
  418-434. doi:10.1111/j.1461-0248.2007.01028.x
\item
  Hahlbrock, K., Scheel, D., Logemann, E., Nurnberger, T., Parniske, M.,
  Reinold, S., . . . Schmelzer, E. (1995). Oligopeptide
  elicitor-mediated defense gene activation in cultured parsley
  cells.\emph{Proceedings of the National Academy of Sciences of the
  United States of America, 92} (10), 4150-4157.
  doi:10.1073/pnas.92.10.4150
\item
  Hall, S. A., Allen, R. L., Baumber, R. E., Baxter, L. A., Fisher, K.,
  Bittner-Eddy, P. D., . . . Beynon, J. L. (2009). Maintenance of
  genetic variation in plants and pathogens involves complex networks of
  gene-for-gene interactions. \emph{Molecular Plant Pathology, 10} (4),
  449-457. doi:10.1111/j.1364-3703.2009.00544.x
\item
  Hamilton, W. D. (1980). Sex versus Non-Sex versus
  Parasite.\emph{Oikos, 35} (2), 282-290. doi:10.2307/3544435
\item
  Hockerstedt, L. M., Siren, J. P., \& Laine, A. L. (2018). Effect of
  spatial connectivity on host resistance in a highly fragmented natural
  pathosystem. \emph{Journal of Evolutionary Biology, 31} (6), 844-852.
  doi:10.1111/jeb.13268
\item
  Hoeksema, J. D., \& Forde, S. E. (2008). A meta-analysis of factors
  affecting local adaptation between interacting species. \emph{The
  American Naturalist, 171} (3), 275-290. doi:10.1086/527496
\item
  Hua, L. (2009). Dissection of salicylic acid-mediated defense
  signaling networks. \emph{Plant Signaling and Behavior, 4} (8),
  713-717. doi:10.4161/psb.4.8.9173
\item
  Jacob, F., Vernaldi, S., \& Maekawa, T. (2013). Evolution and
  Conservation of Plant NLR Functions. \emph{Frontiers in Immunology, 4}
  , 297. doi:10.3389/fimmu.2013.00297
\item
  Jaenike, J. (1978). Host Selection by Mycophagous
  Drosophila.\emph{ECOLOGY, 59} (6), 1286-1288. doi:10.2307/1938245
\item
  Jensen, M. K., Hagedorn, P. H., de Torres-Zabala, M., Grant, M. R.,
  Rung, J. H., Collinge, D. B., \& Lyngkjaer, M. F. (2008).
  Transcriptional regulation by an NAC (NAM-ATAF1,2-CUC2) transcription
  factor attenuates ABA signalling for efficient basal defence towards
  Blumeria graminis f. sp. hordei in Arabidopsis. \emph{The Plant
  Journal, 56} (6), 867-880. doi:10.1111/j.1365-313X.2008.03646.x
\item
  Jones, J. D., \& Dangl, J. L. (2006). The plant immune
  system.\emph{Nature, 444} (7117), 323-329. doi:10.1038/nature05286
\item
  Jones, P., Binns, D., Chang, H. Y., Fraser, M., Li, W., McAnulla, C.,
  . . . Hunter, S. (2014). InterProScan 5: genome-scale protein function
  classification. \emph{Bioinformatics, 30} (9), 1236-1240.
  doi:10.1093/bioinformatics/btu031
\item
  Jousimo, Tack, A. J., Ovaskainen, O., Mononen, T., Susi, H.,
  Tollenaere, C., \& Laine, A. L. (2014a). Ecological and evolutionary
  effects of fragmentation on infectious disease dynamics.\emph{Science,
  344} (6189), 1289-1293. doi:10.1126/science.1253621
\item
  Jousimo, J., Tack, A. J. M., Ovaskainen, O., Mononen, T., Susi, H.,
  Tollenaere, C., \& Laine, A.-L. (2014b). Ecological and evolutionary
  effects of fragmentation on infectious disease dynamics.\emph{344}
  (6189), 1289-1293. doi:10.1126/science.1253621 \%J Science
\item
  Kachroo, A., \& Kachroo, P. (2009). Fatty Acid-derived signals in
  plant defense. \emph{Annual Review of Phytopathology, 47} , 153-176.
  doi:10.1146/annurev-phyto-080508-081820
\item
  Karim, S., Holmstrom, K. O., Mandal, A., Dahl, P., Hohmann, S.,
  Brader, G., . . . Pirhonen, M. (2007). AtPTR3, a wound-induced peptide
  transporter needed for defence against virulent bacterial pathogens in
  Arabidopsis. \emph{Planta, 225} (6), 1431-1445.
  doi:10.1007/s00425-006-0451-5
\item
  Klopfenstein, D. V., Zhang, L., Pedersen, B. S., Ramirez, F., Warwick
  Vesztrocy, A., Naldi, A., . . . Tang, H. (2018). GOATOOLS: A Python
  library for Gene Ontology analyses. \emph{Scientific Reports, 8} (1),
  10872. doi:10.1038/s41598-018-28948-z
\item
  Koff, R. S. (1992). Infectious diseases of humans: Dynamics and
  control. By R.M. Anderson and R.M. May, 757 pp. Oxford: Oxford
  University Press, 1991. \$95.00. \emph{Hepatology, 15} (1), 169-169.
  doi:10.1002/hep.1840150131
\item
  Korneliussen, T. S., Albrechtsen, A., \& Nielsen, R. (2014). ANGSD:
  Analysis of Next Generation Sequencing Data. \emph{BMC Bioinformatics,
  15} (1), 356. doi:10.1186/s12859-014-0356-4
\item
  Kryazhimskiy, S., \& Plotkin, J. B. (2008). The population genetics of
  dN/dS. \emph{PLoS Genetics, 4} (12), e1000304.
  doi:10.1371/journal.pgen.1000304
\item
  Kuniak, E., \& Urbanek, H. (2000). The involvement of hydrogen
  peroxide in plant responses to stresses. \emph{Acta Physiologiae
  Plantarum, 22} (2), 95-203.
  doi:https://doi.org/10.1007/s11738-000-0076-4
\item
  Kursar, T. A., Dexter, K. G., Lokvam, J., Pennington, R. T.,
  Richardson, J. E., Weber, M. G., . . . Coley, P. D. (2009). The
  evolution of antiherbivore defenses and their contribution to species
  coexistence in the tropical tree genus
  \&lt;em\&gt;Inga\&lt;/em\&gt.\emph{Proceedings of the National Academy
  of Sciences, 106} (43), 18073. doi:10.1073/pnas.0904786106
\item
  Lai, Y., \& Eulgem, T. (2018). Transcript-level expression control of
  plant NLR genes. \emph{Molecular Plant Pathology, 19} (5), 1267-1281.
  doi:10.1111/mpp.12607
\item
  Laine, A. L. (2004). Resistance variation within and among host
  populations in a plant--pathogen metapopulation: implications for
  regional pathogen dynamics. \emph{Journal of Ecology, 92} (6),
  990-1000. doi:doi:10.1111/j.0022-0477.2004.00925.x
\item
  Laine, A. L. (2005). Spatial scale of local adaptation in a
  plant-pathogen metapopulation. \emph{Journal of Evolutionary Biology,
  18} (4), 930-938. doi:10.1111/j.1420-9101.2005.00933.x
\item
  Laine, A. L. (2006). Evolution of host resistance: looking for
  coevolutionary hotspots at small spatial scales. \emph{Proceedings of
  the Royal Society B: Biological Sciences, 273} (1584), 267-273.
  doi:10.1098/rspb.2005.3303
\item
  Laine, A. L. (2007). Pathogen fitness components and genotypes differ
  in their sensitivity to nutrient and temperature variation in a wild
  plant-pathogen association. \emph{Journal of Evolutionary Biology, 20}
  (6), 2371-2378. doi:10.1111/j.1420-9101.2007.01406.x
\item
  Laine, A. L. (2008). Temperature-mediated patterns of local adaptation
  in a natural plant-pathogen metapopulation. \emph{Ecology Letters, 11}
  (4), 327-337. doi:10.1111/j.1461-0248.2007.01146.x
\item
  Laine, A. L., Burdon, J. J., Dodds, P. N., \& Thrall, P. H. (2011).
  Spatial variation in disease resistance: from molecules to
  metapopulations. \emph{Journal of Ecology, 99} (1), 96-112.
  doi:10.1111/j.1365-2745.2010.01738.x
\item
  Le Roy, J., Huss, B., Creach, A., Hawkins, S., \& Neutelings, G.
  (2016). Glycosylation Is a Major Regulator of Phenylpropanoid
  Availability and Biological Activity in Plants. \emph{Frontiers in
  Plant Science, 7} , 735. doi:10.3389/fpls.2016.00735
\item
  Leonard, K. J. (1977). Selection pressure and plant
  pathogens.\emph{Annals of the New York Academy of Sciences, 287} (1),
  207-222. doi:10.1111/j.1749-6632.1977.tb34240.x
\item
  Li, B., Ruotti, V., Stewart, R. M., Thomson, J. A., \& Dewey, C. N.
  (2010). RNA-Seq gene expression estimation with read mapping
  uncertainty. \emph{Bioinformatics, 26} (4), 493-500.
  doi:10.1093/bioinformatics/btp692
\item
  Li, G., Xu, X., Bai, G., Carver, B. F., Hunger, R., \& Bonman, J. M.
  (2016). Identification of Novel Powdery Mildew Resistance Sources in
  Wheat. \emph{Acess DL, 56} (4), 1817-1830.
  doi:10.2135/cropsci2015.09.0551
\item
  Li, H., Dong, Z., Ma, C., Tian, X., Xiang, Z., Xia, Q., . . . Liu, W.
  (2019). Discovery of powdery mildew resistance gene candidates from
  Aegilops biuncialis chromosome 2Mb based on transcriptome
  sequencing.\emph{PLoS One, 14} (11), e0220089.
  doi:10.1371/journal.pone.0220089
\item
  Li, M., Zhang, D., Gao, Q., Luo, Y., Zhang, H., Ma, B., . . . Xue, Y.
  (2019). Genome structure and evolution of Antirrhinum majus
  L.\emph{Nature Plants, 5} (2), 174-183. doi:10.1038/s41477-018-0349-9
\item
  Liang, P., Liu, S., Xu, F., Jiang, S., Yan, J., He, Q., . . . Miao, W.
  (2018). Powdery Mildews Are Characterized by Contracted Carbohydrate
  Metabolism and Diverse Effectors to Adapt to Obligate Biotrophic
  Lifestyle. \emph{Frontiers in Microbiology, 9} , 3160.
  doi:10.3389/fmicb.2018.03160
\item
  Love, M. I., Huber, W., \& Anders, S. (2014). Moderated estimation of
  fold change and dispersion for RNA-seq data with DESeq2. \emph{Genome
  Biology, 15} (12), 550. doi:10.1186/s13059-014-0550-8
\item
  McDowell, J. M., \& Simon, S. A. (2006). Recent insights into R gene
  evolution. \emph{Molecular Plant Pathology, 7} (5), 437-448.
  doi:10.1111/j.1364-3703.2006.00342.x
\item
  Meunier, E., \& Broz, P. (2017). Evolutionary Convergence and
  Divergence in NLR Function and Structure. \emph{Trends in Immunology,
  38} (10), 744-757. doi:10.1016/j.it.2017.04.005
\item
  Meyers, B. C. (2003). Genome-Wide Analysis of NBS-LRR-Encoding Genes
  in Arabidopsis. \emph{The Plant Cell Online, 15} (4), 809-834.
  doi:10.1105/tpc.009308
\item
  Michelmore, R. W., \& Meyers, B. C. (1998). Clusters of Resistance
  Genes in Plants Evolve by Divergent Selection and a Birth-and-Death
  Process. \emph{GENOME RESEARCH} (8), 1113-1130.
  doi:10.1101/gr.8.11.1113
\item
  Miedes, E., Vanholme, R., Boerjan, W., \& Molina, A. (2014). The role
  of the secondary cell wall in plant resistance to
  pathogens.\emph{Frontiers in Plant Science, 5} , 358.
  doi:10.3389/fpls.2014.00358
\item
  Mitchell, A. L., Attwood, T. K., Babbitt, P. C., Blum, M., Bork, P.,
  Bridge, A., . . . Finn, R. D. (2019). InterPro in 2019: improving
  coverage, classification and access to protein sequence
  annotations.\emph{Nucleic Acids Research, 47} (D1), D351-D360.
  doi:10.1093/nar/gky1100
\item
  Mock, H. P., Keetman, U., Kruse, E., Rank, B., \& Grimm, B. (1998).
  Defense Responses to Tetrapyrrole-Induced Oxidative Stress in
  Transgenic Plants with Reduced Uroporphyrinogen Decarboxylase or
  Coproporphyrinogen Oxidase Activity. \emph{Plant Physiology, 116} (1),
  107. doi:10.1104/pp.116.1.107
\item
  Mohr, T. J., Mammarella, N. D., Hoff, T., Woffenden, B. J., Jelesko,
  J. G., \& McDowell, J. M. (2010). The Arabidopsis Downy Mildew
  Resistance Gene RPP8 Is Induced by Pathogens and Salicylic Acid and Is
  Regulated by W Box cis Elements. \emph{Molecular Plant-Microbe
  Interactions(r), 23} (10), 1303-1315. doi:10.1094/MPMI-01-10-0022
\item
  Monteiro, F., \& Nishimura, M. T. (2018). Structural, Functional, and
  Genomic Diversity of Plant NLR Proteins: An Evolved Resource for
  Rational Engineering of Plant Immunity. \emph{Annual Review of
  Phytopathology, 56} , 243-267. doi:10.1146/annurev-phyto-080417-045817
\item
  Muller, N. A., Kersten, B., Fladung, M., \& Schroeder, H. (2019).
  RNA-seq of eight different poplar clones reveals conserved
  up-regulation of gene expression in response to insect
  herbivory.\emph{BMC Genomics, 20} (1), 673.
  doi:10.1186/s12864-019-6048-8
\item
  Mundt, C. C. (2002a). Use of multiline cultivars and cultivar mixtures
  fo disease management. \emph{Annual Review of Phytopathology, 40} (1),
  381-410. doi:10.1146/annurev.phyto.40.011402.113723
\item
  Mundt, C. C. (2002b). Use of multiline cultivars and cultivar mixtures
  for disease management. \emph{40} (1), 381-410.
  doi:10.1146/annurev.phyto.40.011402.113723
\item
  Mundt, C. C. (2014). Durable resistance: a key to sustainable
  management of pathogens and pests. \emph{Infection, Genetics and
  Evolution, 27} , 446-455. doi:10.1016/j.meegid.2014.01.011
\item
  Mur, L. A. J., Simpson, C., Kumari, A., Gupta, A. K., \& Gupta, K. J.
  (2017). Moving nitrogen to the centre of plant defence against
  pathogens. \emph{Annals of Botany, 119} (5), 703-709.
  doi:10.1093/aob/mcw179
\item
  Niu, L., \& Liao, W. (2016). Hydrogen Peroxide Signaling in Plant
  Development and Abiotic Responses: Crosstalk with Nitric Oxide and
  Calcium. \emph{Frontiers in Plant Science, 7} , 230.
  doi:10.3389/fpls.2016.00230
\item
  Nurnberger, T., Nennstiel, D., Jabs, T., Sacks, W. R., Hahlbrock, K.,
  \& Scheel, D. (1994). High affinity binding of a fungal oligopeptide
  elicitor to parsley plasma membranes triggers multiple defense
  responses. \emph{Cell, 78} (3), 449-460.
  doi:https://doi.org/10.1016/0092-8674(94)90423-5
\item
  Oksanen, J., Blanchet, F. G., Friendly, M., Kindt, R., Legendre, P.,
  McGlinn, D., . . . Wagner, H. (2018). vegan: Community Ecology
  Package. \emph{R package, version 2.5-3.}
\item
  Perez-Bueno, M. L., Pineda, M., \& Baron, M. (2019). Phenotyping Plant
  Responses to Biotic Stress by Chlorophyll Fluorescence
  Imaging.\emph{Frontiers in Plant Science, 10} , 1135.
  doi:10.3389/fpls.2019.01135
\item
  Phadnis, N., \& Malik, H. S. (2014). Speciation via autoimmunity: a
  dangerous mix. \emph{Cell, 159} (6), 1247-1249.
  doi:10.1016/j.cell.2014.11.028
\item
  Pike, S., Gao, F., Kim, M. J., Kim, S. H., Schachtman, D. P., \&
  Gassmann, W. (2014). Members of the NPF3 transporter subfamily encode
  pathogen-inducible nitrate/nitrite transporters in grapevine and
  Arabidopsis. \emph{Plant and Cell Physiology, 55} (1), 162-170.
  doi:10.1093/pcp/pct167
\item
  Plomion, C., Aury, J. M., Amselem, J., Leroy, T., Murat, F.,
  Duplessis, S., . . . Salse, J. (2018). Oak genome reveals facets of
  long lifespan. \emph{Nature Plants, 4} (7), 440-452.
  doi:10.1038/s41477-018-0172-3
\item
  Polonio, A., Pineda, M., Bautista, R., Martinez-Cruz, J., Perez-Bueno,
  M. L., Baron, M., \& Perez-Garcia, A. (2019). RNA-seq analysis and
  fluorescence imaging of melon powdery mildew disease reveal an
  orchestrated reprogramming of host physiology. \emph{Scientific
  Reports, 9} (1). doi:10.1038/s41598-019-44443-5
\item
  Price, M. N., Dehal, P. S., \& Arkin, A. P. (2010). FastTree
  2--approximately maximum-likelihood trees for large
  alignments.\emph{PLoS One, 5} (3), e9490.
  doi:10.1371/journal.pone.0009490
\item
  Pruitt, K. D., Tatusova, T., \& Maglott, D. R. (2007). NCBI reference
  sequences (RefSeq): a curated non-redundant sequence database of
  genomes, transcripts and proteins. \emph{Nucleic Acids Research, 35}
  (Database issue), D61-65. doi:10.1093/nar/gkl842
\item
  Ragonnet-Cronin, M., Hodcroft, E., Hue, S., Fearnhill, E., Delpech,
  V., Leigh Brown, A. J., \& Lycett, S. (2013). Automated analysis of
  phylogenetic clusters. \emph{BMC Bioinformatics, 14} (317).
  doi:10.1186/1471-2105-14-317
\item
  Rentel, M. C., Leonelli, L., Dahlbeck, D., Zhao, B., \& Staskawicz, B.
  J. (2008). Recognition of the Hyaloperonospora parasitica effector
  ATR13 triggers resistance against oomycete, bacterial, and viral
  pathogens. \emph{Proceedings of the National Academy of Sciences of
  the United States of America, 105} (3), 1091-1096.
  doi:10.1073/pnas.0711215105
\item
  Revers, F., Guiraud, T., Houvenaghel, M. C., Mauduit, T., Le Gall, O.,
  \& Candresse, T. (2003). Multiple Resistance Phenotypes to Lettuce
  mosaic virus Among Arabidopsis thaliana Accessions. \emph{Molecular
  Plant-Microbe Interactions, 16} (7), 608-616.
  doi:10.1094/MPMI.2003.16.7.608
\item
  Ribot, C., Hirsch, J., Balzergue, S., Tharreau, D., Notteghem, J. L.,
  Lebrun, M. H., \& Morel, J. B. (2008). Susceptibility of rice to the
  blast fungus, Magnaporthe grisea. \emph{Journal of Plant Physiology,
  165} (1), 114-124. doi:10.1016/j.jplph.2007.06.013
\item
  Rose, L. E., Bittner-Eddy, P. D., Langley, C. H., Holub, E. B.,
  Michelmore, R. W., \& Beynon, J. L. (2004). The Maintenance of Extreme
  Amino Acid Diversity at the Disease Resistance Gene,
  \textless{}em\textgreater{}RPP13\textless{}/em\textgreater{}, in
  \textless{}em\textgreater{}Arabidopsis
  thaliana\textless{}/em\textgreater{}. \emph{Genetics, 166} (3),
  1517-1527. doi:10.1534/genetics.166.3.1517 \%J Genetics
\item
  Roux, F., \& Bergelson, J. (2016). Chapter Four - The Genetics
  Underlying Natural Variation in the Biotic Interactions of Arabidopsis
  thaliana: The Challenges of Linking Evolutionary Genetics and
  Community Ecology. In V. Orgogozo (Ed.), \emph{Current Topics in
  Developmental Biology} (Vol. 119, pp. 111-156): Academic Press.
\item
  Rozen, S., \& Skaletsky, H. (1999). Primer3 on the WWW for General
  Users and for Biologist Programmers. In S. Misener \& S. A. Krawetz
  (Eds.), \emph{Bioinformatics Methods and Protocols} (pp. 365-386).
  Totowa, NJ: Humana Press.
\item
  Sagar, G. R., \& Harper, J. L. (1964). Plantago Major L., P. Media L.
  and P. Lanceolata L. \emph{Journal of Ecology, 52} (1), 189-221.
  doi:10.2307/2257792
\item
  Salvaudon, L., Giraud, T., \& Shykoff, J. A. (2008). Genetic diversity
  in natural populations: a fundamental component of plant-microbe
  interactions. \emph{Current Opinion in Plant Biology, 11} (2),
  135-143. doi:10.1016/j.pbi.2008.02.002
\item
  Sari, E., Bhadauria, V., Ramsay, L., Borhan, M. H., Lichtenzveig, J.,
  Bett, K. E., . . . Banniza, S. (2018). Defense responses of lentil
  (Lens culinaris) genotypes carrying non-allelic ascochyta blight
  resistance genes to Ascochyta lentis infection. \emph{PLoS One, 13}
  (9), e0204124. doi:10.1371/journal.pone.0204124
\item
  Sari, E., Bhadauria, V., Vandenberg, A., \& Banniza, S. (2017).
  Genotype-Dependent Interaction of Lentil Lines with Ascochyta
  lentis.\emph{Frontiers in Plant Science, 8} , 764.
  doi:10.3389/fpls.2017.00764
\item
  Sarris, P. F., Cevik, V., Dagdas, G., Jones, J. D., \& Krasileva, K.
  V. (2016). Comparative analysis of plant immune receptor architectures
  uncovers host proteins likely targeted by pathogens. \emph{BMC
  Biology, 14} , 8. doi:10.1186/s12915-016-0228-7
\item
  Schulz, M. H., Zerbino, D. R., Vingron, M., \& Birney, E. (2012).
  Oases: robust de novo RNA-seq assembly across the dynamic range of
  expression levels. \emph{Bioinformatics, 28} (8), 1086-1092.
  doi:10.1093/bioinformatics/bts094
\item
  Seppey, M., Manni, M., \& Zdobnov, E. M. (2019). BUSCO: Assessing
  Genome Assembly and Annotation Completeness. In M. Kollmar
  (Ed.),\emph{Gene Prediction: Methods and Protocols} (pp. 227-245). New
  York, NY: Springer New York.
\item
  Shang, X., Cao, Y., \& Ma, L. (2017). Alternative Splicing in Plant
  Genes: A Means of Regulating the Environmental Fitness of
  Plants.\emph{International Journal of Molecular Sciences, 18} (2).
  doi:10.3390/ijms18020432
\item
  Shao, Z. Q., Xue, J. Y., Wu, P., Zhang, Y. M., Wu, Y., Hang, Y. Y., .
  . . Chen, J. Q. (2016). Large-Scale Analyses of Angiosperm
  Nucleotide-Binding Site-Leucine-Rich Repeat Genes Reveal Three
  Anciently Diverged Classes with Distinct Evolutionary
  Patterns.\emph{Plant Physiology, 170} (4), 2095-2109.
  doi:10.1104/pp.15.01487
\item
  Sheshukova, E. V., Komarova, T. V., Pozdyshev, D. V., Ershova, N. M.,
  Shindyapina, A. V., Tashlitsky, V. N., . . . Dorokhov, Y. L. (2017).
  The Intergenic Interplay between Aldose 1-Epimerase-Like Protein and
  Pectin Methylesterase in Abiotic and Biotic Stress
  Control.\emph{Frontiers in Plant Science, 8} .
  doi:10.3389/fpls.2017.01646
\item
  Soneson, C., Love, M. I., \& Robinson, M. D. (2015). Differential
  analyses for RNA-seq: transcript-level estimates improve gene-level
  inferences. \emph{F1000Research, 4} , 1521.
  doi:10.12688/f1000research.7563.2
\item
  Srivastava, A., Sarkar, H., Malik, L., \& Patro, R. (2016). Accurate,
  Fast and Lightweight Clustering of de novo Transcriptomes using
  Fragment Equivalence Classes. \emph{eprint arXiv:1604.03250} .
\item
  Staiger, D., Korneli, C., Lummer, M., \& Navarro, L. (2013). Emerging
  role for RNA-based regulation in plant immunity. \emph{New
  Phytologist, 197} (2), 394-404. doi:10.1111/nph.12022
\item
  Stam, R., Silva-Arias, G. A., Nosenko, T., Scheikl, D., Horger, A. C.,
  Stephan, W., . . . Tellier, A. (2017). A small subset of NLR genes
  drives local adaptation to pathogens in wild tomato. \emph{bioRxiv} .
  doi:10.1101/210559
\item
  Steuernagel, B., Jupe, F., Witek, K., Jones, J. D., \& Wulff, B. B.
  (2015). NLR-parser: rapid annotation of plant NLR
  complements.\emph{Bioinformatics, 31} (10), 1665-1667.
  doi:10.1093/bioinformatics/btv005
\item
  Thompson, J. N., \& Burdon, J. J. (1992). Gene-for-gene coevolution
  between plants and parasites. \emph{Nature, 360} (6400), 121-125.
  doi:10.1038/360121a0
\item
  Thrall, P. H., Laine, A. L., Ravensdale, M., Nemri, A., Dodds, P. N.,
  Barrett, L. G., \& Burdon, J. J. (2012). Rapid genetic change
  underpins antagonistic coevolution in a natural host-pathogen
  metapopulation. \emph{Ecology Letters, 15} (5), 425-435.
  doi:10.1111/j.1461-0248.2012.01749.x
\item
  Tucker, A. E., Ackerman, M. S., Eads, B. D., Xu, S., \& Lynch, M.
  (2013). Population-genomic insights into the evolutionary origin and
  fate of obligately asexual Daphnia pulex. \emph{Proceedings of the
  National Academy of Sciences, 110} (39), 15740-15745.
  doi:10.1073/pnas.1313388110
\item
  Upson, J. L., Zess, E. K., Bialas, A., Wu, C. H., \& Kamoun, S.
  (2018). The coming of age of EvoMPMI: evolutionary molecular
  plant-microbe interactions across multiple timescales. \emph{Current
  Opinion in Plant Biology, 44} , 108-116. doi:10.1016/j.pbi.2018.03.003
\item
  Wallstrom, S. V., Florez-Sarasa, I., Araujo, W. L., Aidemark, M.,
  Fernandez-Fernandez, M., Fernie, A. R., . . . Rasmusson, A. G. (2014).
  Suppression of the external mitochondrial NADPH dehydrogenase, NDB1,
  in Arabidopsis thaliana affects central metabolism and vegetative
  growth. \emph{Molecular Plant, 7} (2), 356-368. doi:10.1093/mp/sst115
\item
  Walters, D. R., \& McRoberts, N. (2006). Plants and biotrophs: a
  pivotal role for cytokinins? \emph{Trends Plant Science, 11} (12),
  581-586. doi:10.1016/j.tplants.2006.10.003
\item
  Van de Weyer, A. L., Monteiro, F., Furzer, O. J., Nishimura, M. T.,
  Cevik, V., Witek, K., . . . Bemm, F. (2019). A Species-Wide Inventory
  of NLR Genes and Alleles in Arabidopsis thaliana. \emph{Cell, 178}
  (5), 1260-1272.e1214. doi:10.1016/j.cell.2019.07.038
\item
  Van Ooijen, G., Mayr, G., Kasiem, M. M. A., Albrecht, M., Cornelissen,
  B. J. C., \& Takken, F. L. W. (2008). Structure--function analysis of
  the NB-ARC domain of plant disease resistance proteins. \emph{Journal
  of Experimental Botany, 59} (6), 1383-1397. doi:10.1093/jxb/ern045
\item
  Varemo, L., Nielsen, J., \& Nookaew, I. (2013). Enriching the gene set
  analysis of genome-wide data by incorporating directionality of gene
  expression and combining statistical hypotheses and
  methods.\emph{Nucleic Acids Research, 41} (8), 4378-4391.
  doi:10.1093/nar/gkt111
\item
  Wickham H. (2016). ggplot2: Elegant Graphics for Data
  Analysis.\emph{Springer-Verlag New York.}
\item
  Wiese, J., Kranz, T., \& Schubert, S. (2004). Induction of pathogen
  resistance in barley by abiotic stress. \emph{Plant Biology, 6} (5),
  529-536. doi:10.1055/s-2004-821176
\item
  Woolhouse, M. E. J., Taylor, L. H., \& Haydon, D. T. (2001).
  Population Biology of Multihost Pathogens. \emph{Science, 292} (5519),
  1109-1112. doi:10.1126/science.1059026 \%J Science
\item
  Xie, Y., Wu, G., Tang, J., Luo, R., Patterson, J., Liu, S., . . .
  Wang, J. (2014). SOAPdenovo-Trans: de novo transcriptome assembly with
  short RNA-Seq reads. \emph{Bioinformatics, 30} (12), 1660-1666.
  doi:10.1093/bioinformatics/btu077
\item
  Yang, J., Duan, G., Li, C., Liu, L., Han, G., Zhang, Y., \& Wang, C.
  (2019). The Crosstalks Between Jasmonic Acid and Other Plant Hormone
  Signaling Highlight the Involvement of Jasmonic Acid as a Core
  Component in Plant Response to Biotic and Abiotic
  Stresses.\emph{Frontiers in Plant Science, 10} .
  doi:10.3389/fpls.2019.01349
\item
  Yang, Z. (2007). PAML 4: phylogenetic analysis by maximum
  likelihood.\emph{Molecular Biology and Evolution, 24} (8), 1586-1591.
  doi:10.1093/molbev/msm088
\item
  Yu, G., Smith, D. K., Zhu, H., Guan, Y., Lam, T. T., \& McInerny, G.
  (2017). ggtree: anrpackage for visualization and annotation of
  phylogenetic trees with their covariates and other associated
  data.\emph{Methods in Ecology and Evolution, 8} (1), 28-36.
  doi:10.1111/2041-210x.12628
\item
  Zhou, J., Sun, A., \& Xing, D. (2013). Modulation of cellular redox
  status by thiamine-activated NADPH oxidase confers Arabidopsis
  resistance to Sclerotinia sclerotiorum. \emph{Journal of Experimental
  Botany, 64} (11), 3261-3272. doi:10.1093/jxb/ert166
\item
  Zhu, Y., Chen, H., Fan, J., Wang, Y., Li, Y., Chen, J., . . . Mundt,
  C. C. (2000). Genetic diversity and disease control in
  rice.\emph{Nature, 406} (6797), 718-722. doi:10.1038/35021046
\end{enumerate}

\section*{Data Accessibility}

{\label{data-accessibility}}

Upon acceptance, all data associated with this study will be deposited
into NCBI (SUB6931705).

\section*{Author Contribution}

{\label{author-contribution}}

A-LL, JS, PS and LH conceived the ideas and JS, PS, LH, MB and A-LL
designed experiment. LH conducted the experimental work and PS, JS and
LH analyzed the data. PS, JS and A-LL led the writing of the manuscript;
all the authors contributed to the drafts and gave final approval for
publication.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Safdari-et-al-Tables---Figures/Safdari-et-al-Tables---Figures}

\end{center}
\end{figure}

\selectlanguage{english}
\FloatBarrier
\end{document}

