Transcriptome analysis reveals that exogenous ethylene activates immune and defense responses in a high late blight resistant potato genotype

Ethylene (ET) is one of the many important signaling hormones that functions in regulating defense responses in plants. Gene expression profiling was conducted under exogenous ET application in the high late blight resistant potato genotype SD20 and the specific transcriptional responses to exogenous ET in SD20 were revealed. Analysis of differentially expressed genes (DEGs) generated a total of 1226 ET-specific DEGs, among which transcription factors, kinases, defense enzymes and disease resistance-related genes were significantly differentially expressed. GO enrichment and KEGG metabolic pathway analysis also revealed that numerous defense regulation-related genes and defense pathways were significantly enriched. These results were consistent with the interaction of SD20 and Phytophthora infestans in our previous study, indicating that exogenous ET stimulated the defense response and initiated a similar defense pathway compared to pathogen infection in SD20. Moreover, multiple signaling pathways including ET, salicylic acid, jasmonic acid, abscisic acid, auxin, cytokinin and gibberellin were involved in the response to exogenous ET, which indicates that many plant hormones work together to form a complex network to resist external stimuli in SD20. ET-induced gene expression profiling provides insights into the ET signaling transduction pathway and its potential mechanisms in disease defense systems in potato.

Plant hormones are important signal molecules that function in regulating biotic and abiotic stress responses 1,2 . Ethylene (ET) is a stress-responsive hormone in addition to regulating plant growth and development. Infestation with Magnaporthe oryzae activates biosynthesis of ET in rice, and the ET content accumulated in the resistant rice varieties is higher than that in susceptible rice varieties 3 . Cochliobolus miyabeanus inhibits host defense by hijacking the ET signaling pathway, which interferes with the synthesis and/or activation of ET, and significantly enhances rice resistance to brown spot disease 4 . After infection by the super race isolate CN152 of Phytophthora infestans (Pi) in the tetraploid potato genotype SD20, multi-signaling pathways of salicylic acid (SA), jasmonic acid (JA), and ET were involved in resistance and defense against Pi, and one ET synthesis pathway enzyme, 1-aminocyclopropane-1-carboxylate oxidase (ACO), was enhanced in Pi infected potato 5 . In addition, the pathogenesis-related (PR) protein 5 Osmotin gene OSML15 (PGSC0003DMG400003057) was down-regulated in Pi 90128-and CN152-infected potato samples 6 .
Exogenous application of plant hormones can improve plant resistance. ET-induced transgenic rice shows broad-spectrum resistance to fungal pathogens M. oryzae and Rhizoctonia solani 7 . The signaling pathway regulated by ET can induce the expression of defensive genes to enable plants to produce local or systemic defenses 8,9 . The defensive gene StOsmotin2 was specifically induced by ET in potato cultivar Desiree 9 . ET Differentially expressed genes in the response to exogenous ET. The levels of gene expression were estimated by FPKM (Fragments per kilobase of transcript per million mapped reads). Differential expression analysis was performed using DESeq. Genes with a false discovery rate (FDR) < 0.05 and a |log 2 (fold change) |≥ 1 were set as differentially expressed. The differentially expressed genes (DEGs) were detected by performing pairwise comparisons at each time point within a single treatment. In total, 7748 and 5063 DEGs were identified for the H 2 O and ET treatments, respectively (Fig. 1A). Further comparison of DEGs in the two treatments, 3911 genes were specifically differentially expressed between any two time points only in H 2 O treatment, and 1226 In the control, 3911 DEGs were categorized into 498 functional GO terms. In biological processes, 384 GO terms were significantly enriched. The top GO terms were cell cycle (GO:0007049, 128 DEGs), cell cycle process (GO:0022402, 92 DEGs), and cell division (GO:0051301, 76 DEGs) ( Fig. 2A). The most DEGs were enriched in metabolic process (1237 DEGs), followed by cellular process (1217 DEGs) and cellular metabolic process (984 DEGs). In molecular function, a total of 32 GO terms were significantly enriched, including structural constituent of ribosome (GO:0003735, 112 DEGs), structural molecule activity (GO:0005198, 122 DEGs), and microtubule motor activity (GO:0003777, 26 DEGs) ( Fig. 2A). The most DEGs were enriched in nucleic acid binding (382 DEGs), DNA binding (263 DEGs), and protein binding (170 DEGs).
In the ET treatment, the 1226 unique DEGs were specifically (p value < 0.05) annotated into 92 GO terms, including 68 biological process, 22 molecular function, and two cellular component categories. The top 20 GO terms in biological process and molecular function are shown in Fig. 2B. Sixty-eight biological processes were divided into categories including metabolic process, response to stimulus, immunity and defense, regulation of transcription, regulation of biosynthetic process, transport, and others (Fig. 3). Metabolic process was the largest category, containing 19 terms (28%), including cell wall macromolecule (GO:0044036, 12 members), lipid (GO:0006629, 55 members), indole derivative (GO:0042434, five members), and phosphate (GO:0006796, 81 members). The second category was GO terms with response to stimulus (15 terms, 22%), including cellular response to endogenous stimulus (GO:0071495), cellular response to hormone stimulus (GO:0032870), response to auxin stimulus (GO:0009733), and response to ET stimulus (GO:0009723). The third category was related to immunity and defense (eight terms, 12%), including regulation of immune system process (GO:0002682, eight DEGs), regulation of innate immune response (GO:0045088, eight DEGs), regulation of defense response (GO:0031347, 12 DEGs), regulation of plant-type hypersensitive response (GO:0010363, six DEGs), host programmed cell death induced by symbiont (GO:0034050, seven), plant-type hypersensitive response (GO:0009626, seven DEGs), regulation of response to stress (GO:0080134, 12 DEGs), and regulation of immune response (GO:0050776, eight DEGs).
In the molecular function category, the top GO terms were transcription regulator activity (GO:0030528, 62 members), transcription factor activity (GO:0003700, 56 members), protein serine/threonine kinase activity (GO:0004674, 54 members), oxidoreductase activity (GO:0016705, 47 members), and calcium ion binding  www.nature.com/scientificreports/  www.nature.com/scientificreports/ (GO:0005509, 19 members) (Fig. 2B). In the cellular component category, DEGs were enriched in terms of peroxisome (GO:0005777, seven members) and microbody (GO:0042579, seven members) (Fig. 2B), which indicated the enriched genes were mainly located on the oxidases and microbodies. As above, specific DEGs in H 2 O treatment were mainly plant growth and development-related genes, and genes related to disease-resistance and regulation were not identified. ET treatment first initiated the response to stimulus, then mobilized the biological synthesis of various biological macromolecules, and finally regulated the expression of related genes through ion transport and signal transduction, which induced the immune defense process and may be related to disease resistance.

KEGG metabolic pathway analysis reveals DEGs involved in multi-hormone signaling pathways and defense pathways in the response to ET. To better understand the interactions among
DEGs in the potato SD20 response to ET, DEGs were compared to the KEGG database. The 1226 unique DEGs were annotated to 101 KEGG metabolic pathways, ten of which achieved significant differences with p value < 0.05 (Table 3), including plant hormone signal transduction (sot04075, 27 DEGs), plant-pathogen interaction (sot04626, 18 DEGs), and glutathione metabolism (sot00480, ten DEGs) ( Table 3). The biosynthesis of secondary metabolites pathway had the most DEGs at 55. Nine DEGs were involved in the cysteine and methionine metabolism pathway, which was associated with the biosynthesis of ET substances.
A total of ten DEGs were enriched in the glutathione metabolism pathway (sot00480) and up-regulated. Among the ten DEGs, seven genes encoded glutathione S-transferase (GST), which is the key enzyme of the glutathione metabolism pathway, two genes encoded 6-phosphogluconate dehydrogenase (6-PGDH), and one gene encoded isocitrate dehydrogenase [NADP] (IDH), which indicated the important role of the glutathione metabolism pathway in response to ethylene stimulus in SD20.
Exogenous ET induced the differential expression of genes in the endogenous ET synthesis pathway in SD20. Cysteine and methionine metabolism (sot00270) contained nine DEGs, of which three genes encoded proteins closely related to ET biosynthesis: one S-adenosylmethionine synthase 2 (METK2) gene (PGSC0003DMG400037250) and two amino-cyclopropanecarboxylic acid (ACC) oxidase genes (PGSC0003DMG400013894 and PGSC0003DMG400016714), which were all up-regulated (Fig. 5). In addition, in the circulating metabolism of methionine, the ET precursor substance, the homocysteine S-methyltransferase gene (PGSC0003DMG400032243) and methylthioribose kinase (MTRK) gene (PGSC0003DMG400025861) were up-regulated, and the tyrosine aminotransferase-like gene (PGSC0003DMG400009374) was down-regulated (Fig. 5). This indicated that exogenous ET induced differential expression of key genes of the endogenous ET synthesis pathway in SD20.

Differential expression of transcription factors response to ET.
Out of the 1226 ET-specific DEGs, functional annotation showed that 81 genes encoded transcription factors (TFs), which belonged to 13 TF families (Fig. 6A). The largest group was zinc finger protein TFs with 33 members, followed by ET response factor (ERF, 13 members), WRKY (11 members), MYB (11 members), bHLH (nine members), bZIP (five members), and so on (Fig. 6B). Up to 72.8% of 81 TF genes were up-regulated with prolonged ET treatment time, indicating that most TF genes played a positive regulatory role in responding to exogenous ET.
In addition to ERF10 (PGSC0003DMG400010285), which is involved in the plant hormone signaling transduction pathway, there were 12 ERFs including ERF1, ERF3, ERF5, and ERF-like genes (Table 4). Among these, ERF003-like gene PGSC0003DMG400014127, ERF3 gene PGSC0003DMG400022305, and ERF034-like gene PGSC0003DMG400026035 were down-regulated, and the other ten ERF coding genes were all up-regulated. ERF2-like gene PGSC0003DMG400002272 and ERF017 gene PGSC0003DMG400002899 had high fold-changes of 28.3 and 319, respectively, at 12 h of ET treatment compared with 0 h. In addition, PGSC0003DMG400011633 and PGSC0003DMG400016769 encoded StWRKY9 and STWRKY8, respectively, which were enriched in the plant-pathogen interaction pathway with up-regulated expression (Table 4). www.nature.com/scientificreports/ Differential expression of protein kinases response to ET. Out of the 1226 ET-specific DEGs, functional annotation showed that 82 genes encoded protein kinases (PKs) (Fig. 6A), which belonged to ten PK families. The largest group was serine/threonine protein kinase (STKs) with 32 members, followed by leucinerich repeat receptor-like PKs (LRR-RLKs, 15 members), putative receptor protein kinases (PRPKs, four members), lectin-domain receptor protein kinases (three members), mitogen-activated protein kinases (MAPKs, two members), and calcium-dependent protein kinase (CDPK, one member) (Fig. 6B). Among 82 PK genes, 67% were up-regulated and two STK genes (PGSC0003DMG400004928 and PGSC0003DMG400025509) and one LRR-RLK gene (newGene_18308) were significantly up-regulated at 3 h; 37% of PK genes were not up-regulated until 6 h and 26% were not up-regulated until 12 h. The expression of LRR receptor-like serine/threonineprotein kinase FLS2 gene PGSC0003DMG400008296 continued to increase by 2-5 times at 3 h, 6 h, and 12 h of ET treatment. Two Avr9/Cf-9 induced kinase genes PGSC0003DMG400002327 and PGSC0003DMG400027893 were up-regulated and increased 7.8 and 8.7 times at 12 h, respectively.
Differential expression of genes related to oxidoreductase activity and calcium ion signal transduction response to ET. In the GO enrichment analysis, a total of 92 DEGs were enriched in oxidoreductase activity (Fig. 6A), including oxygenase (34 members), oxidase (13 members), reductase (12 members), and dehydrogenase (nine members) (Fig. 6B). Out of 34 oxygenase genes, 23 genes (67.6%) encoded cytochrome P450 monooxygenase (CYP450), which contained CYP71, CYP76, CYP78, CYP88, CYP704,  www.nature.com/scientificreports/ CYP736, CYP705, and CYP801 family members. Among 13 oxidase genes, five genes encoded ACO, the key enzyme in the ET synthesis pathway; four ACO genes were up-regulated. A number of biosynthetic key enzyme DEGs were also identified that encoded gibberellin 20-oxidase 4 and ent-kaurenoic acid oxidase (KAO), the key enzymes in GA synthesis, and protein ECERIFERUM 3/WAX2 (CER3/WAX2), cinnamoyl-CoA reductase (CCR), and 3-dehydroquinate dehydratase/shikimate dehydrogenase isoform 2 (DHD/SHD). The last three enzymes are key enzymes in the biosynthesis of plant epidermal wax, lignin, and aromatic amino acids, respectively, and all three coding genes were up-regulated (Table 5). Pheophorbide a oxygenase (PaO) is the key enzyme in chlorophyll degradation and is encoded by the accelerated cell death 1 gene (ACD1) in Arabidopsis. In this study, PaO coding gene PGSC0003DMG400002967 was up-regulated after exogenous ET treatment. The key tropane alkaloids biosynthetic enzymes littorine mutase and alcohol dehydrogenase (ADH) genes were down-regulated (Table 5).
In GO enrichment analysis, there were 19 DEGs in the enrichment of calcium ion binding, and 13 of them encoded calcium-binding protein, two genes encoded calmodulin, one gene encoded calcium-dependent protein kinase (CDPK), and the other three genes encoded thylakoid lumenal 29.8 kDa protein, annexin p35, and respiratory burst oxidase homolog protein C. Fourteen out of 19 DEGs were up-regulated.
Verification of differential gene expression. In order to verify the accuracy of RNA-seq data, we randomly selected ten DEGs for qRT-PCR verification. The results showed that except for NADPH (down-regulated), the other nine genes were up-regulated: zinc metalloprotease EGY2, NAC domain protein, 1-aminocyclopropane-1-carboxylate oxidase, RAV, WRKY1, and PERK1 kinase coding genes were up-regulated at 3 h of ET treatment; glucan endo-1,3-beta-glucosidase gene was up-regulated at 6 h of ET treatment; while Avr9/Cf-9 rapidly elicited protein and ACC oxidase 2 genes were up-regulated at 12 h of ET treatment. This is consistent with the results of the transcripts determined by RNA-seq analysis (Fig. 7), indicating that our transcriptome analysis data is reliable.

Discussion
Exogenous ET initiates the immune defense response pathway in SD20. In this study, 1226 genes were specifically differentially expressed under ET treatment in tetraploid potato SD20. The GO enrichment analysis showed that 1226 DEGs were enriched in 68 biological process terms, which were mainly related to the metabolism processes of lipid, indole derivative, and phosphate, and responses to exogenous stimulation, immunity, and defense. These results indicated that ET first initiated the response to stimulus, then mobilized the biological synthesis of various biological macromolecules, and finally regulated the expression of related www.nature.com/scientificreports/ genes through ion transport and signal transduction, which obviously induced the immune defense process, and may be related to disease resistance. Further KEGG pathway analysis showed that SD20 was first stimulated by hormones, which initiated plant hormone signal transduction; second, hormone signal transduction initiated expression of immune defenserelated genes, and enhanced the plant's defense capabilities; and finally, a series of metabolic processes that produced glutathione, nitrogen, sulfur, cysteine, and methionine were induced to maintain the normal growth and reproduction of SD20 plants. During this process, biosynthesis of secondary metabolites such as steroids, alkaloids, and plant toxins coordinated and promoted the adaptability of plants under stress stimuli. Secondary  www.nature.com/scientificreports/ metabolites are not necessary for cellular life or plant growth and development, but the accumulation of secondary metabolites is related to plant resistance. In this study, the phosphatidylinositol signaling system (sot04070), which was critical in plant cell responses to environmental stress, was significantly enriched, indicating that ET can induce SD20 to improve resistance to stimulus and stress including pathogen infection.

Exogenous ET stimulates differential expression of TFs, PKs, defense, and oxidoreductase-related genes involved in immunity and defense in SD20.
In this study, genes of TFs, PKs, and a large number of immune defense genes were differentially expressed. It was speculated that exogenous ET was essential in enhancing resistance to stimulus or disease in SD20.
TFs are key factors in plant innate immunity. In this study, TFs DEGs mainly encoded zinc finger protein, ERF, and WRKY, including zinc finger protein Zat10, ERF1, ERF2, ERF5, WRKY6, WRKY11, and WRKY53. Many of them have been reported to function directly in plant defense 13 . The plant tolerance to drought, salt, and oxidative stress is improved by elevating the transcript level of Zat10 and regulating the expression of reactive oxygen removal genes 14,15 . ERF1 is co-activated by JA and ET signaling pathways to regulate the expression of downstream pathogen response genes 16 . ERF2 and ERF5 respond to different pathogens by coordinating several chitinase and other defensive pathways 17,18 . OsWRKY6 and OsWRKY53 are positive regulators for defense responses, and several defense-related genes, including PR genes such as OsPR10a and PBZ1, were activated in OsWRKY6-or OsWRKY53-overexpressed transgenic rice plants, finally resulting in enhanced resistance to Magnaporthe grisea 19,20 . Tobacco WRKY11 participated in plant PAMP-triggered immunity (PTI) and effectortriggered immunity (ETI) processes by interacting with NADPH oxidase RBOHB 21 . As above, the TFs found in this study should play a positive regulatory role in the response to exogenous ET treatment in SD20.
PKs can be used as pattern recognition receptors (PRRs) to identify specific molecules produced by pathogens or damaged plant tissues 22,23 , and are thus involved in signal transduction during defensive responses 24 . PKs found in this study included STKs, RLKs, and MAPKs. The LRR-RLK FLS2 is a widely available PRR in plants and initiates PTI through recognizing the conserved N-terminal 22-amino acid sequence (flg22) of bacterial flagellin 25,26 . FLS2 coding genes were up-regulated in our study indicating their positive regulatory role in SD20 response to ET.
Plant cells can produce numerous defense enzymes including glutathione S-transferases (GSTs) after they are stimulated by a variety of biotic and abiotic stresses 27 . GSTs protect plant cells from damage by detoxifying and antioxidant effects 28,29 . GST genes also respond to infection with fungi and other pathogens, and GSTU13 is an integral part of the immune pathway that produces the defensive indole glucosinolates. In A. thaliana, lack of functional GSTU13 leads to increased susceptibility to several fungal pathogens 30 . In this study, the expression of seven GST genes was induced by ET, which is consistent with a previous study 31 , indicating that GSTs may play a protective role in improving resistance to external stimuli and antioxidants in SD20.
Multiple signaling pathways participate in the response of SD20 to ET. The plant hormone signaling network is an important participant in the plant immune regulation network 32 , which combines antagonism and synergy to form a complex defense system 33 . The JA/ET pathway is primarily responsible for protecting against the infection of necrotrophic pathogens, while SA-mediated defense signaling pathways are activated after infection with biotrophic pathogens 34,35 . Crosstalk of these three hormones, along with other hormones such as ABA, IAA, BR, CTK, and melatonin, are also essential in plant disease prevention 36 . In this study, the top KEGG metabolic pathway was plant hormone signal transduction (sot04075) including IAA, CTK, ABA, ET, JA, and SA. The involved genes encoded hormone synthesis and signal transduction key genes including Aux/IAA, AP2/ERF TF, ethylene receptor, and EIN3/EILs. AP2/ERF TFs are the integration factors of the hormone signaling pathway and are essential in the SA-JA-ET signaling network 37 . EIN3/EILs are the nodes of ET, JA, and SA signaling pathways 38 . This indicated that ET not only activates the ERF1-induced defensive reaction in the ET signal pathway, but also activates the defense reaction of JA and SA signal-mediation in SD20.
The plant defensin 1.2 (PDF1.2) gene is specifically induced by JA or ET and the marker genes for the JA/ET dependent signaling pathway. In this study, one defensin-like protein 1 gene (PGSC0003DMG400032549) was significantly up-regulated. We also found that SA signal marker genes PR-2 (five DEGs), PR-5 (one DEG), and NPR1 (one DEG), JA signal pathway marker gene PR-3 (three DEGs), and JA synthesis key enzyme lipoxygenase (LOX) genes (two DEGs) were differentially expressed. Taken as a whole, ET, SA, and JA signal pathways are involved in the response to exogenous ET treatment and may be related to improved resistance in SD20. This is consistent with the results of multiple signals involved in disease resistance to pathogen infection in SD20 5 , and further supports the important role of ET in mediating resistance and defense in potato SD20. DEGs involved in IAA, GA, ABA, and other signal pathways have different signal transduction and regulation roles in the growth and development of SD20.
Exogenous ET induces internal ethylene synthesis and defensive signal conduction in SD20. ET biosynthesis in plants begins with methionine (Met) and is first catalyzed by S-adenosylmethionine synthetase to produce SAM. Then, ET is formed by oxidation cracking of the important intermediate metabolite ACC. ACS and ACO are the key enzymes in ET biosynthesis. In our study, one SAM synthetase gene METK2 and two ACO genes were differentially up-regulated, which indicated that internal ET was synthesized and accumulated in SD20.
Normally the concentration of ET in the air is very low, and under these circumstances, ethylene receptor ETR1 (ETHYLENE RESPONSE 1) and CTR1 (CONSTITUTIVE TRIPLE RESPONSE 1) are combined together and then combined with phosphate downstream signal component EIN2, then EIN2 in the phosphorylation state is degraded by the ubiquinone/26S protease pathway, inhibiting the transmission of signals downstream.  39 . In this study, three ETR genes and one EIN3 gene had up-regulated expression, which indicated an important positive regulatory role in response to ET-induced signal transduction in SD20 through directly or indirectly activating defensive response, hormone regulation, and ubiquitin-proteasome-mediated degradation.
ET may also be involved in resistance and defense to biotrophic action of pathogens during early infection stages in SD20. SA signaling is generally considered to be involved in resistance to biotrophic pathogens, whereas ET/JA signaling is essential to stop infection by necrotrophic pathogens 40,41 . For example, ET treatment enhanced plant resistance to the necrotrophic fungus Botrytis cinerea 42 and the necrotrophic bacterium Erwinia carotovora 43 . Other studies reported that Arabidopsis EIN2 and EIN3 mutants decreased the susceptibility to biotrophic pathogen Heterodera schachtii and hemibiotrophic pathogen Pseudomonas syringae pv. tomato 44,45 . Overexpression of tomato ERF transcription factors, Pti4, Pti5, and Pti6 enhanced the plant defense to P. syringae in tomato and biotrophic fungus Erysiphe orontii in Arabidopsis 46,47 . These results indicated that the effect of ET on plant disease resistance depends on the interactions of specific plants with pathogens. Phytophthora infestans is a hemibiotrophic pathogen with biotrophic action during early infection and necrotrophic action in the later stages of colonization. A previous study reported that the susceptibility of young Nicotiana benthamiana plants to P. infestans is due to a lack of induction of SA signaling, while the resistance of mature plants against this pathogen requires both SA-regulated appropriate induction of cell death and ET-induced production of phytoalexin 48 . We compared 1226 DEGs induced by ET with the DEGs that were specific to P. infestans race CN152 at 24 h, 48 h and 72 h in potato SD20 5 . The results showed that the DEGs with high foldchange under pathogen infection were also differentially expressed under exogenous ET treatment and had the exact same expression trends, which included the up-regulated PR-2 gene PGSC0003DMG400029830 and fatty acid desaturase gene PGSC0003DMG400036004 at 24 h; up-regulated p-coumaroyl quinate/shikimate 3′-hydroxylase gene PGSC0003DMG400007179 and SlTCP3 gene PGSC0003DMG400015377; and down-regulated cytochrome P450 gene PGSC0003DMG400026523 at 72 h. Also, two LOX genes (PGSC0003DMG400010859 and PGSC0003DMG400024693), and one ACO gene (PGSC0003DMG400013894) were up-regulated at 24 h after CN152 infection in SD20. This indicated that ET signaling may also be involved in resistance to the biotrophic action of P. infestans during early infection stages in SD20.

Conclusion
Here, we performed transcriptome analysis of exogenous ET-treated potato leaves in the high late blight resistant potato genotype SD20 using RNA-seq technology, and the main objectives were to reveal the role of exogenous ET in the defense pathway and the specific responses of SD20 to ET. We found ET stimulated the defense response and initiated a similar defense pathway compared to pathogen infection in SD20. Multiple signaling pathways including ET, SA, and JA were involved in the response to exogenous ET in the host SD20. Our results lay a solid foundation for further understanding of the ET signaling transduction pathway and its mechanisms in disease defense systems in potato. Exogenous ethylene treatment. Three-week-old aboveground seedlings of SD20 were sprayed with 0.2 mmol·L −1 ethephon (C 2 H 6 ClO 3 P, Cat # E8021, Solarbio), an ethylene releasing compound or with an equal amount of sterile water as a control. Each treatment was repeated three times and aboveground seedlings were harvested at 0 h, 3 h, 6 h, and 12 h. A total of 21 samples were frozen in liquid nitrogen immediately for RNA extraction and analysis.

Materials and methods
Library construction and RNA sequencing. These methods were based on the published articles with some modifications 11,49 . Total RNA was extracted using the RNAprep Pure Plant Kit (Tiangen Biotech, Beijing, China) according to the protocol provided and treated with DNase I to generate 21 RNA libraries. RNA integrity was assessed with the 2100 BioAnalyzer instrument using the Agilent RNA 6000 Nano Kit (Agilent Technologies, Santa Clara, CA, USA). One μg of RNA was used to generate a sequencing library with the NEBNext Ultra RNA Library Prep Kit (NEB, USA) according to the manufacturer's recommendations, and the library quality was checked on an Agilent 2100 BioAnalyzer system. Libraries were sequenced using the HiSeq X Ten system (Illumina Inc., USA) with the 150-cycle paired-end sequencing protocol.
Analysis of RNA-Seq datasets. Raw data were filtered and reads containing adapters, poly-N sequence, or low-quality reads were removed to obtain clean reads. All downstream analyses were based on high quality Scientific Reports | (2020) 10:21294 | https://doi.org/10.1038/s41598-020-78027-5 www.nature.com/scientificreports/ clean reads and Q20, Q30, and GC content of clean data were calculated. These clean reads were then aligned to the potato reference genome sequence using HISAT2 software 50 . The mapped reads were assembled using StringTie 51 . Alternative splicing (AS) types and their corresponding expression quantities were generated by ASprofile software 52 . FPKM (Fragments Per Kilobase of transcript per Million fragments mapped) was used as an indicator for measuring the level of transcript or gene expression 52 . 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 (ftp://ftp.ensem blgen omes.org/pub/plant s/relea se-34/fasta /solan um_tuber osum/dna/) 53 . Raw data in FASTQ format are available in the Genome Sequence Archive (GSA) in the BIG DATA Center, Beijing Institute of Genomics (BIG), Chinese Academy of Sciences, under the accession number CRA002320.
Identification of differentially expressed genes. DESeq was used to normalize expression levels of genes and perform differential expression analysis based on the negative binomial distribution 54 . Genes with normalized expression fold-change greater than 2 and false discovery rate (FDR) less than 0.05 were considered to be differentially expressed. Differentially expressed genes (DEGs) were annotated based on the functional annotation information of ENSEMBL release Solanum tuberosum SolTub_3.0 and the potato ortholog Arabidopsis genes.
GO and KEGG enrichment analysis of differentially expressed genes.

Validation of RNA-Seq data by real-time quantitative PCR (qRT-PCR).
To validate the results of the RNA-seq data, ten DEGs were selected randomly for qRT-PCR. The RNA samples were identical to the samples used for RNA-seq analysis. Two μg of total RNA was used per 20 μL reaction for reverse transcription. The PCR system was 20 μL containing 10 μL of SYBR Premix Ex Taq (Takara, Japan), 0.5 μL of forward and reverse primers, 7 μL of double-distilled water, and 2 μL of the cDNA gene. All qRT-PCR reactions were performed with an annealing temperature of 60 °C and a total of 40 amplification cycles with three replicates for each cDNA sample. The expression level of each gene was calculated using the 2 −ΔΔCt method with GADPH as an internal reference gene 56 . Primers (Supplementary Table S1) were designed using PRIMER 5 57 .

Data availability
The datasets generated and/or analyzed during the current study are available from the corresponding author on reasonable request. Raw RNA-seq data from the study have been deposited in the Genome Sequence Archive (GSA) under accession number CRA002320.