Comparison of defense responses of transgenic potato lines expressing three different Rpi genes to specific Phytophthora infestans races based on transcriptome profiling

Plant materials

Tetraploid potato cultivar Desiree and its three transgenic lines with the genes R1, R3a and R3b were used. Sterile culture of in vitro plants was carried out in culture bottles (10 cm height × 6.5 cm diameter) containing MS medium by single segmental regeneration and asexual propagation. All samples were cultured for 4 weeks under 24 °C 16 h light/8 h dark conditions with 10 seedlings per bottle.

Inoculation with P. infestans isolates

P. infestans isolates 89148-9 (race 0) and CN152 (race 1, 3b, 4, 5, 6, 7, 8, 9, 10, 11) (Li et al., 2017) were used in this study. Sporangium formation was induced by activation of P. infestans on oat agar medium, and cultured in the dark at 18 °C for 7–10 days. Pre-chilled sterile water was added to the culture dish and placed at 4 °C for 3 h to release zoospores. The spore concentration was counted using a Fuchs-Rosenthal counting chamber. The final spore concentration was adjusted to 5 × 10⁴ spores/mL. At least one week before infection, four-week-old plants were transferred to an infection chamber with 100% humidity and a 10:14 light:dark cycle. Plants were sprayed with an encysted zoospore suspension from P. infestans isolates 89148 and CN152 until the leaf surfaces were fully saturated with the zoospore suspension.

RNA extraction and sequencing

Sampling was carried out at 24-hours after inoculation of P. infestans, and 10 seedlings of each treatment were collected at the same time. The roots were removed and all aboveground tissues of 10 seedlings were mixed as one sample. A total of 8 samples were collected. The samples were frozen with liquid nitrogen and stored at −80 °C for subsequent RNA extraction and sequencing. Total RNA of samples was extracted using TRIzol reagent (Invitrogen, Carlsbad, CA, USA), digested by TURBO DNase I (Ambion, Austin, TX, USA), precipitated and dissolved by ethanol, and detected using β-actin as a reference gene. RNA concentration and OD260/OD280 were observed using an ultraviolet spectrophotometer and 1% agarose gel electrophoresis. Non-degraded RNAs with OD260/OD280 ∼2.0 were used to generate the cDNA sequencing library with the NEBNext Ultra™ RNA Library Preparation Kit (NEB, USA) according to the manufacturer’s recommendations. The RNA quality was checked using an Agilent BioAnalyzer 2100 system and cDNA libraries were sequenced on the Illumina HiSeq 2500 platform.

Raw data quality assessment

Raw data were processed by removing reads containing adapter, reads containing poly-N and low-quality reads, and then the clean data (clean reads) were retained. At the same time, Q20, Q30, GC-content and sequence duplication level of the clean data were calculated. All the downstream analyses were based on high-quality clean data.

Comparison of reference genomes and functional annotations

Clean reads were then mapped to the potato reference genome sequence. The reference accession, the doubled haploid S. tuberosum Group Phureja clone DM1-3 516R44 (hereafter referred to as DM) genome sequence (SolTub 3.0) and annotation files were downloaded from the ENSEMBL plants database (http://plants.ensembl.org/Solanum_tuberosum/Info/Index) (Bolser et al., 2017). Only reads with a perfect match or one mismatch were further analyzed and annotated based on the reference genome. TopHat2 tools soft were used to align RNA-seq reads against the reference genome (Kim et al., 2013).

Identification of differentially expressed genes

Gene expression levels were estimated by fragments per kilobase of transcript per million fragments mapped (FPKM) and normalized using HTseq-count and EBSeq (Florea, Song & Salzberg, 2013; Leng et al., 2013; Anders, Pyl & Huber, 2015). Differential expression analysis was performed using EBSeq based on the negative binomial distribution (Anders & Huber, 2010). Normalized FPKM of genes expressed in the transgenic R1, R3a, and R3b plants after inoculation was compared with FPKM in the wild-type Desiree, and the DEGs were identified. Also, the DEGs involved in incompatible/compatible reactions were analyzed by comparing transcriptomic data in each of the transgenic lines infected by the two tested P. infestans strains with different virulences. The resulting p-values were adjusted using the Benjamini–Hochberg approach for controlling the false discovery rate (FDR). Genes with normalized expression fold-change greater than 2, and FDR < 0.05 were assigned as differentially expressed. The DEGs were annotated based on the functional annotation information of ENSEMBL plants release S. tuberosum SolTub_3.0 genes.

DEGs identified in 89148 treated samples were represented by 8_TR1, 8_TR3a and 8_TR3b, and DEGs identified in CN152 treatments samples were represented by CN_TR1, CN_TR3a and CN_TR3b.

GO and KEGG enrichment

Using the potato SolTub_3.0 database as a reference, the gene ontology GO (Gene Ontology) enrichment analysis of DEGs was carried out by agriGO (http://bioinfo.cau.edu.cn/agriGO/analysis.php). KOBAS software (http://kobas.cbi.pku.edu.cn/index.php) was used to analyze the significant enrichment of DEGs in KEGG pathways (Mao et al., 2005). The significance of enrichment of each GO and KEGG term was assessed by p-value <0.05 or FDR < 0.05.

Validation of DEGs by qRT-PCR

Sixteen DEGs were randomly selected for quantitative real-time PCR (qRT-PCR) to verify the RNA-seq results. Specific primers were designed for qRT-PCR analysis using Primer 5 and listed in Table S1. The RNA samples for qRT-PCR were prepared by the same methods with RNA-seq and reverse-transcribed to cDNAs using PrimeScript™ RT reagent kit with gDNA Eraser (TaKaRa). qRT-PCR was performed using SYBR® Green PCR master mix reagents with 10 µL SYBR Premix Ex Taq (TaKaRa, Japan), 0.5 µL of both forward and reverse primers, 7 µL of double-distilled water and 2 µL (40 ng/µL) cDNA in one reaction system. Then qRT-RCR reactions were performed on the Bio-Rad CFX96 C1000™ instrument with the cycle steps of pre-denaturation at 95 °C, 30 s, 1 cycle; 95 °C, 5 s, 60 °C, 30 s, 40 cycles. Gene expression was analyzed by CFX-manager software (CFX96 Real-Time System; Bio-Rad). The relative expression level of each gene was calculated by the 2^−ΔΔCt method (Livak & Schmittgen, 2001) using the GAPDH gene as an internal reference. All assays were performed in triplicate under identical conditions and correlations between qRT-PCR and RNA-seq data were evaluated using Pearson correlation coefficients.

Clustering of late blight resistance genes

A total of 20 cloned resistance (R) genes to potato late blight, including R1, R3a and R3b were used in this study. The CDS of these genes were downloaded from the website of the National Center for Biotechnology Information (https://www.ncbi.nlm.nih.gov/). The phylogenetic tree of these 20 R genes was constructed using the neighbor-joining (NJ) method in MEGA 5 (Tamura et al., 2011).

Transcriptome sequencing and mapping

The transgenic R1, R3a, R3b lines and wild-type Desiree were infected with two P. infestans isolates, 89148 (universal race 0) and CN152 (super-race 1, 3b, 4, 5, 6, 7, 8, 9, 10, 11). At 24 h post inoculation (hpi), all three transgenic lines showed resistance to 89148 (Figs. 1A, 1B, and 1C), however, all three transgenic lines expressed susceptibility to CN152 (Figs. 1E, 1F, and 1G). Wild-type Desiree was susceptible to both P. infestans isolates while did not show obvious susceptible symptom until 72 hpi under 89148 infection (Figs. 1D and 1H, Yang et al., 2018).

To unravel the role of different R genes in the early molecular response to P. infestans both in incompatible and compatible interactions, the present study included transcriptional profiling of wild-type (WT) and the transgenic plants R1, R3a, and R3b (hereafter referred as TR1, TR3a, and TR3b) 24 h after infection of 89148 and CN152 by RNA-seq. In total, eight RNA libraries derived from seedlings sampled at 24 hpi were sequenced using the Illumina HiSeq 2500 system with the 125-cycle paired-end sequencing protocol.

After data quality assessment, a total of 184,632,455 clean paired-end reads (46.52 Gb) were generated with an average of 23,079,057 paired end reads (5.81 Gb) for each sample. Over 92.99% of these reads were ≥Q30 and the average GC content was 43.85% (Table 1). The clean data were then aligned to the potato DM reference genome, and 74.43%–76.52% of the reads per sample mapped to the reference genome, while 72.28%–74.51% of the reads per sample uniquely aligned to the reference genome (Table 1).

Differentially expressed gene analysis assessment

DEGs were identified by comparing the transgenic materials to wild-type Desiree (Table 2). Under infection by the isolate 89148, the number of DEGs in transgenic R3b line was the highest (1763), followed by R1 (1318), and the number of DEGs in transgenic R3a was the lowest (1171). More down-regulated DEGs (728, 627) were detected than up-regulated DEGs (590, 544) in TR1 and TR3a. Almost the same numbers of up- and down-regulated DEGs were present in TR3b. The same trend was observed under CN152 inoculation with the highest number of DEGs in TR3b (1857) and the lowest number of DEGs in TR3a (1096), including 408, 456, and 929 up-regulated DEGs and 718, 640, and 928 down-regulated DEGs in the TR1, TR3a, and TR3b lines, respectively.

The Venn diagram showed 408, 278, and 779 unique DEGs in TR1, TR3a, and TR3b, respectively, when inoculated with 89148 (Fig. 2A), and 267, 192, and 864 DEGs when inoculated with CN152, respectively (Fig. 2D). The most unique DEGs were identified in transgenic R3b lines, followed by transgenic R1 and R3a lines. Under infection with the same pathogen, more DEGs were involved in the resistance to pathogen infection in transgenic R3b lines, while transgenic R3a lines can defend against infection by mobilizing fewer genes. The tendency of up- and down-regulated DEGs is the same as that of total DEGs in transgenic lines (Figs. 2B, 2C, 2E and 2F). This indicates the weakest resistance in R3b and the strongest resistance in R3a, though both genes exist within the same R gene cluster.

When inoculated with 89148, TR1 and TR3a shared the fewest DEGs, with 134 genes overlapping. Approximately the same number of DEGs were shared by TR1 and TR3b, TR3a and TR3b with 225 and 208, respectively (Fig. 2A). The same trend was evident under CN152 inoculation. The common DEGs shared between TR1 and TR3a, TR1 and TR3b, TR3a and TR3b, were 118, 207, and 252, respectively (Fig. 2D). These results suggest that R1 and R3a have more similar resistance levels than R1 vs. R3b, and R3a vs. R3b.

Enriched KEGG pathways in DEGs

KEGG metabolic pathways were analyzed across all DEGs. Under the P. infestans 89148 treatments, the DEGs derived from R1, R3a, and R3b lines were mapped to 89, 88, and 98 canonical reference pathways, respectively (Fig. 3A). Under CN152 treatments, DEGs were enriched in 86, 77 and 102 canonical reference pathways, respectively (Fig. 3B). More differentially expressed pathways were identified in TR3b than in TR3a and TR3b. As TR3b is the most susceptible genotype, more pathways may be involved to achieve the purpose of disease defense. Eight, 12 and 11 reference canonical pathways were significant (p<0.05) for 89148 isolate (Fig. 3C). For infection by CN152, 14, 14 and 15 reference canonical pathways were significant with p-value <0.05 (Fig. 3D).

When inoculated with 89148, five common pathways in KEGG were significantly enriched (p <0.05) in the three transgenic lines (Fig. 3C). These include protein processing in the endoplasmic reticulum (sot04141), phenylpropanoid biosynthesis (sot00940), amino sugar and nucleotide sugar metabolism (sot00520), starch and sucrose metabolism (sot00500), as well as biosynthesis of secondary metabolites (sot01110). The biosynthesis of secondary metabolites has the most DEGs with 58, 74, and 91 genes enriched for TR1, TR3a, and TR3b (Fig. 4A). Under inoculation with CN152, eight common pathways in KEGG were significantly enriched (p < 0.05) in the three transgenic lines (Fig. 3D), including protein processing in the endoplasmic reticulum (sot04141), phenylpropanoid biosynthesis (sot00940), biosynthesis of secondary metabolites (sot01110), metabolic pathways (sot01100), starch and sucrose metabolism (sot00500), zeatin biosynthesis (sot00908), alpha-linolenic acid metabolism (sot00592), and sesquiterpenoid and triterpenoid biosynthesis (sot00909) (Fig. 4B). Metabolic pathways (sot01100) had the most abundant DEGs, followed by biosynthesis of secondary metabolites (64, 75, 100). Amino sugar and nucleotide sugar metabolism (sot00520) is the specific pathway involved in the incompatible action by 89148 with 13, 10, and 14 DEGs for TR1, TR3a, and TR3b (Table S2); while metabolic pathways (sot01100), zeatin biosynthesis (sot00908), alpha-linolenic acid metabolism (sot00592) and sesquiterpenoid and triterpenoid biosynthesis (sot00909) are the specific pathways involved in the compatible action by CN152.

Specific enriched pathways in different transgenic R1, R3a, and R3b lines

Venn diagram results showed there were two, four, and two significant (p < 0.05) KEGG pathways specific to TR1, TR3a and TR3b, respectively, under inoculation with 89148 (Fig. 3C), and one, two, and two pathways, respectively, specific under CN152 infection (Fig. 3D). The top enriched pathway is unique to each transgenic line.

Under 89148 infection, DNA replication (sot03030), plant-pathogen interaction (sot04626), and pentose and glucuronate interconversions (sot00040) are the top significantly enriched pathways in TR1, TR3a, and TR3b, respectively (Table 3). In DNA replication pathway, nine genes were enriched in TR1 lines, which contained four minichromosome maintenance (MCM) protein genes MCM4, MCM5, MCM6, and MCM7. These MCM coding genes were all down-regulated. The other two genes of unknown function (PGSC0003DMG400011569 and PGSC0003DMG402011204) had highly increased expression of up to 25-fold.

In the pathway of plant-pathogen interaction, a total of 16 genes were differentially expressed in TR3a lines, which encoded pathogenesis-related protein 1 (PR-1, two members), heat shock protein (HSP, three members), respiratory burst oxidase homolog protein (RBOH, two members), receptor protein kinase (two members), calcium-dependent protein kinase (CDPK, two members), and CC-NBS-LRR resistance protein (one member). One PR-1 gene (PGSC0003DMG400005109) had the highest fold change of 27.6. In the pathway of pentose and glucuronate interconversions for TR3b lines, a total of 14 DEGs were enriched and encoded pectinesterase (five members) and pectase lyase (nine members). Except for one pectinesterase coding gene (PGSC0003DMG400000293) that was up-regulated, the other genes were down-regulated in TR3b lines (Table S3). This indicated that most down-regulated DEGs played were essential in the response to 89148 infection.

Under CN152 infection, only the pathway of cutin, suberine, and wax biosynthesis (sot000735) was significant for TR1 lines including five DEGs that encoded feruloyl transferase (PGSC0003DMG400031731), Acyl CoA reductase (PGSC0003DMG400007405), CYP86A33 fatty acid omega-hydroxylase (PGSC0003DMG400002111), Fatty acyl-CoA reductase 3 (PGSC0003DMG400007113) and Fatty acyl-CoA reductase 2 (PGSC0003DMG400008885). For diterpenoid biosynthesis (sot00904) in TR3a lines, one copalyl diphosphate synthase gene (PGSC0003DMG400001948), one GA 3-oxidase gene (PGSC0003DMG400005698), and three gibberellin (GA) 2-oxidase genes (PGSC0003DMG400027632, PGSC0003DMG400021095 and PGSC0003DMG400027631) were enriched. For TR3b lines, the top pathway under CN152 was pentose and glucuronate interconversions (sot00040), which was the same as that under 89148. All 14 DEGs were down-regulated with five genes encoding pectinesterase, eight encoding pectase lyase and one encoding polygalacturonase 7 (PGSC0003DMG400020372) (Table S4). This indicated most down-regulated DEGs functioned directly in response to CN152 infection.

To further investigate the regulation mechanism of resistance-related genes, the up- or down-regulated DEGs separately in the KEGG enrichment were separately analyzed. For TR1, TR3a, and TR3b by 89148 infection, the up-regulated DEGs were significantly enriched in four, ten and eight KEGG metabolic pathways, respectively (Fig. 5A). The top pathways in each were protein processing in the endoplasmic reticulum (sot04141), photosynthesis-antenna proteins (sot00196) and protein processing in the endoplasmic reticulum (sot04141), respectively. The down-regulated DEGs of TR1, TR3a, and TR3b were significantly enriched in ten, eleven, and 13 KEGG metabolic pathways, respectively (Fig. 5B). The top two pathways for all three were phenylpropanoid biosynthesis (sot00940) and MAPK signaling pathway-plant (sot04016). The pathway of plant-pathogen interaction (sot04626) was enriched in the down-regulated DEGs of TR1 (10 members), TR3a (12 members), and TR3b (11 members) (Table S5). This indicated that the down-regulation of DEGs was essential in the incompatible interaction between the transgenic plants and P. infestans isolate 89148.

For TR1, TR3a, and TR3b with CN152 infection, the up-regulated DEGs were significantly enriched in nine, seven, and ten KEGG metabolic pathways, respectively (Fig. 5C). The top pathway was protein processing in the endoplasmic reticulum (sot04141) with 29, 21, and 38 DEGs, respectively. The down-regulated DEGs of TR1, TR3a, and TR3b were significantly enriched in 12, 14, and 17 KEGG metabolic pathways, respectively (Fig. 5D). The top pathways were phenylpropanoid biosynthesis (sot00940), metabolic pathways (sot01100), and MAPK signaling pathway-plant (sot04016). The pathway of plant-pathogen interaction (sot04626) was enriched in the down-regulated DEGs of TR1 (13), TR3a (9), and TR3b (14) (Table S5). This indicated that the down-regulation of DEGs are required in the compatible interaction between the transgenic plants and P. infestans isolate CN152.

Gene ontology classification of DEGs

To get an overview of the functional category of the genes that participated in the P. infestans infection response, the DEGs were subjected to Gene Ontology (GO) enrichment analysis. Under infection by 89148, 57.4% of the total DEGs in transgenic R1 lines (TR1) were enriched for different GO terms, including 53 significantly enriched terms (FDR < 0.05). From DEGs in TR3a and TR3b, 60.6% and 57.1% of the total genes were enriched for different GO terms, and 69 and 68 terms were significantly enriched (FDR < 0.05) (Table 4). After inoculation with CN152, among the DEGs generated from transgenic R1, R3a, and R3b lines, 59.8%, 60.0% and 57.0% of the genes were annotated with GO terms. Of these GO terms, 55, 62, and 74 were significant at FDR < 0.05 (Table 4).

A Venn diagram analysis showed that there were 13, 21, and 15 GO terms specific to TR1, TR3a, and TR3b, respectively, under inoculation with 89148, and eight, one and 17 terms were specific under CN152 infection. In both inoculation conditions, 31 common GO terms were shared across TR1, TR3a, and TR3b (Figs. 6A and 6B). TR3a and TR3b shared the largest number of enriched GO terms (15 for 89148 and 20 for CN152) (Figs. 6A and 6B). TR1 and TR3a shared the smallest number of GO terms under 89148 infection (two common terms), and under CN152 infection the smallest number of GO terms were shared across TR1 and TR3b (six common terms) (Figs. 6A and 6B). This indicated that the defense pathways mediated by R3a and R3b are more similar than those shared between R1 and R3a, or R1 and R3b in both incompatible and compatible interactions of potato and P. infestans.

GO enrichment analysis is broken down into three categories: biological process (BP), molecular function (MF), and cellular component (CC). Among the 15 common GO terms shared between TR3a and TR3b under 89148 infection (Fig. 7A), a total of 11 terms were enriched in the BP category, including cellular amino acid derivative catabolic process (GO:0042219), oxidation reduction (GO:0055114), and glucuronoxylan metabolic process (GO:0010413). The biological process of oxidation reduction (GO:0055114) had the most DEGs, at 129 in TR3a and 150 in TR3b (Table S6). The terms of tetrapyrrole binding (GO:0046906), oxidoreductase activity (GO:0016491), and lyase activity (GO:0016829) were under the MF category. Oxidoreductase activity (GO:0016491) had the most DEGs, at 124 in TR3a and 151 in TR3b.

Among the 20 GO terms shared by transgenic R3a and R3b lines inoculated with CN152, 16 terms belonged to the BP category, three were under MF, and one was enriched in the CC category (Fig. 7B). For the BP category, aromatic compound catabolic process (GO:0019439), polysaccharide catabolic process (GO:0000272), and response to stress (GO:0006950) were among the enriched terms. Response to stress (GO:0006950) had the most DEGs at 96 in TR3a and 132 in TR3b (Table S7). The shared terms enriched in MF are copper ion binding (GO:0005507), carbon-oxygen lyase activity, acting on polysaccharides (GO:0016837), and pectate lyase activity (GO:0030570).

Specific enriched terms in different transgenic R1, R3a, and R3b lines

Among the specific GO terms in different transgenic lines, the top enriched GO term was unique to each transgenic line. Under 89148 infection, DNA unwinding during replication (GO:0006268, 5 DEGs) (Fig. 8A), protein-chromophore linkage (GO:0018298, 15 DEGs) (Fig. 8B), and cell wall organization or biogenesis (GO:0071554, 50 DEGs) (Fig. 8C) were the top significantly enriched terms in BP category in TR1, TR3a, and TR3b, respectively. In MF, hydrolase activity (GO:0016787) was the top term for TR1 with 156 DEGs (Fig. 8A). Further, the top terms in TR3a and TR3b were chlorophyll binding (GO:0016168, 15 DEGs) (Fig. 8B) and carbon-oxygen lyase activity, acting on polysaccharides (GO:0016837, 9 DEGs) (Fig. 8C), respectively.

When inoculated with CN152, the top significantly enriched biological process and molecular function in TR1 were aminoglycan metabolic process (GO: 0006022) and monooxygenase activity (GO: 0004497) with seven and 36 DEGs involved, respectively (Fig. 8D). In TR3a, only one molecular function GO term, UDP-glucosyltransferase activity (GO:0035251), was significantly enriched with specificity to TR3a (Fig. 8E). This term included eight DEGs encoding sucrose synthase, cellulose synthase and amylase synthase. Among these genes, Arabidopsis UDP-glycosyltransferase superfamily protein homologous gene PGSC0003DMG400012111 had the highest expression level across all genotypes, which was down-regulated in the transgenic lines compared to wild-type Desiree. Overexpression of the UDP-glucosyltransferase gene (Ta-UGT 3) can enhance resistance to Fusarium Head Blight in wheat (Xing et al., 2018). In TR3b, more GO terms were significantly enriched compared to TR1 and TR3a. The top three biological process terms for TR3b were pectin catabolism (GO:0045490, 15 DEGs), cellular polysaccharide catabolism (GO:0044247, 15 DEGs), and regulation of catalytic activity (GO:0050790, 43 DEGs). Among the molecular function GO terms, only cellulose activity (GO:0008810) was enriched with 7 DEGs (Fig. 8F).

To further investigate the regulation mechanism of resistance-related genes, the up- or down-regulated DEGs in the GO enrichment were separately analyzed. Under 89148 inoculation, the 590, 544, and 880 up-regulated DEGs of the three transgenic lines were enriched in 20, 42, and 28 GO terms (Table 5), including 16, 20, and 20 biological process (BP) terms, respectively. The top 20 enriched BP terms are in Fig. 9. As with the enriched BP terms in up-regulated genes under 89148 infection, the main enriched BPs were positive or negative regulation of hydrolase activity, peptidase activity, and endopeptidase activity. The response to abiotic stimulus (GO: 0009628) had the most DEGs at 27, 25, and 32 for TR1, TR3a, and TR3b, respectively (Fig. 9A). More GO terms were involved in down-regulated genes than in up-regulated genes. The 728, 627, and 883 down-regulated DEGs of the three transgenic lines were enriched in 92, 51, and 73 GO terms (Table 5) which included 43, 36, and 48 BP terms, respectively. The BP terms of the three transgenic lines all contained response to stress (GO:000695) with 72, 61, and 74 DEGs, respectively (Fig. 9B, Table S8). In addition, 40 DEGs were enriched in defense response (GO: 0006952) in TR1. The highest numbers of DEGs enriched in TR3a and TR3a were oxidation reduction (GO: 0055114) at 83 and 103 DEGs. As above, after 89148 inoculation, the DEGs were mainly down-regulated and mediated the resistance of plants to P. infestans 89148 through response to stress (GO:000695) biological process.

Under CN152 infection, the 408, 456, and 929 up-regulated DEGs of the three transgenic lines were enriched in 30, 15, and 23 GO terms (Table 5), including 24, 12, and 17 BP terms, respectively. The top 10–20 enriched BP terms are in Fig. 10. As with the enriched BP terms in up-regulated genes under CN152 infection, the main enriched BPs were also the positive or negative regulation of hydrolase activity, peptidase activity, and endopeptidase activity. The highest number of DEGs were also enriched in response to abiotic stimulus (GO: 0009628) with 20, 22, and 34 DEGs for TR1, TR3a, and TR3b, respectively (Fig. 10A). More GO terms were involved in down-regulated genes than in up-regulated genes. The 718, 640, and 928 down-regulated DEGs of the three transgenic lines were enriched in 35, 69, and 78 GO terms (Table 5), which included 11, 31, and 40 BP terms, respectively. The highest number of DEGs enriched in all three transgenic lines were oxidation reduction (GO: 0055114) at 92, 95, and 116 DEGs, respectively (Fig. 10B, Table S9). As above, the DEGs were mainly down-regulated and mediated the resistance of the plants to P. infestans strain CN152 through oxidation reduction biological process.

Analysis of DEGs involved in incompatible/compatible reactions in transgenic lines in response to different P. infestans strains

The genes expressed under 89148 infection were compared with those expressed under CN152 infection to identify the specific DEGs involved in incompatible/compatible reactions in the same transgenic lines. TR3b had the most DEGs with 342, followed by TR1 (329 DEGs). TR3a had the least DEGs at 270 (Fig. 11A). This indicated that R3a had the strongest resistance and R3b had the weakest resistance; the trend was the same as that generated by gene comparison between transgenic lines and wild-type Desiree under the same P. infestans strain infection. The number of up-regulated genes (198,174, 231 members, respectively) was higher than that of down-regulated genes (131, 96, 111 members, respectively) in all transgenic lines.

Venn diagrams showed that R1R3a, R1R3b, and R3aR3b shared 39, 40, and 52 genes, respectively. The unique DEGs were 224, 152, and 225 for TR1, TR3a, and TR3b (Fig. 11B). According to a value of log2FC >5 or < −5, a total of 14, nine, and 13 genes were highly differentially expressed in TR1, TR3a, and TR3b, respectively (Fig. 12). The 14 highly expressed DEGs for TR1 encoded anthocyanin 5-aromatic acyltransferase, disease resistance protein, and zinc finger protein 4. The nine highly expressed DEGs for TR3a encoded receptor kinase THESEUS 1, NBS-LRR type disease resistance protein, and polygalacturonase-1 non-catalytic subunit beta. The 13 highly expressed DEGs for TR3b encoded fatty acid desaturase, verticillium wilt disease resistance protein, and bHLH transcription factor. These highly differentially expressed genes may be essential in the incompatible/compatible reactions in transgenic lines response to different P. infestans strains 89148 and CN152.

The 224, 152, and 225 unique DEGs were respectively enriched in 125, 61, and 113 GO terms for TR1, TR3a, and TR3b. A total of 31 GO items were significantly (p <0.05) enriched in TR1 including the biological processes of defense response (GO:0006952, 14 DEGs), response to stress (GO:0006950, 18 DEGs), and response to stimulus (GO:0050896, 21 DEGs). No GO terms were significantly enriched in TR3a, which may be related to TR3a having the fewest DEGs. The top biological process in TR3a was oxidation reduction (GO:0055114). A total of 12 GO terms were significantly enriched in TR3b including response to heat (GO:0009408, 5 DEGs) and response to stress (GO:0006950, 18 DEGs) (Fig. 13). The top 10 biological process GO terms were listed in Fig. 13. Among the genes involved in response to stress biological process, 18 DEGs for TR1 mainly encoded NBS-LRR resistance proteins and disease resistance proteins; five DEGs for TR3a encoded NBS-LRR type disease resistance protein and ethylene-responsive late embryogenesis; 18 DEGs for TR3b mainly encoded pentatricopeptide repeat protein, p-coumarate 3-hydroxylase, and ATP binding protein (Table 6). These specific DEGs may function directly in the response to biotic stress by different P. infestans strains.

Clustering of cloned late blight resistance genes

Cluster analysis of 20 cloned late blight resistance genes including R1, R3a and R3b was performed based on the complete CDS sequences using MEGA5 software. The cluster diagram showed that R3a and R3b were clustered together and supported with a 100% bootstrap value, while R1 was distant from R3a and R3b, and clustered with Rpi-blb2 gene in a separate branch with a 100% bootstrap value (Fig. 14). The clustering relationships of cloned late blight resistance genes were consistent with the results of gene expression profiling.

Validation of RNA-seq data by real-time quantitative PCR (qRT-PCR)

To validate the results of RNA-seq data, we randomly selected 16 DEGs in potato R1-, R3a-, and R3b-expressing transgenic plants 24 h after infection by P. infestans isolates 89148 and CN152, including ERF transcription factor 5, zinc finger protein and leucine-rich repeat family protein coding genes (Table S3). Quantitative RT-PCR was performed on these 16 DEGs. Gene expression was altered in each of the transgenic lines relative to WT Desiree plants. The fold changes of PGSC0003DMG400008364, PGSC0003DMG400000417, PGSC0003DMG400010128, and PGSC0003DMG400010139 genes were over 10 with log2 >3.4 (Table S4). The measured relative expression levels of these 16 DEGs were similar and a high correlation (R²>0.92) was observed between qRT-PCR and RNA-seq data (Fig. 15), confirming the reliability of the transcriptome data. The comparison of RNA-seq to log2 qRT-PCR levels is shown in Table S10.