Potato leafroll virus reduces Buchnera aphidocola titer and alters vector transcriptome responses

Viruses in the Luteoviridae family, such as Potato leafroll virus (PLRV), are transmitted by aphids in a circulative and nonpropagative mode. This means the virions enter the aphid body through the gut when they feed from infected plants and then the virions circulate through the hemolymph to enter the salivary glands before being released into the saliva. Although these viruses do not replicate in their insect vectors, previous studies have demonstrated viruliferous aphid behavior is altered and the obligate symbiont of aphids, Buchnera aphidocola, may be involved in transmission. Here we provide the transcriptome of green peach aphids (Myzus persicae) carrying PLRV and virus-free control aphids using Illumina sequencing. Over 150 million paired-end reads were obtained through Illumina sequencing, with an average of 19 million reads per library. The comparative analysis identified 134 differentially expressed genes (DEGs) between the M. persicae transcriptomes, including 64 and 70 genes that were up- and down-regulated in aphids carrying PLRV, respectively. Using functional classification in the GO databases, 80 of the DEGs were assigned to 391 functional subcategories at category level 2. The most highly up-regulated genes in aphids carrying PLRV were cytochrome p450s, genes related to cuticle production, and genes related to development, while genes related to heat shock proteins, histones, and histone modification were the most down-regulated. PLRV aphids had reduced Buchnera titer and lower abundance of several Buchnera transcripts related to stress responses and metabolism. These results suggest carrying PLRV may reduce both aphid and Buchnera genes in response to stress. This work provides valuable basis for further investigation into the complicated mechanisms of circulative and nonpropagative transmission.


Methods
Plant and insect growth conditions. Solanum tuberosum were propagated using leaf-bud cutting from cv. Désirée 39 in laboratory experiments. Plants were grown in growth chambers under controlled conditions (25/23 °C day/night with a photoperiod of 16/8 h day/night). Non-viruliferous and viruliferous aphid clones of a potato-adapted red strain of Myzus persicae were reared under controlled conditions (25/23 °C day/night with a photoperiod of 16/8 h day/night) on healthy potato. We confirmed our colony was free of secondary symbionts using universal primers for the bacteria 16S-23S ribosomal RNA intergenic spacer region (Supplementary Table S1). All experiments were conducted in the same environmental chambers and conditions, so there were no environmental differences in treatments (25/23 °C day/night with a photoperiod of 16/8 h day/night). 40 was grown at 28 °C in LB broth (+ 10 mM MgSO 4 ), with kanamycin (50 µg/mL), carbenicillin (100 µg/mL) and rifampicin (50 µg/mL) for selection. After 24 h, bacteria were centrifuged to concentrate and resuspended in 10 mM MgCl 2 . One-week-old S. tuberosum were inoculated at an optical density (OD) of 0.70. Three weeks post infection, tissue was collected from all plants, RNA was extracted using the SV Total Isolation Kit as per manufacturer's instructions (Promega, Madison, WI, USA), and cDNA was synthesized using 1500 ng of total RNA and random hexamers (20 ng/µL) with the SMART ® MMLV as per manufacturer's instructions (Takara Bio USA, Mountain view, CA, USA). cDNAs were used in PCRs with PLRV specific primers (F-5' ATG AGT ACG GTC GTG GTT -3' and R-5 'CTA TTT GGG GTT TTG CAA AGC-3'). A set of uninfected potato cuttings were grown at the same time as the plants above to serve as controls. After systemic plant infection was verified plants were immediately used in experiments.

Pathogen infection. Agrobacterium tumefaciens (LBA4404) containing the infectious clone of PLRV
RNAseq, qRT-PCR, and qPCR aphid experiments. One week after infection was verified, five adult aphids were placed on the first fully expanded leaflet of three infected and three healthy plants. After 24 h, all adults were removed and 20 larvae were left to develop. Seven days later all aphids were at the same developmental stage. Ten adult aphids were collected into a tube from each plant (N = 3 plants with 10 aphids per plant) and immediately frozen in liquid nitrogen until use in RNAseq experiments. The entire experiment was repeated a second time for confirmation of RNAseq results using qRT-PCR and to examine the titer of the bacterial symbiont, Buchnera, in aphids. For this experiment 5 aphids were collected for RNA extraction and 5 aphids were collected for DNA extractions from each plant. We also prepared 6-7 plants for each treatment instead of 3 (N = 6-7 with 5 aphids per plant), however all other methods were the same.
RNA and DNA isolation from aphids. RNA was extracted from aphid tissue collected in the first two experiments as described above. The RNA concentration and purity were measured using a NanoDrop. The integrity of RNA was confirmed using the Bioanalyzer 2100 system (Agilent, Santa Clara, CA, USA). DNA was extracted from aphid tissue collected in the second experiment using cetyl trimethylammonium bromide. The integrity of DNA was confirmed using an agarose gel. The DNA concentration and purity were measured using a NanoDrop (Thermo Fisher Scientific, Waltham, MA, USA).
Library preparation, and sequencing. Sequencing libraries were prepared using a multiplexing library protocol 41 . Briefly, oligo(dT) 25 Dynabeads were used to purify mRNA, which was then fragmented, and the first-strand cDNA was synthesized using random primers, dNTP, and reverse transcriptase. The second-strand was synthesized using a dUTP mix, DNA Polymerase I, and RNase H, ends repaired, and adenylated. The cDNA fragments were ligated to adapters, selectively enriched by PCR, and purified using the AMPure XP beads. The library quality was assessed using the Agilent Bioanalyzer 2100 system and sequenced using an Illumina HiSeq 2000 instrument.  Table S1). qRT-PCR reactions were carried out using SYBR green PCR master mix (Applied Biosystems, Carlsbad, CA, USA), in an CFX384 instrument (Bio-Rad, Hercules, CA, USA). Three technical replicates were performed for each individual sample, and a digital pipette was used for all pipetting. Relative transcript abundance was calculated utilizing a standard curve produced from tenfold series dilution of cDNA synthesized from 1000 ng/μL of total RNA according to the standard curve method (Applied Biosystems, Carlsbad, CA, USA). Technical replicates of raw CT values were averaged and transcripts of interest were normalized to the house-keeping transcript ribosomal protein L7 rpl7 for aphids or the 50S ribosomal subunit gene rpIN for Buchnera, as previously described 51-53 . Buchnera titer. Buchnera titer here is defined as the ratio of Buchnera single copy genes to aphid single copy genes. To determine Buchnera titer in whole aphid bodies we used qPCR and measured the ratio of a single copy Buchnera gene (rplN) to a single copy aphid gene (RPL7). qPCR reactions were carried out using SYBR green PCR master mix (Applied Biosystems, Carlsbad, CA, USA), in the CFX384 instrument (Bio-Rad, Hercules, CA, USA). Reactions were performed in triplicate for each sample, and the average was used for analysis. Relative abundance was calculated utilizing a standard curve produced from tenfold serial dilution of DNA.

Statistical analyses.
RNAseq data analyses were performed as described above. All statistical analyses for qRT-PCR and qPCR were determined using a P < 0.05. Analysis of variance (ANOVA) was used to determine significant difference in transcript abundance. For Buchnera titer single factor ANOVA was used to determine difference in relative abundance. The statistical analyses were performed using JMP 8 software (SAS Institute, Cary, NC, USA).
Compliance statement. All research in this paper complied with relevant institutional, national, and international guidelines and legislation.

Results
Differential gene expression in the presence of PLRV. Over 150 million paired-end reads were obtained through Illumina sequencing, with an average of 125,815, 247 100-bp reads per library (Table 1, Fig. 1A). A higher number of average reads were sequenced in control 1 compared to other samples. Nevertheless, reads mapped to the target genome in the same relative proportions across treatments (Fig. 1B) revealing that the sample was not contaminated with reads other than for the target genome. Further, read counts were normalized before differential gene analyses (see methods) accounting for variation in library sizes among treatments. About 73% of the reads mapped to the M. persicae reference genome, with around only 27% being uniquely mapped (Table 1, Fig. 1B). In order to examine biological variability, a principal component analysis (PCA) of the normalized count data was performed (Fig. 1C). The first component of variance separated samples by treatments and accounted for 54% of the variance. Hierarchical clustering confirmed PCA results in visual representation of DEG expression (Fig. 1D). The transcriptome of viruliferous M. persicae was compared to the www.nature.com/scientificreports/  Functional roles of differentially expressed genes. Gene ontology (GO) enrichment analyses were performed with the DEGs from each treatment to identify functions and pathways disturbed in aphids carrying PLRV. One or more gene ontology terms were assigned to each transcript from biological processes, molecular functions, and cellular compartments term using Blast2GO functional gene annotation 49 . The 134 DEGs were assigned to functional GO terms within the three categories, including 125 biological processes, 118 molecular functions, and 148 cellular compartments. Of the 134 DEGs, 53 (39.55% of total DEGs) were classified as "uncharacterized proteins. " The majority of DEGs assigned to biological processes were categorized as metabolic processes (41%), cellular-protein processes (11%), and oxidation-reduction processes (11%) (Fig. 2A). As for DEGs assigned to molecular functions, almost half were associated with catalytic activity (33%) or nucleic acid binding (21%) (Fig. 2B). Within the cellular component category, 24% were related to the membrane and 24% were related to intracellular locations (Fig. 2C).
Next each DEGs was annotated using a single Blast2GO consensus description. Many of the genes upregulated in PLRV aphids were related to cuticle formation and development (16%), and catalytic activity (16%), however the majority of up-regulated transcripts were uncharacterized (31%; Fig. 3A). The largest groups of down-regulated genes in PLRV aphids were related to histones (10%), catalytic activity (10%), transmembrane transport (9%), proteolysis or protein ubiquitination (7%), and nucleic acid binding and metabolic processes (7%; Fig. 3B). A significant proportion of the down-regulated transcripts in PLRV aphids were also uncategorized (47%). The most highly expressed DEGs included transcripts related to cuticle formation and development, and 4C1-like cytochrome P450s (Table 2, Supplemental Table 2). The most down-regulated transcripts in PLRV aphids were related to histones and histone modifying proteins ( Table 3).

Validation of select aphid transcripts via RT-qPCR.
To validate the RNAseq we conducted qRT-PCR analysis of select DEGs in a separate experiment using the same experimental design. We selected one of the most up-regulated genes (MYZPE13164_G006_v1.0_ 000133030.2 (Mpcp5-like)) and one of the most downregulated gene (MYZPE13164_G006_v1.0_ 000070430.1 (Hsp68-like)) and measured transcript abundance using qRT-PCR (bolded genes in Tables 2 and 3). Mpcp5-like is related to cuticle formation and development and Hsp68-like encodes a heat shock protein, which are related to immunity and stress responses. Consistent with the RNA-seq data, abundance of the Mpcp5-like transcript was significantly higher in viruliferous M. persicae compared to virus-free controls (10.045, 2.017, relative expression respectively; p = 0.019; Fig. 4A). Abundance of the Hsp68-like transcript was significantly lower in viruliferous M. persicae when compared to virus-free controls (1.44, 5.51, relative expression respectively; p < 0.01) (Fig. 4A,B).
The impact of PLRV on Buchnera aphidicola titer. Buchnera has been previously implicated in transmission of PLRV and other luteoviruses 23,25,27,31,54 , however Buchnera titer and coding sequence transcripts have not been examined in aphids carrying PLRV. From our experiments, Buchnera titer was ~ 1.5 times higher for virus-free aphids compared to aphids carrying PLRV (ratios 6.42, 4.20 respectively; p = 0.037; Fig. 5A). To investigate the potential mechanisms mediating decreases in Buchnera titer we measured abundance of two transcripts related to stress, dnaK 55,56 and groEL 56,57 , and one transcript related to metabolism, argE 58 . Abundance of all three transcripts were reduced in aphids carrying PLRV compared to controls. Viruliferous aphids had 63.32% less argE transcripts (p = 0.026), 83.33% less groEL transcripts (p = 0.024), and 81.23% less dnaK transcripts (p = 0.046) compared to that of the virus free aphids (Fig. 5B-D).

Discussion
The main focus of this paper was to examine the effect that Potato leafroll virus has on the transcriptome of M. persicae, and their primary endosymbiont Buchnera aphidicola. The largest category of known up-regulated transcripts in viruliferous aphids compared to controls were related to the cuticle and cuticle development. Insect cuticles are largely composed of a protein matrix embedded with chitin filaments 59 . Cuticle proteins (CPs) have been shown to be involved in general development, molting, transmission of non-persistent viruses, and insecticide resistance through changes in cuticle permeability 30,[60][61][62][63] . In Acyrthosiphon pisum, 19 CPs were found to be regulated by photoperiodism and suspected to be involved in the transition from asexual to sexual production 64 . Further, cuticle proteins have been implicated as potentially facilitating transmission of Barley yellow dwarf virus (BYDV-GPV), Cereal yellow dwarf virus (CYDV-RPV), and Turnip yellows virus (TuYV), three related Luteoviridae viruses [65][66][67] . Whilst we cannot know the function of changes in CP transcripts in PLRV-aphid interactions from these experiments, these genes represent promising targets for further investigation.
In addition to many cuticle related proteins, five cytochrome P450s genes were significantly up-regulated in viruliferous aphids compared to controls. Cytochrome P450s play important roles in hormone and pheromone metabolism but are more famous for their roles in the metabolism of insecticides and host plant chemicals. Polyphagous insects, like M. persicae, encounter many different hosts and tend to have a higher number of P450s related to the metabolism of allelochemicals compared to more specialized aphids 68 . Previous work has shown that a cytochrome 450 gene (CYP6CY3) was found to increase nicotine tolerance and aphid host adaptation 69 www.nature.com/scientificreports/ It has been previously hypothesized that up-regulation of p450s could help insect vectors tolerate less desirable hosts which could be beneficial to the virus 71 .
Transcripts encoding a heat shock protein (HSP68-like) was among the most down-regulated in viruliferous aphids compared to controls. HSP68 is a member of the HSP70 family, which are important chaperone proteins that are known to be up-regulated in response to stress. One study found that the Hsp70 from Bemisia tabaci is up-regulated after acquisition of Tomato yellow leaf curl virus (TYLCV) 72 . They went on to show that HSP70 www.nature.com/scientificreports/ protein can directly interact with TYLCV using in vitro studies and that they co-localize together in insect midgut cells using in situ hybridization. The authors suggest HSP70 may play an inhibitory role in virus transmission, as transmission was increased when whiteflies were fed HSP70 antibodies. Given Hsp68-like transcripts were down-regulated in aphids carrying PLRV in our study, it would be interesting to investigate if this has any impact on PLRV transmission. Porras et al. 73 demonstrated that BYDV-PAV, a strain that is only transmitted by Rhopalosiphum padi (bird-cherry oat aphid), up-regulated the abundance of three Hsp70 transcripts in the aphid vector. The authors found BYDV infection increases plant surface temperature and aphid heat tolerance, suggesting a protective role of HSP70 proteins in virus-aphid-plant interactions 73 . Although it is not known if PLRV increases plant surface temperature and vector heat tolerance, it has been shown that potato plants kept at higher temperatures are more susceptible to PLRV than compared to lower temperatures 74 . Also aphid www.nature.com/scientificreports/ acquisition and transmission at higher temperatures resulted in higher transmission rates compared to lower temperatures 75 , however at very high temperatures differences were reduced 76 . It is not known how decreases in Hsp68-like transcripts in aphids carrying PLRV may alter aphid heat tolerance. In this study there was a significant reduction of Buchnera titer and Buchnera gene expression of three genes (dnaK, groEL, and argE) in aphids carrying PLRV compared to control aphids. In general, gene regulation at the mRNA level in Buchnera is thought to be minimal because Buchnera transcription factors are reduced 58 and very few transcriptional responses had been observed previously 77 . Only two transcription initiation factors ( σ 32 and σ 70), the heat shock and housekeeping transcription factors, respectively, remain in Buchnera Myzus's genomes 78 similar to other Buchnera taxa 79,80 . The housekeeping sigma factor ( σ 70) initiates transcription of argE which is regulated by the repressor ArgR when bound to arginine in Escherichia coli 81 . Similar to other Buchnera taxa, Buchnera Myzus's genome 78 has lost the repressor ArgR so it is unclear how this gene is down-regulated in virus-infected aphids compared to un-infected aphids. The other two Buchnera genes (dnaK and groEL) that were down-regulated in this study in aphids carrying PLRV compared to control aphids are associated with the heat shock regulon 80 . Moreover, these Buchnera genes still retain recognizable σ 32 promoter sites up-stream of dnaK and groEL in the Myzus Buchnera G006 genome (NCBI Reference Sequence: NZ_MJNC01000001; Supplemental Table 3) similar to other Buchnera taxa 80 . The σ 32 heat shock response is highly conserved in bacteria and is initiated in response to stress, such as heat shock or other environmental stressors that destabilize proteins [78][79][80][81] . In this study it is unclear how PLRV is either directly or indirectly dampening Buchnera's expression of dnaK and groEL and if it is through a similar mechanism that is also down-regulating the aphid's stress response genes including Hsp70.
A decrease in Buchnera titer has previously been associated with different aphid clones 83 , plant diets 13 , increasing aphid nymphal age 82,84 , and heat shock 81,84,86 . Most obligate pathogens and symbionts, including Buchnera, overexpress the protein GroEL during non-heat shock conditions to rescue misfolded proteins 87 . We hypothesize that PLRV is reducing Buchnera's ability to up-regulate genes that are associated with the heat shock regulon (Fig. 5C,D) and this may lead to increased stress, lysing of Buchnera cells, and ultimately a reduction of Buchnera titer (Fig. 5A). Other insect-plant pathogen systems are known to modulate obligate symbiont titer. For example, in whiteflies Portiera titer is modulated by the co-occurrence of its facultative symbiont Rickettsia and TYLCV 88 . Alternatively, as PLRV-infected plants have higher concentrations of free amino acids 38 , the change in host plant diet may have influenced Buchnera titer similar to Zhang et al. 13 . Because aphids carrying PLRV www.nature.com/scientificreports/  A predicted heat shock protein (Hsp68-like) was significantly down-regulated in expression in individuals with PLRV compared to controls. Transcripts were measured relative to a housekeeping gene RPL7. Significant differences were calculated using an ANOVA (*P < 0.05; Error bars represent ± SEM). www.nature.com/scientificreports/ may obtain higher levels of essential amino acids from virus-infected plants, Buchnera genes that are involved in arginine biosynthesis, such as argE, may be down-regulated compared to aphids feeding on un-infected plants with lower amounts of essential amino acids. Parasites of plants and animals can modify host behavior to improve their own transmission and survival 6,7,89,90 . This work explores the complex relationships that exist between hosts, viruses, vectors, and endosymbionts, and opens up more questions regarding the complexity and depth of these relationships. Aphids and bacterial endosymbionts may benefit from relationships with plant-infecting viruses indirectly or directly but additional studies are needed. Although it is known that Buchnera titer and gene expression responses vary with aphid linages 91 , it is not known how this is impacted by long term associations with plant-infecting viruses. In regions where virus pressure is high or where poor hosts dominate, aphids may more often be associated with plant infecting viruses. This study also expands on previous work that given the mounting evidence of virus manipulation of insect vectors, this could have lasting impacts on the population structures of these insect vectors and their obligate endosymbiont.  www.nature.com/scientificreports/