Investigation of P1/HC-Pro-Mediated ABA/Calcium Signaling Responses via Gene Silencing through High- and Low-Throughput RNA-seq Approaches

The P1/HC-Pro viral suppressor of potyvirus suppresses posttranscriptional gene silencing (PTGS). The fusion protein of P1/HC-Pro can be cleaved into P1 and HC-Pro through the P1 self-cleavage activity, and P1 is necessary and sufficient to enhance PTGS suppression of HC-Pro. To address the modulation of gene regulatory relationships induced by turnip mosaic virus (TuMV) P1/HC-Pro (P1/HC-ProTu), a comparative transcriptome analysis of three types of transgenic plants (P1Tu, HC-ProTu, and P1/HC-ProTu) were conducted using both high-throughput (HTP) and low-throughput (LTP) RNA-Seq strategies. The results showed that P1/HC-ProTu disturbed the endogenous abscisic acid (ABA) accumulation and genes in the signaling pathway. Additionally, the integrated responses of stress-related genes, in particular to drought stress, cold stress, senescence, and stomatal dynamics, altered the expressions by the ABA/calcium signaling. Crosstalk among the ABA, jasmonic acid, and salicylic acid pathways might simultaneously modulate the stress responses triggered by P1/HC-ProTu. Furthermore, the LTP network analysis revealed crucial genes in common with those identified by the HTP network in this study, demonstrating the effectiveness of the miniaturization of the HTP profile. Overall, our findings indicate that P1/HC-ProTu-mediated suppression in RNA silencing altered the ABA/calcium signaling and a wide range of stress responses.


Introduction
P1/HC-Pro is the first identified viral suppressor of potyvirus and can trigger the suppression of RNA silencing in the microRNA (miRNA) and short-interfering RNA (siRNA) regulatory pathways [1][2][3]. Our previous studies indicate that the FRNK motif of HC-Pro plays an essential role in the suppression of the miRNA pathway but still suppresses 40% of the siRNA pathway [2]. Moreover, Hu et al. (2020) demonstrated that various potyviral species of P1/HC-Pro, i.e., turnip mosaic virus (TuMV), zucchini yellow mosaic virus (ZYMV), and tobacco etch virus (TEV), have the same function in RNA silencing suppression [1]. However, the P1/HC-Pro of TuMV (P1/HC-Pro Tu ) triggers ARGONAUTE1 (AGO1) degradation, whereas those of ZYMV (P1/HC-Pro Zy ) and TEV (P1/HC-Pro Te ) do not cause AGO1 degradation, which suggests that viral P1/HC-Pros exhibit functional diversity. Moreover, Sanobar et al. (2021) demonstrated that HC-Pro Tu at the High Throughput Genomics Core of Academia Sinica. The clean reads were trimmed to filter out the adapters and the low-quality reads by using CLC Genomic Workbench (Qiagen, Hilden, Germany). The sequences of the coding regions in the whole Arabidopsis transcriptome were then mapped, and the reads were counted by using Bowtie 2 version 2.2.5 [7] and eXpress version 1.1.5 [8]. The raw RNA-Seq data of the HTP and LTP profiles are available at the National Center for Biotechnology Information Short Reads Archive (NCBI SRA) accession number SRR16916514-SRR16916533 and SRR16916443-SRR16916454 at the following URL: https://www.ncbi.nlm.nih.gov/Traces/study/?acc=PRJNA779609 and https://www.ncbi.nlm.nih.gov/Traces/study/?acc=PRJNA779601 respectively (accessed on 11 November 2021).

Differential Gene Expression Analysis and Functional Annotation
The transcriptome was analyzed using the ContigViews system (www.contigviews. bioagri.ntu.edu.tw, accessed on 11 November 2021) of the NGS core of National Taiwan University. The transcript abundances were based on read counts normalized to FPKM (fragments per kilobase per million). For the HTP expressional network analysis, the differentially expressed genes (DEGs) between the comparative Col-0 and P1/HC-Pro Turelated datasets were identified based on an 80% passing rate, and genes with twofold log 10 FPKM values less than 1.14 were filtered out. Due to the difference in the sequencing coverages achieved with HTP and LTP, different parameters were used in the ContigViews analysis. The LTP datasets were analyzed using appropriately adjusted parameter settings for the correlation thresholds during the construction of the expressional correlational networks. At least 10 samples from the Col-0, P1 Tu , HC-Pro Tu , and P1/HC-Pro Tu profiles were selected to calculate the Pearson correlation based on a 0.975 threshold for a positive correlation and a 0.925 threshold for a negative correlation. For the LTP expressional network analysis, DEGs were identified using an 80% passing rate and a fold-change of 2, and the LTP networks were then constructed using thresholds of 0.95 and 0.90 for positive and negative correlations, respectively.

Quantification of Endogenous ABA and ABA Sensitivity Assay
Ten-day-old Arabidopsis seedlings (20-25 seedlings, approximately 60 mg for each extraction) were freshly collected and ground into fine powder with a tissue grinder pestle in a tube with liquid nitrogen. Fifty microliters of working solution (methanol) containing 0.5 ng of d 6 -ABA were added as a standard to each tube. After the addition of 500 µL of extraction solvent (2-propanol/H 2 O/concentrated HCl, 2:1:0.002), the tube was shaken at 100 r.p.m. and 4 • C for 30 min. Subsequently, 1 mL of dichloromethane was added to each tube, and the tube was shaken for another 30 min and centrifuged at 13,000× g and 4 • C for 5 min. Nine hundred microliters of the lower phase were transferred to a new tube, and the solvent was concentrated using a rotary evaporator (EYELA CVE-3110). An ultra-performance LC/ESI-qMS/MS analysis was conducted by the Metabolomic Core Facility at the Agricultural Biotechnology Research Center, Academia Sinica, Taiwan. A HALO C18 (Advanced Materials Technology, Inc., Wilmington, DE, USA) column (inner diameter, 2.1 mm; column length, 75 mm; particle size, 2.7 µm) was used, and gradient elution was performed with water and 0.05% glacial acetic acid (solvent A) and acetonitrile with 0.05% glacial acetic acid (solvent B) at a constant flow rate of 0.6 mL min −1 . The following gradient profile was applied: t (min), % A): (0, 99), (2.20, 0), (2.40, 0), (2.60, 99), (3,99). The MS and MS/MS experiments were performed with an API 3000 triple quadrupole mass spectrometer (PE Sciex, Concord, Ont., Canada) with the following parameters: temperature of 400 • C, nebulizer gas (N 2 ) 10 (arbitrary units), curtain gas (N 2 ) 12 (arbitrary units), collision gas (N 2 ) 4 (arbitrary units), and the capillary voltage of −3.5 kV. The mass spectrometer was operated in multiple reaction mode (MRM). To germinate seeds for the ABA sensitivity assay, Arabidopsis seeds were sterilized and spread on the MS medium with or without 0.11 µM ABA, and the growth conditions were observed at 7 days after germination.

Expression-Based Heatmaps and Principal Component Analysis (PCA)
The identified DEGs were functionally annotated based on their sequence similarities with known protein annotations in the public TAIR database (www.arabidopsis.org, accessed on 11 November 2021). Heatmaps of DEG expressions were generated with the ClustVis web tool [10]. PCA was performed using the PCA method of the decomposition module in the scikit-learn package of Python. The principal component (PC) scores were plotted with the Matplotlib package. The PCA analysis was performed using the FPKM values for all the transcripts obtained from individual replicates. Thus, the PCA-based comparison of gene expression was performed using 24 libraries of the 4 Col-o, P1 Tu , HC-Pro Tu , and P1/HC-Pro Tu samples from the HTP and LTP datasets.

P1/HC-Pro Tu Suppressor Triggers Plant Defense Responses
The comparative analyses of the HTP RNA-Seq profiles of Col-0 vs. P1/HC-Pro Tu , Col-0 vs P1 Tu , and Col-0 vs. HC-Pro Tu sets revealed 1601, 559, and 777 DEGs, respectively (Table 1). These DEGs were then used for a network analysis using the ContigViews system, which revealed 662, 106, and 162 genes in the networks, respectively (Table 1). A Venn diagram showed that 29 of the network genes were found in all three comparative sets ( Figure 1A). The P1/HC-Pro Tu -only section contained a markedly higher number of network genes (553), whereas the P1 Tu -only and HC-Pro Tu -only sections contained only 18 and 24 genes, respectively ( Figure 1A).
To elucidate the functions of the unique genes in the P1/HC-Pro Tu -only, P1 Tu -only, and HC-Pro Tu -only sections further, we performed gene annotation and functional classification. Unique genes with functions related to stress-responses, cell growth, and plant development were predominantly enriched in the P1/HC-Pro Tu -only section in addition to unknown or unclassified functional genes ( Figure 1B). Similar gene functional classification results were found for the P1 Tu -only and HC-Pro Tu -only sections, although only 18 and 24 genes were identified in these two sections, respectively ( Figure 1A). A high abundance of network genes in the P1/HC-Pro Tu -only section were classified as being involved in ABA phytohormone and calcium signaling pathways ( Figure 1A). Table 1. Identification of the HTP and LTP network genes in the Col-0 vs. P1/HC-Pro Tu , Col-0 vs. P1 Tu , and Col-0 vs. HC-Pro Tu comparative datasets in the ContigViews system.

P1/HC-Pro Tu Alters ABA-Induced Immune Responses
ABA plays pivotal roles in seed germination and stress tolerance [11]. Several studies demonstrate that abiotic stresses increase ABA and Ca 2+ concentrations in the cytosol and induce signaling to modulate common target proteins (Figure 2A) [12]. We found 41 genes belonging to the ABA signaling pathway that could be further classified into the ABA homeostasis, ABA signaling regulation, and ABA response categories ( Figure 2A and Table 2). Among these genes, PMI1/MEE31 (AT3G02570) and SULTR3;1 (AT3G51895) that belong to the ABA positive biosynthesis regulator were respectively repressed and induced in the P1/HC-Pro Tu plants (Figure 2A(panel i),B). The expression of both CYP707A3 (AT5G45340), involved in ABA catabolism, and DTX50 (AT5G52050), involved in ABA transport, was induced in the P1/HC-Pro Tu plants (Figure 2A(panel ii and iii),B). As for the ABA signaling regulation, RCAR/PYLs, PP2C, and SRNK2, the core components of ABA signaling transduction, were not differentially expressed in our study, even though several negative regulators (e.g., ATL27/ATARRE (AT5G66070), PUB22 (AT3G52450), RDUF1 (AT3G46620), RDUF2 (AT5G59550), and HB12 (AT3G61890)), which inhibit RCAR/PYL and PP2C expressions either transcriptionally or posttranscriptionally, had an induced expression pattern in the P1/HC-Pro Tu plants (    ABA response genes accounted for the highest proportion of the genes involved in the ABA signaling pathway found in the P1/HC-Pro Tu -only section and could be divided into the biotic stress response, development, drought stress response, cold stress response, and senescence subcategories (Figure 2A(panel vi and vii) and Table 2). Almost all of these ABA response genes were expressed at a higher level in the P1/HC-Pro Tu plants than in Col-0, P1 Tu , and HC-Pro Tu plants ( Figure 2B). Induced expressions of these ABA response genes might change plant resistance to drought stress, cold stress, and leaf senescence in response to P1/HC-Pro Tu . These results revealed that the regulation of the ABA signaling pathway was disrupted, and its responses might be severely interfered with overexpressing P1/HC-Pro Tu .

Quantification of Endogenous ABA and ABA Sensitivity Assay
To examine ABA accumulation in the P1/HC-Pro Tu plants, the 10-day-old seedlings were extracted and measured by MS/MS analysis. A significantly lower ABA amount was detected in the P1/HC-Pro Tu seedlings than in the Col-0 ( Figure 3A). To test the effects of ABA on the P1/HC-Pro Tu seedlings further, an ABA sensitivity assay was carried out to observe the phenotypical changes with seed germination. P1/HC-Pro Tu plays a major role in PTGS suppression by triggering AGO1 degradation [1]. Therefore, ago1-27 mutant was also used for the ABA sensitivity assay. Without ABA treatment, the seeds of Col-0 were germinated and developed into true-leaf seedlings at phase IV, while seeds of the P1/HC-Pro Tu plants and ago1-27 mutant were more late-germinated or delayed-growth at phase I, II, and III ( Figure 3B), suggesting that the germination rate and post-germination growth were greatly delayed in the P1/HC-Pro Tu plants and ago1-27 mutant. With exogenous ABA treatment, the delayed-germination phenotype became more severe in the P1/HC-Pro Tu plants and ago1-27 mutant than in the Col-0 ( Figure 3B). Specifically, more than half of the P1/HC-Pro Tu seeds remained at phases I and II, which indicated that the P1/HC-Pro Tu plants exhibited a higher sensitivity to ABA during seed germination. The mean values ± SD were obtained from three biological repeats. Comparisons between two groups were performed using a Student's t test. *** p < 0.001. (B) ABA sensitivity assay. The germination phenotypes were classified into four phases: I, rupture of the seed coat; II, radicle protrusion; III, fully-opened cotyledons; and IV, true leaf development.

P1/HC-Pro Tu Triggers Immune Responses in a Calcium-Dependent Manner
We found that the genes in the P1/HC-Pro Tu -only section encode various Ca 2+ transporters, including Ca 2+ channels [CMCU (AT5G66650), CNGC2 (AT5G15410), and CNGC14 (AT2G24610)], Ca 2+ co-transporter [CCX2 (AT5G17850)], and Ca 2+ pump [ACA1 (AT1G27770)] ( Figure 4A(panel i) and Table 3). We also identified several Ca 2+ sensors, including Ca 2+ /calmodulin (CaM), CaM-like (CaML) ( Figure 4A(panel ii)), calcium binding proteins ( Figure 4A(panel iii)), calmodulin-binding proteins ( Figure 4A(panel iv)), and CDPKs/ CIPKs ( Figure 4A(panel v) and Table 3). The members of the various families of Ca 2+ sensors identified in the P1/HC-Pro Tu -only section were expected to contribute to the con-version of Ca 2+ signals into either cellular stress responses or developmental processes ( Figure 4A(panel vi and vii)). For instance, the expressions of genes encoded calmodulinbinding proteins, including CAMBP25 (AT2G41010) and IQM4 (AT2G26190), whose gene expressions are induced by drought and salt/osmotic stress. CPK32 (AT3G57530) and CPK28 (AT5G66210) encode calcium-dependent protein kinases that regulate plant growth in addition to resetting pathogen-associated molecular pattern (PAMP)-induced defense signaling. The majority of genes involved in the Ca 2+ signaling pathway were expressed at higher levels in the P1/HC-Pro Tu plants than in the Col-0, P1 Tu , and HC-Pro Tu plants ( Figure 4B) except for CNGC2 (AT5G15410), which suggests that the induction of the calcium signaling pathway depends on ectopic-expressing P1/HC-Pro Tu . In summary, the results indicate that P1/HC-Pro Tu might trigger various stress responses and developmental processes through the calcium signaling pathway.

Validation of DEGs in the ABA and Ca 2+ Pathways
The expression patterns of several candidate genes in the ABA and Ca 2+ signaling pathways were validated by qPCR. In addition to the Col-0 and P1/HC-Pro Tu plants, P1/HC-Pro Zy plants were also included in gene expressional validation. Compared with the ovate leaves of the Col-0 plants, the serrated-leaf phenotype was observed in the P1/HC-Pro Tu and P1/HC-Pro Zy seedlings, which implies common PTGS suppression effects induced by the heterogeneous P1/HC-Pros. ( Figure 5). Four ABA response genes (ABF4, MYB44, MYB96, and OZF1) exhibited similar expression patterns with higher transcripts in the P1/HC-Pro Tu plants, consistent with the HTP RNA-Seq profiles. Moreover, the expression levels of the selected genes were also highly induced in the P1/HC-Pro Zy plants ( Figure 6A-D). Additionally, CaLB, IQM4, and CPK28 involved in the Ca 2+ signaling pathway showed similar qRT-PCR results to those of the ABA response genes ( Figure 6E-G). The qRT-PCR results indicate that the overall expressions of genes in the ABA and Ca 2+ signaling pathways were consistent in the P1/HC-Pro Tu and P1/HC-Pro Zy plants.  The mean values ± SD were from three biological repeats. Comparisons between two groups were performed with Student's t test. * p < 0.05, ** p < 0.01, *** p < 0.001.

P1/HC-Pro Tu Triggers Drought Response and Stomatal Closure
Our data mining of the P1/HC-Pro Tu -only section revealed 66 drought stress-related genes, e.g., MYB96 (AT5G62470), MYB44 (AT5G67300), NF-YA5 (AT1G54160), ANAC029 (AT1G69490), and TZF1 (AT2G25900), that were highly enriched in terms annotated to drought stress responses (Table 4). Among them, 21 and 4 genes were found to be involved in either the ABA or Ca 2+ signaling pathway, respectively (Table 4). Notably, four regulatory modules composed of 15 drought response genes that function in controlling stomatal guard cell dynamics were identified ( Figure 7A and Table 4). Genes in these regulatory modules could induce ABA/Ca 2+ -mediated stomatal closure, salicylic acid (SA)-or jasmonic acid (JA)-mediated stomatal opening, starch degradation-mediated rapid stomatal reopening, and guard cell division to influence stomatal development ( Figure 7A). The mechanisms through which these genes are controlled and the involvement of phytohormones and environmental stimuli that are involved in the regulation of stomatal dynamics are explained below.
In the ABA/Ca 2+ -mediated stomatal closure modulation ( Figure 7A(panel i)), genes such as WRKY46 (AT2G46400) and RZPF34/CHYR1 (AT5G22920) function in the response to ABA and to water deficit stress ( Figure 5A). However, several genes, such as JAZ2, MYC2, and ANAC019 (AT1G74950, AT1G32640, and AT1G52890), modulate stomatal reopening after microbe-associated molecular pattern (MAMP)-mediated stomatal closure and activate the JA pathway upon bacterial infection ( Figure 7A(panel ii)). Similarly, stomatal reopening is modulated by either the SA-mediated pathogen infection signaling pathway ( Figure 7A(panel ii)) or light-induced starch degradation by the glucan hydrolase β-AMYLASE1, which is encoded by BAM1 (AT3G23920), to promote rapid stomatal opening ( Figure 7A(panel iii)). All the DEGs in the three modulations showed greater induced expression levels in the P1/HC-Pro Tu plants than in the Col-0, P1 Tu , and HC-Pro Tu plants ( Figure 7B).
Moreover, the stress responses induced by P1/HC-Pro Tu could affect the regulatory mechanisms of stomatal development and lead to stable long-term adaptations to stress ( Figure 7A(panel iv)). Induced expression of the negative regulator EPF2 (AT1G34245) and the positive regulator TMM (AT1G80080) might affect the asymmetric divisions of guard mother cells ( Figure 7B). An accumulation of AGL16 (AT3G57230) with silent mutations (AGL16m) in the miR824 recognition site reportedly promotes the development of higherorder stomatal complexes by increasing the number of additional divisions in meristemoid cells [13]. With the exceptions of EPE2 and TMM, all of these genes were expressed at higher levels in the P1/HC-Pro Tu plants than in the Col-0, P1 Tu , and HC-Pro Tu plants ( Figure 7B). Overall, the results indicate that the integration of ABA, calcium, and other hormone signals could simultaneously trigger dynamic closure and opening mechanisms during drought stress or biotic stress in the P1/HC-Pro Tu plants.

P1/HC-Pro Tu Stimulates Cold Response and Leaf Senescence
Twenty cold response genes were identified in the P1/HC-Pro Tu -only section obtained from of the HTP profiles ( Table 5). All of these genes showed higher expression levels in the P1/HC-Pro Tu plants that were different from those found in the Col-0, P1 Tu , and HC-Pro Tu plants ( Figure 8A). Because a transient increase in cytosolic Ca 2+ levels is detected within seconds of cold shock [14], several annotated cold-related genes in our dataset, including CBF2 (AT4G25470), CEJ1/DEAR1 (AT3G50260), ZAT10 (AT1G27730), and ZAT12 (AT5G59820), were found to be involved in the cold stress-induced Ca 2+ signature (Table 5). Moreover, seven of the cold-response genes were involved in the drought stress response (Table 5), which indicates that the gain of P1/HC-Pro Tu function may suggest overlaps between the gene regulatory mechanisms underlying enhancement to cold and drought tolerance. However, there were still many cold-response genes that were independent of the ABA/Ca 2+ signaling pathways, suggesting that P1/HC-Pro Tu could stimulate a cold response via other regulating mechanisms.
Although leaf senescence is a general developmental program, it can be induced by abiotic and biotic stresses through the transcriptional control of senescence-related transcription factors [15]. Twenty-five genes annotated as being involved in leaf senescence, cell death, and chlorophyll degradation were identified in the P1/HC-Pro Tu -only section ( Table 6). Among these genes, the WRKY22 (AT4G01250) and ERF014 (AT1G44830) transcription factor-encoding genes function in mediating dark-induced senescence [16,17]. Chlorophyll degradation is one of the most visually apparent phenomena that occurs during leaf senescence [18,19]. The leaf degreening process appears to be promoted by genes involved in the ABA signaling pathway [ABF4 (AT4G24390)] as well as in the JA [ANAC019 (AT1G52890)] and GA [ANAC029 (AT1G69490)] signaling pathways. Developmentally controlled programmed cell death occurs at the terminal stage of leaf senescence [18] The identified NUDT7 (AT4G12720) and CAD1 (AT1G29690) genes might play distinct roles in plant immunity related to the regulation of programmed cell death [20,21]. With the exceptions of APD7 (AT5G02760) and AT2G44670, all of these genes were expressed at higher levels in the P1/HC-Pro Tu plants than in the Col-0, P1 Tu and HC-Pro Tu plants ( Figure 8B). Taken together, our results suggest that cold tolerance, chloroplast degradation, and leaf senescence are likely to serve as key responses in the P1/HC-Pro Tu plants.

Comparison of the HTP and LTP Profiles
We then compared the data mining efficiency of the LTP profile (~2 M reads) with that of the HTP profile (~20 M reads). A total of 21.94 to 24.66 M clean reads were generated from all the samples for the HTP profiles, and these exhibited an average mapping rate of 79.2% to the coding DNA sequence (CDS) ( Table 7). Approximately one-tenth of the total sequencing depth was used to construct the LTP profiles; thus, the LTP profiles contained 1.89 to 2.96 M clean reads obtained from 1.97 to 3.10 M raw reads ( Table 7). The average mapping rate of the LTP profiles was 78.5%, which was close to that found for the HTP profiles ( Table 7). The equivalent mapping rates obtained for the HTP and LTP libraries indicate that the mapping ability of the RNA-Seq reads does not depend on the RNA-Seq depth.
We then performed a PCA of the HTP and LTP profile data ( Figure 9A). The top two PCs explained 75.7% of all differences among the three varieties, and PC1 accounted for 63.0%, which suggested that PC1 can distinguish between the HTP and LTP profiles ( Figure 9A). We also noted that biological replicates of the HTP profiles were more consistent than those of the LTP profiles ( Figure 9A). Furthermore, the PCA clustering of the HTP data corresponded to the morphological phenotypes: the Col-0 and P1 Tu plants had identical normal developmental phenotypes, whereas the HC-Pro Tu and P1/HC-Pro Tu plants had a serrated leaf phenotype. In contrast, PC2 explained only 12.7% of the overall differences but was most likely to distinguish the P1/HC-Pro Tu samples from the other samples ( Figure 9A). In addition, based on PC2, the clustering of the P1/HC-Pro Tu samples distinctly differed from that of the other samples, and this finding was obtained for both the HTP and LTP profiles.
We compared the Col-0 vs. P1/HC-Pro Tu plant samples, and the results revealed 75 common genes, which were shown in the intersection area of the networks obtained with the HTP and that obtained with the LTP profiles ( Figure 9B and Table 8). These genes were characterized as being involved in ABA/Ca 2+ signaling pathways, drought or cold stress responses, senescence, and gene silencing and RNA regulation (Table 8). We also found that the 75 common genes were located at identical positions in the HTP and LTP networks for comparison ( Figure 9C,D). Moreover, the HTP and LTP profilebased networks of the 75 common genes revealed 132 and 159 gene-gene correlations for the HTP and LTP profiles, respectively ( Figure 9C,D). However, we observed that connections associated with the positive and negative correlations were not 100% identical between the HTP and LTP profiles ( Figure 9C,D). Twenty-six correlations (19.7%), including 25 positive connections and one negative connection, among the 30 common genes in the HTP network remained conserved in the LTP network. Furthermore, the heatmaps of the 75 common genes in the HTP and LTP profiles exhibited similar expression patterns, and the expressions of these genes were upregulated in the P1/HC-Pro Tu plants ( Figure 10).

Functional Classification of the DEGs Identified from the LTP Profile
A total of 188, 234, and 193 network genes were identified in the Col-0 vs. P1/HC-Pro Tu , Col-0 vs. P1 Tu , and Col-0 vs. HC-Pro Tu LTP comparison sets, respectively, whereas the corresponding HTP comparison sets contained 553, 18, and 24 network genes, respectively ( Table 1). The LTP dataset revealed equivalent gene numbers among the three comparison sets, whereas the HTP dataset showed a higher abundance of network genes in the Col-0 vs. P1/HC-Pro Tu comparison. A Venn diagram was generated to determine the unique and shared genes among the Col-0 vs. P1/HC-Pro Tu , Col-0 vs. P1 Tu , and Col-0 vs. HC-Pro Tu comparison sets. Sixty-nine shared network genes were identified from the three comparison sets of the LTP profiles ( Figure 11A). Adequate gene numbers were also obtained in the P1/HC-Pro Tu -only (96 genes), P1 Tu -only (121 genes), and HC-Pro Tu -only (79 genes) sections ( Figure 11A). Moreover, functional characterization revealed that genes involved in stress responses, plant development processes, and the calcium signaling pathway were abundant in the P1/HC-Pro Tu -only section obtained with the LTP profiles, which are similar to the results obtained from the functional characterization of genes in the P1/HC-Pro Tu -only section based on the HTP profiles ( Figures 1B and 11B). Notably, the P1 Tu -only and HC-Pro Tu -only sections obtained from the HTP and LTP profiles were not significantly identical (Figures 1B and 11B).

P1/HC-Pro Tu Alters ABA and the Other Hormones Accumulations
Multiple plant hormones are reported to respond to P1/HC-Pros [1,5]. Endogenous ethylene is maintained at a higher level in the P1/HC-Pro Tu plants, and the comparative network of Col-0 vs. P1/HC-Pro Tu also highlighted critical genes in various hormone signalings (e.g., JA, ethylene, and ABA) [1]. Hu et al. (2020) also proposed that the serrated leaf phenotype of the P1/HC-Pro Tu plants might relate to the endogenous auxin accumulation [1]. These studies implied a comprehensive alternation among different hormone pathways that occurred in response to P1/HC-Pros. Consequently, the coordinated modulations or crosstalk of hormone responses could possibly be interfered by P1/HC-Pros and cause changes in growth and immunity responses.
In this study, the ABA signaling pathway was fundamentally changed in the P1/HC-Pro Tu plants. For example, P1/HC-Pro Tu triggered the ABA negative regulator up-regulation and interfered with ABA positive regulator expressions for the ABA homeostasis and signaling regulation, resulting in low abundant endogenous ABA in the P1/HC-Pro Tu plants. Surprisingly, the ABA response genes were mostly induced in the P1/HC-Pro Tu plants (Figure 2), implying that the PTGS suppression might alter these gene expressions.
Indeed, the endogenous AGO1 was degraded in the P1/HC-Pro Tu plants [1], which also showed ABA-sensitivity in seed germination as ago1-27 mutants, suggesting that AGO1 deficiency might disrupt ABA sensing and ABA responses. However, not P1/HC-Pros of all viral species have the same effect in ABA pathway regulation. Pasin et al. (2020) show that P1 Pp increases the amounts of ABA [5], which is conflicted with the ABA profile of P1/HC-Pro Tu . Our explanation is the sequence divergence of P1 Pp and P1 Tu , which showed only 19.35% amino acid identity, resulting in the difference in endogenous ABA accumulation.

P1/HC-Pro Tu Might Alter ABA and Calcium Signaling Crosstalk during Stomatal Closure and Drought Stress
Abiotic stress and biotic stress can initiate ABA signaling pathways that lead to many molecular and cellular responses [22,23]. Among the ABA-induced stress response genes, those controlling stomatal closure and opening are important during drought conditions [24], and stomatal immunity plays an important role in the restriction of pathogen entry [24]. Thus, stomatal movement serves as a platform for crosstalk between biotic and abiotic stress responses involving ABA action. Studies reporting Ca 2+ oscillations elicited by stimulating ABA and the temporal dynamics of Ca 2+ in ABA signaling provide strong evidence showing that Ca 2+ is a crucial element in the ABA signaling network [25]. The integration of ABA and calcium signalings could govern PP2C-type phosphatase regulators in the responses to abiotic stresses via the modulation of common targets [12]. Our comparative transcriptome profiles demonstrated that upregulated Ca 2+ -related genes in the P1/HC-Pro Tu plants were associated with ABA signaling (Figures 2 and 3; and Tables 4 and 5). The integration of both ABA and Ca 2+ signaling processes might occur and simultaneously induce stress responses to P1/HC-Pro Tu during exposure to drought/cold stress responses and particularly stomatal dynamics. Interestingly, an antagonistic regulatory mechanism controls stomatal movement via crosstalk among ABA, JA, and SA when pathogen effectors, i.e., P1/HC-Pro Tu , ingress into host tissues to induce a rapid defense response ( Figure 5A). In addition, the P1/HC-Pro Tu plants exhibited upregulated genes associated with SA and brassinosteroid (BR) signaling, including NPR3 (AT5G45110), which is involved in the negative regulation of defense responses against bacteria [26], and BAR1 (AT5G18360) and BKI1 (AT5G42750), which have Ca 2+ -dependent functions [27]. These results provided possible links among ABA, other phytohormones, and the secondary messenger calcium in stimulus-response reactions of the P1/HC-Pro Tu plants.

The LTP NGS Strategy Enables the Collection of a Miniature of the HTP Sequencing Data
RNA silencing in plants prevents virus accumulation [28,29], and accordingly, viruses have evolved various strategies to counteract this defense. Viral silencing of suppressor proteins blocks the production of siRNAs or the ability of siRNAs to reach their targets [30,31]. Previous studies proposed models for interfering with viral suppressors in endogenous silencing that contribute to viral symptom development [32,33]; however, the link between plant physiology and the underlying molecular mechanisms remains unclear. NGS technology enables an understanding of the roles of viral suppressors in stress responses and gene silencing mechanisms. Recently, a HTP transcriptome was applied in studies of the mechanism of virus-infected plant cells, and this approach enables a more accurate determination of plant-virus interactions [34,35]. Using a combination of microarray and RNA-Seq data, researchers can identify the molecular mechanisms and physiological alterations that might contribute to viral symptom development during acute infection [36]. Pasin et al. (2020) demonstrated that ABA can control RNA biology, including RNA stability, turnover, maturation, and translation, in viral suppressor transgenic plants [5].
Compared with the HTP NGS, LTP NGS represents a more economical approach for transcriptome studies due to the pooling of additional sample libraries. In this study, the DEGs identified using the LTP profiles could reflect the essence of the HTP profiles because we performed sequencing at only one-tenth of the sequencing depth used in the HTP analysis. The important genes correlated with drought stress, cold stress, senescence, and the RNA silencing pathway were accurately identified with the LTP dataset. Seventy-five network genes were identified with both the HTP and the LTP profiles, and these network genes included gene encoding PTGS components (i.e., AGO1 and AGO2) and miRNA targets (i.e., SPL3, SPL13A; CSD1, CSD2, CCS, and TOE2). These results imply that the LTP strategy has a certain level of accuracy with respect to transcriptome analysis. However, the Pearson correlation coefficients were not 100% identical between the HTP and LTP profiles. A low number of genes might be excluded in the LTP analysis, which would lead to differences in matrix genes and thus differences in the Pearson correlation coefficients between the HTP and LTP profiles. Therefore, we recommend that the LTP approach can be used in preliminary or primary studies to identify critical genes and to identify a backbone network before a more in-depth HTP approach is considered.

Conclusions
By using comparative networks from the HTP RNA-Seq profiles, it revealed that P1/HC-Pro Tu can trigger multiple phytohormone-mediated stress responses, especially the regulatory signaling pathways controlled by ABA/Ca 2+ . Upon stimulation by P1/HC-Pro Tu , expression changes of genes involved in various cellular processes, such as drought and cold stress responses, stomatal dynamics, senescence, cell death, and chlorophyll degradation, were linked to viral suppressor induced responses. The P1/HC-Pro Tu -mediated ABA disturbances via PTGS suppression might alter these gene expressions. Using a LTP NGS approach, we can simulate the HTP network for studying the P1/HC-Pro Tu -mediated PTGS suppression.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.