Proteomic and Transcriptomic Analysis for Identification of Endosymbiotic Bacteria Associated with BYDV Transmission Efficiency by Sitobion miscanthi

Sitobion miscanthi, an important viral vector of barley yellow dwarf virus (BYDV), is also symbiotically associated with endosymbionts, but little is known about the interactions between endosymbionts, aphid and BYDV. Therefore, two aphids’ geographic populations, differing in their BYDV transmission efficiency, after characterizing their endosymbionts, were treated with antibiotics to investigate how changes in the composition of their endosymbiont population affected BYDV transmission efficiency. After antibiotic treatment, Rickettsia was eliminated from two geographic populations. BYDV transmission efficiency by STY geographic population dropped significantly, by −44.2% with ampicillin and −25.01% with rifampicin, but HDZ geographic population decreased by only 14.19% with ampicillin and 23.88% with rifampicin. Transcriptomic analysis showed that the number of DEGs related to the immune system, carbohydrate metabolism and lipid metabolism did increase in the STY rifampicin treatment, while replication and repair, glycan biosynthesis and metabolism increased in the STY ampicillin treatment. Proteomic analysis showed that the abundance of symbionin symL, nascent polypeptide−associated complex subunit alpha and proteasome differed significantly between the two geographic populations. We found that the endosymbionts can mediate vector viral transmission. They should therefore be included in investigations into aphid–virus interactions and plant disease epidemiology. Our findings should also help with the development of strategies to prevent virus transmission.


Introduction
Wheat, Triticum aestivum (L.), the third largest food crop in China, is severely attacked by aphids, including Sitobion miscanthi (Fabricus), one of the most economically important insect pests. This aphid directly pierces wheat plants and sucks the phloem sap, thus indirectly acting as the main vector for barley yellow dwarf virus (BYDV). BYDV is transmitted in a persistent and circulative pattern, causing wheat yellow dwarf virus disease. BYDV is a major phytovirus of the genus Luteovirus (family Luteoviridae), which can adversely affect almost all members of the Gramineae, causing severe crop losses worldwide [1].
Five strains of BYDV, based on their primary aphid vectors, have been identified [1]. Each strain is only transmitted efficiently by its corresponding aphid species [2]. A virus

Viral Transmission
The effect of the antibiotics (ampicillin and rifampicin) on vector transmission of BYDV−PAV isolates (CN and BE) was compared for each S. miscanthi geographic population (STY and HDZ) with the control ( Table 1). The efficiencies of BYDV−PAV transmission by S. miscanthi tested were reduced when the aphids were previously treated with antibiotics. When the STY geographic population was infected with BYDV−PAV−CN isolate, the inhibition rates ranged from 25.0% to 44.2% after antibiotic treatment, corresponding to significant difference with antibiotic−free treatments (t = 7.93 and p < 0.001). For the HDZ geographic population, the inhibition rates ranged from 14.2% to 23.9%, corresponding to a significant reduction in the virus transmission rate (t = 4.37 and p < 0.001). The highest inhibition rate of virus transmission occurred in the STY geographic population treated with ampicillin.
As to the Belgian virus isolate, BYDV−PAV−BE, the observed virus transmission inhibition rates were low after treatment with two antibiotics for the HDZ aphid geographic population (about 3.5%) and not significantly different from the control (t = 0.20; p = 0.842). For the STY geographic population, the inhibition rate was higher (from 21.4% to 25.84%) and similar to the one obtained with the Chinese virus isolate. The inhibition rates of BYDV−PAV−BE isolate transmission were significantly higher with STY than HDZ geographic populations whatever the considered antibiotic (t = 10.18 and t = 7.19 for both p < 0.001).
Whether aphids were infected with BYDV−PAV−CN or BYDV−PAV−BE isolates, the percentage of virus transmission of STY geographic population aphids that were treated with ampicillin was higher than aphids treated with rifampicin; in contrast, the percentage of virus transmission of HDZ geographic population aphids that were treated with ampicillin was lower than aphids treated with rifampicin.

Symbiotic Population Screening
As expected, the primary symbiont, Buchnera aphidicola, was detected in the two aphids' geographic population. However, the composition of the S−symbionts differed in the two geographic populations. PASS1, PASS2, PAUS, Rickettsia2 and Wolbachia were not detected in any sample. After antibiotic treatment, only some symbionts were eliminated from aphid geographic populations. Rickettsia1 was eliminated from both the STY geographic population and the HDZ geographic population (Table S1).
The relative abundance of the endosymbiont genes in the two aphids' geographic populations with different treatments (ampicillin and rifampicin) was analyzed using a comparative ∆∆Ct method ( Figure 1). The abundance of endosymbiont was significantly higher than those aphids fed with antibiotics diets, except Buchnera aphidicola from HDZ clone was treated with rifampicin and Spiroplasma from STY clone treated with ampicillin. Comparison between the two aphids' geographic population in one endosymbiont; ** significantly different (Student's t−test, p < 0.01). Comparison among the expression profiles of endosymbiont treated with different antibiotic in the same aphid geographic population; "abc" significantly different (p < 0.01); "ns" no significantly different (p < 0.01). Comparison between the two aphids' geographic population in one endosymbiont; ** significantly different (Student's t−test, p < 0.01). Comparison among the expression profiles of endosymbiont treated with different antibiotic in the same aphid geographic population; "abc" significantly different (p < 0.01); "ns" no significantly different (p < 0.01).

Transcriptome Overview
The transcriptomes of S. miscanthi feeding on free, antibiotic and BYDV−PAV were sequenced and compared. A total of 129.33 Gb of clean data was obtained from the 18 treatments, and each of these samples contained ≥ 5.4 Gb of data with Q30 quality scores ≥ 92.55% (Table S2), and 56,196 unigenes were identified with 34,941 unigenes having annotation information (Table S3).
The gene expression levels were used to conduct a PCA for each of the biological replicates. Each replicate from the same group was clustered closely together, which suggested that the repeatability of each treatment was satisfactory, and the samples from different antibiotics of S. miscanthi reared with BYDV were clustered far from each other and the control groups, which indicated that aphids feeding on antibiotics induced significant changes in gene expression (Figure 2A). The p value ≤ 0.01 (false discovery rate [FDR] adjusted) and Log2−fold change (Log2FC) ≥ 1 or ≤−1 were set as thresholds for DEGs in aphids at different treatments. Then, these identified DEGs were used for further analysis. Up− and downregulated DEGs were identified between different treatments, respectively ( Figure 2B). The distributions of up− and downregulated genes were calculated for rifampicin or ampicillin and are presented in a Venn diagram ( Figure 2C,D). GO analysis was used for the functional classification of the DEGs in aphids after rearing with antibiotics. The top 30 enriched GO terms of all DEGs are shown in Table 2. Among the STY−Free−vs.−STY−Vir, STY−Free−vs.−STY−Rif + −Vir and STY−Free−vs.−STY− Amp + −Vir, the top 10 upregulated DEGs, three genes (CRC, CRA1, adhesive plaque matrix protein−like) were annotated in three group, and one gene (CRB) was annotated in two antibiotic treatments. Among the top 10 downregulated DEGs, three genes (uncharacter-ized protein LOC100158873 precursor, RNA−binding protein 14, integumentary mucin C.1) were annotated in two antibiotic treatments.  Compared to the rifampicin-S. miscanthi and ampicillin-S. miscanthi, fed with BYDV for 48 h, rifampicin-S. miscanthi had more immune system−, lipid metabolism− and carbohydrate metabolism−related DEGs upregulated, but ampicillin-S. miscanthi had more replication and repaired−related and glycan biosynthesis and metabolism−related DEGs upregulated ( Figure 3).

Protein Identification
A proteomic work was conducted by 2D−DIGE to monitor the different protein expression from two geographic populations: STY geographic population and HDZ geographic population. More than 250 spots were generated but 86 proteins were selected for identification ( Figure 4), mainly (66.0%) with homology with proteins from Acyrthosiphon pisum (which is actually the only aphid species for which the entire genome has been sequenced) ( Figure 5), and classified into 12 functional categories based on their functions ( Figure 6 and Table 3). Compared to the rifampicin-S. miscanthi and ampicillin-S. miscanthi, fed with BYDV for 48 h, rifampicin-S. miscanthi had more immune system−, lipid metabolism− and carbohydrate metabolism−related DEGs upregulated, but ampicillin-S. miscanthi had more replication and repaired−related and glycan biosynthesis and metabolism−related DEGs upregulated ( Figure 3).

Protein Identification
A proteomic work was conducted by 2D−DIGE to monitor the different protein expression from two geographic populations: STY geographic population and HDZ geographic population. More than 250 spots were generated but 86 proteins were selected for identification ( Figure 4), mainly (66.0%) with homology with proteins from Acyrthosiphon pisum (which is actually the only aphid species for which the entire genome has been sequenced) ( Figure 5), and classified into 12 functional categories based on their functions ( Figure 6 and Table 3).
From the variation in 86 proteins, only 14 proteins were upregulated for the inefficient vector against 63 proteins upregulated for the efficient vector.     From the variation in 86 proteins, only 14 proteins were upregulated for the inefficient vector against 63 proteins upregulated for the efficient vector.     Numbers in cells correspond to the spot number on the 2D−DIGE gel. Red represents the downregulated proteins and green represents the upregulated ones of Sitobion miscanthi. The darker the color, the greater the change in protein expression (1− to 5−fold ratio for both geographic populations).

Discussion
For the biological function of an individual symbiont in such complex systems to be understood, a moderate rifampicin treatment of A. pisum and S. miscanthi has been shown to selectively eliminate Buchnera aphidicola, and ampicillin selectively eliminated Regiella and Serratia [12][13][14]. However, in this study, Buchnera aphidicola was found in all S. miscanthi geographic populations after treating with rifampicin, but its concentration was reduced. We speculate that rifampicin treatment might reduce symbiont density but not completely remove Buchnera aphidicola. When S. miscanthi was fed an ampicillin or rifampicin diet, Rickettsia was systematically eliminated in the present study; the Rickettsia symbiont, like other γ−proteobacteria symbionts identified in secondary mycetocytes and sheath cells from A. pisum, was more exposed to antibiotics and thus eliminated [15]. Many studies illustrated that PABS was localized not only in secondary mycetocytes and sheath cells, but also in the hemolymph [4,12], so its concentration was reduced by antibiotics. Arsenophonus and Spiroplasma were successfully eliminated after treatment with rifampicin, but not with ampicillin. This result is similar to a study on Bemisia tabaci where rifampicin inactivated a higher percentage of Arsenophonus than rifampicin [16].
As expected, virus transmission was reduced following the antibiotic treatment; the endosymbionts were presumably killed or inhibited, decreasing the efficiency of BYDV transmission. Since Rickettsia was the only S−symbiont in the HDZ geographic population, Rickettsia might be an important factor in the facilitation of BYDV transmission. Similarly, Kliot et al. [17] showed that a B. tabaci strain infected with Rickettsia acquired more tomato yellow leaf curl virus (TYLCY) from infected plants, retained the virus longer and exhibited nearly double the transmission efficiency than a non−infected strain that had the same genetic background. When TYLCV was acquired, it induced massive activation of gene expression in the Rickettsia uninfected population, whereas in the Rickettsia−infected population, the virus induced massive downregulation of gene expression. Fitness and choice experiments revealed that Rickettsia−infected whiteflies are always more attracted to TYLCV−infected plants [18]. When Sakurai et al. [15] investigated a Rickettsia symbiont using electron microscopy, virus−like particles were sometimes observed in association with Rickettsia cells. So, Rickettsia could play a crucial role in BYDV transmission. We applied the model that could calculate insect symbionts and insect vector contributions to pathogen transmission by insects, proposed by Patricia et al. [19], to test whether Rickettsia is involved in BYDV−PAV transmission. The fraction of the transmission efficiency provided by Rickettsia is equal to 0.14 (ampicillin) and 0.24 (rifampicin); these data indicate that Rickettsia contributes substantially to the BYDV−PAV transmission efficiency, but not as much as the insect vector contribution. In the HDZ geographic population, Buchnera aphidicola density was reduced by rifampicin, and Rickettsia was removed; rifampicin was more effective than ampicillin at reducing virus transmission, providing evidence that Rickettsia may act in concert with Buchnera aphidicola to influence the BYDV transmission of S. miscanthi.
The circulative transmission pathway through an aphid vector involves complex interactions between viral proteins and vector−associated compounds [8]. Using the proteomic and transcriptomic analysis, we identified differentially expressed proteins of the S. miscanthi STY geographic population.

Cell Signaling
The proteasome is a protein−destroying apparatus involved in many essential cellular functions. The 26S proteasome is a large, multi−subunit proteolytic machine found in the nucleus and cytoplasm of mammalian cells. It comprises a 20S cylindrical catalytic core and two 19S regulatory caps. The 20S core contains four heptameric rings, two of which contain seven alpha subunits and two that contain seven bate subunits [20]. The proteasome, protein ubiquitination machinery or both (Ubiquitin/26S proteasome (UPS) pathway) are the major types of proteolytic machinery found in eukaryotes and are associated with immune responses to pathogen invasion, linked to the activation and subcellular localization of virus replication or movement protein complexes [21]. The turnip yellow mosaic virus (TYMV) movement protein is degraded by the proteasome; UPS regulates the accumulation of TYMV during viral infection and therefore decreases viral replication [22]. UPS could protect against viral infection by regulating the proliferation and transport of viruses in host cells via targeting the degradation of many viral proteins [21]. Laodelphax striatellus 26S proteasome played a defensive role against RBSDV infection by regulating RBSDV accumulation [23]. The proteasome of R. padi is strongly implicated as an antiviral immune response against the movement process of BYDV−GPV in the body of its aphid vectors [24]. We found that most proteasomes were upregulated in highly BYDV−PAV transmission−efficient vectors; we inferred that the proteasome may enhance the BYDV−PAV transmission efficiency in S. miscanthi.

Membrane Transport
The nascent polypeptide−associated complex (NAC) is a key regulator of proteostasis to provide the cell with a regulatory feedback mechanism in which translational activity is also controlled by the folding state of the cellular proteome and the cellular response to stress [25]. The alpha subunit is one of two subunits (alpha and beta subunit) of the NAC and contributes to the prevention of inappropriate interactions. The NAC subunit alpha of Sogatella furcifera, which strongly interacted with southern rice black−streaked dwarf virus, is a major outer capsid protein [26]. The relative strengths of the interactions between the BYDV−GPV CP and NAC subunit alpha were greater than the negative control [24]. The NAC domain protein was originally characterized as the first ribosome−associated protein to contact the emerging viral polypeptide chain. Liu et al. [27] found that the NAC of small brown planthopper was confirmed in an interaction with rice stripe virus (RSV) nucleocapsid (pc3), and they proposed that NAC binding to RSV pc3 may play an important role in viral replication. The NAC domain protein can also enhance replication of tomato leaf curl virus by binding the viral replication accessory protein [28]. The NAC subunit alpha was upregulated in the STY geographic population, so the NAC subunit alpha perhaps binds with BYDV and plays an important role in viral replication.

Stress Tolerance
Another well−known protein family related to various stress responses varying between the two geographic populations was that of heat shock proteins (Hsps). In citrus tristeza virus (CTV), the protein P65 (the homologue of Hsp70) was essential for virion assembly and acted to restrict encapsidation by the minor coat protein to the 5 end of the virion [29], and P65 was found have a role in the aphid transmission of the CTV process [30]. The members of the Hsp70 family were upregulated in the STY geographic population; thus, we hypothesize that Hsp70 may be involved in the aphid transmission of BYDV.
Symbionin is abundantly synthesized by endosymbiotic bacteria Buchnera aphidicola harbored in the bacteriocyte cells and is unlikely to be exported into the aphid hemolymph [31]. Symbionin−like molecules are found in major aphid species (including BYDV vectors), except those belonging to Phylloxeridae [24]. The interaction of a coat protein-read−through protein with symbionin was considered an essential factor to stabilize virions in the hostile environment of the aphid hemolymph. Symbionin has been shown to bind to purified luteoviruses in vitro or to a recombinant luteovirus read−through polypeptide [32][33][34][35]. However, the interaction's contribution to transmission is controversial because luteoviruses bind symbionins of both vector and non−vector aphids [35], and recent studies on localization in vivo of the chaperone question its availability for interaction [36,37]. When aphids were cured of endosymbionts by treatment with antibiotics, their ability to transmit the virus was significantly reduced and the amount of coat protein was diminished. Strangely, the amount of read−through protein was not affected [32,33]. After the aphids were treated with rifampicin, the BYDV−PAV transmission efficiency was decreased by a quarter or so. The results of these experiments must be interpreted carefully-the destruction of the endosymbionts is likely to have dramatic effects on the metabolism and physiology of the aphids, and these changes may be directly or indirectly responsible for the effects on luteovirus protein detection and virus transmission. So, we propose that Buchnera aphidicola is involved in virus movement within the aphids, but we do not specify whether the effect of Buchnera aphidicola on transmitting viruses is direct or indirect.

Immune System
Insects rely on their immune system to fight against pathogens [38]. As shown in our results, whether aphids feed with or without antibiotics, after feeding on BYDV−PAV, the DEGs related to immunity in S. miscanthi were upregulated, including the MAPK signaling pathway, lysosomes, antigen processing and presentation, ubiquitin−mediated proteolysis, insect hormone biosynthesis and peroxisomes [39,40]. These results suggest that decreased bacteria Buchnera aphidicola has more of an effect on the immune system than secondary endosymbiont. The proteins involved in the cytoskeleton were also differentially expressed, which may be related to the immune response [41]. There have been previous studies showing that viruses can interact with and reorganize host cytoskeleton components for intercellular trafficking and infection processes [42]. In addition, the cytoskeleton is also commonly involved in the intracellular transport of viruses [43][44][45].
Similarly, the two geographic populations of S. miscanthi were collected from different regions, which differed in the prevalence of wheat yellow dwarf disease. STY was from northwestern China where BYDV disease is severe; HDZ was collected from the Huang−Huai region of China, where BYDV disease is less severe [46]. On the other hand, the STY geographic population has a higher diversity of symbionts than HDZ does, which suggests that the aphid's viral transmission efficiency results from increased fitness to different levels of stress posed by BYDV in the wheat−growing areas and that the symbionts may mediate the evolution of aphid fitness. Such speculation awaits further experimental evidence.

Conclusions
Whether Buchnera aphidicola density was reduced or S−symbiont was removed, BYDV transmission efficiencies of S. miscanthi were all reduced, results which suggest that endosymbiotic bacteria take part in BYDV transmission. When only S−symbiont Rickettrsia was removed, BYDV transmission was reduced significantly, suggesting Rickettsia could play a crucial role in BYDV transmission, but the function of the other S−symbionts needs deeper research. Upon further analysis, we found that the number of DEGs related to the immune system, carbohydrate metabolism and lipid metabolism were increased when Buchnera aphidicola density was reduced, but replication and repair, glycan biosynthesis and metabolism were increased when S−symbionts were eliminated. This result will contribute to further studies on exploring the immune response of S. miscanthi to viruses. As the reports on endosymbionts mediating the interaction of vector and virus transmission are scarce, our research may provide insight into the relationship between endosymbiont and luteovirus transmission. Work on virus transmission efficiencies of aphids as affected by endosymbionts should be promoted to better understand the pathway of the virus in the aphid and to develop new tools to prevent virus transmission. Indeed, identification of molecular receptors in aphids should help discover competitors that prevent binding of the virus and reduce viral transmission.

Aphids and Virus
Two S. miscanthi geographic populations were collected from winter wheat fields in Taiyuan−Shanxi Province (STY) and Dengzhou−Henan Province (HDZ). These two geographic populations were selected from a previous study [10] in which STY was the aphid geographic population that was the most efficient for the transmission of BYDV, contrary to the HDZ geographic population, which had very low efficiency. So that the risk of collecting the same genotype in multiple sampling times was reduced, individual aphids were collected from plants growing at least 10 m apart.
Two geographic populations were reared separately on potted seedlings of wheat cv. Toison d'Or (susceptible to aphids) in the second leaf stage. Each pot was isolated in a transparent, plastic, ventilated, cylindrical cage (10 × 30 cm) covered with gauze on the top. Aphids and plants were maintained in a greenhouse compartment (22 ± 1 • C, 60 ± 5% RH and 16:8 h l:d).

Antibiotic Treatment and Viral Transmission
To selectively eliminate Buchnera aphidicola or S−symbiotic, first−instar (or 24 h old) nymphs of the two geographic populations (STY and HDZ) were fed an artificial diet (15% w/v sucrose solution with and without 50 µg mL −1 rifampicin or ampicillin (Sigma, St. Louis, MO, USA)) confined between two stretched Parafilm ® membranes on an opaque cylinder for 48 h. Aphids were then transferred to the typical virus−acquisition diet (BYDV−infected wheat tissue grinded in a 15% w/v sucrose solution) for 48 h of virus acquisition. After acquiring the virus, aphids were transferred to a 7−day−old healthy wheat seedling (one aphid per test plant) protected by a plastic cage on the pot. After a 5−day inoculation access period, aphids were removed and plants were grown for 15 days in a greenhouse before testing the presence of the virus by DAS−ELISA according to the manufacturer's instructions (DSMZ, Braunschweig, Germany). The artificial diet without antibiotics ("antibiotic−free") was used as a control. Fifty wheats were formed for one biological sample; three biological replicates were performed for each treatment.
The inhibition rate of virus transmission was calculated as: ((Transmission efficiency for treated samples − Transmission efficiency for control samples)/Transmission efficiency for controls) × 100.

DNA Extraction
Aphids were soaked with 70% ethanol and sterile water several times to remove bacteria from their surface. Total DNA was extracted from 50 aphids of each S. miscanthi geographic population (STY and HDZ) following the manufacturer's instructions (DNeasy Tissue Kit, QIAGEN, Frankfurt, Germany). The quantity and purity of extracted DNA were evaluated using a spectrophotometer NanoDrop 1000 (Thermo Fisher Scientific, Pittsburgh, PA, USA). Samples were then diluted to 500 ng µL −1 .

Symbiotic Population Screening
To identify respective endosymbiotic bacteria, DNA from the samples was amplified using the specific primers of Tsuchida et al. [47] and Fukatsu et al. [48]. Amplifications were performed in a reaction volume of 20 µL including 2 µL DNA, 10 µL 2 × Taq PCR MasterMix (Invitrogen, Carlsbad, CA, USA), 1 µL forward primer (10 mM), 1 µL reverse primer (10 mM) and 6 µL ddH 2 O. The PCR cycling conditions were as follows: 95 • C for 4 min, 40 cycles at 95 • C for 30 s, 55 • C for 30 s, 72 • C for 30 s and final extension at 72 • C for 5 min. The amplified product was separated in 2% agarose gel and stained with ethidium bromide (Thermo Scientific, Waltham, MA, USA).
The relative abundance of Buchnera aphidicola and S−symbiont before/after antibiotic control was assessed using quantitative real−time PCR (qPCR). Specific primer pairs for qPCR were designed with Primer 3 (Table S4) The qPCR cycling parameters were 95 • C for 30 s, followed by 40 cycles of 95 • C for 15 s and 60 • C for 30 s. Next, the PCR products were heated to 95 • C for 15 s, cooled to 60 • C for 1 min and 95 • C for 15 s to measure the dissociation curves. qPCR reaction for each sample was carried out with three technical replicates and three biological replicates. Standard curves for reference genes and candidate genes were generated by gradient dilution to identify proper primers with 95-110% amplification efficiency and without nonspecific amplification. The relative abundance of aphid endosymbiont was normalized to the aphid housekeeping gene NADH and calculated using the comparative Ct method according to Vandesompele's method (2 −∆∆Ct ) (2002) [49].

RNA Extraction, Library Construction, and RNA Sequencing
The first−instar nymphs of STY geographic population S. miscanthi were reared on 15% w/v sucrose solution, 50 µg mL −1 rifampicin or ampicillin for 48 h, then one part of aphids transferred to feed with BYTV for 48 h. For each treatment (STY−free, STY−Vir, STY−Rif + , STY−Amp + , STY−Rif + −Vir, STY−Amp + −Vir), three experimental replicates were used. For each replicate sampling, 30 individual aphids were collected and then flash−frozen using liquid nitrogen and stored at −80 • C. Total RNA was extracted using a Trizol reagent kit (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's protocol. RNA quality was assessed on an Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA) and checked using RNase free agarose gel electrophoresis. After total RNA was extracted, eukaryotic mRNA was enriched by Oligo(dT) beads, while prokaryotic mRNA was enriched by removing rRNA by Ribo−ZeroTM Magnetic Kit (Epicentre, Madison, WI, USA). Then, the enriched mRNA was fragmented into short fragments using fragmentation buffer and reverse−transcribed into cDNA with random primers. Second−strand cDNA was synthesized by DNA polymerase I, RNase H, dNTP and buffer. Then, the cDNA fragments were purified with the QiaQuick PCR extraction kit (Qiagen, Venlo, The Netherlands), end−repaired, A base−added and ligated to Illumina sequencing adapters. The ligation products were size−selected by agarose gel electrophoresis, PCR−amplified and sequenced using Illumina novaseq 6000 by Gene Denovo Biotechnology Co. (Guangzhou, China).

RNA−Seq Data Analysis
To obtain high−quality reads, the reads containing adaptor sequences, more than 10% of unknown nucleotides (N), and low−quality (Q−value ≤ 20) bases were removed [50]. Transcriptome de novo assembly was carried out with the short reads assembling program Trinity [51]. The unigene expression was calculated and normalized to RPKM (reads per kb per million reads) [52]. Principal component analysis (PCA) was performed with R package models (http://www.r-project.org/) accessed on 10 February 2022 in this experience. RNA differential expression analysis was performed by DESeq2 [53] software between two different groups (and by edgeR (6) between two samples). The genes with a false discovery rate (FDR) below 0.05 and absolute fold change ≥ 2 were considered differentially expressed genes. Basic annotation of unigenes includes protein functional annotation, pathway annotation, COG/KOG functional annotation and Gene Ontology (GO) annotation. To annotate the unigenes, we used BLASTx program (http://www.ncbi.nlm.nih. gov/BLAST/, accessed on 10 February 2022) with an E−value threshold of 1 × 10 −5 to the NCBI non−redundant protein (Nr) database (http://www.ncbi.nlm.nih.gov, accessed on 10 February 2022), the Swiss−Prot protein database (http://www.expasy.ch/sprot, accessed on 10 February 2022), the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (http://www.genome.jp/kegg, accessed on 10 February 2022) and the COG/KOG database (http://www.ncbi.nlm.nih.gov/COG) on 10 February 2022. Protein functional annotations could then be obtained according to the best alignment results.
The protein extracts (25 µg) were labeled with cyanine dye (Cy2, Cy3, Cy5) following the standard protocol (Lumiprobe, Hannover, Germany). Before labeling, the pH of samples was adjusted to 8.5 with NaOH (100 mM). Two samples (STY or HDZ) labeled either with Cy3 or Cy5 were mixed with an internal reference standard protein mixture (which was pooled from 12.5 µg STY and 12.5 µg HDZ) labeled with Cy2. A conventional dye swap for DIGE was performed by labeling two replicates from each treatment group with one dye (Cy3 or Cy5) and the third replicate with the other two cyanine dyes. A non−labeled 500 µg sample of aphid protein mixture was added on the preparative gel for protein picking. Each mix of labeled proteins was diluted in UT−Tris buffer to obtain a volume of 225 µL. This volume was then adjusted to 450 µL with 225 µL IPG/DTT (4 µL 100× BioLyte ® 3/10 Ampholyte (Bio−Rad), 2 mg DTT (Sigma Aldrich) and 219 µL UT buffer).

2−D DIGE and Gel Analysis
The mix of labeled samples was deposited on 24 cm ReadyStrip™ IPG Strips pH 3-10 NL (Bio−Rad) for the first−dimensional isoelectric focusing (IEF) (Protean ® i12 IEF Cell, Bio−Rad) for 9 h at 50 V and 15 • C. Then, the IEF was carried out at 200 V for 2 h, 10,000 V for 1 h and 10,000 for 4 h 30 min. In an IEF unit, the current was settled at 50 µA/strip.
Before starting the second−dimensional electrophoresis, strips were reduced for 15 min in a buffer containing 30% (w/v) urea, 83% (v/v) equilibration buffer and 0.83% (w/v) dithiothreitol (DTT), and then for a further 15 min in the same buffer but in which DTT was replaced with 2% (w/v) iodoacetamide (IAA). Strips were laid down on 2D HPE TM Large Gels NF 12.5% acrylamide (Serva Electrophoresis GmbH, Heidelberg, Germany) and the second−dimensional electrophoresis was performed with the HPE FlatTop Tower (Serva) according to the manufacturer's instructions. Then, the preparative gel was placed overnight in a fixation buffer (10% acetic acid, 30% ethanol and 60% H 2 O) and stirred. The scan of gels was performed at wavelengths corresponding to each cyanine dye with a Typhoon Ettan DIGE Imager (GE Healthcare, Freiburg, Germany). Gel images were analyzed using Nonlinear Progenesis Samespots (Nonlinear Dynamics, Newcastle Upon Tyne, United Kingdom), and protein spots were excised from the gel using an Ettan spotpicker robot (GE Healthcare). Selected gel pieces were processed as described by Bauwens et al., 2013 [54].

Protein Identification
Protein identification was possible thanks to the NCBI Database (restricted to Arthropoda) and a homemade aphid symbiont database. Searches were treated on the Mascot server 2.2.06 with BioTools TM 3.2 (Bruker Daltonics). Proteins were retained only when their score was at least 45 and matched at least four peptides with error values < 100 ppm. The identified proteins were categorized according to metabolic function using the Kegg pathway database (http://www.genome.jp/kegg/pathway.html, accessed on 10 February 2022) and Expasy Proteomic tools (http://www.expasy.org/tools/, accessed on 10 February 2022), particularly the Biochemical-Metabolic pathway sections on 10 February 2022.

Statistical Analysis
For the viral transmission, an analysis of variance (ANOVA) was performed on the percentage of virus transmission of infected plants in different treatments using the GLM procedure in the SAS 9.1 program. Data were analyzed with Student's t−test. For the qPCR, differences in transcript expression of same endosymbiont among different treatments were statistically analyzed with a one−way ANOVA using SAS 9.1 followed by Duncan's Multiple Range Test. Differences in transcript expression of same endosymbiont with the same treatment between STY geographic population and HDZ geographic population were analyzed with Student's t−test.
Quantitative differences in spot intensity among the two groups were analyzed by analysis of variance implemented in SAMESPOT, version 3.5. Differential regulation of proteins was compared by a log 2 −fold change approach. A Pearson's chi−squared independence test implemented in R software (R−Core−Team, 2014) was used to test the association between groups (STY and HDZ geographic populations) and protein regulation (up− and downregulation). A heatmap was elaborated using Excel (Microsoft Corp., Redmond, Washington, DC, USA) to visualize proteins displaying increased and decreased expression.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/plants11233352/s1, Table S1: Endosymbiont detected in the Shanxi Taiyuan (STY) geographic population and Henan Dengzhou (HDZ) geographic population of Sitobion miscanthi with and without antibiotic treatment; Table S2: Summary of transcriptome data; Table S3: List of unigenes annotated by Nr, KEGG, KOG and SwissPort; Table S4: Specific primers used in this study.