A Detailed Spatial Expression Analysis of Wing Phenotypes Reveals Novel Patterns of Odorant Binding Proteins in the Soybean Aphid, Aphis glycines

The wide range of insect niches has led to a rapid expansion of chemosensory gene families as well as their relatively independent evolution and a high variation. Previous studies have revealed some functions for odorant-binding proteins (OBPs) in processes beyond olfaction, such as gustation and reproduction. In this study, a comparative transcriptomic analysis strategy was applied for the soybean aphid, Aphis glycines, focusing on various functional tissues and organs of winged aphids, including the antenna, head, leg, wing, thorax, cauda, and cornicle. Detailed spatial OBP expression patterns in winged and wingless parthenogenetic aphids were detected by RT-qPCR. Twelve OBPs were identified, and three new OBPs in A. glycines are first reported. All OBPs showed comparatively higher expression in sensory organs and tissues, such as the antenna, head, or leg. Additionally, we found some novel expression patterns for aphid OBPs (Beckendorf et al., 2008). Five OBPs exhibited high-expression levels in the cauda and four in the cornicle (Biasio et al., 2015). Three genes (OBP2/3/15) were highly expressed in the wing (Calvello et al., 2003). Two (OBP3/15) were significantly more highly expressed in the wingless thorax than in the winged thorax with the wings removed, and these transcripts were significantly enriched in the removed wings. More details regarding OBP spatial expression were revealed under our strategy. These findings supported the existence of carrier transport functions other than for foreign chemicals and therefore broader ligand ranges of aphid OBPs. It is important for understanding how insect OBPs function in chemical perception as well as their other potential physiological functions.


INTRODUCTION
The olfactory system plays an important role in directing insect behaviors, such as foraging, mating, oviposition, and predation. Similar to other insects, aphids, especially those with a migratory phenotype (winged morph), rely heavily on chemical signals, including plant volatiles and speciesspecific pheromones, to locate hosts, find mates, and avoid natural enemies (Pickett, 2009). In general, odorant-binding proteins (OBPs) are involved in the first step of olfactory recognition (e.g., Pelosi et al., 2018). They bind and transport external odors through the hemolymph to activate corresponding olfactory receptors, which are responsible for transmitting environmental chemicals into electrophysiological signals. As one of the most important groups of chemo-reception proteins in insects, OBPs have been studied since 1981 (Vogt and Riddiford, 1981). Regarding aphids, OBPs have been widely reported for Acyrthosiphon pisum (Ji et al., 2016) genes encoding complete OBPs with a signal peptide, Zhou et al., 2010), Myzus persicae (Ji et al., 2016), Megoura viciae (Iovinella et al., 2011;Daniele et al., 2018), Sitobion avenae (Kim et al., 2013;Xue et al., 2016), Aphis gossypii (Grabherr et al., 2011;Gu et al., 2013), and Aphis glycines (Grabherr et al., 2011;Wang et al., 2019). Consistent with the simple facultative parasitic lifestyle of aphids, the aphid OBP family has few members, and they are generally more highly conserved than in other insects. The spatial expression profiles of these proteins have also been broadly investigated. OBP expression is not limited to the chemosensory system (e.g., Xue et al., 2016;Gao et al., 2018;Wang et al., 2019), but also occurs in non-sensory tissues and organs, such as the wings (Calvello et al., 2003;Pelosi et al., 2005;Wang et al., 2020), reproductive organs (Li et al., 2008;Sun Y. F. et al., 2012), mandibular glands (Iovinella et al., 2011), and salivary glands .
The soybean aphid, A. glycines, is an important phytophagous pest that feeds on plants by sucking sap from leaves, stems, and pods, significantly reducing soybean yield and quality (Wang et al., 1996;Beckendorf et al., 2008;Wu et al., 2009). Moreover, in soybean plants heavily infested with aphids, sugary excretions ("honeydew") produced by aphids indirectly damage plants by reducing photosynthesis . Plant viruses can be transmitted during aphid infestations (Clark and Perry, 2002). Accordingly, A. glycines is now used as a model for studying the evolution of biotypic virulence (Wenger et al., 2017). In a recent study, we attempted to identify the OBP family of soybean aphids based on a homologous cloning strategy (Wang et al., 2019). In this study, a more detailed comparative transcriptomic analysis strategy focusing on various functional tissues and organs of winged aphids, including the antenna, head, leg, wing, thorax, cauda, and cornicle, was successfully applied to A. glycines. In addition to identifying more OBPs, detailed spatial expression patterns of both winged and wingless parthenogenetic aphids were analyzed by RT-qPCR and the findings were discussed.

Insects and Sampling
Aphis glycines was reared from parthenogenetic aphids initially collected on soybean plants at the Minzhu Experimental Station, Heilongjiang Academy of Agricultural Sciences, China, and cultured in an air-conditioned insectary [24 ± 1 • C, 75 ± 5% relative humidity (RH), 16-h light:8-h dark photoperiod]. Newborn aphids (0-12 h) were transferred to new plants for synchronization of developmental stages. The insects were reared for 7 days at different densities (20 aphids/plant and 100 aphis/plant) and winged and wingless aphids were separately collected in the last developmental stage (adult). For transcriptome sequencing, antenna (500), head (200), wing (500), leg (500), thorax (200), cornicle (200), and cauda (500) specimens of winged aphids were carefully dissected under the microscope. For RT-qPCR analysis, the same tissues or organs of adult wingless and winged morphs were collected. Each experiment was carried out in biological triplicate. Samples were flash-frozen in liquid nitrogen and stored at -80 • C until RNA extraction.

Total RNA Extraction and Illumina HiSeq 4000 Sequencing
Total RNA was isolated from the above samples using TRIzol (Invitrogen, Carlsbad, CA, United States), and DNA fragments were removed with RNase-free DNase I (Takara, Shiga, Japan). An Agilent 2100 Bioanalyzer (Agilent Technologies, Carlsbad, CA, United States) was used to determine the concentrations, integrity, and 28S/18S values of the RNA samples, and a NanoDrop 2000 spectrophotometer (NanoDrop products, Wilmington, DE, United States) was used to access purity. mRNA was then enriched using oligo (dT) beads (Agilent Technologies) and fragmented using fragmentation buffer (Agilent Technologies) and then used for the synthesis of firststrand cDNA. After purification and the repair of cohesive ends, the DNA samples were ligated to adapters, and fragment selection and PCR amplification were conducted. The final quality assessment was performed using an Agilent 2100 Bioanalyzer (Agilent Technologies). Three DNA libraries were examined using the Illumina HiSeq 4000 sequencing platform.
Fragmentation involved the use of divalent cations under elevated temperature in NEBNext and first-strand synthesis reaction buffer (5×). Single-stranded (ss) cDNA was synthesized using a random hexamer primer, M-MuLV reverse transcriptase, DNA polymerase I, and RNase H (NEB, United States). After the adenylation of the 3 ends of the fragments, NEBNext adaptors with a hairpin loop structure were ligated for hybridization. The library fragments were purified using the AMPure XP system (Beckman Coulter, United States), selecting cDNA fragments 150-200 bp long. Then, 3 mL of USER enzyme (NEB, United States) was applied to the size-selected, adaptor-ligated cDNA, and the reaction was incubated at 37 • C for 15 min, followed by 5 min at 95 • C before PCR. PCR was then performed using Phusion high-fidelity DNA polymerase, universal PCR primers, and an index (X) primer. The products were purified (AMPure XP system), and library quality was assessed using an Agilent Bioanalyzer 2100 system (Agilent Technologies, United States). Clustering of the index-coded samples was performed with a cBot cluster generation system using TruSeq PE Cluster Kit v3-cBot-HS (Illumina, China) according to the manufacturer's instructions. The library preparations were sequenced on the Illumina HiSeq 4000 platform, and pairedend reads (PE125 sequencing strategy) were generated after cluster generation.

RNA-Seq Data Generation and de novo Transcriptome Assembly
After sequencing, the raw reads were processed by NGS-QC to remove low-quality sequences (≥ 15% bases with Q ≤ 19), excess adaptors (≥ 5 bp adaptor bases in reads), and reads with a high content of unknown bases (≥ 5%; CASAVA FASTQ files). The clean reads were then assembled into unigenes using Trinity r20140413p1 with min_kmer_cov:2 and the other parameters set to the default values (Langmead and Salzberg, 2012). Gene expression levels in each sample were estimated by RSEM (Li and Dewey, 2011): (I) Clean data were mapped to the transcript sequence, and (II) the read count for each gene and isoform was obtained from the mapping results. The fragments per kilo base per million (FPKM)-mapped reads value of each gene was calculated based on gene length and the mapped read number using HTSeq v0.5.4p3 and Cufflinks v2.2.1 (Mortazavi et al., 2008).

Differentially Expressed Genes and Annotation of OBP-Encoding Transcripts
Reads for the A. glycines transcriptomes from seven different tissues (antennae, head, wing, leg, thorax, cornicle, and cauda), with three replications for each tissue, were produced based on next-generation sequencing (NGS) results. Expression analysis was performed using TopHat and Cufflinks (Trapnell et al., 2012;Kim et al., 2013). Differential expression analyses comparing each tissue to the antenna were separately performed using the DESeq R package (version 1.10.1), which provides statistical routines for determining differential expression using a model based on the negative binomial distribution. To control the false discovery rate, the resulting p-values were adjusted using Benjamini and Hochberg's approach. Genes with a fold change (FC) > 2 and an adjusted p-value < 0.05 according to DESeq analysis were considered DEGs. The log2 (fold change) values and p values are shown in a volcano plot.
We used the BLASTx program of the National Center for Biotechnology Information (NCBI, 1 ) to predict genesencoding OBPs.
The basis of the annotation was a hand-curated database of OBPs containing known aphid candidate OBP sequences. The assembled sequences were compared with the reference dataset using BLASTx. All sequences that generated a hit were further analyzed by a motif search program based on a 5-6 conserved OBP cysteine pattern consisting of C1-X 15−39 -C 2 -X 3 -C 3 -X 21−44 -C 4 -X 7−12 -C 5 -X 8 -C 6 for OBPs (Zhou et al., 2008).

Real-Time Quantitative PCR
Real-time quantitative PCR (RT-qPCR) experiments were carried out using a 7,500 Fast Real-Time PCR System (Applied Biosystems-Life Technologies, Carlsbad, CA, United States) and the cDNA samples prepared from winged and wingless aphid antennae, heads, legs, cornicles, caudae, and wings (only for winged morphs). Two reference genes, 18S rRNA and GAPDH dehydrogenase, were used for normalizing target gene expression and correction for sample-to-sample variation (Wang et al., 2019). Specific primers were designed for each A. glycines OBP gene and for the two reference genes using Primer Premier v5.0 software; the primer information is listed in Supplementary Data 1. PCR amplification was conducted in a volume of 20 µL containing 10 µL of 2 × SYBR Mix, 1 µL of diluted cDNA template, 7.8 µL of PCR-grade water, and 0.6 µL of each primer at 10 µM. The PCR conditions were as follows: 95 • C for 30 s, 40 cycles of 95 • C for 15 s, 60 • C for 30 s, and 72 • C for 45 s. The OBP expression status was calculated using the 2 − Ct comparative CT method (Livak and Schmittgen, 2001), and a CT value greater than 35 was considered no expression. The fold changes of OBPs in the tissues of both winged and wingless morphs are reported relative to the antennal transcript levels of OBP3 in the wingless morph (wingless antennal OBP3). Means and standard deviations were calculated using data from experiments performed in triplicate, and the results were presented as n-fold differences in expression. Differences in transcriptional characteristics among various OBPs in different tissues were analyzed using SPSS 16.0. Statistical significance was determined using one-way ANOVA and post hoc Duncan multiple range tests. Significance was established at p < 0.05. External reference genes were randomly chosen from among OBPs to first perform a preliminary assessment, after which, we defined those with broader expression profiles, such as wingless OBP2 and OBP3, as candidate external genes. Wingless OBP3 was ultimately chosen as the external gene.

Differential Expression Analysis
First, antenna-specific genes were compared among different groups and further analyzed using Venn diagrams. Although genes were found to exhibit antenna specific expression, there were no OBPs (Figure 2A; gene lists see Supplementary Data 3).
Next, genes specifically expressed in each tissue were screened by the same strategy, and the numbers showing specific expression in the heads, legs, wings, thoraxes, caudae, and cornicles were 226, 2,005, 580, 1,735, 1,667, and 1,741, respectively, also with no OBPs included ( Figure 2B, see Supplementary Data 3 for gene lists).

OBP Prediction and Functional Enrichment
Odorant-binding proteins genes were predicted by homology comparison, the relative expression levels (FPKM) of target genes in different tissues were visualized by a heat map, and the tissue or organ specificity of expression was preliminarily analyzed.
Similar to peach aphids (Ji et al., 2016), OBP1 was neither predicted nor identified in A. glycines. Twelve OBPs FIGURE 1 | Volcano plots for differentially expressed genes between antennae and each of the other six tissues (heads, legs, caudae, cornicles, thoraxes, and wings). Frontiers in Physiology | www.frontiersin.org were identified (GenBank accession numbers MW924836-MW924846, and MW930727) using the NCBI BLASTX program and named according to their ortholog names in other aphids, including three newly reported OBPs: OBP13, OBP14, and OBP15. The OBPs included in the heatmap and the sequences of these genes are listed in Supplementary Data 4; their FKPM values are provided in Supplementary Data 5. The heat map in Figure 3 illustrates that most of the OBPs, such as OBP2-OBP10, OBP13, and OBP14, were found to mainly be expressed in sensory organs and tissues (i.e., antennae, heads and legs) but that OBP15 showed relatively low expression in all specimen types. Our results showed that OBPs are also expressed in organs such as caudae and cornicles, which are not chemical sensory organs. As OBP3 and OBP7 were assembled into one transcript (DN1170_c0_g1_i8, Supplementary Data 4) following the TRINITY instructions (Grabherr et al., 2011), their gene expression values were ultimately quantified as equal and then were corrected by RT-qPCR.
In addition to its remarkably higher expression in winged antennae, OBP2 was found to be systemically expressed in all tissues of both morphs, including the antenna, head, leg, wing, thorax, cornicle, and cauda ( Figure 5A). OBP3 was also systemically expressed, with significantly higher expression in the wingless head, thorax, and cauda, and therefore showed a winged aphid-specific expression pattern in those tissues ( Figure 5B, p < 0.05, N = 3). OBP4 was expressed at quite a low level but still showed a phenotypic correlation with wingless aphids, with comparatively high expression in the antenna, leg, and cauda. OBP5 was leg specific in both morphs, whereas OBP15 was highly expressed in the wings and legs of both winged and wingless morphs (Figure 5).  In addition to OBP2, OBP6 maintained high expression levels in the head of winged and wingless aphids ( Figure 4B). Specifically, OBP2 was more highly expressed in the winged aphids head and OBP6 in the wingless aphid head.
Similar to OBP3, OBP15 was more highly expressed in the wingless aphids thorax. OBP2 was highly expressed in both winged and wingless morphs.
In the leg, the expression of most OBPs was quite low and showed no phenotypic correlation. Among them, OBP2/5/6 displayed comparatively high expression, with OBP5 being significantly more highly expressed in the leg than in other tissues ( Figure 4C).
Surprisingly, we found that three OBPs, OBP2, OBP3, and OBP15, were expressed at much higher levels in the wing than other OBPs (Figure 4D).
The results for the cornicle indicated that five OBPs (OBP2/3/6/9/14) were highly expressed; among them, OBP2/6 showed differential expression between morphs, though they all presented significantly elevated expression in the winged morph. Nevertheless, the expression of the other three OBPs (OBP3/9/14) did not differ between morphs. In addition to OBP3 and OBP9, the gene coding the other EBF-binding protein, OBP7, was expressed in the cornicle at relatively low expression levels, with no significant difference between winged and wingless morphs.
OBP3/5/7 generally appeared to be specific to the wingless morph cauda. OBP6 was highly expressed in wingless aphids; OBP2 was also highly expressed, but with no significant difference between winged and wingless aphids.

DISCUSSION
In this study, three newly reported OBPs were identified based on TRINITY, which has been demonstrated to recover more full-length transcripts across a broad range of expression levels, with a sensitivity similar to the methods that rely on genome alignments (Grabherr et al., 2011). Our differential expression analysis (Figure 1) and enrichment results (Figure 3) showed that OBPs are widely expressed in the soybean aphid, although none was found to be antenna specific or specific to any of seven organs/tissues (Figure 2). The results derived from the aforesaid transcriptome data prompted us to further carry out a detailed FIGURE 5 | RT-qPCR detection of each OBP expression pattern in seven tissues; WL, wingless; WD, winged; An, antennae; Hd, head; Th, thorax; Ci, cornicle; Lg, leg; Ca, cauda; Wg, wing. Bars represent the standard deviation of the mean for three independent experiments. (A-L) are the detection results of each OBP expression pattern in seven tissues and the sequence is from OBP2 to OBP10 and OBP13 to OBP15, respectively. The letters above bars (a-e) are the result of a multicomparison, which indicated significant differences from other samples with different letters (p < 0.05).
investigation and analysis of OBP expression levels among different wing phenotypes and different tissue parts based on RT-qPCR technology. The qPCR results show that relatively high expression of most OBPs in antennae, the main olfaction organ of aphids, of both phenotypes, and further, higher expression in the winged phenotype (Figure 4) are consistent with the fact that winged aphids have more developed olfaction (Pickett, 2009) and support that OBPs play key roles in aphids' olfaction.
Our work provides insight into the potential functions of OBPs correlating with their spatial expression among seven body parts, including various functional organs and tissues, including the antenna, head, leg, thorax, wing, cauda, and cornicle. We further report the breakthrough of the acquisition of aphid cauda and cornicle transcriptomes and their publication. The insect OBP family is believed to participate in chemosensory perception due to their high abundance in chemosensory organs, such as antennae, heads, and legs (e.g., Vogt and Riddiford, 1981;Sun et al., 2013). Without exception, all 12 OBPs identified in this study exhibited relatively high expression in the antenna, head, or leg of both A. glycine wing morphs. Most OBPs showed relatively high expression in the antennae, the main olfaction organ of aphids, in both phenotypes, with higher expression in the winged phenotype (Figure 4), consistent with the fact that winged aphids have more developed olfaction (Pickett, 2009).
OBP3 (Qiao et al., 2009), OBP7 Zhong et al., 2012;Fan et al., 2017), and OBP9 (Qin et al., 2020) in aphids are known for their high affinities for E-β-farnesene, the key component of the aphid alarm pheromone. In this study, the genes encoding the three EBF-binding proteins (OBP3/7/9) showed different antennal expression patterns from each other (Figure 5), providing a new perspective for understanding relationships among them. OBP7 was significantly highly expressed in the antenna of both phenotypes, with higher levels in winged aphids. In contrast, there was no difference in the expression of OBP3 or OBP9 between the two phenotypes. Further analysis showed that OBP3 was systemically expressed; significantly higher expression in the head, thorax, and cauda of wingless aphids was detected. However, OBP9 presented with high expression in the antenna, leg, and cornicle. The higher expression level of OBP7 in the antennae of winged aphids suggests that it may contribute more to the EBF sensitivity of winged aphids. Notably, the genes encoding these three reported EBF-binding proteins were all expressed in the cornicle ( Figure 4G). Cornicles comprising a pair of tubular tissues involved are the alarm pheromone (E-β-farnesene, EBF) storage and release organ, and OBP3, OBP7, and OBP9 may be involved in alarm pheromone activity preservation, release or biosynthesis, and metabolism by binding to and releasing EBF.
Insect OBPs have been reported to act as carrier proteins in the male reproductive apparatus of mosquitoes (Li et al., 2008). After matting, the OBP expressed by male moths is found on the surface of fertilized eggs, which helped the larvae to avoid cannibalistic behaviors . In this study, caudae were dissected with dorsal segments of the distal and abdominal segments, anus, and gonapophysis. Therefore, the high expression of OBPs observed in this tissue suggests potential functions in reproduction or excretion. In addition, carrier proteins function in the binding of or transfer of foreign chemicals or signal ligands.
Odor-binding proteins have also been found to be expressed in wings, such as in Polistes dominula (Calvello et al., 2003), Vespa crabro (Pelosi et al., 2005), and Helicoverpa armigera . Wang et al. further demonstrated lipid binding by OBPs indicating roles beyond their typical functions in the olfactory system to support insect flight activity. In this study, two OBPs (OBP3 and OBP15) were found to be expressed in the thorax of the wingless phenotype and were significantly downregulated in the winged thorax with wings removed (p < 0.05, Figures 4D,E, 5B,L). In addition, the removed wings expressed significantly high levels of both. Hence, there is a possibility that OBP3 and OBP15 were enriched from the thorax to the wings and that they may also function in other ways, such as lipid-binding proteins in the energy supply of flight or carrier proteins which we discussed above in the section on caudae.
Although OBP3/7/9 all exhibit an affinity for EBF, they showed differential expression patterns in this study. OBP3 was extensively expressed throughout the aphid body, OBP7 was antenna-specifically expressed, and OBP9 was highly expressed in the antenna, leg, and cauda. As the cornicle is the alarm pheromone (E-β-farnesene, EBF) storage and release organ, it was not surprising to find that all previously reported EBFbinding proteins were expressed in the cauda (OBP3, OBP9, and low expression of OBP7).
SaveOBP2, SaveOBP4, and SaveOBP5 have been reported to have a limited affinity for wheat volatile benzaldehyde (Zhong et al., 2012). However, no potential ligand has yet been reported for OBP6, one of the most highly expressed OBPs, which suggests that the ligand spectrum for insect OBPs may be far greater than our expectations.
More details regarding OBP spatial expression were revealed under our strategy. These findings supported the existence of carrier transport functions other than for foreign chemicals and therefore broader ligand ranges of aphid OBPs. It is important for understanding how insect OBPs function in chemical perception as well as other physiological functions of OBPs.

DATA AVAILABILITY STATEMENT
The transcriptomic datasets presented in this study can be found in online repositories (http://www.ncbi.nlm.nih.gov/bioproject/ 744282) with temporary submission ID: SUB9958208.