Parthenogenetic adults of P. bambucicola used in this study were collected from Fuzhou, China in 2020. Paired salivary glands and guts of P. bambucicola were dissected in ice-cold phosphate-buffered saline solution (PBS, 10 mM NaH₂PO4, 1.8 mM KH₂PO4, 140 mM NaCl and 2.7 mM KCl, pH=7.4) using fine tweezers. The dissected tissues were quickly washed twice in PBS solution and immediately snap-frozen in liquid nitrogen, and then stored in a -80°C freezer. Approximately 200 pairs of salivary glands were used for transcriptome and proteome sequencing, respectively, and 200 pairs of salivary glands and 150 guts were also prepared for RT-qPCR for comparison of gene expression between samples. Each sample consisted of three biological replicates. Most previous studies used artificial diets to collect saliva and identified aphid salivary proteins successfully (Harmel et al., 2008; Carolan et al., 2009; Rao et al., 2013; Yang et al., 2018). However, it is difficult to rear P. bambucicola with traditional artificial diet as other aphids due to its unique feeding habitat. Therefore, we turned to identify injected salivary proteins by comparative proteomic analysis of fed and unfed bamboos with aphid salivary gland transcriptome data as the search database. For identification of injected salivary proteins by P. bambucicola, the Bambusa multiplex stems that had been fed continuously by over thousands of P. bambucicola individuals for three days were collected, with stems that had not been fed using for control. Several samples were unqualified and thus discard during the sequencing process, two bamboo samples fed and one sample unfed by P. bambucicola were used finally. Transcriptome and proteome sequencing were performed by Sangon Biotech (Shanghai, China) and Jingjie PTM BioLab (Hangzhou, Zhejiang, China), respectively.

For identifying all PCWDEs in P. bambucicola whole body, individuals of different morphs and developmental stages, including newborn 1st instar normal nymphs, newborn 1st instar soldiers, older 1st instar normal nymphs, older 1st instar soldiers, medium instar normal nymphs, viviparous adult females producing soldiers and viviparous adult females producing normal nymphs, were also collected and subjected for transcriptome sequencing.

The P. bambucicola samples across different morphs and different developmental stages and approximately 200 pairs of dissected salivary glands of P. bambucicola were used for RNA extraction using TRIzol Reagent (Qiagen, CA). The RNA concentration and quality were assessed by a NanoDrop spectrophotometer, gel electrophoresis and an Agilent 2100 Bioanalyzer system (Agilent Technologies, CA, USA). Qualified RNA was then used for cDNA library construction. The generated libraries were sequenced using the DNBSEQ sequencing platform. The obtained raw data was subjected to removing adapters, low quality sequences and ambiguous nucleotides (reads with more than 5% N bases). The obtained clean data was used for de novo assembly with the Trinity (Grabherr et al., 2011) to obtain final unigenes. Bowtie (Bowtie, RRID: SCR_005476) (Langmead and Salzberg, 2012) was used for aligning clean reads to the unigene library, and then RSEM (RSEM, RRID: SCR_013027) (Li and Dewey, 2011) was used to calculate the gene expression level of unigenes. The relative abundance of unigenes was measured by FPKM, which represents fragments per kilobase of transcript per million mapped reads. For functional annotation, all predicted unigenes were run blast against multiple public databases, including non-redundant protein (Nr) database, Nt, Swiss-Prot, Kyoto Encyclopedia of Genes and Genomes (KEGG), euKaryotic Ortholog Groups (KOG), Pfam and Gene Ontology (GO) databases.

Insects can adapt to plant defense responses by utilizing effectors from a variety of sources, such as salivary proteins, intestinal proteins and symbiotic microorganism derived functional proteins (Zhao et al., 2019). A total of 11 genes associated with aphid feeding, including five genes related to digestion and six genes related to defense, were randomly selected for detection of gene expression levels between salivary glands and guts, with the HSP70A1 (heat shock protein 70 A1-like) and MGST1 (microsomal glutathione S-transferase 1-like) used as reference genes to normalize selected genes’ expression ( Table S1 ). Primer Premier 5.0 (Premier Biosoft, CA, USA) was used to design RT-qPCR specific primers for selected and reference genes, as shown in Table S1 . cDNA was synthesized using FastKing gDNA Dispelling RT SuperMix (Tiangen, Beijing, China). RT-qPCR was performed with Green qPCR SuperMix Kit (TransGen Biotech, Beijing, China) following the manufacturer’s instructions. Three biological replicates were performed on each salivary gland and intestinal tract sample, and each biological replicate was run in three technical replicates. All data were analyzed by Graphpad (GraphPad, RRID: SCR_000306) (http://graphpad.com/) version 9.0 software with unpaired t-test (P < 0.05).

The label-free LC-MS/MS quantitative proteomic analysis was performed by the Jingjie PTM BioLab. The salivary gland samples were grinded with liquid nitrogen into cell powder and transferred to 5 ml centrifuge tube. After adding four times the volume of lysis buffer (including 1% SDS and 1% protease inhibitor cocktail), the cell powder samples were boiled with a metal bath at 95°C for 10min, and were sonicated with a high intensity ultrasonic processor (Scientz, Ningbo, China). To remove cell debris, the protein solution was spun for 10 min (12000 g at 4°C) and the supernatant was pipetted into clean tubes. The protein concentration was determined using BCA Protein Assay kit (Beyotime, Shanghai, China) following the manufacturer’s instructions.

For digestion, an equal amount of protein for each sample was used and lysis buffer was added to adjust to the same volume. After adding dithiothreitol (DTT) to a final concentration of 5 mM, the protein solution was incubated at 56°C for 30 min, followed by adding iodoacetamide (IAA) to 11 mM final concentration and incubating 15 min at room temperature in the dark to alkylate cysteines. The alkylated protein samples were transferred to ultrafiltration tubes, centrifuged at 12000 g for 20 min at room temperature. The protein was re-suspended in 8 M urea (Sigma) for 3 times, and then urea was also re-suspended with 100mM ammonium bicarbonate solution for 3 times. Trypsin was added for a final trypsin:protein ratio of 1:50 (w/w) and incubated overnight. The peptides were recovered by centrifugation at 12000 g for 10 min at room temperature, and then recovered again with ultrapure water. The two peptide solutions were then combined.

The peptides were dissolved with solvent A (0.1% formic acid and 2% acetonitrile in water) and then separated using the Easy-nLC 1200 ultra-high-performance liquid system. The separated peptides were ionized by injection into an NSI ion source and then analyzed by Orbitrap Exploris™ 480 mass spectrometers (Exploris 480, Thermo Fisher Scientific, USA). The electrospray ionization voltage was set to 2.3 kV, and a high-resolution Orbitrap was used to detect and analyze the peptide parent ions and their secondary fragments. The primary mass spectrum range was 400-1200 m/z with the scanning resolution was set to 60000. The fixed start point of the secondary mass spectrum scan range was 110 m/z with the scanning resolution of 15000, and TurboTMT was set to Off. A data dependent scanning (DDA) program based on Cycle time was used as the data acquisition mode. Specifically, within a 1.0-s cycle period, the parent ions of the peptide were selected according to the sequence of the signal intensity from high to low, and then entered the HCD collision pool to fragment with 27% of the fragmentation energy. The secondary mass spectrometry analysis was also performed sequentially. To improve the efficient utilization of MS, the automatic gain control (AGC) was set to 100%, the signal threshold was set to 5E4 ions/s, the maximum injection time was set as Auto, and the dynamic exclusion time of tandem MS scanning was set to 20 s to avoid repeated parent ion scanning.

The resulted MS/MS data were analyzed using Maxquant search engine (version v1.6.15.0) (Prianichnikov et al., 2020) with the protein sets (22,597 sequences) of salivary gland transcriptome using as the retrieval database and an inverse decoy library used to calculate the false positive rate (FDR). The cleavage enzyme was Trypsin/P allowing up to 2 cleavages; the minimum peptide length was 7 amino acid residues; the maximum number of modifications of the peptide was set as 5; the mass tolerance of precursor ions was 20 ppm in the first search and 4.5 ppm in the main search, respectively. The mass tolerance of fragment ions is 20 ppm. Carbamidomethyl on cys was set as a fixed modification with oxidation on Met, acetylation on protein N-terminal and decarboxamidation as variable modifications. The false discovery rate (FDR) for both protein identification and peptide−spectrum matches (PSMs) identification was 1%.

For protein extraction of bamboo tissues, samples of bamboo fed and unfed by P. bambucicola were grinded with liquid nitrogen. The powder samples were sonicated with a high intensity ultrasonic processor (Scientz, Ningbo, China) after adding four times the volume of lysis buffer (including 10 mM dithiothreitol and 1% protease inhibitor cocktail). An equal volume of Tris-saturated phenol was then added and centrifuged for 10 min (5500 g at 4°C). The supernatant was collected in clean centrifuge tubes and five times the volume of 0.1 M ammonium acetate/methanol were added and incubated at -20°C overnight. After centrifugation at 4°C for 10 min, the supernatant was removed, and the precipitate was washed with cold methanol once and cold acetone for three times, respectively. The precipitate was redissolved with 8 M urea (Sigma), and the protein concentration was determined using BCA Protein Assay kit (Beyotime, Shanghai, China) following the manufacturer’s instructions. For enzymatic digestion, an equal amount of protein for each sample was taken and adjusted to the same volume with lysis buffer. TCA was added slowly to a final concentration of 20% TCA, mixed by vortex, and precipitated for 2h at 4°C. After centrifugation at 4500g for 5min, the supernatant was discarded and the precipitate was washed with precooled acetone for two to three times. After the precipitation was dried, TEAB was added to a final concentration of 200 mM, and the precipitation was broken up by ultrasound. Trypsin was then added at 1:50 (trypsin: protein, m/m) ratio for digestion overnight. Dithiothreitol (DTT) was added to make the final concentration of 5 mM and reduced at 56°C for 30 min. Then iodoacetamide (IAA) was added to 11 mM final concentration and incubated for 15 min at room temperature under in the dark. The subsequent LC-MS/MS analysis were then performed as described above. And the obtained MS/MS spectra data were searched separately against the salivary gland transcriptomic database and Bambusa protein database (including 14361 proteins) download from Nr protein database of the NCBI (accessed on April 19, 2022).

The identified salivary gland proteins, saliva proteins and bamboo proteins was annotated with multiple public databases, including Nr, KEGG, Swiss-Prot, Pfam, GO and KOG databases. Functional enrichment analysis was then conducted with the functions phyper in R software, with FDR adjust P-value (Qvalue) < 0.05 as the threshold.

Aphid effectors are likely secreted proteins delivered into the saliva secreted by salivary glands to mediate plant defenses (Bos et al., 2010). For identification of candidate effectors, signal peptides were predicted from the amino acid sequences of dual transcriptomic-proteomic data from salivary glands as well as proteomic data from saliva using SignalP (SignalP, RRID: SCR_015644) (https://dtu.biolib.com/SignalP-6) v6.0 (Teufel et al., 2022), followed by DeepTMHMM (Hallgren et al., 2022) to identify transmembrane domains for proteins containing signal peptides. Proteins containing an N signal peptide but no transmembrane domain were regarded as candidate effectors.

The degradation of plant cell wall components requires a large repertoire of highly specialized carbohydrate-active enzymes (CAZymes) that are produced by the organism itself or its associated symbiotic microbes (Ni and Tokuda, 2013; Scully et al., 2013; Bredon et al., 2019). Firstly, Hmmscan program in the HMMER (Hmmer, RRID : SCR_005305) (http://hmmer.janelia.org/) version 3.1b2 (Eddy, 1998) was used to search amino acid sequences of transcriptome and proteome against the family specific HMM profiles of CAZymes within dbCAN HMMdb v11 to identify CAZymes and assign them to CAZy families, with an e-value cutoff 1e-3 (≤80 aa) or 1e-5 (>80 aa) and coverage above 30% as the filter threshold. CAZy families can be classified into glycoside hydrolases (GHs), glycosyltransferases (GTs), polysaccharide lyases (PLs), carbohydrate esterases (CEs), auxiliary activities (AAs) and carbohydrate-binding modules (CBMs) (Drula et al., 2022). Enzymatic activity of all identified CAZymes were detected using Hotpep (Busk et al., 2017) to determine whether they are candidate plant cell wall degrading enzymes (PCWDEs). The identified candidates were further confirmed as putative PCWDEs by reference to Tokuda (2019). The transcriptome data used for PCWDEs identification in this study include salivary gland transcriptome and whole-body transcriptome of P. bambucicola across different morphs and developmental stages.

An average of 47,326,665 bp raw reads was yielded from transcriptome of P. bambucicola salivary glands. After data filtering, a total of 19.47 Gb clean data was used for de novo assembly, resulting in 48,028 unigenes with an average length of 1,310 bp and N50 of 2,400 bp. There were 28,512 (59.37%), 24,816 (51.67%), 21,093 (43.92%), 22,891 (47.66%), 20,427 (42.53%), 20,570 (42.83%) and 13,385 (27.87%) unigenes homologous to known sequences in the Nr, Nt, SwissProt, KEGG, KOG, PFAM and GO databases, respectively. About 32,111 (66.86%) unigenes were functionally annotated in at least one of the used databases, and many of them could be annotated by multiple databases ( Figure 1A ). The unigenes showed the most similarity with Sipha flava according to the matched species distribution of annotation based on Nr database ( Figure 1B ).

GO function classification was conducted on predicted unigenes and showed that most of them were enriched in cellular anatomical entity, cellular process, binding, catalytic activity and metabolic process ( Figure 1C ). The top 20 highly-expressed unigenes include some genes associated with mitochondrial activity, genes encoding a invertebrate-type lysozyme 6, a odorant-binding protein 2, a prohormone-2, a alpha-N-acetylgalactosaminidase, a putative sheath protein and some genes encoding proteins with unknown functions ( Table S2 ). These highly expressed genes in salivary glands such as the gene encoding putative sheath protein and genes with unknown functions were worth for further study. In addition, many homologous genes encoding salivary proteins that are known to paly an important role in aphid-plant interactions, such as some digestive enzymes, detoxifying enzymes, antioxidant enzymes and some effector proteins modulating plant defenses, were also identified in P. bambucicola salivary glands ( Table S3 ). These results suggest that some salivary components are conserved across different aphids.

In addition to the salivary glands, guts are also important for insect feeding. To understand the relative role of salivary gland and gut in P. bambucicola feeding, expression levels of 11 genes, including five genes related to digestion and six genes involved in detoxification and antioxidant activities, were detected between salivary glands and guts of P. bambucicola by RT-qPCR. All 11 genes expressed in both salivary glands and guts. Among the five digestive-related genes, significant differences were found in the expression of beta-galactosidase-like (GLB1), lysosomal alpha-mannosidase (MAN2B1), carboxypeptidase E-like (CPE) and methionine aminopeptidase 1D, mitochondrial (METAP1D) between salivary glands and guts except the AAEL006169 (lysosomal aspartic protease), among which the GLB1, MAN2B1 and CPE showed evidently higher expression levels in salivary glands ( Figures 2A–E ). For six genes involved in detoxification and antioxidant activities, the expression level of glucose dehydrogenase [FAD, quinone] (CHDH) in salivary glands was remarkably higher than that in guts; and superoxide dismutase [Cu-Zn] 1 (SOD1) also showed higher expression levels in salivary glands. Phospholipid hydroperoxide glutathione peroxidase (GPX4) and glutathione S-transferase-like (GST) were highly expressed in guts compared with salivary glands ( Figures 2F–K ). These results suggest that both salivary glands and guts may play important roles in digestion and detoxification in P. bambucicola plant feeding.

A total of 4793 proteins were detected from the salivary gland proteome. Of them, 3115 proteins attributed to at least one GO term, with the cellular metabolic process and organic substance metabolic process, organelle and cytoplasm, protein binding and hydrolase activity being the two most represented terms in each of the three categories, respectively ( Figure S1A ). Proteins without unique peptides in two of three replicates and those with an average of unique peptides less than two were filtered out, resulting in 2442 candidate proteins. All proteins were annotated with the KEGG pathway database to characterize the general metabolic functions of the salivary gland proteome, and many proteins were classified and associated with global and overview maps, signal transduction, endocrine system, translation, and transport and catabolism pathways ( Figure 3A ), consistent with the biological roles of salivary glands.

The raw MS/MS data of bamboo tissues unfed and fed by aphids were then analyzed and searched against the transcriptomic data of P. bambucicola salivary glands to identify putative salivary proteins in secreted saliva. In total, 3244 proteins were detected, and 1496 of them were expressed only in bamboo tissues after aphid feeding, and were regarded as aphid salivary candidates secreted into host plants. To reduce false positives, 1160 proteins with at least two unique peptides were selected for further metabolic functional analysis. The majority of proteins participated in pathways of global and overview maps, signal transduction, translation, transport and catabolism and immune system, which was generally similar to that of salivary gland proteome ( Figure 3B ).

KEGG pathway enrichment analysis showed that the salivary gland proteins were significantly enriched in pathways of proteasome, protein processing in endoplasmic reticulum, carbon metabolism, citrate cycle (TCA cycle), oxidative phosphorylation and protein export, while ribosome and oxidative phosphorylation were the most representative pathways for putative saliva proteins ( Figure S1B, C ). These results may indicate the important roles of salivary glands in protein secretion and energy metabolism.

From transcriptome data of salivary glands, 1213 putative secretory proteins were predicted. Annotation against NCBI Nr database showed that 812 (66.94%) of all identified secretory proteins were functionally annotated, 359 (29.60%) proteins were annotated with unknown functions, and 42 (3.46%) proteins showed no similarity with all known sequences in the Nr database ( Table S4 ). Functional enrichment analysis of these putative secretory proteins showed that the most enriched GO terms were structural constituent of cuticle, extracellular region and carbohydrate metabolic process, and the most enriched KEGG pathways were RNA polymerase, lysosome and other glycan degradation ( Figure S2 ).

When we identified secreted proteins from mass spectrometry proteins, a total of 196 and 114 putative secreted effector candidates were predicted from the salivary gland and saliva proteome of P. bambucicola, respectively. Many of them were hypothetical proteins with unknown functions or functionally annotated proteins whose roles in aphid-plant interactions are not clear. Some candidate secretory effectors were homologous proteins also characterized in secretome of other aphid species, which had been supposed to play important roles in aphid-host interactions. For example, some detoxifying and antioxidant enzymes including glucose dehydrogenase, glutathione S-transferase, several carboxylesterases and peroxidases, were identified from salivary gland or saliva of P. bambucicola ( Table S5 ). The putative salivary secretory effectors also contained some digestive enzymes such as sugar-degrading enzymes, carboxypeptidase, cathepsins, serine proteases and phospholipase, and some effectors modulating plant immunity and defense such as apolipophorin, odorant binding protein and yellow-like protein. In addition, some salivary glue proteins and cuticle proteins were also identified in both salivary gland and saliva of P. bambucicola, while two sheath proteins were detected only in salivary glands ( Table S5 ). There were 44 putative secretory effectors in both salivary gland and saliva, including three salivary glue protein, a sheath protein (mucin-5AC protein), a venom serine carboxypeptidase, cathepsin L, phospholipase A-2-activating protein, odorant binding protein in addition to the above mentioned common salivary gland protein and several cuticle proteins. These putative secretory effectors may help promote aphid stylet penetration, digestion and detoxification activities, or contribute to suppression or activation of plant defense responses.

The PCWDEs from the transcriptome of whole body and salivary glands were identified ( Table 1 , Table S6 ). A total of eight potential PCWDEs were identified based on salivary gland transcriptome, including four β-glucosidases, one endo-β-1,4-glucanase and three β-mannosidases. However, most of genes encoding PCWDEs showed very low expression levels ( Figure 4A , Table 1 ). We did not identify any transcripts with potential pectinase activity, indicating that the P. bambucicola may loss ability to secrete and degrade pectin. When we detected PCWDEs expressed in translational levels using the proteome data of salivary glands and secreted saliva, only one β-glucosidase and one β-mannosidase were detected at protein level in salivary gland, and no other PCWDEs were detected at protein level in the secreted saliva except for four β-glucosidases with putative cellulolytic activity. Of these, only one β-mannosidase identified from salivary gland proteome and one β-glucosidase from saliva proteome were predicted containing secretory signal. Besides, two of the four β-glucosidase transcripts identified in the proteome of secreted saliva were not full-length and so it is uncertain whether they contain secretory signal or not.

To obtain more complete PCWDEs in P. bambucicola, PCWDEs from whole-body transcriptome across different morphs and developmental stages of P. bambucicola were also identified. A total of three β-glucosidases, one endo-β-1,4-glucanase and one β-mannosidase were identified, among which β-glucosidase 1, the endo-β-1,4-glucanase 1 and beta-mannosidase 1 were found in both the whole-body and salivary gland ( Table 1 , Table S6 ). Except for two glucosidases (β-glucosidase 2 and β-glucosidase 3) found only in whole-body, the three putative PCWDEs identified in both the whole-body and salivary gland all contained secretory signals. All predicted PCWDEs from body transcriptome of P. bambucicola exhibited high expression levels (with an average of FPKM > 16) ( Figure 4B ).

An comparative proteomic analysis of bamboo tissues unfed and fed by aphids may reveal the changes in protein expression and plant cellular process modulated by P. bambucicola. A total of 171 proteins were differentially expressed between two types of bamboo tissues, with 71 of them up-regulating and 100 of them down-regulating in bamboo after being fed respectively ( Table S7 ). The highly expressed proteins were mainly enriched in categories of non-membrane-bounded organelle, nucleolus, cell surface, hydrolase activity and maintenance of protein location in cell, while the downregulated proteins were mainly associated with vesicle, golgi apparatus, cellulose synthase activity and plant cell wall biogenesis ( Figure S3 ).