Mal de Río Cuarto virus infection causes hormone imbalance and sugar accumulation in wheat leaves

Effect of MRCV-infection on the wheat transcriptome

Controlled infection experiments were performed to compare the wheat transcriptome between MRCV-infected and mock-inoculated plants. We were particularly interested in understanding the changes that may contribute to the establishment of a systemic infection rather than the late, pleiotropic effects of the infection. Therefore, the youngest fully expanded systemic leaf of each plant was sampled at 12 and 21 dpi. Plants at 12 dpi are asymptomatic, whereas the latter time point coincides with the appearance of macroscopic symptoms (Additional file 1: Figure S1). For each time point, four infected plants with similar viral loads and four mock-inoculated plants treated with non-viruliferous planthoppers were selected. We extracted total RNA from the individual leaf samples, constructed RNA-seq libraries and sequenced them using Illumina technology. On average, 50.0 million reads per library were obtained, 98.2% of which were retained after adapter and quality trimming (Additional file 1: Figure S2). Reads were mapped to the wheat genome (assembly version TGACv1) and, using only uniquely mapped reads (~ 43% of the total), we identified DATs between infected and control plants at each time point. At 12 dpi, before the appearance of disease symptoms, we identified just two DATs (FDR-adjusted P-value< 0.05, Additional file 2: Table S1). This result indicates that MRCV infection has almost no impact on the wheat transcriptome at this early time point. In contrast, at 21 dpi, following the onset of visual disease symptoms, we identified 3233 DATs, 42% with lower and 58% with higher transcript levels in infected plants compared to mock-inoculated samples (Additional file 2: Table S1). Twelve DATs at 21 dpi were randomly selected to validate the RNA-seq by qRT-PCR. We found consistent results using both technologies (Student’s t Test, P < 0.05), indicating that RNA-seq data is reliable (Additional file 1: Figure S3).

Using a set of annotated wheat transcription factors (TFs) [27], we identified 180 TFs from 33 different families among DATs at 21 dpi (Fig. 1 and Additional file 2: Table S2). Although 58% of the total DATs showed increased accumulation in infected plants, 68.9% of the transcription factor DATs showed reduced levels upon the infection. TFs from the basic helix-loop-helix (bHLH), MYB-related, basic leucine zipper (bZIP), NAC and WRKY families were the most represented. When analyzing the proportion of altered TFs within each family, CONSTANS-like (CO-like), WHIRLY, DNA BINDING WITH ONE FINGER (Dof), TEOSINTE-BRANCHED1/CYCLOIDEA/PCF (TCP) and VASCULAR PLANT ONE-ZINC FINGER (VOZ) families were most affected by MRCV infection (Fig. 1). Some of the TFs found to be altered upon MRCV infection are involved in hormone responses and are discussed below.

Effect of MRCV infection on primary metabolism and sugar allocation

Viral replication demands the biosynthesis of a great amount of new molecules in infected tissues and thus has a strong impact on primary metabolism [23]. At the transcript level, we identified 34 DATs involved in carbohydrate synthesis, degradation and transport, according to their associated GO terms and Ensembl annotations (Additional file 2: Table S3). Four sucrose synthase transcripts were significantly increased, including four transcripts with homology to Arabidopsis sucrose synthase 5 and 1. In addition, 12 transcripts coding for invertases and/or fructosyltransferases (evolutionarily related to vacuolar invertases) were altered at 21 dpi (Additional file 2: Table S3). Vacuolar invertases and fructosyltransferases directly affect the levels of soluble sugars and may thus contribute to the generation of turgor, which further drives cell expansion [28]. Since wheat invertases and fructosyltransferases are not fully annotated in the TGACv1 assembly version of the wheat genome, we built a phylogenetic tree including sequences of previously reported homologs to better classify them [29, 30]. We identified four significantly decreased DATs coding for vacuolar invertases and eight significantly increased DATs encoding fructosyltransferases (1 sucrose:sucrose 1-fructosyltransferase, 1-SST; 4 fructan:fructan 1-fructosyltransferases, 1-FFT; and 3 sucrose:fructan 6-fructosyltransferases, 6-SFT). In line with a possible accumulation of fructans in the infected plants, a fructan 6-exohydrolase involved in fructan degradation was decreased upon the infection (Additional file 1: Figure S4, Additional file 2: Table S3). In addition, 12 transcripts encoding homologs of enzymes involved in starch metabolism (2 glucoamylases, 2 granule-bound starch synthases, 2 branching enzymes, and 6 alpha-glucan phosphorylases) were significantly increased, while four β-amylase transcripts were significantly decreased and one was increased at 21 dpi (Additional file 2: Table S3). We also found evidence of alterations in cell wall remodeling transcripts. The GO categories of cellulose biosynthetic and catabolic processes were over-represented among transcripts with increased levels at 21 dpi (Table 1). These categories include several cellulose synthase, endoglucanase and xyloglucan endotransglucosylase/hydrolase transcripts (Additional file 2: Table S3).

The changes observed in transcripts coding for enzymes involved in primary metabolism prompted us to quantify soluble sugars and starch. MRCV-infected plants showed a significant increase in glucose, fructose, sucrose and starch at 21 dpi (Fig. 2, P < 0.05). In agreement with the gene expression analysis, we found no significant differences for any of these metabolites between MRCV-infected and control plants at 12 dpi. To further investigate the metabolic status of MRCV-infected plants, we used liquid chromatography coupled to tandem mass spectrometry (LC-MS/MS) to determine the levels of phosphorylated intermediates and organic acids. No differences were observed for any of these metabolites at 12 dpi; conversely, 10 out of 23 metabolites were significantly accumulated in MRCV-infected leaves at 21 dpi (Fig. 3, P < 0.05).

Interestingly, the signal metabolite Tre6P was 5-fold more abundant in MRCV-infected leaves at 21 dpi than in controls (Fig. 3). Tre6P is both a signal and a negative regulator of sucrose levels in plants, thus linking growth and development to carbon status [31, 32]. In illuminated Arabidopsis leaves, Tre6P promotes the synthesis of organic and amino acids by post-translational activation of phosphoenolpyruvate (PEP) carboxylase and nitrate reductase [33]. Similarly, the levels of Tre6P and certain amino acids concomitantly increased in axillary buds from decapitated pea plants [34]. In line with these findings, we observed accumulation of several tricarboxylic acid (TCA) cycle intermediates (citrate, iso-citrate, 2-oxoglutarate and malate), PEP, shikimate and total amino acids in MRCV-infected plants at 21 dpi (Figs. 2 and 3). Total protein levels were also increased, whereas nitrate content was not significantly impacted by MRCV infection (Fig. 2).

The GO category trehalose biosynthesis was over-represented among transcripts with decreased levels at 21 dpi (Table 1). Previous phylogenetic analyses of plant TPS proteins showed that they cluster in two distinct sub-families [35, 36]. Class I TPS proteins are catalytically active, whereas the function of class II TPS proteins remains unknown [32, 37, 38]. Because wheat TPS genes are not fully annotated in the TGACv1 assembly version of the wheat genome, we built a phylogenetic tree including 21 putative wheat TPS proteins (Additional file 2: Table S4) and previously reported TPS sequences from Arabidopsis, rice, poplar and common bean [35, 36, 39]. Additional file 1: Figure S5 shows that 11 DATs coding for class II TPSs (TaTPS2, TaTPS5, TaTPS6 and TaTPS10) were significantly decreased, whereas two transcripts encoding TaTPS4 (another class II TPS) were increased.

The GO category carbohydrate transport was over-represented among DATs with decreased levels at 21 dpi (Table 1). We considered this interesting since it has been shown that bacterial and fungal plant pathogens are able to hijack sucrose transporters [22, 40, 41]. In Arabidopsis leaves, sucrose efflux from phloem parenchyma cells to the apoplast is driven by AtSWEET11 and 12 [42, 43]. Interestingly, we found 13 altered SWEET transcripts in MRCV-infected leaves at 21 dpi (Additional file 2: Table S3). To further characterize these transcripts, we searched for sequences coding for SWEET transporters in the wheat genome and found 90 putative members for this family (Additional file 2: Table S4). Next, we built a phylogenetic tree including sequences from Arabidopsis and rice [40, 44]. The levels of nine transcripts coding for TaSWEET13 were decreased in MRCV-infected plants. The proteins encoded by these transcripts cluster with AtSWEET11/12 and OsSWEET14 (Additional file 1: Figure S6), which are sucrose exporters located in phloem parenchyma cells [42].

Amino acid metabolism is also modulated by viral infections [23]. Interestingly, we observed a significant accumulation of amino acids and total proteins upon MRCV infection (Fig. 2, P < 0.05). The levels of several transcripts encoding enzymes involved in amino acid synthesis were altered at 21 dpi (22 increased and 6 decreased; Additional file 2: Table S5). The recent identification of wheat amino acid transporters [45] allowed us to detect 18 transcripts significantly altered upon MRCV infection (15 decreased and 3 increased; Additional file 2: Table S5), suggesting that amino acid transport is also affected by MRCV infection. Overall, our results show that MRCV infection has a profound effect on carbon partitioning and amino acid metabolism.

Effect of MRCV infection on phytohormone metabolism and transport

Changes in hormone metabolism, transport and signalling are frequently associated with virus symptoms. Thus, we analyzed transcripts involved in phytohormone synthesis, catabolism, transport and signalling that were significantly altered in MRCV-infected plants at 21 dpi. We also quantified several hormones including their biosynthetic precursors and metabolites by ultra-high performance liquid chromatography–electrospray tandem mass spectrometry (UHPLC–MS/MS). The transcripts associated with the different hormones are listed and classified in Additional file 2: Table S6, while Fig. 4 shows the fold change of hormone levels between infected and mock-inoculated plants. Each hormone is discussed separately in the sections below.

Plant inoculation and RNA sequencing

Six-day-old wheat (T. aestivum cv. ProINTA Federal) seedlings were infected with MRCV as previously described [14]. This cultivar was selected because it is susceptible to MRCV infection [4]. Briefly, second instar Delphacodes kuscheli (delphacid planthopper vector) nymphs fed on MRCV-infected wheat plants were used for 1:1 transmission assays. Nymphs fed on non-infected wheat were used as control treatment. Plants were grown in 750 cm³ pots with soil sterilized by solarization in greenhouses with temperature-controlled conditions (24 ± 3 °C) in 16h light / 8h dark photoperiod and daily irrigation with no fertilization. The plants were rotated regularly within the greenhouse to reduce any positional effects. The youngest fully-expanded leaf was collected at 12 and 21 dpi, frozen in liquid nitrogen and kept at − 80 °C until use. Wheat plants were identified as MRCV infected or non-infected by DAS-ELISA at 50 dpi [4]. Absolute viral load was quantified as previously described [14, 120] and plants with similar viral load (log₁₀(MRCV-S3 molecules)/(μg of total RNA) = 20.60 ± 0.96) for each time point were selected. To avoid leaf oxidative stress caused by ovipositions, only plants infected by male delphacids were used. Total RNA samples from mock-inoculated and MRCV-infected plants were extracted at 12 and 21 dpi (n = 4, a total of 16 individual samples) using a mirVana kit (Thermo Fisher Scientific Inc., Waltham, MA, USA) following the total RNA extraction protocol. cDNA libraries were synthesized using an Illumina TruSeq RNA sample prep kit (Illumina, San Diego, CA, USA). Library quality was checked using a High Sensitivity DNA chip in a Bioanalyzer 2100 kit (Agilent Technologies, Santa Clara, CA, USA). Libraries were sequenced on the HiSeq 3000 platform (Illumina, San Diego, CA, USA) at the UC Davis Genome Center using a single-end 50 bp (SE 50) module. Raw sequences of the 16 libraries were deposited in NCBI Sequence Read Archive (SRA) under the accession SRP160433.

Sequence data analysis

Raw reads were processed as previously described [121]. Briefly, “Scythe” (https://github.com/vsbuffalo/scythe) was used to remove adapter contamination (default options) and “Sickle” (https://github.com/najoshi/sickle) to remove low-quality reads (default options except –l 25 –q 25). Trimmed reads were mapped to the Chinese Spring (CS42) TGACv1 genome assembly [12] using GSNAPl (version 2017-09-11, default parameters except -m 2 -n 1 -N 1 -A sam) [122]. Raw counts for gene features were obtained using ht-seq count (default parameters except -m union -a 30) [123] using the GFF file available in Ensembl plants release 37. A threshold for mapping quality was selected so that only uniquely-mapped reads were retained, ensuring homoeologue-specific expression profiles [124]. DESeq2 version 1.17.39 [125] and R version 3.4.0 were used for the normalization and classification of differentially accumulated transcripts (DATs) between infected and mock-inoculated samples (FDR-adjusted P-value ≤0.05, Additional file 2: Table S1).

Annotation data were retrieved from Ensembl plants release 37 using BioMart [126, 127]. GO annotation was used for functional enrichment analysis using the ‘R’ package TopGO version 2.29.0. “Biological Process” ontology was used and the significance values were calculated using ‘classic’ Fishers’ exact test (P < 0.01). When detailed annotation data were not available for T. aestivum, data from the closest homolog of Aegilops tauschii (goatgrass, the diploid progenitor of the D-genome of hexaploid wheat) was used. A reference-quality genome sequence was recently reported with more than 76% of genes annotated [128].

qRT-PCR validations

One μg of DNase-treated total RNA, SuperScript III Reverse transcriptase (Invitrogen, Carlsbad, CA, USA) and random primers were used to synthesize the first-strand of cDNA, according to the manufacturer’s instructions. Twelve DATs were randomly selected and the respective oligonucleotides (listed in Additional file 2: Table S7) were designed using Primer3 [129]. SYBR green PCR reactions were carried out in 20 μl final volume containing: 3 mM MgCl₂, 200 nM each primer, 0.2 mM dNTPs, 1 U Platinum Taq DNA Polymerase (Invitrogen, Carlsbad, CA, USA), 20 mM Tris-HCl (pH 8.4), 50 mM KCl. Cycling conditions were as follows: 95 °C for 5 min, 40 cycles at 95 °C for 15 s, 60 °C for 60 s and a final extension at 72 °C for 30 s, in an Applied Biosystems Step One Plus Real Time PCR system. qRT-PCR data analysis and primer efficiencies were obtained using LinRegPCR software [130]. The genes coding for GTP-BINDING PROTEIN (GTPB, Genbank ID: JQ673330.1) and ELONGATION FACTOR 1 ALPHA (EF1α, Genbank ID: AK455225.1) were used as reference [131]. Relative expression ratios and statistical analysis were performed using fgStatistics software interface [132]. Six biological replicates (the same four as in the RNA-seq plus two additional replicates) and two technical replicates were performed. Statistical significance between RNA-seq and qPCR Fold Change values was assessed by Two-tailed Student’s T tests.

Extraction and analysis of metabolites

Soluble sugars (glucose, fructose and sucrose) were extracted from the same previously described leaf samples of five individual plants (the same four as in the RNA-seq plus one additional replicate) using ethanol and enzymatically quantified [133]. Total amino acids and nitrate were determined using the same extract, while starch and proteins were assayed in the insoluble material obtained after the ethanolic extraction, following the protocols described by Cross et al. [134]. Tre6P, phosphorylated intermediates and organic acids were extracted with chloroform/methanol, and measured by high-performance anion-exchange LC-MS/MS as described by Lunn et al. [135], with the modifications introduced by Figueroa et al. [33].

Phylogenetic analysis

Protein sequences for TPS, SWEET, invertases and fructosyltransferases from Triticum aestivum (Ta), the progenitors of the A and D subgenomes Triticum urartu (Tu) and Aegilops tauschii (AEGt), Triticum turgidum (Tt), Aegilops searsii (As), Oryza sativa (Os), Arabidopsis thaliana (At), Solanum lycopersicum (Sl) and Populus trichocarpa (Pt) were downloaded from the NCBI database (https://www.ncbi.nlm.nih.gov/). Protein sequences from wheat were selected using the BLASTP algorithm in Ensembl Plants (http://plants.ensembl.org/index.html). Briefly, representative sequences from rice (for TPS, SWEET and invertases) or wheat (for fructosyltransferases) were used to search for homologs encoded in the T. aestivum genome (assembly TGACv1). The resulting sequences were downloaded from UniProt (http://www.uniprot.org/), manually curated and automatically aligned with those obtained for other species using the ClustalW program in GenomeNet (http://www.genome.jp/). Phylogenetic trees were built using the neighbor-joining algorithm (bootstrap of 1000) in Seaview 4.3.0 [136] and figures were prepared with the FigTree 1.3.1 program (http://tree.bio.ed.ac.uk/). Sequences used in this study are listed in Additional file 2: Table S4.

Quantitative analysis of plant hormones

The endogenous phytohormones concentrations were determined in leaf samples by an ultra-high performance liquid chromatography–electrospray tandem mass spectrometry (UHPLC–MS/MS) using stable isotope-labelled internal standards as a reference. The isolation of CKs was performed according to Antoniadi et al. [137] and their levels were quantified by UHPLC–MS/MS [138]. Briefly, samples (20 mg FW) were extracted in 1 ml of modified Bieleski buffer [139] together with a cocktail of stable isotope-labeled internal standards used as a reference (0.25 pmol of CK bases, ribosides, N-glucosides, and 0.5 pmol of CK O-glucosides, nucleotides per sample added). The extracts were purified using the Oasis MCX column (30 mg/1 ml, Waters) and cytokinin levels were determined using the LC-MS/MS system consisting of an ACQUITY UPLC System and a Xevo TQ-S triple quadrupole mass spectrometer (Waters). The acidic phytohormones (JA; jasmonoyl-L-isoleucine, JA-Ile; cis-12-oxo-phytodienoic acid, cis-OPDA; indole-3-acetic acid, IAA; ABA; and SA) were extracted using an aqueous solution of methanol (10% MeOH/H₂O, v/v) [140]. A cocktail of stable isotope-labelled standards was added as follows: 5 pmol of [¹³C₆]IAA, 10 pmol of [²H₆]JA, [²H₂]JA-Ile, and [²H₆]ABA, and 20 pmol of [²H₄]SA and [²H₅]OPDA (all from Olchemim Ltd., Czech Republic) per sample. The extracts were purified using Oasis HLB columns (30 mg/1 ml, Waters), targeted analytes were eluted using 80% MeOH and then analyzed by LC-MS/MS method [140]. The brassinosteroid content in the samples was determined using the previously described UHPLC–MS/MS method [141]. All experiments were performed in five independent biological replicates (the same four as in the RNA-seq plus one additional replicate).