A single viral amino acid shapes the root system architecture of a plant host upon virus infection

Abstract

Background

Grapevine fanleaf virus (GFLV) is one of the most detrimental viral pathogens of grapevines worldwide but no information is available on its effect on the root system architecture (RSA) of plant hosts. We used two wildtype GFLV strains and their single amino acid mutants to assess RSA traits in infected Nicotiana benthamiana and evaluate transcriptomic changes in host root gene expression in replicated time course 3’RNA-Seq experiments. Mutations targeted the multi-functional GFLV-encoded protein 1E^Pol*/Sd, a putative RNA-dependent RNA polymerase and determinant of foliar symptoms in N. benthamiana plants.

Results

Plant infection with wildtype GFLV strain GHu and mutant GFLV strain F13 1E^Pol _G802K, both carrying a lysine in position 802 of protein 1E^Pol*/Sd, resulted in a significantly lower number of root tips (-30%), and a significantly increased average root diameter (+ 20%) at 17 days post inoculation (dpi) in comparison with roots of mock inoculated plants. In contrast, the RSA of plants infected with wildtype GFLV strain F13 and mutant GFLV strain GHu 1E^Pol _K802G, both carrying a glycine in position 802 of protein 1E^Pol*/Sd, resembled that of mock inoculated plants. Modifications of RSA traits were not associated with GFLV titer. Root tissue transcriptome analysis at 17 dpi indicated dysregulation of pattern recognition receptors, plant hormones, RNA silencing, and genes related to the production of reactive oxygen species (ROS). For wildtype GFLV strain GHu, RSA modifications were correlated with an abundant accumulation of ROS in the pericycle of primary roots at 7 dpi and the duration of vein clearing symptom expression in apical leaves. Dysegulation of a hypersensitive response was an overarching gene ontology found through enrichment analyses of 3’RNA-Seq data.

Conclusions

Our findings revealed the causative role of lysine in position 802 of protein 1E^Pol*/Sd in a novel RSA phenotype during viral infection and documented GFLV-N. benthamiana interactions at the root level based on (i) antiviral response, (ii) receptor mediated production of ROS, and (iii) hormone regulation. A correlation between above and below ground symptoms was reported for the first time in plants infected with wildtype GFLV strain GHu. Further work is warranted to test whether the modified RSA of a plant host might impact GFLV acquisition and transmission by the ectoparasitic dagger nematode Xiphinema index.

Introduction

Roots play a crucial role in plant development, contributing largely to overall health and growth. Plants rely on this underground network of tissue to sequester nutrients, acquire water, provide energy sinks for long chain carbohydrates, and act as anchorage to combat adverse environmental stressors. The configuration of a root system in the soil is commonly referred to as root system architecture (RSA). The overall RSA is dynamically influenced by abiotic and biotic factors, including pathogens such as viruses that can negatively modify RSA upon infection of plant hosts [1,2,3].

Some plant viruses are transmitted by soil-borne vectors [4] and are known to manipulate host root systems. For example, tobacco rattle virus (TRV) alters RSA and root volatiles of Nicotiana benthamiana, affecting transmission characteristics by trichodorid nematode vectors [5]. TRV manipulates the production pathway of 2-ethyl-1-hexanol to attract nematodes for feeding on infected tissue, and subsequentially to increase virus acquisition for more efficient transmission. Similarly, the RNA3-encoded p25 protein of beet necrotic yellow vein virus (BNYVV) determines both foliar and root symptoms in sugar beet (Beta vulgaris subsp. vulgaris) [6]. The constitutive expression of BNYVV protein p25 in Arabidopsis thaliana induces hormonal changes and a distinct RSA phenotype [1].

Grapevine fanleaf virus (GFLV, species Nepovirus foliumflabelli, genus Nepovirus, family Secoviridae, order Picornavirales) is a causative agent of fanleaf degeneration disease, and the most detrimental viral pathogen of grapevines worldwide [7]. In the vineyard, GFLV is transmitted by the ectoparasitic dagger nematode Xiphinema index in an exclusive manner, solely dependent upon the viral coat protein [8,9,10]. The GFLV genome is bipartite with two single-stranded RNA molecules (RNA1 and RNA2). GFLV RNA1 encodes proteins involved in polyprotein processing [11], replication [12, 13], and suppression of RNA silencing [14], while GFLV RNA2 encodes proteins involved in replication [15], cell-to-cell movement [16], and encapsidation [17]. Performing experiments with GFLV in the natural host (Vitis spp.) can be laborious and time-consuming; therefore, herbaceous model hosts are commonly utilized for GFLV-plant interaction studies [18]. Among herbaceous hosts, N. benthamiana is frequently used to study virus-plant interactions because it can easily be infected via mechanically inoculation or through Agrobacterium tumefaciens-mediated infiltration of GFLV infectious clones, among other factors [19].

Previous work revealed a single amino acid of GFLV strain GHu is responsible for a transient vein clearing phenotype in apical leaves of N. benthamiana plants [20, 21]. This symptom determinant (Sd) residue is located at position 802 (or position nts 2,404–2,406) of the GFLV RNA1-encoded putative RNA-dependent RNA polymerase (Pol*), refered to as protein 1E^Pol*/Sd [20,21,22]. GFLV strain GHu encodes a lysine (K) at position 802 of protein 1E^Pol*/Sd which is necessary but not sufficient to produce these symptoms in GFLV strain F13 [20, 21]. The viral symptom determinant region is hypothesized to interact with host components that ultimately elicit or repress an immune response, although structural or functional homologues could not be identified by advanced protein modeling techniques [22]. A combined time-resolved transcriptomic and proteomic approach documented a reactive oxygen species (ROS) molecular pathway in addition to the dysregulation of transcripts and proteins of plant immune response, transcriptional regulation, and translational regulation in N. benthamiana leaves upon GFLV infection [23]. Transient leaf symptoms can be observed on leaves of Vitis and Nicotiana species alike. Although experimental evidence of ROS production in N. benthamiana following GFLV infection is lacking, other studies have correlated ROS production with pronounced symptoms in GFLV-infected Vitis cultivars [24] or detected ROS in GFLV-infected N. occidentalis [25].

Studies on GFLV-host interactions have centered on above-ground plant tissue to investigate underpinnings of symptom development in leaves and determine how fruit yield and quality are affected by viral infection [7, 16, 19,20,21, 23,24,25]. No research has investigated the effects of GFLV on the RSA of a plant host. The aim of this study was to explore below ground features of N. benthamiana plants and probe the molecular dysregulation of root tissue upon GFLV infection. We hypothesized that RSA traits follow similar phenotypic trends to foliar symptomology, and that GFLV elicits similar transcriptomic profiles in roots as previously observed with leaves [19]. Here, we document GFLV protein 1E^Pol*/Sd single amino acid-specific modifications of RSA traits, ROS production in roots, and root transcriptome dysregulation upon infection.

Materials and methods

Nicotiana benthamiana plants

Seeds of wildtype N. benthamiana were sown in Cornell LM-2 soil, transplanted ~ 14 days after germination, and maintained in a controlled growth chamber at 25 °C and 80% humidity, as previously documented [19]. Plants were watered daily and fertilized (Scotts Miracle-Gro® All Purpose Plant Food 24–8-16) weekly. Four independent experiments were performed with 75, 74, 107, and 80 plants referred to as cohorts 1–4, respectively.

GFLV strains

Four GFLV strains were selected for this study: wildtype GFLV strains GHu and F13, and their mutants carrying a single amino acid change at position 802 of the RNA1-encoded RNA-dependent RNA polymerase [21]. These strains will be referred herein to as wildtype GHu, wildtype F13, mutant GHu 1E^Pol*/Sd_K802G, and mutant F13 1E^Pol*/Sd_G802K.

Plant mechanical inoculation with GFLV strains and leaf phenotyping

The genome of GFLV strains used in this study consisted of RNA1 of wildtype F13, wildtype GHu, mutant F13-1E^Pol*/Sd_G802K, or mutant GHu-1E^Pol*/Sd_K802G, alongside RNA2 of GFLV-GHu [19]. These combinations of GFLV RNA1 and RNA2 were selected for their high infection rates in N. benthamiana [21]. Infected plant tissue stored at -80 °C was used for mechanical inoculations [20]. The infection status of inoculated plants was confirmed in apical uninoculated leaves by double antibody sandwich enzyme-linked immunoassay (DAS-ELISA) at 14 days post inoculation (dpi) using GFLV specific antibodies (BIOREBA, Switzerland) (19). At 4 dpi, the status of GFLV infection was tested in inoculated leaves through RT-PCR and in roots by both RT-qPCR and 3’RNA-Seq of total RNA.

Leaves of N. benthamiana plants inoculated with wildtype GHu were rated daily for vein clearing symptoms on a binary scale: 0 for asymptomatic, and 1 for symptomatic. No symptoms were observed on the other GFLV treatments or mock control plants, as previously reported [21, 23].

Root system excavation, image acquisition and processing, and root trait analysis

The root excavation pipeline used in this study was previously reported as an open access BioProtocol [26]. Briefly, potted plants were gently released of soil by hand and root systems were dipped in water until all soil was removed. If persistent soil was adhering to roots, a spray bottle of distilled water was used to remove excess debris. The roots were dried on a paper towel and immediately placed on clear plastic imaging trays [26].

Imaging trays were placed on the scanning bed of an Epson Expression XL12000 photo scanner (Model:12000XL-PH) connected to a computer for image acquisition of the roots with the following parameters: ≥ 400 dpi, grayscale, and best quality. Images were analyzed in RhizoVision Explorer version 2.0.3 [27] with the following parameters: analysis mode (Broken roots), pixels converted to physical units at ≥ 400 dpi, image thresholding level = 158, (158 for cohort 1, 186 for cohort 2, 182 for cohort 3) (manually adjusted per image ± 10), filter non-root objects, enable edge smoothing = 2, feature extraction levels by root diameter (in pixel units) for fine (0.00 – 0.50), small (0.50 – 0.75), medium (0.75 – 1.75), and large (more than 1.75) roots. Using dpi and conversion factors of image acquisition, these measurements corresponded to thresholds of 0.00 mm, 0.32 mm, 0.48 mm, and 1.11 mm. Ultra-fine roots were not considered in this study because the excavation pipeline was optimized to analyze crown root tissue only pertaining to the first several orders of the root system including the primary root with secondary, and at least tertiary lateral roots extending off the primary. Output data was transferred to a Microsoft Excel file. All raw data, including images, are available on the GitHub page or within the supporting information.

Measurement of N. benthamiana plant height and biomass

The height of N. benthamiana plants was measured from the soil line barrier to the apical most part of the central stem with a ruler one day prior to root excavation. Shoots of N. benthamiana plants were collected into brown paper bags and thereafter roots were placed in the corresponding numbered bag. Dry biomass was measured after ~ 1 week in a drying oven at 105 °C. The ratio of root to shoot biomass was calculated for select cohort of N. benthamiana plants to assess biomass allocation.

3, 3'-diaminobenzidine (DAB) staining of N. benthamiana primary roots

Primary root samples of N. benthamiana plants were excavated from soil at 7- and 11-dpi for histological staining. For this, five plants per treatment were sacrificed for tissue collection, where at least three or more of the plants tested true to infection status through RT-PCR. Roots were vacuum infiltrated with 1 mg/mL DAB solution for five minutes and kept in the dark at ambient conditions overnight (8–12 h) [28]. De-staining was not performed because there is no chlorophyll to extract, and roots produce high amounts of ROS when exposed to ethanol or other abiotic stressors. Cross-sections of roots were made by hand just above where lateral roots formed (e.g., first branching points of root system) and performed two centimeters from the cut made at the stem of plants during excavation. Root cross-sections were placed on a glass microscope slide (Fisher, Cat No. 12–550-A3) mounted with 0.1 mL of 50:50 v/v glycerol:PBST (137 mM NaCl, 2.7 mM KCl, 10 mM Na₂HPO₄, 1.8 mM KH₂PO₄), and protected with a coverslip (Fisher, Cat No. 12–545-88). Imaging of root cross sections was performed immediately following slide mounting to ensure no supplemental induction of ROS production due to additional stress. Images were acquired on an Olympus SZX12 microscope using the following parameters: Gamma 1.16, Contrast 0.79, Saturation 1.53, Color Temp 4188, Exposure 160–200, Gain 1.

N. benthamiana root tissue collection, nucleic acid extractions, RT-PCR, and RT-qPCR

After root system images were acquired, five N. benthamiana plants that tested true to virus infection status were randomly chosen for GFLV titer analysis. Roots of mock inoculated N. benthamiana plants were used as controls in quantification assays. Approximately half of the root system was sectioned by excising primary, lateral, and fine roots with a sterile razor blade. Samples were placed in 2 mL Eppendorf tubes, flash frozen in liquid nitrogen, and stored at -80 °C until further processing. Total nucleic acid was extracted from root samples using E.Z.N.A.® Plant RNA Kit (SKU:R6827-01) with a final elution volume of 120 µL nuclease-free water. Then, 40 µL of the eluted nucleic acids were used for DNAse I treatment at 20 °C for 15 min followed by RNA clean-up with Zymo Research RNA Clean & Concentrator®-5 (Cat # R1013). RNA integrity and quantity was assessed on an Agilent 2100 Bioanalyzer with the RNA 6000 Nano Kit Assay (No. 5067–1511) according to manufacturer instructions.

Total RNA was used to amplify the GFLV RNA1 nucleotide sequence region covering amino acid in position 802 of protein 1E^Pol*/Sd with the QIAGEN® OneStep RT-PCR Kit (No. 210210) for Sanger sequence confirmation, as previously described [21]. Total RNA was also used for GFLV titer analysis by Luna® Universal One-Step RT-qPCR (Cat. No. E3005) with primers specific to GFLV RNA1 [21] and to the endogenous N. benthamiana control FBOX gene [29] (Supplementary Table 1). Samples were tested in triplicate, and quantification cycle (Cq) values for SYBR green were collected after each cycle, as previously reported [21]. GFLV treatments were compared with mock inoculated control samples to calculate ∆∆Cq values, in which relative expression was calculated as log(2^−∆∆Cq) values [30].

GFLV strain identification

Amplified RT-qPCR products of ~ 250 bp in size obtained with GFLV RNA1 primers from roots of two to four infected plants at 17- and 26 dpi were cleaned up with Zymo DNA Clean & Concentrator®-5 (Cat #D4004). The integrity of the two wildtype and two mutant GFLV strains was tested in roots through Sanger sequencing of the cleaned amplicons. Sequences obtained in this study were compared with GFLV sequences deposited in GenBank (NC003615 for wildtype F13 RNA1 and JN391442 for wildtype GHu RNA1) using SeqMan Ultra™ (v.17.3.1) of the DNASTAR Navigator Package (v17.3.0.57) (DNASTAR, Inc., Madison, Wisconsin, USA).

3 ‘-ribonucleic acid sequencing

The transcriptome of N. benthamiana roots was characterized via 3’RNA-Seq for wildtype GHu, mutant GHu 1E^Pol*/Sd_K802G, and mock control with samples collected at 4-, 17-, and 26 dpi. Three biological replicates for each time point per viral treatment were purified, concentrated, and DNAse I treated to produce a total of 27 samples. Libraries were prepared with the Illumina® TruSeq RNA Library Prep Kit (RS-122-200X) and sequenced via Illumina® NextSeq500/550. Transcript reads were trimmed with Trimmomatic (v0.39) [31] and subsequently mapped to the host genome, as well as the previously curated GFLV library ‘Nepo’ [23] using HiSat2 (v2.2.1) [32]. Sequences were checked for quality control using the FastQC program (v0.11.9) [33].

Bioinformatic analysis

Sequence reads were mapped to the N. benthamiana genome from two different sources: Niben v1.0.1 (Boyce Thompson Institute, Ithaca, NY, USA) [34] with 61,328 gene models and Nbe_v1 [35] with 57,583 high confidence gene models. Through a bioinformatic pipeline, the Niben v1.0.1 genome was reannotated in this study using an updated list of gene ontology identifiers by subjecting all protein encoding sequences to a BLAST search through EggNOG (v5.0) [36] [Supplementary Fig. 1]. After curation of additional annotation terms, the files were parsed and linked through the GO API (Open Biological and Biomedical Ontology Foundry, http://api.geneontology.org/api) to create a .gmt file for use in gene enrichment. Data files were then manually curated in RStudio and R [37, 38].

The DESeq2 package (v1.40.2) [39] was used to normalize, test differential transcript abundance, and prepare data for downstream analysis, as performed previously [23]. Gene set enrichment analysis was performed using gProfiler (v0.2.1) [40] and clusterProfiler (v4.8.2) [41] packages. Additional packages used for visualization and data analysis include the following: enrichplot (v1.20.1) [42], pathview (v1.40.0) [43], dplyr (v1.3.1) [44], ggplot2 (v3.4.3) [45], and GO.db (v3.16.0) [46].

Statistical analyses

Matrices of acquired root traits were subject to Pearson’s correlation analysis [47] with the hierarchal clustering method. RStudio package corrplot (v0.92) [48] was employed to calculate correlation coefficients (-1 to + 1) and p-values at four significance intervals (* p < 0.05, ** p < 0.01, *** p < 0.001, **** p < 0.0001). In the instances of parametric data and normal distribution, root data were analyzed through analysis of variance (ANOVA) followed by student’s t-test for mean testing differences amongst treatments to controls. The p-value for parametric data is reported through the adjustment of Holm’s method. For instances of data with non-parametric or non-normal distributions, Kruskal–Wallis [49] tests followed by a Bonferroni correction [50] were performed. Tukey’s HSD with associated letters of significance were utilized for follow up statistical testing in root phenotype analyses. A Dunn’s test with Holm’s method adjusted p-values was utilized for mean comparison testing in non-parametric data. For the analysis of virus titer by RT-qPCR, logarithmic transformed values were used in a geometric mean comparison [51] and multiple t-test. Further analysis of RNA-Seq data utilized native statistical testing within clusterProfiler, DESeq2, WGCNA, enrichplot, and gProfiler packages.

Days showing symptoms (DSS) on apical leaves of N. benthamiana plants infected with wildtype GHu were calculated by subtracting the first and last day of symptom observation, and this metric was used as a predictor variable for RSA traits. Subsequent linear models were constructed and visualized in RStudio through the intrinsic stats package (v4.1.3). Mixed models were also constructed by using the predictors of day first symptomatic (DFS), day last symptomatic (DLS), and the interaction of these two terms with regards to foliar symptoms of wildtype GHu.

Results

Confirmation of GFLV infection status and strain integrity in N. benthamiana plants

At 4 dpi, GFLV was detected in most inoculated leaf samples of representative infected plants by RT-qPCR (75%, 15/20) or by RNA-Seq of root tissue (100%, 6/6), but was not detectable by DAS-ELISA, as expected for such an early infection stage. At 14 dpi, a high GFLV infection rate was observed in three cohorts of plants [100% (60/60) for cohort #1, 73% (52/71) for cohort #2, and 91% (85/93) for cohort #3] in apical leaves through DAS-ELISA. None (0%, 0/38) of the mock inoculated control plants tested positive for GFLV in DAS-ELISA in any of the replicated experiments, as expected.

A glycine (GCC in RNA1 nt position 6,787–6,789) was found in the progeny of wildtype F13 and mutant GHu 1E^Pol*/Sd_K802G, while a lysine (AAG in RNA1 nt position 6,797–6,799) was observed in the progeny of wildtype GHu and mutant F13 1E^Pol*/Sd_G802K, as anticipated (Supplementary Table 1). These analyses revealed high GFLV infection rates of inoculated N. benthamiana plants, and the genetic integrity of the GFLV progeny (Supplementary Table 2).

RSA trait acquisition and analysis

Roots of GFLV-infected and mock inoculated N. benthamiana plants were excavated to result in 263 crown root samples being scanned at ≥ 400 DPI. A total of 77 images captured a first and second side of the whole RSA per plant after a 180º rotation. The three time points of excavation, i.e., 4-, 17-, and 26 dpi, displayed distinct phenotypes with the most dramatic at 17 dpi for plants infected with wildtype GHu (Figs. 1 and 2, Supplementary Fig. 2).

Fig. 1

Root system architecture traits for Nicotiana benthamiana infected with four strains of grapevine fanleaf virus (i.e., wildtypes GHu and F13, and mutants GHu 1E^Pol*/Sd_K802G and F13 1E.^Pol*/Sd_G802K) at 17 days post-inoculation. a Number of root tips per plant, (b) average root diameter in millimeters, and (c) number of branching points per treatment group. Distribution of each continuous value is plotted along the y-axis as a violin plot with boxplot within, showing additional mean with standard error (central point with error bars). ANOVA significance values are displayed in top left of each panel, as well as students t-test relative to control as significance values (< 0.05 = *, < 0.01 = **, < 0.001 = ***, < 0.0001 = ****). The number of samples per treatment is displayed above panel A. d Representative images of mean phenotype observed between treatments with border color indicative of treatment as in plots A, B, and C. White bar displayed in bottom right corner of each image is 1.0 cm, calculated in ImageJ software (v1.8.0.172)

Fig. 2

Root system architecture traits for Nicotiana benthamiana infected with four strains of grapevine fanleaf virus (i.e., wildtypes GHu and F13, and mutants GHu 1E^Pol*/Sd_K802G and F13 1E.^Pol*/Sd_G802K) at 26 days post-inoculation. a Number of root tips per plant, (b) average diameter in millimeters, and (c) number of branch points plotted per viral treatment group. Distribution of each continuous value is plotted along the y-axis as a violin plot with boxplot within, showing additional mean with standard error (central point with error bars). Pairwise Wilcox test results are displayed as significance values (< 0.05 = *, < 0.01 = **, < 0.001 = ***, < 0.0001 = ****). The number of samples per treatment is displayed above panel A. d Representative images of mean phenotype observed between treatments, border color is indicative of treatment as in plots A, B, and C. White bar displayed in bottom right corner of each image is 1.0 cm, calculated in ImageJ software (v1.8.0.172)

By focusing on ten root traits (number of root tips, number of branch points, total root length, branching frequency, network area, average diameter, median diameter, perimeter, volume, and surface area) of the 32 metrics obtained by Rhizovision analysis, significant differences were found with several co-variable or often highly correlated traits (Fig. 3a). Five root traits (number of root tips, root diameter, total root surface area, and root branching frequency) were selected for further analysis through pairwise treatment testing at each time point.

Fig. 3

The root system architecture (RSA) of Nicotiana benthamiana plants upon grapevine fanleaf virus infection is dependent on the identity of amino acid 802 in protein 1E^Pol*/Sd. a RSA traits across four experimental replicates show highly correlated phenes that are dependent upon one another across all treatment groups. The variable ‘days showing symptoms’ was calculated from observational documentation of vein clearing symptoms in apical leaves for wildtype GHu infected plants only. b A histogram displays the first and last day of symptoms (in teal and purple) before full recovery for apical leaves in wildtype GHu infected N. benthamiana. c A normal distribution appears when calculating the total days symptomatic across three experimental replicates (n = 56). At 17 days post inoculation (d) the number of days showing symptoms had a negative relationship with the number of total root tips and inversely, (e) the days showing symptoms in leaves had a positive relationship with the average diameter in millimeters

At 4 dpi, all GFLV strains indicated some level of initial RSA perturbation (Supplementary Fig. 4d) with an overall 7–23% reduction in root tip number, a 20–25% increase in average diameter, and a 30–48% reduction in branching points, when compared with mock inoculated plants. However, no significant differences were observed for the number of root tips for virus strain, while the average diameter significantly differed between all viral strains compared to mock controls, as shown by Tukey’s test and pairwise analysis (Supplementary Fig. 4a,b). The number of branching points was less conclusive because all viral treatments, except wildtype GHu, were significantly different from mock inoculated plants, and wildtype GHu overlapped with all treatments according to Tukey’s test (Supplementary Fig. 4c). The limited changes in the overall RSA at 4dpi were likely due to younger plants being exposed to the effects of viral infection for a short period.

At 17 dpi, wildtype GHu and mutant F13 1E^Pol*/Sd_G802K caused up to a 33% reduction of the total number of root tips (****p < 0.0001, Fig. 1a), and a 21%-22% increase in average root diameter compared with mock inoculated control plants (***p < 0.001, Fig. 1b, Supplementary Tables 3, 4). These two GFLV strains also decreased root branching frequency, a measure of the number of branch points per total root length by up to 10.1% (**p < 0.01, Fig. 1c, Supplementary Tables 3, 4). Reduced root tip numbers, less fibrous but thicker roots were visually apparent for wildtype GHu and mutant F13 1E^Pol*/Sd_G802K compared with the other treatments (Fig. 1d). Wildtype F13 and mutant GHu 1E^Pol*/Sd_K802G increased the average root diameter by 13% and 10%, respectively, compared with mock inoculated control plants (*p < 0.05, Fig. 1d). Dry root biomass was reduced in plants infected wildtype GHu (21%), wildtype F13 (9%), and mutant GHu 1E^Pol*/Sd_K802G (10%) when compared with mock controls while plants infected with mutant F13 1E^Pol*/Sd_G802K had the same biomass as mock controls (Supplementary Fig. 5). These results documented drastic changes in the RSA of plants infected with wildtype GHu and mutant F13 1E^Pol*/Sd_G802K, both containing K at position 802 of protein 1E^Pol*/Sd, at 17 dpi.

At 26 dpi, plants infected with wildtype GHu exhibited a 25% reduction in the total number of root tips (*p < 0.05) and a 6.0% increase in average root diameter when compared with mock control plants (Fig. 2a, b, Supplementary Tables 3, 4). Plants infected with mutant F13 1E^Pol*/Sd_G802K displayed an increased average root diameter of 13% (***p < 0.001) but little change in the number of root tips (p > 0.05) when compared with mock control plants. Branching frequency was consistent across all treatment groups (Fig. 2c). These characteristic RSA traits were apparent on pictures of representative root systems (Fig. 2d). Dry root and dry shoot biomass decreased by 26% and 20%, respectively, for plants infected with wildtype GHu, while an increase in these metrics was observed with all other virus strains in comparison with mock inoculated plants (Supplementary Table 5). These results documented drastic changes in the RSA of plants infected with wildtype GHu and to a lesser extent with mutant F13 1E^Pol*/Sd_G802K, both containing K at position 802 of protein 1E^Pol*/Sd, at 26 dpi. Additionally, many traits such as perimeter, surface area, root to shoot biomass, and plant height were less pronounced when comparing treatments at all time points (Supplementary Fig. 5).

GFLV titer in root tissue

No significant differences in GFLV titer were found at 4 and 17 dpi, as shown by student’s t-test of RT-qPCR data (Supplementary Fig. 6a). At 26 dpi, the titer of wildtype F13 was significantly higher from both wildtype GHu and mutant GHu 1E^Pol*/Sd_K802G, however no other viral treatments differed through student’s t-test (Supplementary Fig. 6a). The titer of wildtype F13 was significantly different between 17- and 26-dpi, as shown by a Wilcoxon test (*p = 0.016), however no other viral strain had significantly different titer between time points (p > 0.05) (Supplementary Fig. 6a). As expected, mock control plants retained no viral titer. These data suggested that virus titer was unrelated to RSA traits during symptom development.

To complement RT-qPCR data, GFLV RNA1 and RNA2 reads from 3’RNA-Seq data were extracted and analyzed. No significant differences were observed amongst viral strains at 4 dpi or 17 dpi, but highly significant differences in RNA1 (**p = 0.00014) and RNA2 (****p = 0.0000008) were observed when comparing wildtype GHu to mutant GHu 1E^Pol*/Sd_K802G at 26 dpi, as shown by Tukey’s HSD test (Supplementary Fig. 6b,c). This trend was numerically similar but not statistically significant in the RT-qPCR data. Together, these data suggested that GFLV titer is not directly correlated with root phenology at 17 dpi but differed amongst viral strains and corresponding root phenotype later in infection.

Effect of GFLV infection on plant height

At 4 dpi, plants infected with wildtype GHu and F13 had a 2% and 6% decrease in height, respectively, while plants infected with mutant GHu 1E^Pol*/Sd_K802G and F13 1E^Pol*/Sd_G802K displayed a 1% and 9% increase in height, respectively (Supplementary Table 5, Supplementary Fig. 5). At 17 dpi, GFLV-infected plants had increases in height (up to 18%) for asymptomatic strains and symptomatic wildtype GHu (up to 3%) (Supplementary Table 5, Supplementary Fig. 5). At 26 dpi, height was reduced by an average of 16% and 9% for plants infected with wildtype GHu and F13, while mutant GFLV strains displayed no change or a minor increase in height when compared with mock inoculated plants (Supplementary Fig. 5). These data revealed that GFLV infection had overall more drastic effects on plant height at late infection stages rather than during early infection.

Correlation of RSA traits and foliar symptomatology for wildtype GHu

A strong positive correlation was observed for root diameter and days showing symptoms (DSS) in apical leaves. Conversely, a strong negative correlation was found between DSS and root branching traits such as root tip number, number of branching points, and network area (Fig. 3a, Supplementary Fig. 3). At 17 dpi, dry biomass, network area, perimeter, and surface area were significantly reduced for plants infected with wildtype GHu, and to lesser extents across other viral strains several, illustrating modifications above and below ground (Supplementary Table 5). At 26 dpi, additional traits were significantly impacted but not to the level described at 17 dpi (Supplementary Table 5). While plants infected with other viral strains remained close to mock control plants, wildtype GHu induced a lessened amount of root tips, number of branch points, total root length (mm), network area (mm²), perimeter (mm), volume (mm³), and surface area (mm²) (Supplementary Table 5).

Since many RSA traits were significantly and strongly correlated (positively or negatively) with leaf symptomology in N. benthamiana infected with wildtype GHu, general linear models and mixed models were constructed to predict RSA traits. A relatively normal distribution of symptom emergence centered at 7 dpi (average 6.77 dpi) with symptom disappearance around 11 dpi (average 11.37 dpi) and vein clearing symptoms apparent for an average of 5.66 days across all 53 N. benthamiana plants infected by wildtype GHu (Fig. 3b,c, Supplementary Table 6). A linear model constructed using the continuous response variables from acquired RSA traits against the predictor variable of days symptomatic documented an inverse relationship between the number of root tips and the number of days symptomatic in leaf tissue (RSE = 32.35, adjusted R² = 0.254, p = 0.00266, Fig. 3d). The observed variance in average root diameter could be significantly explained by the symptomatic days (RSE = 0.1454, adjusted R² = 0.1517, p = 0.01912, Fig. 3e), albeit explaining less variation compared with root tip number. These results illustrated an association between symptoms above and below ground for N. benthamiana infected with wildtype GHu.

A mixed model incorporating DFS and DLS as date component rather than difference of the two dates performed nearly twice as well for modeling root tip number (RSE = 23.95, adjusted R² = 0.6067, p = 1.105E-5), and nearly as well for median root diameter prediction (RSE = 0.07733, adjusted R² = 0.3649, p = 0.002929). At 26 dpi, correlation analysis was strong although linear modeling returned insignificant results [p > 0.05] (Supplementary Fig. 7). These observations suggested a correlation between apical leaf symptoms and long-term root phenotypes for wildtype GHu infected plants.

Reactive oxygen species are produced in roots during GFLV infection

ROS production was monitored in root tissue of N. benthamiana via histological cross sections of primary roots followed by DAB staining at 7 dpi (n = 25) and 11 dpi (n = 25) (Fig. 4). Filamentous roots were very dark from DAB staining, likely due to the stress of digging up the plants from the soil and the abundance of ROS in emerging root tips (Fig. 4a). At 7 dpi, roots of plants infected with wildtype GHu and mutant F13 1E^Pol*/Sd_G802K contained a ring of dark stain bordering the pericycle and vascular bundle, but not the cortex (Fig. 4b). This was not observed for any GFLV strain at 11 dpi. Roots of plants infected with wildtype F13 and mutant GHu 1E^Pol*/Sd_K802G did not have the same level of contrast but showed stain slightly darker than roots of control plants in the outermost parts of the vascular bundle (Fig. 4b). Mock control plants at 7 dpi (n = 5) and 11 dpi (n = 5) yielded minimal background staining of any vasculature or cortex, but epidermal tissue consistently showed a dark brown or black color. These results suggested that the nature of residue 802 of protein 1E^Pol*/Sd largely influences the amount of ROS produced in the primary roots of N. benthamiana.

Fig. 4

Reactive oxygen species (ROS) transiently produced in Nicotiana benthamiana primary root pericycle by four grapevine fanleaf virus (GFLV) strains (i.e., wildtypes GHu and F13, and mutants GHu 1E^Pol*/Sd_K802G and F13 1E^Pol*/Sd_G802K) at 7- and 11- days post-inoculation (dpi). a Anatomical view of primary root cross sections from N. benthamiana immediately stained with 3′3-diaminobenzidine, dark portions represent tissue with high amounts of ROS reacting with the stain. Fine roots produced high amounts of ROS as expected from rapidly growing root tissue. b Primary root cross sections of virus treatments display high amounts of precipitated 3–3’-diaminobenzidine stain in the pericycle for both wildtype GHu and mutant F13 1E^Pol*/Sd_G802K at 7 dpi and a diminished ROS abundance at 11 dpi. All other treatments have much lower levels of ROS detectable within the vasculature and cortex at both time points. (Epi = root epidermis, C = cortex, V = vascular bundle, LR = lateral root primordium). Images are representative of at least three individuals (n = 3) of five sacrificed plants per treatment

3’RNA-sequencing, and data curation

Twenty-seven N. benthamiana root samples were sequenced on Illumina 500/550 platform, three replicates of each treatment (wildtype GHu, mutant GHu 1E^Pol*/Sd_K802G, mock control) at three time points (4-, 17-, 26 dpi). The average number of total reads per sample was 7,338,802. After identification of reads passing quality control thresholds and proper adapter trimming, an average of 7,271,449 reads remained per sample. Alignment to the Niben v1.0.1 genome [34] resulted in an average of 74.63% reads per sample, after outlier removal of sample mutant GHu 1E^Pol*/Sd_K802G 17 dpi replicate #3 (low alignment rate ~ 6%). Alignment to the Nbe_v1 genome resulted in an average of 5% less aligned reads, therefore, downstream analysis utilized the Niben v1.0.1 alignment. Quality control as reported through FastQC analysis of RNA-sequencing files returned average sequence lengths of 86 base pairs, 36–39% GC content, and Phred quality scores maintaining a majority above 30, with many sequences at 35 (1 in 5000 base pair call inaccuracy), which were all acceptable for downstream analysis (Supplementary Table 7). Re-annotation of the N. benthamiana genome was performed to retrieve 11,525 GO terms that matched to the Niben v1.0.1 protein models (Supplementary Fig. 1).

Initial statistical considerations of the 3’RNA-Seq dataset

FeatureCounts files were parsed to produce a gene count matrix that was utilized with a metadata matrix containing treatment, time, and a group factor of treatment*time for DESeq2 input. The resulting matrix represented 37,784 genes with a threshold of at least 20 counts across all samples. Count filtering, dispersion of counts, hierarchal clustering, and principal component analysis were performed on data prior to removal of mutant GHu 1E^Pol*/Sd_K802G replicate #3 at 17 dpi, because it showed large separation from all other samples (Supplementary Fig. 8a-d). After removal of this outlier, the analysis pipeline continued with the remaining 26 samples. Normalized counts returned roughly ~ 200,000 reads per sample across all time points and treatments (Supplementary Fig. 9a). Further, dispersion estimate curves were largely normal across the DESeq2 experiment, indicating that the normalization of counts and quality of reads were acceptable. A strong grouping by both treatment and time in hierarchal clustering and principal component analysis was observed (Supplementary Fig. 9c-d). All mock control samples grouped away from GFLV-infected samples, and GFLV-infected samples separated by time point through McQuitty method hierarchal clustering (Supplementary Fig. 9c). Further, samples were discriminated by treatment and time through principal component analysis with PC1 of 29.89% and PC2 of 12.42%. Samples at 4 dpi grouped in the middle of the plot, loosely associated, while samples at 17 dpi grouped to the right of the plot. Samples collected at 26 dpi were located on the left of the PCA plot in tight association, suggesting that these samples were differentiated from the other time points based on PC1, but less discriminated by PC2. Based on these results, the dispersion of counts being acceptable after normalization (Supplementary Fig. 9c), and the successful alignment of 201,108 reads per sample, data were deemed acceptable for further analysis.

Differential transcript abundance analysis

Contrasts were generated within the DESeq2 package to examine differences in the transcriptome between treatments. Thresholds of P < 0.05 and |Log2FoldChange|> 1 were used for determining differentially abundant genes (DEGs). When comparing wildtype GHu to mutant GHu 1E^Pol*/Sd_K802G at 4-, 17- and 26 dpi, a total of 704, 761, and 1041 transcripts was found to be differentially abundant, respectively (Table 1, Fig. 5a-c). An increased number of differentially abundant transcripts was obtained when comparing viral treatments to the mock controls at each time point. When comparing wildtype GHu to mock control transcripts, a total of 713, 4351, and 1692 transcripts were differentially abundant at 4-, 17- and 26 dpi, respectively (Table 1, Fig. 5 d-f). Lastly, when comparing mutant GHu 1E^Pol*/Sd_K802G to mock control, a total of 437, 3293, and 2124 DEGs were found at 4-, 17-, and 26 dpi, respectively (Table 1, Fig. 5g-i).

Table 1 Differentially abundant transcripts as compared through DESeq2 analysis of Nicotiana benthamiana plant roots infected with two grapevine fanleaf virus strains, i.e., wildtype GHu and mutant GHu 1E^Pol*/Sd_K802G) at 4-, 7- and 11 days post-inoculation (dpi). DEGs were defined as having both a p-value < 0.05 and an absolute Log2FoldChange > 1

Fig. 5

Volcano plot contrasts of Nicotiana benthamiana plants infected with two grapevine fanleaf virus (GFLV) strains (i.e., wildtypes GHu and mutant GHu 1E^Pol*/Sd_K802G) reveal perturbations of the root transcriptome at three time points; 4-, 17- and 26 days post-inoculation (dpi). Contrasts are divided by treatment in columns such as wildtype GHu vs mutant GHu 1E^Pol*/Sd_K802G (a, b, c), wildtype GHu vs mock control (d, e, f), and mutant GHu 1E^Pol*/Sd_K802G vs mock control (g, h, i) and by time (rows) at 4 dpi (a, d, g), 17 dpi (b, e, h) and 26 dpi (c, f, i). Logarithmic transformation of transcript abundance comparisons is shown on the x-axis. Confidence values as logarithmic transformations of p-value are represented on the y-axis. Genes shown as purple data points meet qualifications of p-value < 0.05 and absolute Log2FoldChange ≥ 1. Genes in gray meet only one qualification, while black data points denote neither qualification was met

The top ten N. benthamiana root genes of differential abundance per treatment contrast and time point were identified (Table 2). By contrasting wildtype GHu and mutant GHu 1E^Pol*/Sd_K802G, host root genes relating to immune response, signaling, disease resistance, hydrolase activity, and components of viral silencing were consistently found to be dysregulated. Transcripts for regulators of transcription, translation, and post-translational processes were most abundant (Table 2). When contrasting wildtype GHu to mock inoculated plants, host genes relating to transcription, mitochondria and signaling patterns, protein synthesis, and immune response were dysregulated in the top 10 DEGs across all time points. Further, glutamine synthase and phosphatidylinositol-related genes appeared in the top 10 DEGs at 17 dpi, consistent with previously identified candidates in apical leaf tissue during GFLV infection and foliar symptom development [23]. Additional regulators of transcription and translation, immune response and signaling genes were dysregulated in the contrast of mutant GHu 1E^Pol*/Sd_K802G and mock control across the three time points. Across all contrasts, genes relating to innate immune response were consistently present, suggesting residue 1E^Pol*/Sd₈₀₂ impacts the level in which immune responses are activated during infection. Notably, five unknown proteins were dysregulated to high magnitudes up to Log2FC of 14.66. These genes could have important functions related to morphology or viral infection not yet characterized in plant root tissue. All differentially abundant transcripts were documented (Supplementary Table 8).

Table 2 The top ten differentially abundant transcripts present in each viral contrast with their respective gene ID, name, Log2FoldChange, and p-value, ranked by p-value

Gene set enrichment analysis

From each of the treatment contrasts, DEGs were extracted and subject to GSEA to detect overrepresented gene ontology terms. An ordered list of DEGs was used to produce enriched terms between each contrast at each respective time point (Supplementary Fig. 10). Few terms were produced during contrasts of wildtype GHu and mutant GHu 1E^Pol*/Sd_K802G root transcriptome profiles to include response to external stimulus (GO:0009605) and pigment binding (GO:0031409) at 4 dpi, negative regulation of plant-type hypersensitive response (GO:0034051) at 17 dpi, and regulation of steroid metabolic process (GO:0019218) at 26 dpi. The relatively low number of DEGs detected between these two treatment contrasts might have prevented further resolution of these data related to the observed changes in RSA.

Many more gene ontologies were detected through GSEA by contrasting wildtype GHu with mock inoculated control plants. At 4 dpi, ontologies related to dephosphorylation (GO:0016311) and phosphatase activity (GO:0016791) amongst other closely related GOs were overrepresented. Out of 59 total terms at 17 dpi, the top three GO’s overrepresented according to p-value were DNA-binding transcription factor (GO:0003700), transcription regulator activity (GO:0140110), and response to stimulus (GO:0050896). Several terms relating to biosynthetic processes, phosphorylation, response to oxygen species, and hormone response were overrepresented at this time point as well (Fig. 6a). Pigment binding (GO:0031409), obsolete anchored component of membrane (GO:0031225), and photosynthesis, light harvesting in photosystem I (GO:0009768) were the top three overrepresented GOs at 26 dpi (Fig. 6a). Further terms related to photosynthetic components and secondary metabolite regulation were found at 26 dpi.

Fig. 6

Downstream overabundance of transcript dysregulation in Nicotiana benthamiana roots infected with wildtype grapevine fanleaf virus strain GHu through gene set enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) reveals several conserved pathways of overabundance. a The dotplot representation of wildtype GHu infected plants against mock control GSEA shows pathways dysregulated at all three time points past inoculation. At early infection, 4 days post inoculation (dpi), categories related to signaling and phosphorylation are displayed. At 17 dpi, many components of metabolism, immune pathways, and signaling are overrepresented. Finally at 26 dpi, components related to physiology, photosynthesis, and signaling were significantly overrepresented. The size of each data point represents the ratio of genes dysregulated in this category related to the total number of genes in the genome. b KEGG analysis shows similar pathways at early infection for specific signaling and innate immune encoding genes. c At 17 dpi unique signaling molecules through gene regulation were overrepresented in addition to structural and transport genes. d Finally, similar genes were found at 26 dpi to 17 dpi, with the addition of differential abundance of plant hormone pathways and a general response factor chitinase

Contrasting mutant GHu 1E^Pol*/Sd_K802G and mock inoculated control plants revealed gene ontologies that were reminiscent of viral infection without all changes occurring to root morphology, as observed through the previously described RSA phenotyping (Supplementary Fig. 10). Thylakoid (GO:0009579), photosynthetic membrane (GO:0034357), and singlet oxygen-mediated programmed cell death (GO:0010343) were among the overrepresented, non-redundant GO’s present at 4 dpi. At 17 dpi, plasma membrane (GO:0005886), cell periphery (GO:0071944), and response to stimulus (GO:0050896) were the top three overrepresented GOs by p-value, with other terms relating to hormone response and phosphorylation activity also being abundant at this time. At 26 dpi, 29 terms were overrepresented with the top three by p-value being regulation of hormone levels (GO:0010817), hormone metabolic process (GO:0042445), and lipid metabolic process (GO:0006629). Other closely related terms of response to stimulus, plasma membrane, and phosphorylation were also present during late infection (Supplementary Fig. 10).

A secondary evaluation of enrichment for differentially abundant transcripts found in roots was performed through Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis. Significant KEGG terms were defined by the thresholds of p < 0.05 and q < 0.05 to evaluate the 3’RNA-Seq differentially abundant genes against the newly curated genome annotation. The enriched KEGG terms confirmed previous GSEA results to display terms of immune response, gene regulation, and signal transduction (Fig. 6b-d, Supplementary Fig. 11). Signaling genes such as calcium-binding protein CML (ko:k13448), phosphatase 2C of the MAP kinase pathway (ko:k14497), and calcium dependent protein kinase (ko:k13412) were present at various contrasts and time points during GFLV infection (Supplementary Fig. 11a,d,h). Interestingly, a new term relating to cell elongation and division appeared as well throughout multiple contrasts and time points: xyloglucan transferase (ko:k08235) (Supplementary Fig. 11b,e,f). The results of the KEGG analysis further documented phytohormone pathways (ethylene, abscisic acid, transcription factor TGA) to be impacted upon viral infection although their relation to the development of root symptoms may not be immediately apparent. Several terms related to protein repair such as chaperone proteins, oxidases, and hormone regulators were found at 26 dpi from KEGG analysis (Supplementary Fig. 11c,f,i).

Dysregulation of host pathways associated with viral infection and root symptom development

The three major host pathways resolved from GSEA, KEGG, and top candidates of differential abundance analysis were: (i) regulation of gene expression, (ii) immune signaling and oxidative response to pathogen invasion, and (iii) production of secondary metabolites (Fig. 6). Substantial overlap of these pathways was previously identified in apical leaf tissue through transcriptomic and proteomic analysis [23]. A fourth pathway of interest was considered based on the recent description of GFLV suppression of RNA silencing (14). From top gene candidates in root RNA-Seq analysis, immune signaling pathways were highly overrepresented, especially innate immune receptors (positive and negative dysregulation), cAMP phosphoprotein downregulated in both viral contrasts at 17 dpi, and proteinase inhibitors negatively regulated for wildtype GHu while upregulated for mutant GHu 1E^Pol*/Sd_G802K when compared against the mock controls (Fig. 7). Several host-encoded gene silencing transcripts, including suppressor of gene silencing 3, SDE 1 and 2, AGO4, HEN1, and DRBP4, were less abundant in GFLV-infected N. benthamiana roots. Other members of the RNA-induced silencing complex such as AGO1 and AGO2 were upregulated, but other transcripts were either too low in number or provided orthologs with conflicting patterns to deduce a trend (Fig. 7).

Fig. 7

The dysregulation of three major pathways in the roots of Nicotiana benthamiana infected by grapevine fanleaf virus (GFLV) revealed by ontology and enrichment analysis. a Immune receptors and innate immune pathways were disrupted at 4- and 17-days post inoculation (dpi) when comparing viral treatments to mock inoculated plants. Further, wildtype GHu infected plants were more dysregulated across multiple types of innate immune pathways that feed into the PI3K cascade. b The host viral silencing pathway had several genes dysregulated throughout the course of GFLV infection, in that, either both viral treatments dysregulated or wildtype GHu altered expression to a greater extent. c Inositol biosynthesis was dysregulated in collected root RNA samples throughout the course of infection, in which the downstream effects of inositol metabolites can have profound signaling capabilities in hormone regulation, inorganic phosphate signaling, and calcium signaling that can amplify the PI3K cascade. Ultimately, these three pathways converge on the hallmarks of hypersensitive response during plant-pathogen interactions

Discussion

In this study, we documented changes in RSA traits in N. benthamiana upon GFLV infection (Figs. 1 and 2). The average number of N. benthamiana root tips decreased by up to a third and the average root diameter increased over 20% upon infection by wildtype GHu and mutant F13 1E^Pol*/Sd_G802K when compared with roots of mock inoculated plants (Fig. 1a,b). The RNA1 of these two GFLV strains differ by 17% at the nucleotide sequence level [20]; however, the RNA1-encoded protein 1E^Pol*/Sd of both wildtype GHu and mutant F13 1E^Pol*/Sd_G802K carry a lysine in position 802 [21]. In contrast, N. benthamiana plants infected with wildtype F13 and mutant GHu 1E^Pol*/Sd_K802G, both carrying a glycine at position 802 of protein 1E^Pol*/Sd [21], displayed RSA traits resembling those of mock control plants. These findings indicate that a single GFLV amino acid is associated with changes in RSA traits rather than other genetic differences of up to 17% between wildtype F13 and GHu. This is the first report of a single amino acid of a secovirus altering RSA traits of a plant host.

Efforts to correlate the vein clearing phenotype in apical leaves upon infection of N. benthamiana with wildtype GHu to RSA traits showed that the duration of foliar symptoms was significantly positively correlated with root diameter traits, and significantly negatively correlated with branching patterns, root tip number, and network metrics (Fig. 3a, Supplementary Fig. 3). Linear models confirmed these correlated above and below ground symptoms (Fig. 3d,e). The other three strains of GFLV used in this study were asymptomatic in N. benthamiana leaves [20, 21, 23], preventing any comparison between above and below ground features.

The root phenotype was correlated to ROS produced in the pericycle of the primary roots (Fig. 4) and the dysregulation of the host transcriptome (Fig. 5). Abundant ROS production was specific to root systems infected with GFLV strains containing K at residue 1E^Pol*/Sd₈₀₂ (Fig. 4) with transcriptomic analysis confirming an active regulation of these transcripts across viral treatments (Fig. 6). There was a direct correlation between the number of days showing symptoms in apical leaves and the severity of changes in RSA traits for plants infected by wildtype GHu (Fig. 3), while a correlation of ROS production was shown in the roots infected with wildtype GHu and mutant F13 1E^Pol*/Sd_G802K (Fig. 4B). These results suggested a conserved reaction across tissue types that might be required for GFLV systemic infection and perhaps an evolutionarily conserved mechanism that elicits symptoms in N. benthamiana. Previous studies found that additional modifications to the amino acids in this region often reduced pathogenicity or virulence, or completely abolished the infectivity of the virus, indicating this protein’s critical involvement in the infection process [21, 52].

While several new host gene candidates were identified as most differentially abundant in root tissue of GFLV-infected N. benthamiana, general trends remained like those previously identified in apical leaves [23]. The enrichment term of singlet oxygen-mediated programmed cell death (GO:0010343) was present in all GFLV-infected plants tested in this study. This result further supports the hypothesis that all GFLV strains elicit a certain level of ROS with symptomatic strains being the most active. Additionally, five genes specified as PROTEIN OF UNKNOWN FUNCTION were among the top differentially abundant transcripts. Many of these candidate genes could benefit from follow up functional validation analysis in knock-out assays, with confirmation of gene classification and identification if resultant functions are specific to root tissue. The continued functional validation of root specific transcripts will aid in systems biology analyses for plant-pathogen interactions and provide genetic improvement resources for plant biologists [53].

Curated transcriptomics data and our updated annotation informed molecular mechanisms of the resultant N. benthamiana root phenotypes upon GFLV infection (Fig. 7). In roots, upregulation of inositol-tetrakisphophate 1/2 kinase and the antagonistic expression patterns of other inositol synthesis genes prevented the full realization of a mechanism, although an imbalance of this signaling molecule could explain similar events occurring for the transient ROS production in root tissue (Fig. 7). Rather than a 1:1 model, we hypothesize synergistic effects of inositol-mediated signaling, immune receptors, and hormones to contribute to the observed root modifications. Interestingly, the top DEGs return to phosphatidylinositol and extended pathways that relate to the PI3 kinase cascade. This central innate immune response regulator is widely involved in host–pathogen interactions. The top gene according to p-value at 17 dpi when contrasting wildtype GFLV-GHu to mock control is the SEC14P-LIKE PHOSPATIDYLINOSITOL TRANSFER FAMILY PROTEIN, with a Log2FoldChange of 1.98 (Table 2). This gene is a key regulator of many downstream biological processes [54], further supporting the proposed hypersensitive response during GFLV infection. As with the innate mammalian immune system [55], the PI3K cascade was identified as a major factor in bacterial plant-pathogen hypersensitive response [56]. It is notable that genes related to this pathway were documented in both shoot [23] and root tissue (this study), making evident that the combination of signaling pathways are utilized to mount the robust production of ROS (Fig. 7). Future work could address the benefit of ROS accumulation to the host and the pathogen, and how this cascade might have influenced the evolutionary adaptation of the virus in plant hosts.

It is known that plant viruses proficiently manipulate host gene regulatory pathways [57]. Generally, transcriptional and translational regulators can be activated or repressed, directly or indirectly, by RNA/protein binding plant hormones, metabolites, or proteins in plant root tissue. Some factors described in the literature to great extents include abscisic acid, indole 3-butyric acid, ethylene, gibberellic acid, brassinosteroids, nitrogen, glutamate, nitrous oxide, phosphate, sulfate, temperature, and water availability [58,59,60]. The manipulation of one or more of these factors can alter the way in which root systems are formed, usually in a dose-dependent manner. The underlying molecular mechanisms behind abiotic factor changes have been identified in many systems from model plants like Arabidopsis thaliana [61, 62] to perennial crop systems like grape, apple, and orange [60]. A summary of historically understood models expands upon gene regulation, hormone pathways, and nutrient availability to determine root structure [59]. The documentation of modified root morphology upon virus infection overlaps considerably with nutrient deficiencies/toxicities such as phosphorus [63], nitrogen [64], calcium [65], and others. As our study was performed in controlled conditions, the impacts of nutrient, light, or temperature differences were limited. Thus, root traits could be ascribed to a single amino acid of the GFLV-encoded protein 1E^Pol*/Sd. We further hypothesize that the severity of symptoms occurs in a similar manifestation to dose-dependent concentrations of stressors (i.e., secondary signaling molecules, nutrient availability, plant hormones, etc.,), but in this instance to be dependent upon the protein encoding factor(s) of the virus exposed to the host machinery.

GFLV-host interactions documented at the root level by 3’RNA-Seq included multiple overlapping pathways of immune response and hormone regulation. The upregulation of oxidoreductases in roots at 17 dpi indicate lessened ROS levels and a diminished effect on plant response thereafter (Fig. 6). Tolerance of oxidative stress and respective signaling patterns to auxin mediate large morphological changes in plants [66], which might explain the root phenotype in GFLV-infected N. benthamiana. Overall, ROS production due to a hypersensitive response (HR) agrees with current understandings of virus-host interactions at the molecular level [67]. The regulation of HR is not linear and can utilize any combination of RNA silencing, pattern triggered immunity, and effector triggered immunity to elicit such responses [67]. Overall, HR, ROS burst, hormone manipulation are documented as common denominators in symptomatic viral infections [68] but the appearance and disappearance of ROS and symptomology above and below ground depending on a single viral amino acid is unique to GFLV.

Other plant systems have documented changes in RSA upon viral infection. For example, single infections of sweet potato chlorotic stunt virus or multiple infections of sweet potato feathery mottle virus, sweet potato virus 2, sweet potato virus C, and sweet potato virus G decrease lateral root development and adventitious roots of sweet potatoes (Ipomoea batatas) [2]. In a similar study, viral co-infections reduced the total marketable yield between 60–95% depending on the level of tolerance or susceptibility of the sweet potato variety [69]. Additionally, tomato (Solanum lycopersicum L.) infected with tomato brown rugose fruit virus differ in multiple RSA traits when compared with uninfected plants with a 15% reduction in root diameter, and a 30% decrease in root tip number at three weeks post inoculation [3]. Cassava brown streak virus (CBSV) causes root necrosis and tuber yield loss [70]. Almost all (99%) of above ground symptomatic cassava plants display necrotic lesions below ground due to the presence of CBSV, while only 16% of asymptomatic infection events above ground also have root necrosis [71]. Cassava mosaic disease is associated with at least ten different viral species and has been shown to have anywhere from minor changes to up to 87% losses in tuber yield in the most aggressive instances [72]. While cassava and sweet potato are tuber crops, significant impacts to above ground crop production can be associated with total root health, and thereby increasing the importance on documenting the most resilient and sturdy roots for prolonged cropping systems such as with perennial crops. In the case of BNYVV, three viral amino acids are crucial to determine a resistant or a susceptible phenotype in beets, with amino acid 68 of p25 being of critical stature for host-genotype interactions [73].

Our findings on GFLV strain-specific alterations in N. benthamiana RSA have direct implications in transmission biology. The ectoparasitic dagger nematode vector of GFLV, X. index, acquires GFLV in a non-circulative, non-propagative manner by feeding on infected plant roots [74,75,76]. We hypothesize that the phenotype of RSA induced by wildtype GHu and mutant F13 1E^Pol*/Sd_G802K deter continued feeding of X. index on N. benthamiana roots. This is because plants with significantly less root tips and increased root diameter are less appetent to X. index, which prefers fibrous over mature roots for sustained feeding and reproduction [75]. Modified RSA traits of plants infected with wildtype GHu or mutant F13 1E^Pol*/Sd_G802K could push nematode vectors away from an already infected host and increase the likelihood of feeding on neighboring healthy plants at a faster rate than unperturbed root systems of healthy plants or plants infected with other GFLV strains such as wildtype F13 or mutant GHu 1E^Pol*/Sd_K802G. Virus manipulation of plants and vectors to enhance acquisition and/or promote transmission to new plants is well documented for hemipteran-transmitted viruses [77]. It would be interesting to test how modified N. benthamiana RSA traits influence X. index-mediated GFLV transmission in controlled assays with different virus strains. Additionally, if modified RSA traits induced by specific GFLV strains are observed in the natural host Vitis spp., and X. index feeding strategies altered by GFLV-induced changes in root morphology are confirmed, new management strategies of GFLV based on changes in the host rhizosphere could be explored to reduce transmission events.

It is unknown if the GFLV-mediated root phenotype in N. benthamiana fully manifests in the natural host, Vitis spp. In the instance of BNYVV, the transient expression of protein p25 induced similar effects in A. thaliana as observed for the natural host, sugar beet [1, 73]. It would be interesting to investigate if distinct GFLV strains, including some that cause severe foliar symptoms in Vitis spp., affect RSA traits of infected grapevines over the vineyard’s lifespan. Similarly, a correlation analysis of above and below-ground traits in grapevines would be beneficial to better understand the mechanism of disease induced by GFLV. Such work could inform new management strategies for growers experiencing the impacts of fanleaf degeneration.

Conclusions

We documented for the first time the impact of GFLV on RSA traits of infected N. benthamiana. The most notable changes in roots were observed in plants infected with wildtype GHu and mutant F13 1E^Pol*/Sd_G802K, both carrying a lysine at the 802nd position of the putative RNA-dependent RNA polymerase (protein 1E^Pol*/Sd). These two GFLV strains caused a significant reduction in the number of root tips (up to 30%), an increased root average diameter (up to a 20%), and an abundant accumulation of ROS in the pericycle of the primary roots when compared with roots of mock inoculated control plants. The severity of root phenotypic changes was associated with ROS production in the root pericycle. Transcriptomic analysis revealed potential host gene candidates related to viral immune response, RNA silencing, and plant hormone pathways that could contribute to the observed root phenotype. More work is needed to understand any potential downstream impacts of root morphology, particularly on GFLV transmission by the ectoparasitic nematode vector X. index. Finally, unknown root transcripts dysregulated due to a single GFLV amino acid should be validated through reverse genetic assays to improve our understanding of how GFLV manipulates root architecture.