Multiple quantitative trait loci contribute to resistance to bacterial canker incited by Pseudomonas syringae pv. actinidiae in kiwifruit (Actinidia chinensis)

Pseudomonas syringae pv. actinidiae (Psa) biovar 3, a virulent, canker-inducing pathogen is an economic threat to the kiwifruit (Actinidia spp.) industry worldwide. The commercially grown diploid (2×) A. chinensis var. chinensis is more susceptible to Psa than tetraploid and hexaploid kiwifruit. However information on the genetic loci modulating Psa resistance in kiwifruit is not available. Here we report mapping of quantitative trait loci (QTLs) regulating resistance to Psa in a diploid kiwifruit population, derived from a cross between an elite Psa-susceptible ‘Hort16A’ and a resistant male breeding parent P1. Using high-density genetic maps and intensive phenotyping, we identified a single QTL for Psa resistance on Linkage Group (LG) 27 of ‘Hort16A’ revealing 16–19% phenotypic variance and candidate alleles for susceptibility and resistance at this loci. In addition, six minor QTLs were identified in P1 on distinct LGs, exerting 4–9% variance. Resistance in the F1 population is improved by additive effects from ‘Hort16A’ and P1 QTLs providing evidence that divergent genetic pathways interact to combat the virulent Psa strain. Two different bioassays further identified new QTLs for tissue-specific responses to Psa. The genetic marker at LG27 QTL was further verified for association with Psa resistance in diploid Actinidia chinensis populations. Transcriptome analysis of Psa-resistant and susceptible genotypes in field revealed hallmarks of basal defense and provided candidate RNA-biomarkers for screening for Psa resistance in greenhouse conditions.


Introduction
Pseudomonas syringae is a hemi-biotrophic bacterial complex 1 that can infect a range of plant species. It comprises pathovars which cause similar symptoms on their host plants and several pathovars can lead to severe crop loss. P. syringae pv. actinidiae (Psa) infects several species of Actinidia (kiwifruit) 2,3 and virulent Psa strains induce a range of symptoms on the main stem of the vine, foliage, floral buds and fruits 4 . Psa pathovar strains can be grouped into five biovars based on their genetic and biological characteristics 4,5 . Strains of biovar 3, previously known as Psa-V (referred to here as Psa), are currently the most aggressive and were responsible for outbreaks from 2008 [6][7][8][9] . Psa has cost the kiwifruit industry billions of dollars worldwide and its incursion in New Zealand in 2010 completely destroyed vines of the Psa-susceptible diploid A. chinensis 'Hort16A' 4,8,10 .
Screening of thousands of Actinidia genotypes from 24 taxa in the breeding program at The New Zealand Institute for Plant & Food Research Limited (PFR) for resistance to natural and artificial Psa infections 25,26 revealed that diploid (2×) A. chinensis var. chinensis are more susceptible to Psa infection than tetraploid (4×) A. chinensis var. chinensis, which in turn are more susceptible than diploid and hexaploid (6×) A. chinensis var. deliciosa [25][26][27] . Many species outside the A. chinensis complex are more resistant to Psa than A. chinensis and the germplasm holds diverse genetic potential for Psa resistance 28 . Information on the genetic markers and molecular mechanisms associated with Psa resistance and resistance in the commercial cultivars producing taxas including A. chinensis and A. arguta is however limited. In this study we provide the first detailed view of the genetic loci modulating Psa resistance and tissue-specific response in diploid A. chinensis, utilizing an intensively phenotyped population of seedlings developed from a cross between Psa-susceptible 'Hort16A' and a resistant breeding parent (P1), as our experimental material for quantitative trait locus (QTL) analysis.

Intensive phenotyping targets diverse developmental stages and environmental conditions
Initially, a pilot population of 53 genotypes from 'Hort16A' × P1 were replicated 3 times and phenotyped following natural field infection with Psa. The purpose of this population was to record and differentiate the types of symptoms elicited in response to natural Psa infection in the field and their segregation with genotype and was reported previously 28 . The response to Psa infection in the expanded population was measured on 236 genotypes of the 'Hort16A' × P1 population, which were clonally replicated~30 times. Phenotyping of the population was performed under field conditions following natural infection, as well as using two bioassays (scheme for phenotyping is laid out in Supplementary Fig. 1a). Multiple phenotypes were recorded in field (Fig. 1a-e) to develop a combined score referred to as Psa_score_Field (Fig. 2a). The mean clonal repeatability for this score was 0.65, while the repeatability of clonal means at 0.8. For the stab assay 26 , various tissue-specific phenotypic responses were recorded, including Stem_necrosis, Leaf_spots, Ooze, Stem_collapse, Tip_death and Wilt ( Fig. 1f-k, Fig.  2b), with repeatability of clonal means for these scores as 0.60, 0.766, 0.64, 0.79, 0.78, and 0.71, respectively. A Psa_score_Stab was also calculated (Fig. 2b) from all phenotypes assessed in the Stab assay (see Experimental Procedures). In the flood bioassay, adapted from previous work 29 , overall health was scored at weekly intervals postinoculation (Flood Assay/FA_Week1 to FA_Week5) (Fig.  1l). The frequency distribution of phenotypes and Psa_scores revealed that most exhibited non-Normal distribution based on the Shapiro-Wilk test (Fig. 2). However, for the stab assay, the majority of the observations displayed normal distributions. (Fig. 2b). As such, the correlation among the phenotypic scores from the field assessments and the bioassays was found to be poor ( Supplementary Fig. 1b). The correlation among different phenotypes within the bioassays was positive and medium to high (between 0.5 and 0.9) as defined in the correlation matrix ( Supplementary Fig. 1b). A 3-dimensional principal components analysis (PCA) on the correlation matrix of the field assessment, stab assay and flood assay displays a high degree of divergence in the rankings of the population for Psa response and resistance when assessed through different approaches (Fig. 2d).
QTL mapping from field phenotype scores confirmed polygenic nature of Psa field resistance A QTL for control of field resistance to Psa, Psa_s-core_Field, was identified in 'Hort16A' on the upper arm of LG27 (Fig. 3a) using multiple models for QTL discovery. At a LOD score of 7.02 (Fig. 3a), the location of the LG27 QTL on the Red5 genome (version 1.69.0) 33 is between~3.4 and 4.6 Mbp. The LG27 QTL was also identified for Psa_score-Field, in 'Hort16A' from the pilot trial (Supplementary Table 1). A SNP marker G9P1 developed from Acc30822, a gene of unknown function underlying the QTL and a multi-allelic Simple Sequence Repeat (SSR) marker SSRLG27_439F4R4, contributed 16% (favorable allele b) and 19% (favorable allele v, band size 428 bp) of the population phenotypic variance, respectively ( Fig. 3g and Supplementary Table 1). The multi-allelic SSR marker revealed the contribution of the favorable 428 bp A. chinensis grandparental allele v, to Psa resistance (Fig. 3g), compared to the other 408 bp allele u which is associated with susceptibility.
From these, the effect of favorable grandparent alleles from P1 on field resistance was verified for 3 QTLs i.e., by analysis of an SSR marker designed in the region underlying the LG22 QTL (SSRLG22_8032664) (Fig. 3g Fig. 5). Screening of SSRLG27_439F4R4 in another set of field-grown 'Hort16A' × P1 progeny confirmed association of the 428 bp allele from 'Hort16A' with  Table 1): this combination identified~80% of the resistant 'Hort16A' × P1 genotypes in the field ( Supplementary Fig. 7).
Validation of the LG27 QTL marker and candidate alleles for susceptibility and resistance Since the QTL identified in 'Hort16A' has the greatest effect, we predicted that this locus might be linked to susceptibility observed in diploid kiwifruit breeding parents. For this purpose, we performed validation of the LG27 QTL in four different field grown A. chinensis populations with a genetic background related to the parents of the mapping population, using the SSR marker SSRLG27_439F4R4. The G9P1 and SSRLG27_439F4R4 markers are tightly linked on the genetic map and located within~300 kb on the physical map ( Supplementary  Fig. 5b).
The first validation population is a cross between a Psa resistant female V1 (which is a sister of P1) and Psa resistant male P2, the son of 'Hort16A'. The 428 bp allele of SSRLG27_4396125F4R4 linked to Psa resistance in P2 contributed 10.3% phenotypic variance ( Supplementary  Fig. 6b). Similarly, the same allele contributed 13.2% phenotypic variance in the population derived from a cross between another diploid A. chinensis female V2 and P2 ( Supplementary Fig. 6c). We demonstrated that the LG27 QTL marker SSRLG27_4396125F4R4 is not associated with Psa resistance in the population derived from a full-sib cross of Psa resistant V1 and P1, where the favorable 428 bp allele is exhibited by neither parent (Supplementary Fig. 6d). In the final population between two Psa resistant parents, V3 × P2, association of the 428 bp allele from P2 with resistance is suggested, but is not statistically significant ( Supplementary Fig. 6e), suggesting interaction with a locus from V3.
The region underlying the QTL on chromosome 27 spans orthologues of genes with putative functions involved in Pathogen-Associated Molecular Pattern (PAMPs)-triggered immunity (PTI). In this region, genes exhibiting non-synonymous substitutions in the coding region have predicted functions in cell wall/ carbohydrate metabolism, cold tolerance and plant immune signaling indicating that more than one mechanism may be involved in control of Psa resistance (Supplementary Table 2).
Additional genetic hotspots associated with tissue and environment-specific phenotypic responses to Psa infection identified using bioassays Analysis of Psa resistance in P1 using the stab assay and leaf infection The stab assay which targets the vascular system enabled a range of different phenotypes to be scored following Psa infection (Fig. 1f, k). P1 appeared to be relatively resistant in comparison with 'Hort16A' in the stab assay, as in the field and grouped close to Psaresistant A. arguta and A. chinensis var. deliciosa for the Stem_necrosis response to infection ( Supplementary Fig.  8). Consistent with this, 'Hort16A' hosted significant growth of endophytic populations of Psa in the leaves, 10 days post-inoculation ( Supplementary Fig. 9), compared with P1 and a Psa-resistant tetraploid A. chinensis genotype both of which did not support endophytic growth of Psa over the same time period ( Supplementary  Fig. 9).
QTLs for control of stem necrosis and collapse, tip death and Psa score determined from the stab bioassay Multiple interval mapping methods identified QTLs for control of Stem_necrosis on LG13 in 'Hort16A' at three positions; S13_6915810 and S13_10678547, with a LOD score ranging between 4.5 and 9 (Fig. 3c) and S13_13629983 (Supplementary Table 3). Moreover, QTLs were detected in the same region on LG13 for control of Stem_collapse and Psa_score_Stab, indicating these were genetic hotspots for host-pathogen interaction in vascular tissues. Interestingly, QTLs for the control of Stem_necrosis in P1 were identified on different chromosomes from those of 'Hort16A', namely the upper arm of LG16 and lower arm of LG23 (Supplementary Table 4 and Supplementary Fig. 10). As for 'Hort16A', QTLs from P1 coincided with those for other phenotypes including Tip_death and Psa_score_Stab. A significant QTL for control of Psa_score_Stab was also detected on LG1 of P1 (Fig. 3f, Supplementary Fig. 11). It was noticeable that the (see figure on previous page) Fig. 2 Distribution of phenotypes in the 'Hort16A' × P1 mapping population of 236 genotypes. a Least squares mean (LSM) of Psa_score_Field after 15 months in the field. The x-axis displays the progression of susceptibility from left to right, while the y-axis represents frequency in the population for the distribution of the trait on the x -axis. b LSM of phenotypes from the stab assay, including Stem_necrosis, Stem_collapse, Tip_death, Psa_score_Stab, Ooze, Leaf_spot, and Wilt. c Means of the health score from Flood bioassays (FA_Week 1 to FA_Week5). The WSTATISTIC is from the Shapiro-Wilks test for the null hypothesis that the distribution is normal. Phenotypic scores with P < 0.001 are rejected for the hypothesis that these distributions are normal. d Principal components analysis on the correlation matrix of the field assessment, flood assay and stab assay measures. Genotypes are shown as points and measurements are shown as vectors (lines pointing from the origin) defined by their correlation with the three principal components Tip_death phenotype generated multiple putative QTLs from both 'Hort16A' and P1 (Supplementary Fig. 12 and 13).

Oozing as a symptom of Psa infection
Oozing of a bacterial exudate was observed following Psa infection and QTLs for control of this phenotype were identified on LGs 2, 13 and LG15 (Supplementary Table 3 Table 4) and LG13. QTLs for control of the Ooze phenotype detected on LG13 and 27 overlapped QTLs detected in 'Hort16A' for the Stem_necrosis phenotype, as well as Psa_score_Field. Other QTLs identified in 'Hort16A' and P1 using KW analysis for the Ooze phenotype are listed in Supplementary Tables 3 and 4, respectively.

Leaf spots and Wilt
We observed symptomatic responses to Psa infection in leaf tissues distant from the point of inoculation in the stem. In 'Hort16A', QTLs for Leaf_spots (Supplementary  Table 3 Table 4 and Supplementary  Fig. 11) were detected on LGs 1 and 5. A significant QTL was detected on LG10 of P1 for Wilt ( Fig. 3e and Supplementary Table 4).

Phenotypic resistance to Psa exposure in tissue culture
When the 'Hort16A' × P1 population grown aseptically in tissue culture was challenged with Psa, multiple QTLs were detected for a health score at each weekly time-point (FA_Week1 to FA_Week5) (Supplementary Table 5). For 'Hort16A', K values were significant on LG15 at the third and fourth weeks following infection. A QTL on LG27 with lower significance overlapped the major QTL on LG27 identified in 'Hort16A' for Psa_score_Field. For P1, a significant QTL identified on the upper arm of LG13 for 3 and 4 weeks post-infection and overlapped the QTL region identified from phenotypes in the stab assay. Plant phenotypes changed dramatically during the period post-infection and additional QTLs were identified for health score at different time points (Supplementary Table 5).
The coordinates for all the QTLs in the Red5 genome versions 1.69.0 33 and 1.68.5 31 are provided in Dataset 1.
Patterns of innate immunity revealed by RNA-seq of 'Hort16A', P1 and F1 genotypes exhibiting Psa resistance or susceptibility in the field RNA-seq performed on healthy young leaf tissues from three groups of 'Hort16A' × P1 F1 genotypes differing in field resistance to Psa demonstrated clear differences in gene expression. The first group included three relatively resistant-to medium-resistant genotypes, including P1 (Psa-RMR), while the second group included three fully susceptible genotypes, including 'Hort16A' (Psa-Sus). At the same time, samples were harvested from the three most resistant 'Hort16A'xP1 genotypes, which had shown resistance for four years in the field (Psa-FR). Heat maps and PCA plots of expression data from the pair-wise comparison between the three groups demonstrated extreme variation between the susceptible (Psa-Sus) and two resistant groups (Psa-RMR and Psa-FR) (Fig. 4a, d). Differential gene expression analysis conducted between the groups of resistant and susceptible genotypes at α < 0.005 with adjusted p values (padj) < 0.1 revealed that from 31,588 genes, 23 (0.076%) were upregulated and 88 (0.28%) were downregulated in Psa-RMR compared with Psa-Sus (Dataset 2). Psa-FR genotypes exhibited 712 differentially expressed genes (DEGs) when compared with Psa-Sus. Of these, 172 (0.59%) were upregulated and 539 (1.9%) were downregulated in Psa-FR genotypes compared to Psa-Sus (Dataset 2). Seventy-seven genes (0.24%) were differentially expressed in common among resistant genotypes of the Psa-FR and Psa-RMR groups when each was compared with Psa-Sus ( LG22 of P1, both for Psa_score_Field. From stab assay phenotypes major QTLs on: (c) LG13 in 'Hort16A' for Stem_necrosis, and in P1 on (d) the upper arm of LG27 for Ooze, (e) LG10 for Wilt and (f) LG1 for Psa_score_Stab. SNPs at peaks are indicated and a key for the QTL mapping models is provided. The red line at LOD 3 represents a minimum threshold level for a candidate QTL and the dashed line at LOD 4.5 represents a statistically significant level for a QTL using a genome-wide permutation test at alpha 0.05. A geom_smooth, a function of ggplot2, is added to the graph in form of a blue line with pink background to represent a smoothed conditional mean with span = 1.1. g shows dot-plot analysis of the allelotypes of markers underlying quantitative trait loci derived from the Psa_score_Field for the population 'Hort 16A' × P1' Fig. 4 RNA-seq analysis of F1 'Hort16A' × P1 genotypes, exhibiting resistance or susceptibility to Psa in the field. RNA-seq was performed on young healthy leaf tissues of field-grown plants belonging to three groups based on relative resistance/susceptibility. The first group included three relatively Psa-resistant plants (Resistant to medium resistant /Psa-RMR). The second group included three fully Psa-susceptible genotypes, including 'Hort16A' (Psa-Sus) and had been exposed to the natural levels of Psa infection in field for 2 years. The third group represents the three most resistant genotypes over 4 years in the field (Psa-FR). All genotypes were from the QTL mapping study (see Methods). a and b show heat-maps for the genome-wide differential expression (DE) analysis in Psa-RMR vs. Psa-Sus and Psa-FR vs. Psa-Sus, respectively. c and d are plots of Principal component analysis for the DE in Psa-RMR vs. Psa-Sus and Psa-FR vs. Psa-Sus, respectively. e shows volcano plots for the DE, with significantly (Log10 adjusted p-value (padj), Log2 fold-change) upregulated and downregulated genes highlighted in red and blue respectively in the two comparisons, Psa-RMR vs. Psa-Sus and Psa-FR vs. Psa-Sus. The gray dots indicate non-significantly expressed genes, whereas the green dots highlight the genes that are differentially expressed in common between Psa-RMR vs. Psa-Sus and Psa-FR vs. Psa-Sus. f shows the DE genes in common or unique, respectively, between Psa-RMR vs. Psa-Sus and Psa-FR vs. Psa-Sus (Table 1). Psa-FR genotypes exhibit upregulation of a high number of genes with functions related to defense. The genes significantly downregulated in common in both resistant genotypes, Psa-RMR and Psa-FR compared with Psa-Sus, are orthologues of protein coding genes involved in chromatin modulation such as histone encoding proteins, auxin efflux, and abiotic and biotic defense ( Table 1).
Furthermore, we explored the expression of these genes in leaf tissues of 'Hort16A' and P1 plants, inoculated in the glasshouse with Psa for bacterial growth assessments ( Supplementary Fig. 9), at 0 and 24 h post-infection. We found that Acc16485.1 (alpha-glucan phosphorylase) and Acc03527.1 (AGAMOUS-like) were significantly upregulated in P1 at 0 and 24 h post-infection compared with 'Hort16A' suggesting that their expression is naturally higher in the resistant parent or suppressed in the susceptible parent and is not induced by pathogen infection (Supplementary Fig. 17). The remaining candidate genes were not differentially expressed between the two parents at either time point, except for Acc08664.1 (Ammonium transporter) which was significantly upregulated in P1 within 24 h post-infection compared with 'Hort16A' (Supplementary Fig. 17).

Discussion
This study provides the first information about genetic loci involved in the host-pathogen relationship between A. chinensis and Psa. Although a genetic map of the chromosomal location of basal defense and R-genes has been reported 34 , there has been no previous genetic mapping of Psa resistance. Our study employed natural field and artificial infection data in three environments over multiple years, combined with genetic and transcriptomic experiments in a segregating population resulting from a cross between Psa-susceptible 'Hort16A' and a resistant male P1, to develop an understanding of the genetic factors underpinning quantitative resistance to Psa in diploid A. chinensis.
QTL mapping of the field phenotypic data following natural infection demonstrated the polygenic nature of this field resistance, with a single major-effect QTL for resistance identified on LG27 in 'Hort16A' and six minoreffect QTLs on LGs 3, 14, 15, 22, 24, and 28 of P1. In addition, we demonstrated the interaction of four of the QTLs (LGs 27, 14, 22, 28), accounting for 30 to 40% of the total variance. Our results are consistent with reports of quantitative resistance against sub-species of Pseudomonas syringae [35][36][37][38] in other hosts and reinforce the longstanding view that no single genetic model can account for incomplete or partial resistance 39,40 . The major QTL on LG27 of 'Hort16A' (initially identified in the field for control of resistance and expressed as Psa_score_Field) overlaps QTLs for tissue specific responses (Fig. 6). These were for the Ooze phenotype in the stab bioassay in both parents (on LG27.1 S27_4621046 in P1 and LG27 4358305 in 'Hort16A') and for the FA_Week3 phenotype in 'Hort16A' (on LG27, S27_4853516). In addition, a number of other QTLs identified from the stab bioassay overlapped in the genomic regions S13_6915810 and S13_10678547 on LG13 (Ooze, Tip_death, Stem_necrosis and Psa_score_Stab) (Fig. 6). As stem necrosis leads to collapse of the vascular structure, we suggest that oozing, together with stem necrosis, is not only an important phenotype for assessing resistance to Psa, but also possibly points towards diverse mechanisms providing field resistance in A. chinensis, that might involve cell wall strengthening and basal defense.
Validation of SSR markers underlying the QTL on LG27 in an independent population of the same cross, as well as three other diploid A. chinensis populations supports association of this region with Psa resistance. Genetic analysis of the polymorphism under the LG27 QTL region in 'Hort16A' × P1 and other populations indicated that resistance to Psa is recessive and there is likely a susceptibility gene(s) in this region of diploid A. chinensis. Further investigation in the kiwifruit germplasm for resistance-associated haplotypes in this region will aid in fine mapping and the search for candidate gene(s) for Psa resistance.
Pyramiding of pest and disease resistance loci to enhance durability is an important focus of most crop breeding programs 40,41 . Marker-Assisted Selection (MAS) has been recognized as a useful tool in breeding perennial fruit crops for major traits such as disease resistance, flowering, ripening [42][43][44][45] and is the most efficient route to pyramiding of resistance loci. The first step towards using MAS to improve the efficiency of breeding new Psaresistant A. chinensis cultivars is the identification of key genetic loci controlling field resistance to Psa. The moderate-high to high resistance to Psa identified in diploid A. chinensis seedlings in PFR breeding populations was reported to be under polygenic control 25 and our study has identified a number of genetic loci associated with field resistance and tissue-specific responses to Psa.
The polygenic nature of resistance to the pathogen is both an advantage and a disadvantage for breeders. Quantitative resistances that aggregate small effects from multiple genes are relatively durable in comparison to qualitative resistances, as virulent pathovars can more readily evade single Resistance (R) gene-based resistance 46,47 . Furthermore, quantitative resistances can also improve the durability of R-gene mediated resistances 48 . However, validation of genetic markers for multiple QTLs in the populations of different ploidy levels that exist in A.  chinensis can be a challenge. As multiple sources of resistance to Psa from a range of species exist in New Zealand kiwifruit germplasm 25,27,49 , resistance pyramiding based on multiple QTLs is a sustainable first approach in a kiwifruit breeding program and can be strengthened in future with yet unidentified R gene resistances against Psa.
The polygenic resistance to Psa in A. chinensis that we have described provides a framework that could lead to the development of durably Psa-resistant cultivars. Pathovars of P. syringae have a complex relationship with their hosts 50 and develop a range of phenotypes in annual or perennial plant species 51 . Additional QTLs were identified associated with tissue-specific responses of A. chinensis to Psa in the stab and flood bioassays and some of these overlapped. For example, QTLs for phenotypes in vascular tissues including Stem_necrosis, Stem_collapse and Ooze were adjacent or overlapped on LGs 13 and 16, but QTLs for leaf-associated phenotypes in the stab assay including Wilt, Leaf_spots and Tip death and overall health score recorded in the flood assay (FA_Week1-5) were located on LGs 3, 5, 7, 10, and 18 (Fig. 5). This is consistent with a previous finding where distinct quantitative genetic variation underlies leafspecific and stem-specific phenotypic responses to a pathogen 52 . As the QTLs located using bioassays were not identified for field Psa resistance, it appears probable that different genetic mechanisms regulate the response to Psa infection in different environments and in different tissues. Many environmental factors differ in greenhouse and in in vitro growth conditions compared to the field, so they might contribute to the plasticity of plant phenotypic responses. This includes factors such as temperature [53][54][55] , humidity [56][57][58][59][60] , other microbial communities 61 in the field, as well as physiological changes during the growth and aging of A. chinensis vines may have an effect. In the future, elucidation of the role of the genetic loci regulating the observed tissue-specific responses to Psa infection will be helpful in determining the dynamics of the hostpathogen relationship in the disease triangle of the A. chinensis/Psa patho-system 62 . Remarkably, a number of the QTLs identified in the bioassays overlie differentially expressed genes, identified from RNA-seq data from field resistant and susceptible genotypes (Fig. 6).
In general, the association of genes determining quantitative resistance with a range of mechanisms of innate immunity or PTI enables them to act effectively to counter the virulence strategies of pathogens during different stages of plant development 63 . In A. chinensis, the genome assembly has demonstrated that more genes are associated with PTI, than with R gene based Effector-Triggered Immunity (ETI), implying a strong selective pressure on the expansion of genes involved in PTI 32 . Further evidence for this idea comes from studies exploring the transcriptome of the kiwifruit-Psa interaction in the period directly following inoculation [64][65][66] . Data obtained from our study have provided a list of classes of gene families underlying the QTLs that might be directly or indirectly involved in the innate immune response of Actinidia, as well as its host-pathogen relationship with Psa over the longer term in the field.
A number of genes in the region underlying the most significant QTL on chromosome 27 are associated with plant defense (Supplementary Table 2). A gene encoding a putative cell wall protein Acc15766 (Acc15766.1), located under the P1 LG14 QTL for field resistance, was employed to design SNP marker E6P3. Two QTLs on LG13 of 'Hort16A' were repeatedly identified in association with control of stem necrosis and health, and Psa score in bioassays, as well as in field screens. Underlying these QTLs were two genes, one an orthologue of Ethylene production protein 1/ETO1 (Acc14810.1) that is intricately linked with a plant's susceptibility to pathogens 67 , and the other a Protein ENHANCED DOWNY MILDEW 2/EDM2 (Acc14938.1), which is involved in DNA methylation, transcriptional regulation and plant resistance to an oomycete pathogen 68 .
In the present study we performed RNA-seq on different groups of F1 genotypes from a single population exhibiting extreme variation in field resistance and susceptibility to natural Psa levels for at least 3 years to explore genes that are associated with Psa resistance and susceptibility in field over an extended time period. A putative orthologue of UGT72B1, which is highly expressed in resistant Psa_RMR and Psa_FR genotypes is localized within 2-LOD interval of Psa_score_Field QTL on LG 27. Association of UGT72B1 with non-host resistance against a fungal pathogen has been suggested, as it encodes an enzyme of the phenylpropanoid pathway 69 . RT-qRT-PCR analysis on samples from controlled inoculation further showed that this gene is induced 24 h post-Psa infection in both 'Hort16A' and P1, however whether this gene is directly responsible for Psa resistance needs to be validated. EDM2 is significantly upregulated in the field-resistant Psa-FR genotypes and co-localizes with the QTL on LG13 associated to stem necrosis and collpase. A gene encoding putative cellulose synthase (Acc15562.1), located on the upper arm of LG14, was upregulated in field-resistant genotypes and might play a role in strengthening the vascular system. On LG24, an orthologue of a Histone protein coding gene (Acc27699.1) that was downregulated in field-resistant genotypes (Psa-FR and Psa-RMR) underlies a P1 QTL that is associated with field resistance.
Furthermore, we verified the expression of the candidate genes associated with plant immunity in Psa-RMR, Psa-FR and Psa-Sus genotypes using gene-specific primers. Consistent with the RNAseq data, we found these genes to be significantly differentially expressed in the resistant genotypes compared to susceptible genotypes. Specifically, Acc16485.1 (Alpha-glucan phosphorylase), Acc03527.1 (AGAMOUS-like) and Acc08664.1 (Ammonium transporter) genes were confirmed to be significantly induced in P1 in greenhouse and field resistant genotypes. Acc03527.1 (AGAMOUS-like) is located very close to the QTL on LG3 in P1 for Psa_score_Field and an ammonium transporter gene has been recently shown to be involved in stem rust resistance in wheat 78 . Our study therefore provides new resource for candidate RNAbiomarkers for predicting resistance in kiwifruit field breeding nurseries, that might lead to improvement of the speed of breeding for multi-genic traits 79 .
A Circos plot of all the QTLs and the DEGs anchored on the Red5 genome 1.69.0, highlighted a number of DEGs that co-localized with the QTL regions (Fig. 6). Circos diagrams for individual phenotypes are presented in Supplementary Figs. 14 to 16 for phenotypes from the field, Stab bioassay and Flood bioassay, respectively.
Expansion of the pathogenic P. syringae strains and their divergence with respect to virulence factors and toxins, as well as antimicrobial compounds 5,80,81 , indicate that the capabilities of this pathogen in suppressing plant defense are remarkable and likely based on targeting multiple host proteins involving diverse post-translational modifications 82 . These modifications have origins in genetic permutations and provide a good target for future breeding strategies 82 . Advances in the genomics of both A. chinensis and Psa make them a powerful plant-pathogen model system in the context of perennial host species. Results from this study will be utilized to develop MAS for Psa resistance in diploid breeding populations and to elucidate the molecular mechanisms to combat the virulent strain of Psa.

Plant material
The two populations for genetic mapping of resistance to Psa were each progeny of a cross between Psasusceptible 'Hort16A' (female) and resistant P1 (male). The first, a pilot population, comprised 53 genotypes that were clonally propagated 3 to 5 times through cuttings,

Phenotyping
The pilot field population was phenotyped monthly for symptoms arising from natural Psa infection between 2013 and 2015 and the data used to develop the phenotypic scoring for the expanded population, which was monthly from February 2017 to September 2018. Traits scored included cane death, ooze, shoot death and tip death (Fig. 1). Presence/absence of leaf spots was not recorded, as scores in the pilot study exhibited high between-plant variability. A cumulative Psa score (Psa_score_Field) was calculated as follows. Under field conditions seedlings developed leaf spots, shoot death, tip death, cane death and oozing cankers, with no ordinal progress for those symptoms. Major secondary symptoms such as a large percentage of cane deaths (two or more cane deaths in a small plant) and oozing requires immediate plant removal under the New Zealand Biosecurity Act (1993). Seedlings displaying those symptoms were given a maximum score of 4 with sub-identifiers such as 4c for oozing cankers and 4d for large cane death compared to a resistant plant with a score of 0. Seedlings with minor symptoms such as tip death and shoot death were given a score of 2; seedlings with leaf spots scored 1. By tracking monthly monitoring, a cumulative score was given to each individual seedling. Least Square Means (LSMs) were calculated for each genotype based on the performance of its clonal replicates at different points in time, to register the progression of the disease. The validation populations were scored for surviving and dead plants after three to seven years of exposure to natural Psa infections in the field.
The bioassays were performed in controlled environments, with the stab bioassay 26 , being performed between September to December and February to April, in 2016, 2017, and 2018. Inoculations were performed with 10627 SmR, a naturally occurring streptomycin-resistant isolate of Psa biovar 3 84,85 , in the greenhouse with temperatures of 22 to 30°C. In total, 200 genotypes were phenotyped using the stab bioassay, with 35 batches phenotyped across three years. Details are in Supplemental Methods S1. The flood bioassay 29 was performed by flooding six biological replicates of each genotype with Psa, that had been grown on tissue-culture media in an aseptic growth medium in a tub for 4 to 6 weeks. Details are provided in Supplemental Methods S1.
Bacterial inoculations for assessment of growth curve in resistant vs. susceptible plants Assessment of the growth curve for Psa in 'Hort16A' and P1 was performed using multiple biological replicates in the greenhouse, as described for the stab test bioassay. Young potted kiwifruit plants, grown under standard growth conditions, were inoculated with Psa, on 8 to 10 biological replicates of each genotype in February, 2018. Further details are provided in Supplemental Methods S1.
Genotyping, genetic maps and QTL mapping DNA was extracted from freeze-dried leaves using the Cetyl trimethylammonium bromide (CTAB) method 86 . GBS libraries were prepared for 53 individuals from the pilot population and 236 individuals from the expanded population, as well as the two parents, using a previously described method 87 , modified from the standard GBS protocol 30 . The individual and pooled libraries were checked for quality with a Fragment Analyzer (Advanced Analytical) and pooled libraries with satisfactory values from quality checks were dried down and dispatched to the Australian Genome Research Facility (AGRF) for single-end sequencing on an Illumina® HiSeq™ platform. The sequencing reads were demultiplexed based on GBS library preparation barcodes using the ea-utils.1.1.2-537 package and those reads starting with the approved barcode immediately followed by the remnant of the BamHI cut site sequence were retained for further analysis. Variant calling and genotyping was performed using TASSEL v3.0 and 5.0 and~60,000 and 80,000 SNP calls were generated for the individuals in the two populations, respectively. SNP calling was performed using an early version of the Red5 genome (1.68.5), which preceded the 1.69.0 version 33 , and the 'Hongyang' genome 32 as references. Genome coordinates for the Red5 version 1.68.5 were converted to those of the published version using in-house PERL scripts (available on request from Ross Crowhurst, PFR). The coordinates for SNPs associated with the QTL peaks in the published Red5 genome are listed in Dataset 1. The SNP calls represented 70% coverage of the expanded population. In our data sets (Dataset 1), Red5 markers begin with S, whereas markers generated from 'Hongyang' begin with HY, followed by the number of the linkage group and the position of the marker on the respective physical genome (for example S1_10661198 or HY10_1385907). The SNP data were subsequently filtered to obtain 9875 and 9327 SNP markers polymorphic between 'Hort16A' and P1, respectively (3364 for P1 in pilot study). JoinMap v 5.0 88 was used to develop genetic linkage maps for both the parents, at a LOD score between 15 and 22. QTL mapping was performed using the rQTL package 89 and MapQTL5 software 90 . Multiple QTL models, including Maximum likelihood (EM), Haley-Knott regression, multiple imputation and Non-parametric/ Kruskal-Wallis analysis (KW) were employed for single QTL scans. All the QTLs > LOD 3 were candidates. The statistical significance of each QTL was validated by performing a genome-wide permutation test. The permutations were performed 1000 times and the genome-wide LOD score (for all 29 chromosomes) was identified at a significance threshold of p < 0.05. The QTL interval was then calculated to identify potential causative SNPs.

RNA-seq of Psa-resistant and -susceptible field-grown plants
The RNA-seq study was performed on leaves harvested from nine different genotypes falling into three different groups of F1 genotypes of the 'Hort16A' × P1 mapping population. These groups are made based on variation observed in resistance response to natural level of Psa infection in the field, over several years. Each group has three different genotypes. The first two groups included genotypes that were grown over two years. These are defined, respectively as: (1) Psa-resistant/medium resistant (Psa-RMR) group, consisting of three relatively Psaresistant genotypes including P1 (genotype IDs 77,78,79) and (2) Psa-susceptible (Psa-Sus) group, comprising three fully Psa-susceptible genotypes including 'Hort16A' (genotypes IDs 82,83,84). These genotypes were part of the large QTL mapping population of 230 genotypes used in this study and resistance and susceptibility in these genotypes was assessed based on observations of 8 to 11 biological replicates. The third group consists of F1 genotypes of "Hort16A" × P1 that belonged to the pilot QTL mapping population (comprising 53 genotypes) and exhibited full resistance to natural levels of Psa infection in the field over four years. This groups is referred to as Psa-FR (Genotypes IDs 73,74,76). Genotypes in all the groups were planted in the same orchard in Te Puke, under standard growth conditions.
Total RNA was extracted from healthy young leaves. Soft green leaves at the sixth to ninth position from the apical leaf on the shoot were harvested at the same time point from all genotypes during Feb, 2018. Leaves samples were snap frozen in liquid nitrogen. RNA extraction was performed using the Spectrum Total Plant RNA kit (Sigma-Aldrich, Auckland, New Zealand) and QC was performed with the Fragment Analyzer to select RNA with RNA Integrity Number (RIN) of 7.1-8.2.
Each sample for a genotype in Psa-RMR and Psa-FR group, consists of a pool of three clonal replicates of the respective genotype. However, the sample for each genotype in Psa-FR is derived from a single clonal replicate as these genotypes are from the pilot trial for which only a single replicate each genotype was retained in the orchard. Library preparation at the Australian Genome Research Facility used the TruSeq Stranded kit and subsequent paired-end Illumina® sequencing employed the NovaSeq6000 platform. An average of~19 million, 150 bp paired-end reads were retrieved for each sample (~6 Gb) and read sequences of low-quality, ribosomal RNA, as well as adapters were filtered out using Trimmomatic 91 and SortMeRna 92 . RNA-seq reads were aligned to the Red5 reference gene models using STAR and differential expression analysis was performed, using DESeq2 93 .
Differential analysis was performed firstly in between the Psa-RMR and Psa-Sus group and secondly in between the Psa-FR and Psa-Sus group. The differentially expressed genes (DEGs), at α < 0.005 with adjusted p values (padj) < 0.1 and log2 fold-change, that were in common between the two comparisons, were considered candidates for control of Psa resistance since they are differentially expressed between Psa resistant and susceptible genotypes. The expression profile for these genes was than verified using gene-specific primers with RT-qRT-PCR in representative samples from all the groups. The details of the RNAseq target file, read statistics and DEGs are provided in Dataset 2.

RT-qRT-PCR for the DEGs
To validate the expression of the DEG genes in the RNAseq study, RT-qRT-PCR was performed on the RNA used for the RNA-seq study using gene-specific primers. This included samples from 3-7 replicates per genotype in the Psa-RMR and Psa-Sus group and 5 different genotypes in the Psa-FR group. Results for the three representative genotypes are presented. Total RNA (~2 μg) was treated with DNase I (Roche Applied Sciences) and used for cDNA synthesis using SuperScript IV Reverse Transcriptase (Life Technologies-Invitrogen). The cDNA was diluted 20-fold and used for qRT-PCR employing a LightCycler ® 480 SYBR Green 1 Master PCR labeling kit (Roche Applied Sciences) and RotorGene 3000 Real time PCR machine (Corbett Research, Sydney, Australia). Relative transcript abundance was determined by normalizing to the global mean of the expression of the two house-keeping genes, Actin and Ubiquitin in the same sample. Comparative quantification was performed using the mathematical model for relative quantification 94 . Primers used for gene amplification are provided in Supplementary Table 6.

RT-qRT-PCR for candidate genes in plants when challenged with Psa artificially
Due to strict Kiwifruit Vine Health (KVH) regulations in New Zealand, artificial Psa infection in field was not possible. The expression of the candidate genes, identified from the RNA-seq study, was therefore tested in 'Hort16A' and P1 plants challenged with Psa under standard greenhouse conditions. Details for bacterial infection and harvesting of samples for gene expression is provided in Supplemental Methods 1. In brief, a 10 mm leaf disc was harvested for RNA extraction from the region infected with Psa, at 0 h (before infection) and at 6, 24, 48 h time points post-infection, and frozen in liquid nitrogen. One leaf disc from inoculated and one from sterile water-treated area was harvested per time point and three biological replicates were harvested at each time point. RNA extraction, cDNA synthesis and relative transcript abundance was determined as described above.

SSR and SNP marker design and screening
Repeats were identified manually in the genome sequence underlying the QTLs. PCR primers for SSR markers were designed using Primer3 and employed to screen DNA extracted from the populations 95 . Analysis and scoring of the alleles in the amplicons was performed on a Hitachi ABI3500 Applied Biosystems genetic analyzer. Primers were also designed around SNPs in the genes identified in the genomic sequence of Red5 underlying the QTLs. The SNP markers were screened using real-time High Resolution Melting analysis 96 . All primer sequences are provided in Supplementary Table 6.