Several benchmarks were used to assess the throughput and quality of sequencing reads including the total raw read count for each sample (library), the proportion of clean reads in the raw reads, their base-wise accuracy or sequencing error rate, and GC content. These metrics are summarized in Table 1 .

In terms of overall throughput, 596.43 million raw reads were generated from the four libraries. After filtering, 586.56 million reads were retained as clean reads, which averaged over 40 million per replicate or sample library ( Table 1 ). The mean single base error rate was lower than 1%, and GC content per read averaged 42%–43%. The Q20 scores (percentage of bases whose base call accuracy exceeds 99%) averaged 97%, while the Q30 scores (percentage of bases whose base call accuracy exceeds 99.9%) averaged 92%.

In addition to the overall alignment rates, mapped reads were annotated in terms of the proportions of reads that aligned uniquely or to multiple loci, as well as the proportions of complete reads versus reads composed of spliced (or truncated) RNA. The results of the alignment are presented in Figure 1A .

Across all treatments, the alignments achieved read coverage rates of over 90%, with unmapped reads averaging 4.4%–6.1% of clean reads. Meanwhile, 52%–55% of reads mapped uniquely across treatments, whereas 40%–43% mapped to more than one locus. By comparison, the proportion of spliced reads averaged 20%–22% except for PVY°, which was significantly higher (43%). As a corollary, the fraction of complete reads for PVY° was lower than the 30% observed for PVY^N-Wi and mock treatment, although a similarly lower annotation of complete reads was also observed for PVY^NTN.

The aggregate distribution of mapped reads over exonic, intronic, and intergenic features of the reference genome was another mapping parameter generated (Oshlack et al., 2010). The resulting distribution statistics are shown in Figure 1B . As shown, exons are the highest fraction of mapped reads (88%–91%), followed by intergenic regions (5.7%–6.6%), with 3.4%–4.6% of the reads falling within introns. No significant differences in reads mapping to these genome features were found among the inoculation treatments.

To validate the differential gene expression patterns related to PVY infection as determined by RNA-Seq, the expression levels of selected candidate differentially expressed genes (DEGs) were determined by quantitative reverse transcription-PCR (qRT-PCR). The expression trends obtained by qRT-PCR were consistent with those of mRNA-Seq, further supporting the accuracy and reliability of the sequencing data ( Figure 2 ).

Gene co-expression across the PVY treatments was assessed based on the per-gene read counts [fragments per kilobase of transcript per million mapped reads (FPKM)] obtained from the FeatureCounts program as input data. The resulting pattern is shown in Figure 3 .

As Figure 3 illustrates, 569, 530, and 479 virus-responsive genes were found to be co-expressed between PVY° and PVY^N-Wi, PVY° and PVY^NTN, and PVY^NTN and PVY^N-Wi, respectively. All three PVY strains shared an overall total of 845 co-expressed genes. The number of uniquely expressed genes was the largest for the mock inoculation (1,623), followed by PVY^N-Wi (1,360), PVY° (1,068), and PVY^NTN infection (816). In the context of susceptibility, these overlapping and strain-specific gene sets may be factors in the differences in strain virulence displayed by PVY.

In addition to co-expression analysis on virus-responsive genes, differential analysis with the read count matrix was employed to identify genes whose expression levels were significantly altered by PVY infection. Gene expression levels were considered significant at false discovery rate (FDR) ≤0.05 and log2 ratio ≥1. Figure 4 is a graphic representation of the DEG counts.

Total DEG counts varied among strains, being the highest under PVY^N-Wi treatment (2,778), followed by PVY^NTN (1,744), and the lowest under PVY° infection (1,549) ( Figure 4A ). Across the three PVY-mock comparisons, the number of DEGs under downregulation was greater than that under upregulation. Under PVY° infection, for example, 913 DEGs were downregulated and 636 were upregulated. In the case of PVY^NTN infection, 1,083 DEGs were downregulated and 661 were upregulated. Of the 2,778 DEGs induced by PVY^N-Wi, 1,546 were downregulated and 1,232 were upregulated. On aggregate, the three PVY strains yielded 6,071 DEGs. Among these, 306 overlapped between PVY^N-Wi and PVY^NTN only, with 242 being downregulated and 112 upregulated.

Among the three PVY strains, the total DEG count ( Figure 4B ) was the highest between PVY^NTN and PVY^N-Wi (1,175), followed by PVY° and PVY^N-Wi (794), and PVY° versus PVY^NTN (128). In terms of the direction of differential expression, the number of DEGs downregulated was greater only between PVY^NTN and PVY^N-Wi. Between PVY° and PVY^N-Wi, fewer DEGs were downregulated. In the case of PVY° and PVY^NTN, the number of up- and downregulated DEGs was comparable.

To characterize the biological significance of the DEGs induced by PVY infection, enrichment or over-representation analysis of Gene Ontology (GO) terms was first performed on the DEG sets for each strain. A total of 357 GO terms were annotated across the three PVY treatments, and their distribution is given in Figure 5 .

As shown in Figure 5 , a total of 100 GO terms were enriched for the DEGs related to PVY°. Among these, 74 GO terms were enriched for the downregulated DEGs and comprised 38 biological processes and 36 molecular functions ( Figure 5A ). The remaining 26 GO terms were enriched for the upregulated DEGs and included 16 biological processes, six cellular components, and four molecular functions ( Figure 5B ).

For the DEGs associated with PVY^NTN infection, 85 GO terms were annotated. Of these, 60 GO terms were enriched for the downregulated DEGs and included 26 biological processes and 34 molecular functions ( Figure 5A ). Further, 25 GO terms were enriched for the upregulated DEGs and included 13 biological processes, five cellular components, and seven molecular functions ( Figure 5B ).

Among the three PVY strains, 172 GO terms were enriched in response to PVY^N-Wi inoculation. The majority of these (141 GO terms) were enriched for the downregulated DEGs and included 87 biological processes, four cellular components, and 50 molecular functions ( Figure 5A ). By contrast, 31 GO terms were enriched for the upregulated DEGs and included 12 biological processes, 13 cellular components, and six molecular functions ( Figure 5B ).

Whereas the overall GO term representation differed among the three PVY strains, the analysis also identified a subset of 36 common GO terms ( Table 2 ). Among them, 31 GO terms were enriched for downregulated DEGs, with 15 being biological processes and the other 16 molecular functions.

Of the 15 biological processes, five were related to lipid metabolism, including the lipid biosynthetic process (GO:0008610), steroid biosynthetic process (GO:0006694), steroid metabolic process (GO:0008202), sterol metabolic process (GO:0016125), and sterol biosynthetic process (GO:0016126). Also detected were isoprenoid biosynthetic process (GO:0008299), isoprenoid metabolic process (GO:0006720), terpenoid metabolic process (GO:0006721), and terpenoid biosynthetic process (GO:0016114), all of which participate in terpenoid metabolism. Carbohydrate catabolic process (GO:0016052) and polysaccharide metabolic process (GO:0005976) were shared GO terms related to carbohydrate metabolism. Likewise, organic hydroxy compound biosynthetic process (GO:1901617) and organic hydroxy compound metabolic process (GO:1901615) were the only GO terms related to the metabolism of organic hydroxy compounds. The other shared biological processes were cofactor metabolic process (GO:0051186) and coenzyme metabolic process (GO:0006732), which are necessary for the proper functioning of enzymes.

Of the 16 shared molecular functions, four had transferase activity annotation, including glucosyltransferase activity (GO:0046527) transferase activity, transferring hexosyl groups (GO:0016758) transferase activity, transferring nitrogenous groups (GO:0016769), O-methyltransferase activity (GO:0008171), and NAD+ ADP-ribosyl transferase activity (GO:0003950). A further four were annotated with steroid dehydrogenase activity, including steroid dehydrogenase activity (GO:0016229), 3-beta-hydroxy-delta5-steroid dehydrogenase activity (GO:0003854), and steroid dehydrogenase activity, acting on the CH-OH group of donors, NAD or NADP as acceptor (GO:0033764). Three molecular functions had oxidoreductase activity, including oxidoreductase activity, acting on single donors with incorporation of molecular oxygen, incorporation of two atoms of oxygen (GO:0016702); oxidoreductase activity, acting on the CH-CH group of donors, NAD or NADP as acceptor (GO:0016628); and oxidoreductase activity, acting on single donors with incorporation of molecular oxygen (GO:0016701). The other overlapping molecular functions were intramolecular lyase activity (GO:0016872), dioxygenase activity (GO:0051213), transaminase activity (GO:0008483), xenobiotic transmembrane transporter activity (GO:0042910), and starch binding (GO:2001070).

As shown in Table 2 , the five GO terms that were enriched for upregulated DEGs were all cellular components that included nucleosome (GO:0000786), DNA packaging complex (GO:0044815), protein–DNA complex (GO:0032993), chromatin (GO:0000785), and chromosomal part (GO:0044427).

In contrast to these GO terms that enriched in all three treatments, there were 10 GO terms that overlapped uniquely between PVY^N-Wi and PVY^NTN only ( Table 3 ). Of these, six were molecular functions that all enriched for downregulated DEGs and included 1-deoxy-d-xylulose-5-phosphate synthase activity (GO:0008661); transferase activity, transferring aldehyde or ketonic groups (GO:0016744); lyase activity (GO:0016829); pyridoxal phosphate binding (GO:0030170); vitamin B6 binding (GO:0070279); and vitamin binding (GO:0019842).

The other four GO terms were enriched for upregulated DEGs and comprised three biological processes and one molecular function ( Table 2 ). Two of the three biological processes, trehalose biosynthetic process (GO:0005992) and trehalose metabolic process (GO:0005991), are related to trehalose metabolism. The remaining biological process, disaccharide biosynthetic process (GO:0046351), participates in the formation of disaccharides. The only molecular function, sequence-specific DNA binding (GO:0043565), is related to the regulation of gene expression.

Diverse sets of overrepresented GO terms were also observed among the PVY strains and are presented in Figure 6 . For the DEGs between PVY° and PVY^NTN, the downregulated DEGs yielded no GO term. In the case of upregulated genes, the analysis yielded one significant biological process (copper ion binding, GO:0005507). For the DEGs between PVY° and PVY^N-Wi, no GO term was annotated for the downregulated DEGs. However, 41 GO terms were enriched for the upregulated DEGs, of which 29 terms were biological processes and 12 were molecular functions. For the DEGs between PVY^NTN and PVY^N-Wi, the downregulated DEGs were involved in 36 GO terms, of which 20 terms were biological processes and 16 were molecular functions. However, two significant GO terms were enriched for the upregulated genes, both of which were molecular functions (GO:0016799, hydrolase activity, hydrolyzing N-glycosyl compounds; GO:0019104, DNA N-glycosylase activity).

None of the three PVY treatments resulted in the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment for upregulated DEGs.

By contrast, annotation of the downregulated DEGs yielded a combined total of 22 different significant KEGG pathways. A comparison of these pathways among the PVY treatments is shown in Figure 7 .

As shown in Figure 7 , 14, 14, and 17 virus-responsive pathways were annotated for PVY°, PVY^NTN, and PVY^N-Wi, respectively. Among them, nine pathways were shared by all three PVY strains. This suggests that these pathways are active under PVY-related stress. Of these, four were related to biosynthesis of other secondary metabolites, including flavonoid biosynthesis (KEGGID: sot00941), flavone and flavanol biosynthesis (KEGGID: sot00944), phenylpropanoid biosynthesis (KEGGID: sot00940), and stilbenoid, diarylheptanoid and gingerol biosynthesis (KEGGID: sot00945). Steroid biosynthesis (KEGGID: sot00100) and linoleic acid metabolism (KEGGID: sot00591) were two pathways related to lipid metabolism. The remaining pathways across the three strains that overlapped were zeatin biosynthesis (KEGGID: sot00908), terpenoid backbone biosynthesis (KEGGID: sot00900), and starch and sucrose metabolism (KEGGID: sot00500). Only one pathway, carotenoid biosynthesis (KEGGID: sot00906), was shared between PVY^NTN and PVY°. Likewise, PVY° and PVY^N-Wi shared only the pentose phosphate pathway (KEGGID: sot00030). By contrast, PVY^N-Wi and PVY^NTN shared three pathways: amino sugar and nucleotide sugar metabolism (KEGGID: sot00520), alpha-linolenic acid metabolism (KEGGID: sot00592), and ubiquinone and other terpenoid-quinone biosynthesis (KEGGID: sot00130).

In addition to these overlaps, seven pathways were enriched in a strain-specific manner. Three of these were enriched for PVY°: butanoate metabolism (KEGGID: sot00650) related to carbohydrate metabolism, synthesis and degradation of ketone bodies (KEGGID: sot00072) linked to lipid metabolism, and circadian rhythm—plant (KEGGID: sot04712) linked to environmental adaptation. One carbon pool by folate (KEGGID: sot00670) was enriched exclusively for PVY^NTN, while PVY^N-Wi impacted only cutin, suberine and wax biosynthesis (KEGGID: sot00073), monoterpenoid biosynthesis (KEGGID: sot00902), phenylalanine metabolism (KEGGID: sot00360), and phenylalanine, tyrosine and tryptophan biosynthesis (KEGGID: sot00400).

In addition to the pathways enriched between PVY and mock control, we evaluated pairwise enrichment between the necrotic strains and PVY°. No significant pathways were enriched for DEGs between PVY° and PVY^NTN. By contrast, five KEGG pathways for the DEGs detected between PVY° and PVY^N-Wi were statistically significant. Fructose and mannose metabolism (KEGGID: sot00051) and pentose phosphate pathway (KEGGID: sot00030) were enriched for the downregulated DEGs. Conversely, the plant MAPK signaling pathway (KEGGID: sot04016), amino sugar and nucleotide sugar metabolism (KEGGID: sot00520), and phenylpropanoid biosynthesis (KEGGID: sot00940) were three pathways enriched for the upregulated DEGs.

As a rule, plant–pathogen interactions are associated with quantifiable changes in the host transcriptome, proteome, metabolome, and epigenome (Garcia-Seco et al., 2017). Thus, omics studies have been widely adopted in researching disease pathogenesis, as they can reveal not only the underlying mechanisms of plant resistance but also the strategies used by the pathogen leading to host susceptibility and the subsequent establishment of infection.

In the present study, the response of PVY-susceptible cv. Russet Burbank to infection by strains PVY^NTN, PVY^N-Wi, and PVY° was characterized by mRNA-Seq. Unlike their common progenitor PVY°, both PVY^NTN and PVY^N-Wi are major triggers of tuber necrosis in susceptible potato cultivars. However, the molecular mechanisms inducing necrosis remain poorly understood.

Bioinformatics analysis of the mRNA-Seq data was carried out using the recently published autotetraploid genomes (Hoopes et al., 2022) to explore key aspects of the host–virus interaction, including the differential gene expression outcomes of PVY infection, and their functional significance in inducing host susceptibility and necrosis in the potato host. The throughput and quality of the sequencing reads obtained were within the range considered suitable for high-quality mRNA/whole transcriptome profiling. For example, over 90% of the total clean reads are mapped to the reference genome. The proportion of clean reads obtained, which averaged over 40 million per sample, was also within the 40–60 million range for paired-end reads required for higher-accuracy sequencing projects, including studies of alternate splicing (Conesa et al., 2016).

Whereas over 90% of the reads aligned to the reference genome, unmapped reads averaging 4.4%–6.1% of the clean reads from the four libraries were also detected across all treatments. The significance of these unmapped reads warrants further exploration, as their elucidation could reveal new and meaningful biological information (Laine et al., 2019), including sequences related to the studied PVY strains and genes missing in the reference genome (Isakov et al., 2011; Samuels et al., 2013). For this, de novo assemblies could be generated from the unmapped reads, followed by BLAST alignment of the generated sequence contigs to the National Center for Biotechnology Information (NCBI) non-redundant nucleotide database, which can reveal known matching sequences (Laine et al., 2019).

Although most of the clean reads mapped uniquely to the two genomes, a significant fraction (40%–43%) also mapped to multiple positions. The relatively high multi-mapping rate can be attributed to the autotetraploid nature of the reference genome used (Hoopes et al., 2022), as well as the high rates of gene-level copy number variation generally encountered in tuber-bearing Solanum plants (Hardigan et al., 2017).

A significant fraction (20%–22%) of mapped reads were annotated as spliced chimeric reads. This is not surprising, as eukaryotic genomes generally contain large numbers of introns interspersed among the coding sequences (exons) of their genes. In plants, nearly 60% of such genes are subject to alternative splicing (Filichkin et al., 2010; Mandadi and Scholthof, 2015), where some intronic sequences present in pre-mRNA can be retained in, or removed, from mature mRNA. Consequently, many mRNA variants are produced from one gene (Ben-Dov et al., 2008), which makes alternative splicing a major source of cellular transcriptome and proteome diversity (Pan et al., 2008; Nilsen and Graveley, 2010). In the present study, the proportion of spliced reads obtained from PVY°-treated plants was noticeably higher than observed for PVY^NTN and PVY^N-Wi. Further investigation is required to determine whether these splicing patterns contribute to the added virulence of PVY^NTN and PVY^N-Wi over the PVY° strain.

Other mapping parameters that were generated included the aggregate distribution of mapped reads over exonic, intronic, and intergenic features of the reference genome (Oshlack et al., 2010); this is essential to understanding the gene expression behavior. Further, it indicates the quality of a genome assembly. Exon-mapped reads were the most abundant type of reads, indicating that both reference genomes are well-annotated (Parekh et al., 2018).

The goal of this study was to determine the likely transcriptomic and functional features of PVY susceptibility and tuber necrosis in potatoes. Bioinformatics analyses were carried out to systematically analyze the patterns of gene co-expression, differential gene expression, gene ontology, and KEGG pathway enrichment, from which we infer the following conclusions.

Through co-expression analysis, 845 co-expressed genes induced by the three PVY strains were revealed in this study. As such, these genes may be delineated as a common gene expression signature of PVY infection in this cultivar (Goyer et al., 2015). The results of differential gene expression analysis on the quantification results indicate that the DEG set elicited by PVY^N-Wi inoculation was significantly larger (2,778) than the DEG counts under PVY^NTN (1,744) and PVY° (1,549). On this basis, PVY^N-Wi can be seen as the most aggressive (or virulent) of the three strains. Interestingly, downregulated genes dominated the DEG counts for all three PVY strains. Thus, downregulated gene expression may largely sustain PVY compatibility.

A common feature of such transcriptional modifications is the reported suppression of host defense responses, which are essential for virus replication and symptom development (Whitham et al., 2006). A previous study of the potato–PVY pathosystem (Ross et al., 2022) revealed widespread downregulation during the early response in Russet Burbank after PVY^N-Wi inoculation. Key features of this reported downregulation were genes involved in plant immunity, including plant immune signaling and viral suppression of the antiviral RNAi pathway (Amari et al., 2012).

Among the potyviral proteins, the helper component proteinase (HCPro) is a well-characterized suppressor of antiviral RNAi (Ivanov et al., 2016). Thus, HCPro-mediated suppression of RNAi may be an important molecular explanation for the downregulated gene expression induced by PVY infection of this cultivar. Additionally, viruses can also modify phytohormone levels to attenuate defense or development signaling and change the cellular conditions to favor their replication and spread (Mishra et al., 2020). The differential expression analysis also uncovered 306 genes that achieved DEG significance under both PVY^N-Wi and PVY^NTN, which made them potential candidate genes for tuber necrosis in potatoes.

When enrichment of GO terms of the DEG sets from each PVY-mock comparison were analyzed, the GO term profiles generally mirrored the strain behavior depicted by DEG analysis, where the GO term count elicited by PVY^N-Wi inoculation was also significantly greater (172) than the counts under inoculation with PVY^NTN (85) and PVY° (100). Analysis for similarity and differences between PVY treatments revealed 36 common GO terms across the three strains, of which 31 were enriched for downregulated DEGs and five were enriched for upregulated DEGs. Together, these overlapping GO terms could be considered essential for the susceptibility of potatoes to PVY.

Among the 31 shared GO terms enriched for downregulated DEGs, 15 were biological processes and the other 16 were molecular functions. Five of the 15 biological processes were related to lipid and steroid metabolism, including lipid biosynthetic process, steroid biosynthetic process, steroid metabolic process, sterol metabolic process, and sterol biosynthetic process. In general, plant steroids, including brassinosteroids (BRs) and their precursors (phytosterols), are known to play important roles in growth tolerance to different abiotic and biotic stresses (Divi & Krishna, 2009), and their role in inducing resistance to plant viruses has been reported previously (Yu et al., 2018). Thus, the downregulation of these processes may suggest that BR repression is a key factor in compatibility between potato and PVY.

Four biological processes, isoprenoid biosynthetic process, isoprenoid metabolic process, terpenoid metabolic process, and terpenoid biosynthetic process, participate in terpenoid metabolism. Terpenoids mediate plant responses to biotic and abiotic factors (Tholl, 2006; Nagegowda, 2010), including fine-tuning host–vector–virus interactions to facilitate vector transmission (Wamonje et al., 2020), and the cellular environment for viral replication (Wu et al., 2019). Thus, the observed repression of these terpenoid-related processes could modulate aphid transmission and pathology of PVY in a similar fashion.

On the basis of the foliar transcriptome data, a model for the necrotic behavior of PVY^N-Wi and PVY^NTN is proposed. Based on plant morphology, at least two classes of programmed cell death can be distinguished, including vacuolar cell death and necrosis (Van Doorn et al., 2011). Vacuolar cell death is common during tissue and organ formation and elimination and occurs via the removal of cell contents by a combination of an autophagy-like process and the release of hydrolases from collapsed lytic vacuoles (Van Doorn et al., 2011).

Necrosis occurs mainly under abiotic stress and is distinguished from vacuolar cell death by mitochondrial dysfunction and early loss of protoplast integrity (Van Doorn et al., 2011; Minina et al., 2013). Mitochondrial dysfunction is characterized by the initial permeabilization of the mitochondrial membrane. This implies the formation of pores or channels across the membrane whose inner layer, at physiological homeostasis, acts as a negatively charged system. The resulting membrane polarity is essential for the proper functioning of the respiratory chain and the generation of ATP. Classically, permeabilization of the mitochondrial membrane triggers the rapid influx of Ca²⁺ ions and water into the mitochondrial matrix, resulting in loss of Ca²⁺ homeostasis, dissipation of membrane potential, and swelling and rupture of mitochondria (Scott and Logan, 2008). Oxidative stress is known to contribute to mitochondrial dysfunction, as reactive oxygen species (ROS)-mediated oxidation can cause damage to the mitochondrial DNA, membrane, and respiratory chain (Guo et al., 2013), leading to impaired oxygen consumption and a profound drop in intracellular ATP (Clarke, 1990; Leist and Jäättelä, 2001).

Permeabilization of the plasma membrane causes the early loss of protoplast integrity (Gao & Showalter, 1999; Van Doorn et al., 2011). Cell membrane permeabilization can be triggered by factors such as ROS bursts and associated oxidation of membrane lipids (Yadav et al., 2018), excessive osmotic (turgor) stress (Van Doorn et al., 2011), and pore-forming toxins secreted by pathogens (Brito et al., 2019). A link between cell membrane permeabilization and cell wall integrity is also suggested by the fact that an aberrant cell wall architecture could adversely alter cell wall/plasma membrane interactions and organization, including interruptions in material trafficking, cell wall signaling, and cell wall–plasma membrane cross-linking (Liu et al., 2015).

In the permeabilized state, the increased permeability triggers Ca²⁺ flux into the cytosol (Fink and Cookson, 2005), thereby creating an intracellular ionic imbalance that induces water inflow. Accumulation of Ca²⁺ and water to pathological levels leads to a cascade of excessive protoplast swelling and turgor stress (Trump and Berezesky, 1996), lysis of the plasma membrane (Kroemer et al., 2009) and its retraction from the cell wall, and cell death. Based on these characteristics, a proposed model for the development of tuber necrosis induced by PVY is shown in Figure 6 .

Our data suggest that after entering the cell, all three PVY strains induce downregulation of a mix of pathways of primary and secondary metabolism, which appear to attenuate host defense responses and promote/favor virus survival, propagation, and transmission. As such, these pathways can be considered generic mechanisms of potato susceptibility to PVY infection.

Among them, the common pathway of primary metabolism is starch and sucrose metabolism. In both source and sink tissue cells, sucrose is a raw material for the production of energy and carbon skeletons required for the biosynthesis of amino acids, nucleotides, starch, structural carbohydrates, and defense-related metabolites (Stein and Granot, 2019). Thus, we believe that the observed viral suppression of starch and sucrose metabolism weakens the potato host by inducing starvation, impeding tuber starch accumulation, and decreasing carbon supply for induced defense responses (Engelsdorf et al., 2013).

The suppression of starch and sucrose metabolism also correlates with the eight downregulated pathways of secondary metabolism that, collectively, could limit the scale and effectiveness of any host defense response against PVY. We speculate on the key features of attenuated secondary metabolism that could support or be essential for compatible potato–PVY interaction and inferred these to include disorders in brassinosteroid and cytokinin homeostasis, RNAi suppression, vasiRNA-directed host gene silencing, PVY-induced ROS toxicity, attenuated JA-mediated defense signaling, impaired cell membrane integrity, and disrupted primary carbon flow for secondary metabolism.

Our data suggest that in addition to these generic pathways, PVY^N-Wi and PVY^NTN infections also degrade the activities of another eight pathways that are not affected by PVY° inoculation. Three of these pathways were shared between PVY^N-Wi and PVY^NTN and include amino sugar and nucleotide sugar metabolism, α-linolenic acid metabolism, and ubiquinone and other terpenoid-quinone biosynthesis. The other five pathways were enriched in a strain-specific manner, where 1C metabolism was impacted specifically by PVY^NTN infection, while cutin, suberine and wax biosynthesis, phenylalanine metabolism, phenylalanine, tyrosine, and tryptophan metabolism, and monoterpenoid biosynthesis pathway were impacted uniquely by PVY^N-Wi inoculation. A common consequence of these pathway changes is the generation and supply of substrates for cell wall construction. This supports the conclusion that the necrotic phenotypes induced by both strains of PVY are preceded by initial cell wall permeabilization, which satisfies the first criterion for plant necrosis. Since both strains also cause downregulation of starch and sucrose metabolism, interference with carbon flux for cell wall synthesis could be a key factor promoting its permeabilization.

As indicated earlier, we believe that this abnormal cell wall architecture can also abolish the cell membrane’s selective permeability, resulting in the cascade of excessive Ca²⁺ and water intake, cell swelling, organelle and protoplast deformation, and subsequent rupture and cell death, key events to tuber necrosis. At the same time, viral downregulation of defense-related pathways would attenuate the power of redox defense against PVY-induced ROS stress and induce the mitochondrial dysfunction necessary for cell and tuber necrosis.

Therefore, what distinguishes PVY^N-Wi-induced necrosis mechanistically from that of PVY^NTN? As indicated earlier, 1C metabolism was enriched specifically by PVY^NTN infection, while cutin, suberine and wax biosynthesis, phenylalanine metabolism, phenylalanine, tyrosine, and tryptophan metabolism, and monoterpenoid biosynthesis pathway were impacted uniquely by PVY^N-Wi inoculation. This deviation between PVY^N-Wi and PVY^NTN is intriguing, considering that the three aromatic amino acids (tryptophan, phenylalanine, and tyrosine) and folate, which mediates 1C metabolism (Kaiser & Leistner, 1990), are all synthesized from the same precursor chorismate (Knaggs, 2003). Thus, the primary difference in the necrotic behaviors of these two strains could be related to cutin, suberine and wax metabolism.

Thus, we speculate that the downregulation of this pathway causes suboptimal deposition of these cell wall components, resulting in weakened and permeabilized cell walls. Since the cell wall protects cells from lysing by excessive turgor pressure and osmotic disequilibrium resulting from the net influx of water due to salt imbalance, this cell wall configuration can induce excessive water uptake, cell swelling, and loss of cell wall recalcitrance to large turgor pressure, leading to rupture of cells on a tissue-scale. We speculate that this cell rupture mirrors growth cracking, as it also induces suberin deposition (Woolfson et al., 2022), resulting in the suberized tuber cracking symptom associated with PVY^N-Wi (Benedict et al., 2015).

In this model, the absence of manifestations of foliar necrosis would be related to the criteria for defining the source–sink relationship between aerial and underground organs. In the context of water supply, plants extract water through their roots and transport it to their leaves for photosynthesis. As such, roots could be considered sources of water, while leaves and other aerial tissues represent net sinks for water transported from the root system (White et al., 2016). In addition, tuber tissues with defective cell walls would be more vulnerable to cell swelling and the associated cellular cascade that leads to plant necrosis than leaf cells, as foliar tissues can release excess water from their tissues through transpiration.

The experimental design was described in our previous study (Manasseh et al., 2023). Briefly, the experimental units were single potted potato plants of PVY-susceptible cv. Russet Burbank grown under controlled greenhouse conditions (soil, 21° ± 2°C, with a photoperiod of 16 hours and relative humidity of 70%). PVY strain was the sole treatment factor. To delineate possible differences and commonalities between PVY strains, the strain factor was tested for differences in potato response to inoculation with PVY^NTN, PVYN-Wi, or PVY°. Mock inoculation was performed with sodium phosphate buffer as a control.

To offset the effects of unanticipated micro-environment variations within the greenhouse, the potted plants were arranged in BugDorm-2400 Insect Rearing Tents (MegaView Science. Co. Ltd., Taichung, Taiwan) in a randomized complete block design, with one tent as a block. Four weeks after transplanting, two fully expanded compound leaves selected from the medium plant canopy were inoculated as previously described (Goyer et al., 2015). Each treatment (inoculation) had three replicates, the minimum recommended for RNA-Seq experiments (Conesa et al., 2016). At 7 days post‐inoculation (dpi), two leaflets, one from each inoculated compound leaf, were harvested for total RNA extraction and their subsequent use in Illumina sequencing (mRNA-Seq) and RT-qPCR analysis.

Total RNA was extracted with a modified protocol that sequentially combined the use of TRIzol reagent (Life Technologies, Gaithersburg, MD, USA) with phenol-chloroform extraction. Following extraction with TRIzol, the air-dried RNA pellets were first dissolved in 500 μL of diethylpyrocarbonate (DEPC)-treated water by incubation at 55°C for 10 minutes. After incubation, 500 μL of a mixture containing buffer-saturated phenol (VWR Chemicals, Radnor, PA, USA), chloroform, and isoamyl alcohol (VWR International, USA) in the ratio of 25:24:1 (v/v/v) was added to the solution.

The resulting suspensions were then mixed by inverting each tube several times before centrifuging at 4°C and 12,000 g for 20 minutes. Following centrifugation, ~400 μL of each supernatant was transferred into a fresh 1.5-mL microcentrifuge tube. To each supernatant was added 80 μL of 1 M sodium acetate at pH 5.2 followed by 0.7 times the resulting cumulative volume (supernatant and 1 M sodium acetate) of cold 100% isopropyl alcohol. The tubes were inverted several times to mix the contents followed by incubation at −20°C for 30 minutes.

After incubation, the solutions were centrifuged at 4°C and 12,000 g for 20 minutes. The resulting supernatants were discarded, and 1 mL of 3 M sodium acetate (pH = 5.2) was added, followed by a 5-minute centrifugation at 4°C and 10,000 g. The supernatants were again discarded, and the pellets were washed three times by adding 1 mL of 75% ethanol, followed by gentle vortexing to resuspend the pellets, and then centrifugation at 4°C and 7,500 g for 5 minutes. Following the final wash, the ethanol was decanted, and the pellets were air-dried at room temperature for 10 minutes. As a final step, the dried pellets were dissolved in 30 μL of DEPC-treated water at room temperature for 10 minutes.

Prior to the preparation of the complementary DNA (cDNA) libraries, the RNA samples were DNase-treated (Ambion-Life Technologies, Waltham, MA, USA) to minimize the contribution of sequence reads derived from residual genomic DNA in each sample. RNA quality was then assessed using the Agilent Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA), which assigned an RNA integrity number (RIN) to each sample.

The mRNA was purified from 1.5 μg of total RNA (RIN ≥ 7.0) using Oligo(dT) magnetic beads oligo(dT), according to the manufacturer’s instructions. The cDNA library was prepared from mRNA using the NEBNext^® Ultra™ RNA Library Prep Kit for Illumina^® (New England BioLabs, Inc., Ipswich, MA, USA).

Following purification, divalent cations were used to fragment the mRNA at high temperatures in NEBNext^® First Strand Synthesis Reaction Buffer (5X) (New England BioLabs, Inc., USA). After fragmentation, random hexamer primers and M-MuLV Reverse Transcriptase (RNase H-) were used to synthesize the first cDNA strands, with the single-stranded mRNA as a template.

Subsequently, second-strand cDNA was synthesized with dUTP-containing dNTPs, DNA polymerase I, and RNase H. Double-stranded cDNA products with overhangs converted into blunt-ended cDNA using T4 DNA polymerase and Klenow DNA polymerase. After the end-repair, the cDNA fragments were then monoadenylated on the 3′ ends, to which hairpin looped NEB Next Adaptors were ligated via T-overhangs at their 3′ ends. To select adaptor-ligated cDNA fragments of preferentially 370–420-bp length, the ligation products were purified with the AMPure XP system (Beckman Coulter, Beverly, MA, USA). After purification, the selected cDNA fragments were enriched by PCR performed with Phusion High-Fidelity DNA polymerase, universal PCR primers, and Index (X) primer containing an 8-bp index sequence that allows identification of each library. Finally, amplification products were purified (AMPure XP system) to create the final cDNA library, and library quality was assessed on the Agilent Bioanalyzer 2100 system, following the manufacturer’s instructions.

Index-coded cDNA libraries were clustered on a cBot Cluster Generation System using TruSeq PE Cluster Kit v3-cBot-HS (Illumina) according to the manufacturer’s instructions. After cluster generation, library preparations were sequenced on an Illumina NovaSeq platform, and 150-bp paired-end reads were generated. Library preparation, sequencing, and the subsequent bioinformatics analysis were performed at Novogene Corporation (Sacramento, CA, USA).

Analysis of the RNA-Seq data addressed the following questions: 1) How is the potato transcriptome modified by PVY infection? 2) What are the potential functional consequences of the transcriptome modifications associated with PVY infection? 3) To what extent could these transcriptomic and functional adjustments subserve the necrotic behavior of PVY^NTN and PVY^N-Wi? 4) What unique pathways are associated with necrotic ringspot and tuber cracking?

To answer these questions, the analysis followed the basic steps for a reference-based RNA-Seq data analysis as described by Conesa et al. (2016). These included quality check and preprocessing of raw sequence reads, mapping of the processed reads to the potato reference genome, quantification of gene and transcript levels, assessment of the scope of PVY-induced differential gene expression, and functional annotation of the identified DEGs.

The analysis pipeline was run using the recently published tetraploid genome for Castle Russet as a reference. Castle Russet is resistant to PVY (Quick et al., 2020). The genome assembly and the corresponding annotation files were downloaded from Dryad Data (https://datadryad.org/). In this study, the genome assembly with phased pseudomolecules, unphased pseudomolecules, and unplaced scaffolds (CR_v2.0_asm.fasta.gz) was used. The corresponding set of working gene models (cr.working_models.pm.locus_assign.gff3.gz) was used for annotation.

After sequencing, FastQC (Ver. 0.11.9) was used to examine the sequencing quality of the paired-end raw reads from each sample based on sequence quality, GC content, the presence of adaptors, overrepresented k-mers and duplicated reads, PCR artifacts, or contaminations (Conesa et al., 2016). Following the evaluation of the FastQC results, in-house perl scripts were used to clean the raw reads by removing adapter sequences from the reads and low-quality reads with >10% uncertain (N) base calls, as well as those with more than 50% low-quality (Q-value ≤ 20) bases. All the downstream analysis steps were run on the clean reads.

The reference genome was indexed, and clean reads from each sample were mapped to the indexed reference genome using Hisat2 v2.0.5. For assembly of the mapped reads in each sample into full-length transcripts, the alignments of each file were passed to StringTie (v1.3.3b) (Pertea et al., 2016) with the corresponding reference annotation file, which guided the assembly process. Following the initial assembly, the assembled transcripts were merged by StringTie.

To quantify gene expression, FeatureCounts v1.5.0-p3 (Liao et al., 2014) counted the uniquely mapped read numbers mapped to each gene. For this, FPKM values, which consider the read count and length of each gene (Mortazavi et al., 2008), were computed as a measure of gene expression level (Trapnell et al., 2009). Co-expression analysis was then performed to evaluate the variability in the observed gene expression levels (FPKM) between treatments.

To assess the scope of PVY-induced changes in host gene expression, the per-gene counts were used as input for differential gene expression analysis with the DESeq2 R package (1.20.0). In addition to determining the fold changes between PVY and mock treatments, the default Wald test in DESeq2 was used to test the significance of the observed differential expression. The resulting Wald test p-values were adjusted for multiple hypothesis testing using the Benjamini–Hochberg procedure to control the false discovery rate. Genes with an adjusted p-value ≤0.05, FDR < 0.05, and |log2FC| > 1 were assigned as differentially expressed.

To understand the biological significance of the observed differential gene expression profiles, lists of significant DEGs were subjected to GO (Young et al., 2010) and KEGG pathway enrichment analyses (Kanehisa, 2004). Both analyses were performed using the ClusterProfiler R package while correcting for gene set testing biases toward GO categories containing longer genes. Only the GO terms and KEGG pathways enriched with an adjusted p-value <0.05 (Benjamini–Hochberg correction for multiple testing) were considered significantly enriched by DEGs.

To evaluate the reliability of the RNA-Seq data, two putative DEGs, omega-hydroxypalmitate O-feruloyl transferase (Soltu.Cru.03_4G014040; cutin, suberine and wax biosynthesis) and serine hydroxymethyltransferase (Soltu.Cru.05_3G016520; one-carbon metabolism), were reanalyzed using qRT-PCR. These genes were selected because of their statistical significance of repression or upregulation (p-value <0.05), as well as their potential role in PVY-induced tuber necrosis in potatoes.

Aliquots of total RNAs isolated for RNA-Seq were used for qRT-PCR experiments. For each sample, 1 μg of RNA was first processed to eliminate genomic DNA using a TURBO DNA-free™ kit (Invitrogen, Carlsbad, CA, USA) following the manufacturer’s instructions. First-strand cDNA synthesis was then performed using a commercial cDNA synthesis kit from Sangon Biotech (Shanghai, China). A 1:5 dilution of each cDNA was made in nuclease-free water (Thermo Fisher Scientific, Waltham, MA, USA) and stored at −20°C.

Coding sequences of the selected DEGs were obtained from the Potato Genomics Resource database (http://spuddb.uga.edu). Gene-specific primers for each DEG were then designed using the online PrimerQuest Tool (https://eu.idtdna.com/Primerquest/Home/Index) following default parameters and synthesized by Sigma-Aldrich (St. Louis, MO, USA). Prior to the qPCR assays, the annealing temperature and concentration for each primer were optimized using gradient PCR.

RT-qPCR was performed in 96-well plates with a C1000 Touch Thermal Cycler and CFX Real-Time PCR Detection System (Bio-Rad, Hercules, CA, USA). Each reaction (final volume of 10 µL) contained 1.0 μL of cDNA template, 5.0 μL of SsoAdvanced Universal SYBR Green Supermix (Bio-Rad, USA), 0.5 μL of each primer (10 μM), and 3.0 μL of nuclease-free water. The reaction profile had an initial pre-denaturation at 94°C for 3 minutes, followed by 30 cycles of denaturation at 94°C for 10 seconds, and a combined annealing and extension step at 58°C for 30 seconds. Two technical replicates of each of the three biological replicates along with a no template control were included in the qPCRs for each DEG.

The initial cycle threshold (Ct) or quantification cycle (Cq) data generated were then exported to an Excel file and used as input to compute 2^−Δδct values (or fold change estimates) in the expression of each DEG (Livak and Schmittgen, 2001), normalized to the eukaryotic elongation factor 1α (eEF1A), as previously used in similar experimental settings (Nicot et al., 2005; Goyer et al., 2015; Tang et al., 2017). Finally, the 2^−ΔΔCt data were then log-transformed and subjected to t-tests for statistical significance of the observed fold changes at a threshold of p < 0.05.