A Comprehensive Transcriptomic and Proteomics Analysis of Candidate Secretory Proteins in Rose Grain Aphid, Metopolophium dirhodum (Walker)

1. Introduction

Aphids (Hemiptera: Aphididae) are voraciously phytophagous insects that ingest plant phloem sap by their needle-like stylets (mouth parts). They are unusual herbivores because their feeding site is a single phloem cell in the sieve element buried deep within plant tissues; yet, they represent one of the most important insect pests in temperate agriculture Gupta [2,3]. Globally, there are approximately 5558 aphid species belonging to 703 genera and 30 subfamilies [4]. Among the recorded aphid species, around 450 of them feed on crop plants, and 250 aphid species are considered economically important agricultural pests causing significant economic losses [5,6]. Four aphid species cause serious damage to wheat farms in China. The major aphid species infesting wheat in China are the grain aphid (Sitobion avenae Fabricius), green bug (Schizaphis graminum Rondani), bird cherry oat aphid (Rhopalosiphum padi Linnaeus) and rose grain aphid (Metopolophium dirhodum Walker) [7]. Some Phloem feeding insects, such as pea aphids (Acyrthosiphon pisum Harris) have a wide host range often colonizing many leguminous plant species. Unlike the other types of aphids, M. dirhodum are primarily restricted to crops that belong to the Rosaceae (mainly cherries and peaches), but are a secondary pest of cereal crops (wheat and other cereal crops).

Certain plants have the ability to withstand aphid feeding without negative consequences, either through inherent genetic resistance or by adjusting their interactions with aphids. This natural defense mechanism has proven effective in controlling aphids in various crops. The interactions between plants and aphids can be seen as a process of coevolution, where both parties adapt to each other [8]. Aphids extract nutrients from host plants by piercing the phloem to access sap with their stylets, releasing saliva from salivary glands that contain effector proteins to counter plant defenses [9,10]. The composition of watery saliva, containing a diverse blend of enzymes and other substances, varies significantly among aphid species and even within the same species depending on their diet [11,12,13]. The ability of aphids to adapt to different host plants is closely tied to the diversity in watery saliva content. This generally involves the secretion of secretory molecules specifically effectors to target host molecules [14]. Many studies have shown that insects, including aphids, produce and secrete effectors that suppress or induce plant defense responses [15,16,17,18]. These aphid secretory proteins are thought to be produced predominantly inside the head and are secreted along with saliva during probing and feeding [19,20,21]. Nonetheless, limited research has delved into the connection between feeding patterns and the genes expressed in aphid salivary glands during shifts in host preferences.

The recent availability of data on the genome and transcriptome sequences of aphids has enabled the development of an approach to identify potential secretory protein of aphids [9,15,19]. Bioinformatics pipelines to identify putative secreted proteins have been developed and applied to several aphid species [17,22]. In addition, saliva collection methods based on artificial feeding systems combined with mass spectrometry allow for the detection of many proteins in the saliva of several types of aphids released during probing and feeding [15,23]. Significant progress has been made in the identification and functional annotation of potential effectors of economically important aphid species of cereals, such as bird cherry-oat aphids, and green bug aphids, through proteomic and transcriptomic analysis. However, limited gene and/or proteomic data for rose-grain aphids have been documented, despite the fact that rose-grain aphids are the fourth most dominant and destructive aphids affecting wheat production in many temperate regions, including China. Furthermore, the molecular mechanisms underlying these differences have not been documented because there is a lack of sufficient genetic data for rose grain aphids in public databases. Consequently, this study employed a combined transcriptomic and proteomic approach to mark an initial step in identifying candidate secretory proteins involved in plant infestation and environmental adaptation.

2. Materials and Methods

2.1. Aphid Handling and RNA Preparation

Both apterous and alate adult rose grain aphids from a clonal lineage were raised on susceptible wheat (Triticum aestivum L. Var. Zhongmai-175) in a controlled laboratory setting. Aphid heads along with salivary glands were dissected under a stereomicroscope and temporarily stored in liquid nitrogen. RNA was extracted using TRIzol reagent and then purified using a Tianmobio total RNA extraction kit. The degradation and contamination of RNA were examined using a 1% agarose gel. RNA purity was measured with a Nano-drop spectrophotometer (IMPLEN, CA, USA). The quantity and integrity of RNA were assessed using the RNA Nano 6000 Assay Kit (Agilent Technologies, CA, USA).

2.2. cDNA Library Preparation, Cluster Generation and Sequencing

The process of constructing the cDNA library commenced with a 1.5 µg RNA sample. The synthesis of the cDNA library involved fragmenting mRNA using divalent cations in the presence of random hexamer primers and M-MuLV Reverse transcriptase. DNA Polymerase I and RNase H were employed to construct the first and second cDNA strands. The overhangs were then converted into blunt ends through exonuclease/polymerase activity. Three microliters of enzyme (NEB, USA) were used with size-selected, adaptor-ligated cDNA at 37 °C for 15 minutes, followed by 5 minutes at 95 °C. PCR amplification was carried out using Phusion High-Fidelity DNA polymerase, universal PCR primers, and Index Primer. The PCR products were purified using the AMPure XP system (Beckman Coulter, Beverly, USA), the library quality and insert size was assessed on the Agilent Bioanalyzer 2100 system. Furthermore, the library's quantity was measured using a Qubit 2.0 Fluorometer. A paired-end library of 250-300bp size was prepared according to Illumina's protocol/instructions. Subsequently, the samples were clustered with the aid of cBot Cluster Generation System along with the TruSeq PE Cluster Kit v3-cBot-HS from Illumina (Illumina, China).

2.3. Quality Control, Transcriptome Assembly, and Gene Annotation

The first stage of the analysis involved processing of the raw data (raw reads) in FASTQ format through custom Perl scripts. Subsequently, clean reads were acquired by removing intact adapter reads (the proportion of reads with connectors or linkers), ambiguous reads (The proportion of reads with unspecified base information), and low-quality reads (Q Phred ≤ 20 is the number of bases in the total read length of 50% or more of the reads) from the raw data. At the same time, the clean data was evaluated for Q20 (The percentage of bases with a Q_Phred value greater than 20 to the total bases, where Q_Phred=-10log10(e)), Q30 (The percentage of bases with a Q_Phred value greater than 30 to the total bases), GC content, and sequence duplication level. The transcriptome was assembled using Trinity software with min_kmer_cov set to 2, and all other parameters with default settings [24]. Following the assembly, reads were aligned back to the transcriptome assembly using the align_and_estimate_abundance.pl script, with RNA-seq by Expectation Maximization (RSEM) method and Bowtie to estimate abundance and mapping the reads. Transcripts with fewer than 0.5 transcripts per million mapped reads (TPM) or representing less than 10% of the expression value of the dominant isoform for each unigene were excluded from the transcriptome assembly. Protein coding regions were identified using Trans-decoder v.3.0.1, retaining the highest scoring open reading frame (ORF) for each transcript with the single_best_orf option. Finally, transcripts without open reading frames were removed from the assembly.

Functional annotation of the derived genes was done using the Basic Local Alignment Search Tool (BLAST2go) against public databases, including the National Center for Biotechnology Information (NCBI), Non-redundant Nucleotide Acid Database, Non-redundant Protein Database, and BLASTx against the Swiss-Prot Database. The best matched alignment results were considered to annotate all unigenes; if the alignment results among these databases differed, preference was given to the results from the NR database. All genes that matched the Swiss-Prot database entries were further classified according to Gene Ontology terms (matched genes with an E-value < 10^-6). Furthermore, all genes that matched the Swiss-Prot database were classified based on the Eukaryotic Orthologous Groups (KOG)/Clusters of Orthologous Groups of proteins (COG) (matched genes with an e-value < 10⁻³), and the pathways were annotated using the Kyoto Encyclopedia of Genes and Genomes (KEGG). Automatic Annotation Server (KAAS) with a cut-off E-value < 1.0x10⁻¹. Gene homology was established by searching public databases, such as non-redundant proteins (Nr, with genes having an e-value < 10⁻⁵), nucleotide sequences (Nt, with genes having an e-value < 10⁻⁵), Pfam (with genes having an e-value < 0.01), and Swiss-Prot (with genes having an e-value < 1.0x10−5), for all genes. Gene Ontology enrichment analysis of the differentially expressed genes (DEGs) was conducted using the GOseq R package. The KOBAS software [25] was utilized to assess the statistical enrichment of the differentially expressed genes in the KEGG pathways.

2.4. Prediction of Secretory Proteins

To discover potential effectors, the coding sequences of every unigene were translated into peptide sequences and subsequently analyzed using a web-based tool for predicting signal peptides (https://services.healthtech.dtu.dk/services/SignalP-6.0/). Out of 705 peptide sequences, clusters containing a signal peptide at the N-terminal of the polypeptide chain were identified. These sequences with signal peptides were then assessed using a transmembrane domain prediction web-based tool (https://dtu.biolib.com/DeepTMHMM) to identify any transmembrane domains within the protein sequence. Ultimately, protein sequences featuring a signal peptide at the N-terminal and those with one or devoid of transmembrane domains were classified as candidate secretory proteins.

2.5. Saliva Collection, Extraction and Protein Identification

2.5.1. Saliva Collection and Extraction

The rose grain aphids were dislodged from wheat plants by shaking and then placed on sterile diets (consisting of 15 % sucrose, 100 mM L-serine, 100 mM L-methionine, and 100 mM L-aspartic acid, pH 7.2 (KOH)) [21]. The diet was prepared under aseptic conditions and filtered using a 0.22μm syringe filter (Millipore, MA, USA). Approximately 250 aphids were confined in each glass tube, with 1 mL of the diet provided between two layers of Parafilm (Neenah, WI, USA). After 72 hours of feeding, the diets were collected from the space between the two layers of Parafilm using a pipette. Ultrafiltration was performed with a 3-kDa molecular-weight cutoff Amicon Ultra-4 Centrifugal Filter Device (Millipore) at 5000g at 4 °C for 30 minutes. The concentrated samples were then precipitated using a trichloroacetic acid protein precipitation kit (Sangon, Shanghai, China). The pellets were solubilized in 200µl of SDT buffer and Bicinchoninic Acid (BCA) protein determination assay (Bio-Rad, USA). The proteins were then separated on a gel using SDS-PAGE and stained with Coomassie Blue R-250 dye to visualize protein bands. Gel slices were processed further by destaining, reduction, alkylation, digestion with trypsin, and peptide extraction. The resulting peptides were dehydrated and reconstituted in 0.1% formic acid for downstream analysis.

Figure 1. Aphid feeding on artificial diet and saliva collection.

Figure 1. Aphid feeding on artificial diet and saliva collection.

2.5.2. Liquid Chromatography Tandem Mass Spectrometry (LC-MS/MS)

LC-MS/MS analysis was performed on a Q Exactive mass spectrometer (Thermo Scientific) coupled with Easy nano-liquid chromatography (Thermo Fisher Scientific). While the sample eluted from the LC column, thousands of mass spectra were acquired, and subsequently, the mass spectrum of all peptides was measured at that time point. The mass analyzer separated the peptides based on their mass-to-charge ratio, and the detector detected precursor peptide ions. The most abundant precursor peptide ion (300-1800 m/z) was further fragmented by collision with neutral gas molecules after passing it through the filter unit of the mass spectrometer.

2.5.3. Proteomic Data Analysis and Similarity Search

The acquired mass spectrometry (MS) data were analyzed using the MASCOT engine (Matrix Science, London, UK; version 2.4). MS data were searched against the UniProtKB Escherichia coli database. The search followed an enzymatic cleavage rule for trypsin, and allowed a maximum of two missed cleavage sites and a peptide mass tolerance of 20ppm for fragment ions. The cutoff global false discovery rate (FDR) for peptide and protein identification was set to 0.01 or (FDR≤0.01).

BLAST score ratio (BSR) tests were employed to compare the proteome similarities among various aphid species, including M. dirhodum, Aphis craccivora, Aphis glycines, Acyrthosiphon pisum, Aphis gossypii, Schizaphis graminum, Sitobion avenae, Rhopalosiphum padi, Myzus persicae, Melanaphis sacchari, and the yellow sugar cane aphid, Sipha flava. The reference comparisons were conducted using the BLAST software available at http://blast.ncbi.nlm.nih.gov/Blast.cgi. The BSR index was determined by dividing the BLAST query score by the reference score and normalized the result to a scale from 0 to 1. The reference library included A. craccivora, A. glycines, A. pisum, Aphis gossypii, Schizaphis graminum, S. avenae, R. padi, M. persicae, M. sacchari, and S. flava.

$$\mathrm{B}\mathrm{L}\mathrm{A}\mathrm{S}\mathrm{T} \mathrm{s}\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{e} \mathrm{r}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o} 1 \left(\mathrm{B}\mathrm{S}\mathrm{R}1\right)={\displaystyle \frac{\mathrm{Q}\mathrm{u}\mathrm{e}\mathrm{r}\mathrm{y} \mathrm{s}\mathrm{e}\mathrm{q}\mathrm{u}\mathrm{e}\mathrm{n}\mathrm{c}\mathrm{e}}{\mathrm{R}\mathrm{e}\mathrm{f}\mathrm{e}\mathrm{r}\mathrm{e}\mathrm{n}\mathrm{c}\mathrm{e} \mathrm{s}\mathrm{e}\mathrm{q}\mathrm{u}\mathrm{e}\mathrm{n}\mathrm{c}\mathrm{e} 1}}$$

3. Results

Altogether, 47,565,328 bp of raw reads were retrieved from the sequencing of the established library. After removing low-quality reads, ambiguous reads, and those related to adapters, 46,238,772 clean reads were produced using Illumina paired-end RNA-seq technology, achieving a Q20 value of 97.61% in the transcriptomic assembly. Approximately 39,293,370 reads (84.98% of the total clean reads) were successfully aligned to the reference sequence (longest transcript unigene). After filtering out adapter-related reads (609272), ambiguous reads (18194), and low-quality reads (35812), clean reads were retained. The transcriptome of the rose grain aphid head was then de novo assembled using the short-read assembling program (Trinity), resulted in clustering of 31344 transcripts and 18030 genes (Table 1). The transcripts varied in length from 301 to over 2000 base pairs, with an average size of 1532 base pairs. Out of the transcripts, 7830 (24.98%) were in the range of 301 to 500 base pairs, 6510 (20.77%) were between 501 to 1000 base pairs, 8229 (26.25%) ranged from 1001 to 2000 base pairs, and 8775 (28%) were longer than 2000 base pairs. Furthermore, among the genes that were assembled, 4949 (27.45%) had ranged between 301 to 500 base pairs, and 4261 (23.63%) fell within the range of 501 to 1000 base pairs. Among the genes analyzed, 4377 (24.27%) fell within the range of 1001 to 2000 base pairs, while 4443 (24.64%) surpassed 2000 base pair mark (Figure 2). The mean size of the genes derived from the transcripts was 1413 base pairs. In comparison, the transcripts produced by Trinity had an average length of 1532 base pairs (N50=2335bp), demonstrating longer transcripts than the average gene length of 1413 base pairs (N50=2205 bp). All unigene assembled in this study are illustrated in supplementary table 1. The raw sequencing data have been submitted to the Sequence Read Archive (SRA) of the National Center for Biotechnology Information (NCBI) database with the accession number PRJNA1134911 (http://www.ncbi.nlm.nih.gov/bioproject/ PRJNA1134911).

3.1. Functional Annotation of Gene Transcripts

Head tissue libraries were annotated utilizing BLASTX against the NCBI database, and all unigenes underwent homology searches. The functional annotation of all unigenes of M. dirhodum were subjected to BLAST search against seven databases including NCBI's non-redundant protein (Nr) and nucleotide (Nt) sequences, protein family (Pfam), the manually annotated and reviewed Swiss-Prot protein sequence database, Eukaryotic Clusters of Orthologous Groups of proteins (KOG/COG), the Kyoto Encyclopedia of Genes and Genomes (KEGG) Ortholog database, and Gene Ontology (GO). This comprehensive approach was employed to annotate the functions of each gene transcript. Of the total assembled unigenes (18030), 12589 (69.82%), 14253 (79.05%), 5587 (30.98%), 9467 (52.5%), 9314 (51.65%), 9314 (51.65%), and 5850 (32.44%) genes were matched to known proteins in the non-redundant protein sequences, non-redundant nucleotide sequences (Nt), KEGG Orthology (KO), Swiss-Prot, Protein family (Pfam), Gene Ontology (GO), and Eukaryotic Orthologous Groups (KOG) databases, respectively (Table 2). The alignment result of the total unigenes were not similar across all databases and hence, the NR database results were preferentially employed to annotate all unigenes.

The newly de novo assembled transcriptomic sequences, were examined using a BLAST search against the NCBI database to deduce the functional and structural characteristics of the proteins encoded by these transcripts. The comparison involved matching all coding sequences of M. dirhodum genes with established transcriptomic sequences from various aphid species found in recognized protein databases such as Nr, Swiss-Prot, and Pfam. Consequently, the de novo assembly of M. dirodum transcripts resulted in the prediction of 12589, 9467, and 9314 functionally known and unknown gene products (proteins) in the Nr, Swiss-Prot, and Pfam genome databases, respectively. The E-value distribution for the top matches against the nr database indicated that 17.18% of the sequences exhibit strong homology (E-value < 1.0E^-60), with most E-values falling between 0 and 1.0E^-100 (Figure 3A). Conversely, the similarity distribution revealed that 97.18% of the unique sequences with top matches have a similarity greater than 60%, whereas only 2.82% of the hits showed a similarity below 60% (Figure 3B). On top of that, the sequences of assembled unigenes were aligned to gene sequences of other organisms accessible in prominent databases to look for conserved genes (homologous sequence search). Most of the sequences of unigenes (68.2%) were aligned to A. pisum sequences, followed with M. persicae (18%), for which the genome sequences were made available in the database (Figure 3C). However, some gene transcripts (13.7%) revealed less semblance to other economically important aphid species.

3.2. Gene Ontology and Eukaryotic Orthologous Groups classification (KOG)

The gene transcripts characterized based on the Nr and Swiss-Prot databases were further annotated against Gene ontology to determine the number of genes associated with each functional category. Of the 18030 assembled gene clusters, 9314 were assigned to three functional groups. The three most abundant functional categories were biological processes, cellular components, and molecular function with 4493 (48.24%), 2762 (29.65%) and 2059 (22.1%) unigenes associated with each functional category, respectively. Accordingly, a wide distribution and assignment of unigenes were implicated in biological process. Based on the GO enrichment analysis results, the top 54 over-represented genes were cellular processes (5212), metabolic processes (4516), single organism processes (4165), biological regulation (2113), regulation of biological process (1992), and localization (1823) respectively. On top of that among the genes associated with molecular functions, 5389, 3902, and 848 unigenes were involved in binding, catalytic and transport activities respectively. In light of cellular components, most unigenes were associated in cells (2777), cell parts (2777), organelles (1922), and membranes (1914) (Figure 4).

Eukaryotic Orthologous Groups (KOG) were extrapolated from Clusters of Orthologous Groups (COG) to specifically predict the sequence of unknown genes from previously identified and stored in known databases of evolutionarily related sequences of different organisms. Based on KOG, 6575 genes were annotated into known (6126 genes) and unknown (449 genes) functional ortholog groups. Among the known functional groups, the highest percentage of genes (914,13.9%) was annotated to the general function prediction, followed by signal transduction mechanisms (800, 12.17%), post-translational modifications, protein turnover, and chaperones (605, 9.2%) (Figure 5).

3.3. Metabolic Pathway Analysis by Kyoto Encyclopedia of Genes and Genomes (KEGG)

The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway was used to describe the structural elements of the gene products and metabolic pathways in cell. Based on KEGG, 5587 gene transcripts were assembled from the transcriptomic sequences of the heads of M. dirhodum aphids. The transcripts of these genes were annotated to 228 functional pathways.

Based on KEGG pathway annotations, the transcripts of genes associated with the KEGG pathway were grouped into five functional categories: Cellular Processes, Environmental Information Processing, Genetic Information Processing, Metabolism and Organismal Systems. In the cellular process category, the highest number of unigenes were associated with the transport and catabolism (347, 5.74%) and cellular (262, 4.33%) functional pathways, respectively. In the environmental information processing category, significant number of unigenes were implicated in signal transduction (718, 11.88%) while extremely few numbers of genes involved in membrane transport (55), signal molecules and interaction (118). Of these genes associated to genetic information processing category, the highest number of genes (410, 6.78%) were involved in Translation and folding, followed by sorting, and degradation (330, 5.46%) KEGG pathways. In the metabolism group, carbohydrate metabolism pathway was significantly enriched with 296 (4.89%) unigenes. In the organismal systems category, the number of the transcripts of genes associated with endocrine system, digestive system, immune system and nervous system were 402 (6.65%) and 253 (4.18%), 230 (3.8%) and 230 (3.8%) respectively (Figure 6). These findings offer significant insight into the environmental information processing pathways within the head tissue of M. dirhodum. They may serve as a foundation for future research aimed at identifying and characterizing the genes and transcripts involved in key signal transduction pathways.

3.4. Prediction of Secretory Proteins

Salivary proteins are believed to be involved in plant interactions only when they are secreted during aphid probing and feeding. All genes were compared with the priority order of non-redundant protein sequences (NR protein), Protein family (Pfam) and the Swiss-Prot protein library to analyze the presence of signal peptides. Therefore, in this study, 705 putative secretory genes were predicted of all unigenes assembled from the transcriptome of the head tissue. These proteins which had signal peptides associated to the N-terminal of amino-acid sequence along with one and /or zero transmembrane domains, indicating their secretory nature (supplementary files Table S2). Of all clusters of transcripts of secretory genes, 504(71.5%), 487 (69.08%), and 429 (60.85%) genes were characterized as encoding known secretory proteins with known functions in the Nr, Pfam, and Swiss-Prot databases, respectively. However, 193 (27.37%), and 7 (0.99%) transcripts of genes were predicted to encode known proteins with uncharacterized functions in the Nr, and Swiss-Prot databases, respectively. In contrast, the type and function of gene transcripts of some secretory genes were not known in the NR (8 gene), Swiss-Prot (269 genes) and Pfam (218 genes) databases. Most functionally annotated products of putative secretory genes were highly conserved (orthologous) among different aphid species, including pea aphids. Thus, they had similar biological and molecular functions in plant-aphid interactions during feeding.

3.5. Salivary Proteins and Sequence Similarity among Aphid Species

The release of aphid salivary proteins plays indispensable role during aphid-plant interaction. A thorough examination of the saliva from M. dirhodum uncovered the existence of 30 salivary proteins (Table 3). However only 16.67% of the salivary proteins (Glucose dehydrogenase like-protein 2, Heat shock protein, Protein slit, Putative sheath protein, and Vitellogenin domain-containing protein) had signal peptides associated to the N-terminal of amino-acid sequences of each protein. It is not always true that all salivary proteins originate in the head tissue of sap sucking insects. As hemolymph is continually circulated through the salivary glands, it is likely that various biological molecules, including proteins originate from different tissues, are transported into the salivary glands along with the hemolymph or through other, currently unknown, mechanisms.

The identity of the majority of these proteins were determined based on sequence homology search (38.8%) followed by transcript expression levels (13.3%) and sequence prediction (20%). Among the salivary proteins actin had depicted high similarities with 4066 insects and micro-organisms. The sequence of salivary proteins like Elongation factor 1-alpha (1691), 60S acidic ribosomal protein P (1611), Tubulin beta chain (391), Heat shock protein (331) and ACYPI010077 protein (288) had also shown close resemblance with protein sequence of considerable insects and micro-organisms at 50% similarity level (Table 3). The biochemical activity and property of these salivary protein comprises like enzymes, binding proteins, putative effectors, and regulatory proteins (mainly transcription factors). The presence of salivary proteins like glucose dehydrogenase, elongation factors, putative sheath, protein slit, heat shock protein, actin, histone and odorant receptor proteins in the saliva of M. dirhodum, which was widely found in other phloem feeding insects and plant pathogenic micro-organisms, indicating us the significance of salivary proteins in determination of behavioral and functional interaction of aphids with crop plants. Overall, we observed a significant level of similarity across species, typically exceeding 50% identity between the rose grain aphid and the rest of phloem-feeding insects and plant pathogens. This allowed us to anticipate conserved salivary proteins and contribute to a deeper comprehension of the molecular processes involved in the interaction between aphids and host.

The discovery of diverse salivary proteins in the liquid saliva of sap-sucking insects led to the exploration of homologous proteins through sequence comparisons among aphid species. The BLAST score ratio analysis was employed to visualize the extent of salivary protein similarity among different aphid species. The findings revealed that 30% of the salivary proteins were conserved among the three or more aphid species (Table 4). Using the predicted salivary proteins of the rose grain aphid as a query sequence, it was found that 53.3% of the total salivary proteins exhibited close similarity with proteins in the cowpea aphid, Aphis craccivora, while 36.67% of these proteins showed high similarity to those in A. pisum and the cotton aphid, Aphis gossypii. Nevertheless, only 33.3% and 30% of the total salivary proteins exhibited close similarities with sugar cane aphids (Melanaphis sacchari and Sipha flava). Meanwhile, only few numbers of salivary proteins of rose grain aphid, 3.3% and 6.67% displayed close similarity with peach aphid, Myzus persicae and oat aphid, Rhopalosiphum padi respectively.

In general, of the total salivary proteins discovered in the saliva of M. dirhodum, 90% had shown close resemblance to aphids that primarily depend on pulse crops as source food and sustenance, (Aphis craccivora, Aphis glycines, A. pisum, and A. gossypii). In contrast, only 50% of the salivary protein sequences of the rose grain aphid displayed close similarities with cereal aphids, including S. graminum, S. avenae, R. padi, M. persicae, Melanaphis sacchari, and the yellow sugar cane aphid, Sipha flava (Table 4). The close similarity in the biological makeup and composition of salivary proteins between rose grain aphids and pulse crop aphids may indicate that they employ comparable functional and/or behavioral adaptation strategies, setting them apart from cereal aphids.

4. Discussion

In this research, we employed a multidisciplinary strategy to characterize the transcriptome sequence of rose grain aphid, M. dirhodum, and to identify potential secretory proteins. The study of the rose grain aphid's head transcriptome and salivary proteome revealed four distinct trends: (1) a significant number of head transcriptome secretory proteins possess signal peptide sequences, (2) some portion of the head transcriptome proteins associated with secretory proteins have unknown functions, (3) with the exception of five proteins most saliva proteins had no signal peptides associated to their sequence, and (4) Except five proteins, most secretory proteins predicted in the transcriptome of the head tissue could not detected in the saliva of rose grain aphid. This inconsistency is likely attributable to the highly specialized role of aphid saliva to adapt to the diverse and evolving defense mechanisms of plants.

Some of the inconsistencies between the salivary gland transcriptome and proteome might stem from the differences in expression levels between gene transcripts and proteins. Expression profiles of mRNA are highly dynamic, and there is often a lack of direct correlation between mRNA expression and protein levels [26]. Therefore, discrepancies between the expressed transcript and the resulting protein are likely to be common. Nonetheless, the goal of our study was to provide an initial catalogue of potential secreted proteins from the head tissue, utilizing both transcriptomic and proteomic data. The integration of these data sets marks a significant preliminary step in identifying candidate secretory proteins for future aphid research.

4.1. Putative Effectors

At present, there is no functional data to clarify the roles of the aphid salivary proteins identified in this study. Nevertheless, possible effector functions for many of these candidate proteins can be inferred from their homology or similarity to effector proteins involved in pathogenesis or parasitism, as observed in other sap-sucking insects and plant pathogens. Host plants meticulously detect and respond to incoming injuries caused by the release of secretory proteins and mechanical damage from herbivorous insects. However, plant-feeding insects subsequently suppress the immune responses triggered by these elicitors. Consistent with the finding of Huang, et al. [27] putative sheath proteins (encoded after the gene shp-1) were predicted in both saliva and transcriptomic sequence of rose grain aphid heads. According to Huang, et al. [27] putative sheath protein (LsSP1) binds to the salivary sheath via mucin-like protein (LsMLP1) when released into plants during feeding and probing [28]. The protein specifically synthesized in the head tissue of sap sucking insects and release into plant to suppress plant defense system.

Secretory proteins like glucose dehydrogenase protein-1 (secretory proteins encoded after the gene gld-1) and elongation factor 1 alpha (regarded as functional product of EF1A_0 gene) were also detected in the secretome and transcriptome of rose grain aphid heads. Glucose dehydrogenase like protein 1 facilitates oxidation-reduction process [29]. Ubiquitous plant defense responses against invading plant pathogens and insect pests may be degraded by proteins involved in oxidoreductase activity [30]. Thus, it is a potential effector that deters plant defense responses against aphid infestation. Glucose dehydrogenase like protein 1 also plays indispensable role in insect development and suppressing plant defense mechanisms.

Lipid-binding and transporting proteins (apolipophorins), (Cluster-3959.1281), was predicted in the transcripts of secretory genes. This protein is homologous among some aphid species like potato aphid Myzus euphorbiae and M. persicae. The fecundity rate of M. persicae was increased following the overexpression of Me10 in N. benthamiana, indicating species specific ability to deter plant defense responses [31]. Beta-glucosidase transcript (Cluster-9628.0) were predicted from the transcriptome of M. dirhodum head tissue. It is a typical enzyme in termites (Neotermes koshunensis) that aids in degrading hemicellulose [32]. On top of that it may also serve as a signaling pathway that triggers plant defense responses against herbivores. Treatment of cabbage leaves with beta-glucosidase induces the release of volatile compounds like salicylic acid (SA), ethylene, and H2O2 [33].

Gene transcripts specifically, Cluster-3959.2598, Cluster-3959.726, Cluster-1501.0, Cluster-45.0, Cluster-671.0, and Cluster-4581.0 were predicted encoding of glucose dehydrogenase. By the same token, glucose dehydrogenase was also detected in the saliva of rose grain aphid (glucose dehydrogenase-like protein 1 and glucose dehydrogenase-like protein 2). These proteins were reported in the saliva and salivary glands of some aphid species like M. persicae [6], S. avenae [30,34], and A. craccivora [35]. Glucose dehydrogenase is a potential effector that deters pattern triggered immune responses against aphid infestation.

The transcripts of lipase (Cluster-7728.1) and phospholipases (Cluster-3959.2585, Cluster-10963.0 and Cluster-7682.0) were extrapolated from the transcriptomic sequence of M. dirhodum heads with salivary glands. Lipase and phospholipase hydrolyze the linkages of triacylglycerols and phospholipids, respectively [36]. Unlike phospholipase, lipid processing, food digestion and transport has been enhanced by the activity of lipase [33]. The transcript of secretory proteins (Cluster-3959.2585, Cluster-12716.1, Cluster-10963.0, & Cluster-7682.0), phospholipase often involve in lipid biosynthesis, acyltransferase activity (catalysis of the transfer of an acyl group to an oxygen atom), calcium ion and receptor binding activities. The calcium binding activity of the phospholipase protein subsequently led to prevention of calcium-mediated occlusion and obliteration of phloem sieve elements against aphid feeding. Characterization of putative effectors of pea aphid using proteome mass spectrometry indicated both metalloprotease and calcium-binding proteins inactivate plant defense responses and inhibit calcium-mediated blockade of plant sieve elements, respectively [9]. However, both lipase and phospholipase have versatile roles in various biological and industrial applications. The biproduct of phospholipase, phosphatidic acid, plays an important role in signal transduction cascades and lipid metabolic pathways, thus determining plant response to stress. When plants are exposed to stress such as mechanical stress (wounding and frost) and biological stress (pathogen and insect attack), the expression and activity of the phospholipase D protein increases rapidly [37]. Based on the findings of [38], rice infection with Xanthomonas oryzae increased phospholipase accumulation in its plasma membrane. Phospholipases also trigger plant defense responses (release of abscisic acid, ethylene, and nicotinamide adenine dinucleotide phosphate ion (NADPH)) against herbivorous insects and diseases [37]. The transcript of trehalose (Cluster-608.0) facilitates the expression of trehalose transporter Tret1-2 homolog isoform X2, trehalose-phosphatase, and alpha-trehalose-phosphate synthase. In plants, trehalose plays a regulatory role in sugar metabolism, growth, development, and responses to stress induced by biotic and abiotic factors [39]. The overexpression of trehalose biosynthetic gene enacts plant and microbes to receive exogenous trehalose, thereby increase tolerance to stress [40,41]. Trehalose has also been reported in the salivary glands of B. tabaci and S. avenae [42,43]. Based on the findings of [43], trehalose plays an important role in deterring the accumulation of excess trehalose in plant cells to suppress immune response.

4.2. Detoxifying Secretory Proteins

When plants are attacked by herbivorous insects, they release secondary toxic metabolic products (mainly protease inhibitors). Nonetheless, herbivorous insects secrete detoxifying proteins to degrade the toxins produced upon feeding, thus paving the way for adaptation to adverse conditions. Esterase is one of the detoxifying proteins released by aphids manipulate plant secondary metabolites and insecticides [44]. The transcripts of secretory genes (Cluster-6992.0, and Cluster-6189.0) were predicted to encode secretory protein, esterase E4. This protein breaks down the effect of chemical pesticides when farmers sprayed on farm to manage aphids. A prevalent strategy utilized by Myzus persicae aphid to build resistance against insecticide involves boosting the expression of genes that encode for Esterase E4 [45].

Secretory proteins such as glutathione S-transferase 6 and CAAX prenyl protease 1 homolog (ortholog metalloproteases) were also predicted in the transcriptome assembly of the head tissue after gene transcripts of Cluster-11681.0 and Cluster-3959.1466. Glutathione S-transferases allow insects to survive under chemical stress through metabolism of xenobiotics or by providing protection against oxidative stress [46]. Although there is insufficient evidence to support our findings, a report [47] indicated, many types of metalloproteinases are involved in detoxifying plant secondary metabolites.

Clusters of gene transcripts (Cluster-7473.0, Cluster-477.0, Cluster-3959.111, Cluster-3959.1431, Cluster-3959.824, and Cluster-3959.2743) were also predicted encoding secretory proteins, specifically peroxidase and oxidoreductase. These proteins are essential for basic biological functions, such as response against stress and defense mechanisms. Previously, the secretory proteins were discovered in some aphid species, mainly in S. avenae [43], A. pisum [9], Mayetiola destructor [48] and Megoura viciae [47]. Unlike other species of aphids, only superoxide dismutase (Cluster-3959.779) was predicted in the transcriptomic sequences of M. dirhodum. Research conducted on parasitoid wasp (Leptopilina heterotoma), revealed that extracellular superoxidase dismutase was synthesized and released along with the fluid and act as virulence factors to counteract the immune response of Drosophila species [49]. It is also used by pathogenic bacteria, fungi, and pathogenic protozoan parasites as virulence factors against the resistant gene of the host at early infection. This protein initiates plant defense responses when aphids begin to probe and feed on plants. Thus, it is considered to be a signaling pathway in the plant defense system. Peroxidase transcripts (Cluster-3959.111, Cluster-3959.824, and Cluster-7473.0) were predicted in transcriptomic sequences of secretory proteins from M. dirhodum aphids. Based on the findings of [43], the presence of peroxidase in aphid saliva serves as an antioxidant involved in the scavenging of hydrogen peroxides and thus plays significant role as a detoxifying enzyme and in suppressing reactive oxygen species-induced plant defense responses.

4.3. Digestive Secretory Proteins

Probing and feeding activities of aphids are also enhanced by digestive enzyme release into the host. The transcript of serine protease (Cluster-12396.0) was assembled from the transcriptomic sequences of rose-grain aphids. Serine protease is a digestive enzyme secreted by the midguts of phloem-feeding insects to facilitate the digestion of food assimilates while ingesting phloem saps from the plant [50]. Serine protease also acts against host-related serine protease inhibitors and phenol oxidase defense response. Secretory gene transcripts, namely Cluster-1194.0, Cluster-1892.0, Cluster-2182.0, Cluster-2350.0, Cluster-3666.0, Cluster-3959.1850, Cluster-3959.2179, Cluster-3959.964, Cluster-6237.0, & Cluster-7233.0 were assembled and predicted encoding digestive secretory proteins, most likely trypsin. Trypsin is secreted by phloem feeding aphids to facilitate the break down and oral digestion of plant phloem constituents [44].

4.4. Ca2+ Binding Secretory Proteins

Putative calcium-binding protein (regucalcin) was predicted in the transcriptional product of secretory gene (Cluster-11867.0) assembled from the transcriptomic sequence of M. dirhodum heads with salivary glands. Similarly, in this experiment we had predicted sarcalumenin (SAR) from secretory gene transcript, Cluster-3959.1170. Sarcalumenin is associated with calumenin calcium buffer proteins to modulate the uptake of calcium ions and release the excitation and contraction of muscle fibers. It also plays an important role in other physiological functions such as muscle resistance against fatigue, muscle development, sarco-endoplasmic reticulum calcium ATPase stabilization, and store-operated calcium entry mechanisms [51]. Plants release calcium ions into the sieve element lumen to block the outflow of phloem saps during mechanical damage, feeding, and probing. In response, sap-sucking insects release water saliva consisting of calcium ion-binding proteins into the lumen of the sieve element to prevent plugging and clogging of feeding sites for insects [52]. With the aid of transcriptomic and proteomic approaches, this protein was detected in the saliva and salivary glands of some homopterous insects such as A. pisum and Nephotettix cincticeps. These results revealed that the synthesis and secretion of calcium-binding proteins could be a means of suppressing plant defense responses against all phloem sap-sucking insects [53].

4.5. Zn-Binding Secretory Proteins

The transcripts of aminopeptidase N (Cluster-8442.0, Cluster-11978.0, and Cluster-7461.0) were detected in the transcriptomic sequences of M. dirhodum secretory proteins. This protein is implicated in binding zinc ions and also known for metallopeptidase activity during aphid-host interactions. Generally, Aminopeptidase N is a multifunctional protein and expressed at different organs and cells of many insects [54]. This protein catalyzes the breakdown of the amide group of amino acids from the N terminal of proteins or peptides. Other versatile types of aminopeptidases have also been identified in M. dirhodum aphid transcripts. Endoplasmic reticulum aminopeptidase 2 was predicted to be synthesized from the transcripts of secretory proteins (Cluster-3959.627) and was found to play important molecular functions, such as metallopeptidase, zinc binding, and hydrolysis. The protein acts against ester bonds of the amino acid, thereby facilitate the disintegration of the disulfide bonds of polypeptides and/or proteins. The functional product of the transcript of gene (Cluster-3959.1656), that is, thyrotropin-releasing hormone-degrading ectoenzyme, is also an integral part of the aminopeptidase protein typically secreted in the salivary glands of insects and functions as the binding of zinc ions (to the active site of the protein, aminopeptidase) and other metallopeptidase activities (penetration of plant cell wall and membrane while probing their stylet to the layers of a cell’s phloem sieve element) at a specific amino acid location of the protein. Clusters of unigenes transcripts, Cluster-3959.1011, Cluster-3959.1012, and Cluster-7461.0 were assembled from the transcriptomic sequences of aphid heads and predicted to encode metallopeptidase and lysosomal alpha-mannosidase proteins which in turn bind zinc. Zinc, a crucial micronutrient vital for the growth, development, and defense of all living things, plays a significant role in the interaction between hosts and pests in plant systems, where competition for Zn can influence the outcome [55].

4.6. Reproduction and Development

Vitellogenin (encoded after gene vit-6) was detected in the proteome and transcriptome (Cluster-1741.0) of M. dirhodum head tissue. It is synthesized in the fat body, and plays multiple roles in insect reproduction, mainly during oocytes and embryo development. Vitellogenin receptors (VgR) present on the surface of oocytes, are responsible for vitellogenin transportation from haemolymph to oocytes [56]. The vitellogenin genes have been extensively investigated in root-knot nematodes, Caenorhabditis elegans [57]. Clusters of secretory gene transcripts, Cluster-3959.2780, & Cluster-4827.0, were also predicted to encoding secretory proteins (protein slit). On top of that it was detected in the secretome of rose grain aphid. This secretory protein plays essential role as midline repellents, stopping longitudinal axons from crossing the central nervous system's midline in a range of insects, nematodes, and planarians [58].

4.7. Protein Synthesis and Secretion

Disulfide isomerase was identified through transcriptomic analysis of head tissues, but it was not notably present in the proteomic profile of salivary proteins. Protein disulfide isomerases are widespread and versatile enzymes found in eukaryotic organisms. These enzymes belong to the thioredoxin superfamily. Classical disulfide isomerases are highly conserved and include an N-terminal signal peptide [59]. This enzyme is believed to play a role in protein folding regulation. Additionally, disulfide isomerase is found in the salivary secretions of plant-parasitic nematodes and has been linked to an increased production of salivary proteins [60,61].

4.8. Odorant Binding and Chemosensory Proteins

In addition to visual signals, aphids largely depend on certain metabolic products as cues to identify their host within complex and spatially diverse environments [62]. In line with the information of Shih, et al. [62], chemosensory protein was predicted among the secretary genes (Cluster-3959.2546) assembled from the transcriptomic sequence of M. dirhodum aphids. Previously, chemosensory proteins were reported in the transcriptome of salivary glands of N. lugens, M. persicae, M. cerasi, and bird cherry-oat aphid R. padi [63]. According to the review of Shih, et al. [62], multiple gene products typically chemosensory proteins, odorant binding proteins, gustatory receptors, and olfactory receptors are important to locate and exploit host resources following the release of phenolic compounds from the plant. In insects, the transportation of organic compounds that convey specific chemical signals (semio-chemicals) is performed with the help of chemosensory proteins. Aphids communicate with each other and locate their food sources and ovipository sites with the aid of odorant binding and chemosensory proteins expressed and released through their antennae and mouth [64]. However, odorant binding and chemosensory proteins may also be located in other parts of insect tissues or organs, mainly the leg, head, thorax, abdomen, and salivary glands of the head. The function of these proteins may differ depending on the location of their expression, in which they may aid in insect development, regeneration of insect legs, host interaction, and immune responses [65,66]. The product of unigene transcripts (Cluster-3959.2546), i.e., chemosensory proteins, predicted from the transcriptome sequence of M. dirhodum showed high similarity to proteins reported in S. avenae and M. persicae aphid species. According to a functional assay carried out on Nicotiana benthamiana, the overexpression of Mp10 protein in N. benthamiana suppressed bacteria-associated molecular patterns and local cell death that impede the fecundity of M. persicae [67]. The odorant binding and chemosensory proteins of Spodoptera frugiperda recognize chemical, behavioral and/or physiological responses of plants [68].

5. Conclusion and Recommendation

The successful colonization of plants by invading pathogen hinges on a variety of secretory proteins released into the host to either suppress or modulate specific innate defenses. This study aimed to determine whether the rose grain aphid possesses a similar array of potential secretory proteins as other phloem-feeding insects. The analysis of the rose grain aphid head tissue secretome revealed three key characteristics: (1) a significant portion of head tissue proteins possess a peptide secretion signal; (2) six salivary proteins identified in the saliva of rose grain aphid consisting of signal peptides associated to the N-terminal of amino-acid sequences; and (3) there have been notable similarities in secretory protein composition between rose grain aphid head tissue and other aphid species. A transcriptomic analysis of M. dirhodum aphids revealed 18,030 genes from 31,344 transcripts in the head tissue along with salivary glands, with 705 genes encoding secretory proteins. Notably, 28.5% of the total gene transcripts were deem of unknown functions in the Nr database. The proteomic analysis unveiled the discovery 30 salivary proteins, including enzymes, binding proteins, putative effectors, and regulatory proteins, highlighting their role in aphid-plant interactions. Despite lacking signal peptides, salivary proteins such as TTF-type domain-containing protein and ATP synthase subunit were found in the watery saliva, suggesting unknown secretion mechanisms. 30% of these proteins are conserved across three or more aphid species. Furthermore, 53.3% of the total salivary protein showed close similarity to Aphis craccivora, followed by 36.67% of these proteins had close similarity with A. pisum and Aphis gossypii. However, only 3.3% and 6.67% of these proteins were similar to those of Myzus persicae and Rhopalosiphum padi, respectively. 90% of the salivary proteins of rose grain aphid closely matched to aphids targeting leguminous crops. Hence, rose grain aphid may share functional and/or behavioral adaptive strategies with pulse aphids. The study underscores the need for further research on gene expression profiling and functional analyses of salivary proteins to elucidate the nature, expression and their role in aphid-plant interactions.