Insights into wing dimorphism in worldwide agricultural pest and host-alternating aphid Aphis gossypii

The worldwide pest Aphis gossypii has three-winged morphs in its life cycle, namely, winged parthenogenetic female (WPF), winged gynopara (GP), and winged male, which are all produced by a wingless parthenogenetic female (WLPF). Most studies on A. gossypii have focused on WPF, while few have investigated GP and male. The shared molecular mechanism underlying the wing differentiation in the three wing morphs of A. gossypii remains unknown. The wing differentiation of WPF was explored in a previous study. Herein, GP and male were induced indoors. The characters of the body, internal genitals, wing veins, and fecundity of GP and male were compared with those of WPF or WLPF. Compared with WLPF, the shared and separate differentially expressed genes (DEGs) were identified in these three-wing morphs. Newly-born nymphs reared in short photoperiod condition (8 L:16D, 18 °C) exclusively produced gynoparae (GPe) and males in adulthood successively, in which the sex ratio was GP biased. A total of 14 GPe and 9 males were produced by one mother aphid. Compared with WLPF, the three-wing morphs exhibited similar morphology and wing vein patterns but were obviously discriminated in the length of fore- and underwings, reproductive system, and fecundity. A total of 37 090 annotated unigenes were obtained from libraries constructed using the four morphs via RNA sequencing (RNA-Seq). In addition, 10 867 and 19 334 DEGs were identified in the pairwise comparison of GP versus WLPF and male versus WLPF, respectively. Compared with WLPF, the winged morphs demonstrated 2 335 shared DEGs (1 658 upregulated and 677 downregulated). The 1 658 shared upregulated DEGs were enriched in multiple signaling pathways, including insulin, FoxO, MAPK, starch and sucrose metabolism, fatty acid biosynthesis, and degradation, suggesting their key roles in the regulation of wing plasticity in the cotton aphid. Forty-four genes that spanned the range of differential expression were chosen to validate statistical analysis based on RNA-Seq through the reverse transcription quantitative real time polymerase chain reaction (RT-qPCR). The comparison concurred with the expression pattern (either up- or downregulated) and supported the accuracy and reliability of RNA-Seq. Finally, the potential roles of DEGs related to the insulin signaling pathway in wing dimorphism were discussed in the cotton aphid. The present study established an efficiently standardized protocol for GP and male induction in cotton aphid by transferring newly-born nymphs to short photoperiod conditions (8 L:16D, 18 °C). The external morphological characters, especially wing vein patterns, were similar among WPFs, GPe, and males. However, their reproductive organs were strikingly different. Compared with WLPF, shared (2 335) and exclusively (1 470 in WLPF, 2 419 in GP, 10 774 male) expressed genes were identified in the three-wing morphs through RNA-Seq, and several signaling pathways that are potentially involved in their wing differentiation were obtained, including insulin signaling, starch and sucrose metabolism.


Background
Wings are essential organs for insects to find hosts and spawning sites, playing a vital role in maintaining their population continuation. Insects develop alternative wing dimorphism, which includes two categories (i.e., long or short wing and winged or wingless morphologies), to adapt to environmental changes (Campbell and Tregidga 2005;Zhang et al. 2019a). The former is commonly observed in planthoppers or crickets, wherein host plant nutrition affects wing developmental plasticity and the insulin signal pathway participates in regulating wing differentiation (Cisper et al. 2000;Xu et al. 2015;Zhao and Zera 2002). By contrast, the latter mainly exists in aphid species, which consist of fully winged (flight-capable) versus fully wingless (flight-incapable) morphs (Brisson 2010;Ogawa and Miura 2014). However, the molecular mechanism of winged morph selection in aphids is still not well understood.
Considerable attempts have been made to uncover the molecular mechanism underlying winged morph switch in Acyrthosiphon pisum, a non-host-alternating aphid species. This species has two types of wing dimorphisms, including environmental-induced wing polyphenism and genetically determined wing polymorphism (Brisson 2010;Ogawa et al. 2012). The former occurs among parthenogenetic females and is regulated by ecdysone or octopamine (Vellichirammal et al. 2017;Wang et al. 2016). By contrast, the latter is only found in males, in which the sex-linked locus on the X-chromosome named aphicarus controls wing differentiation (Braendle et al. 2005;Caillaud et al. 2002). Studies in pea aphid have considerably advanced the understanding of the molecular mechanism underlying wing phenotypic plasticity in non-host-alternating aphid species. However, only 12 species belong to this category and own two types of wing dimorphisms, which account for less than 2% of the 731 recorded aphid species (Blackman and Eastop 2007). Comparison with non-host-alternating aphids, few studies have focused on the wing plasticity in hostalternating species.
Host-alternating aphids are often more economically important than non-host-alternating species because their secondary hosts are usually herbaceous species and frequently a crop plant. Moreover, the males in these species are exclusively winged (Hardie 2017). Gynoparae (GPe), another exclusively winged morph, capable of giving birth to sexual female offspring, are observed in host-alternating aphids (Campbell and Tregidga 2005;Hales et al. 1989;Hardie 1980;Hong and Boo 1998;Kwon and Kim 2017). Three winged morphs exist in these species, namely, winged parthenogenetic females (WPFs), winged GPe, and winged males, which are different from non-host-alternating aphids. WPFs, GPe, and males have been induced by crowding, poor plant nutrition, shortened photoperiod, and reduced temperature in several aphid species (Hardie 1980;Kwon and Kim 2017;Liu et al. 2014). However, the shared molecular mechanism underlying the wing plasticity of these three-wing morphs in host-alternating aphids remains unclear.
The present study focused on the representative of host-alternating aphids, Aphis gossypii Glover, which is also known as the cotton aphid or the melon aphid. This aphid species is among the top 10 agricultural pests; it is distributed in 171 countries and damages various crops, including those in Cucurbitaceae, Malvaceae, Solanaceae, and Rutaceae (Willis 2017). In this species, the WPFs appearing in spring and summer facilitate population expansion among various host plants and different regions, whereas the GPe and males produced in autumn could fly from second hosts to seek winter hosts and promote gene communication of colonies from different secondary hosts or regions (Hales et al. 1989;Kwon and Kim 2017). Thus, the cotton aphid is a compelling laboratory model for studying the molecular mechanism of wing dimorphism phenomenon in host-alternating aphids. Several strategies for inducing WPFs, GPe, and males in the cotton aphid have been explored (Kwon and Kim 2017;Liu et al. 2014;Shi et al. 2010;Takada 1988). However, in these studies, winged morphs were produced qualitatively not quantitatively, and the methods were complicated and time-consuming.
In a previous study, the effects of postnatal crowding on winged morph induction were identified, and several signaling pathways potentially involved in the wing differentiation of WPF were obtained in cotton aphid (Ji et al. 2019). In the present study, GPe and males were induced in the laboratory to identify the shared important genes and signaling pathways probably involved in the wing differentiation of the three wing-morphs in cotton aphids, and comparative transcriptome analysis was performed to identify differentially expressed genes (DEGs) among GPe, males, and wingless parthenogenetic female (WLPF). Moreover, thousands of shared DEGs among these three-winged morphs and WLPF were obtained through a conjoint analysis with DEGs identified in the comparison of WPF versus WLPF (Ji et al. 2019). The putative roles of these shared DEGs and significantly enriched signaling pathways potentially involved in the winged morph switch of the cotton aphid were discussed. Taken together, the results provide valuable resources for the gene expression profiles in the three-wing morphs of A. gossypii. These results can improve our understanding of the molecular mechanisms underlying the wing mode switch of aphid species.

Induction of GPe and males
A. gossypii colony, collected originally in Anyang, Henan, China, was reared on cotton seedlings at the Institute of Cotton Research of China Academy Agricultural of Sciences under controlled laboratory conditions for more than 50 generations before the start of all experiments (25°C ± 1°C, 75% relative humidity, and 16 L:8 D photoperiod). Winged GPe and winged males were induced by rearing newly-born WLPF nymphs (G1) at short-day (SD) conditions (18°C ± 1°C, 75% relative humidity, and 8 L:16 D photoperiod; Additional file 1: Fig. S1). The offspring (G2) produced by G1 aphids were translated to a separate cotton seedling daily, and their morphs were identified in adulthood. The progeny numbers of G1 mothers in SD condition were counted every day.

Morphological characters and fecundity
WLPFs were obtained by rearing aphids in solitarily condition, while WPFs were generated by rearing aphids in crowding condition at a density of 20 nymphs·cm − 2 , as shown in Additional file 1: Fig. S1 (Ji et al. 2019). Then, the morphological characters of the bodies, wing vein patterns, and internal genitals of one-day-old adult WLPFs, WPFs, GPe, and males were visualized using a SteREO Discovery V8 microscope (Zeiss, Germany). In addition, the length of the body, forewing and underwing of each morph was measured, and the offspring numbers of WLPFs, WPFs, GPe were counted. The internal genital characters of GPe and males captured in the field were used as a reference to identify the GPe and males produced indoors (Additional file 2: Fig. S2).
All eight morphs in the annual life cycle of cotton aphid, including egg, fundatrix, fundatrigenia, WLPF, WPF, GP, sexual female, and male, captured in the field were used as criteria for verifying the induction protocol.
Sample preparation for RNA sequencing (RNA-seq) Adult GPe and males were used for RNA sequencing. The adult GPe were gathered individually after they produced all offspring to eliminate the potential influence of embryos embedded in their ovaries, and males were collected 1 day after their last molting. Adult WLPFs and WPFs were collected separately only after their reproductive cycles ended (Ji et al. 2019). Four biological replicates were conducted for each morph. Each biological replicate contained 50 aphids. The samples were immediately frozen in liquid nitrogen and stored at − 80°C.

Transcriptome assembly and gene annotation
Total RNA was isolated from the pooled whole body of aphid samples by using a TRIzol reagent (Promega, USA) according to the manufacturer's instructions. The purity and integrity of RNA were assessed using an Agilent 2100 Bioanalyzer (Agilent, USA). The cDNA libraries were sequenced in paired-end modes on a BGISEQ-500 system at Beijing Genomics Institute (BGI, Shenzhen, China). Sixteen libraries from GPe, males, WLPFs, WPFs were constructed. Raw reads were obtained after sequencing, and adapter sequences and low-quality reads were subsequently filtered to acquire clean reads. Then, they were assembled de novo using Trinity software (Grabherr et al. 2011) in accordance with a previous study (Ji et al. 2019). Sequences larger than 200 bp were aligned to protein databases, including Nr, Swiss-Prot, KEGG, Pfam, and KOG by blastx, and to the nucleotide Nt by blastn, with a cutoff E-value of 10E-5 (Altschul et al. 1990). A custom-made analysis was used for the comparison between WLPFs and WPFs (Ji et al. 2019). Here, the clean reads from previous studies were used for in-depth data mining.

Identification of DEGs
The gene expression levels in each library were calculated in fragments per kilobase per million reads using Bowtie2 and RESM (Langmead and Salzberg 2012;Li and Dewey 2011). Pairwise comparisons were conducted between GPe versus WLPFs, males versus WLPFs, respectively. The fold change (FC, the relative expression level of a gene in one morph to another) and P-value were used to determine the DEGs between two morphs via DEGseq (Wang et al. 2010). P-values in multiple tests were adjusted using the false discovery rate (Reiner et al. 2003). FC > 2 and adjusted P < 0.001 were set as the thresholds to determine significant differences in gene expression. Gene ontology enrichment and significantly enriched KEGG pathways were identified through hypergeometric tests at P < 0.05 by using phyper function in R software.

Validation of RNA-Seq data by RT-qPCR
The same RNA samples in transcriptome sequencing were used to assess the reliability of sequencing and analysis via RT-qPCR. Single-stranded cDNA was reverse-transcribed and synthesized using 1 μg of RNA from each sample with PrimeScript RT reagent Kit with gDNA Eraser (Takara, Japan) in accordance with the manufacturer's recommendations. RT-qPCR was performed on an ABI StepOnePlus Real-Time PCR System (Thermo, USA) using GoTaq qPCR Master Mix (Promega, USA) and initiated at 95°C for 5 min, followed by 40 cycles of 95°C for 5 s and 60°C for 30 s. Melting curve analysis was conducted to verify the specificity of amplification. The relative expression levels were calculated using the 2 -ΔΔC T method (Livak and Schmittgen 2001) and normalized by the housekeeping gene GAPDH (Gao et al. 2017). The primers used for RT-qPCR validation are shown in Additional file 5: Table S1. Statistical analyses were performed using Student's t-test on SPSS 19.0 software. The significant difference was considered at P < 0.05. Values were reported as mean ± SE.

Production procedure of GPe and males
A. gossypii is a typical host-alternating aphid. Its life cycle is shown in Fig. 1, in which the egg, fundatrix, fundatrigenia, GP, sexual female, and male were captured in pepper, and WLPF and WPF were obtained in cotton. As shown in Fig. 1, the cotton aphid was characterized by three representative winged morphs, namely, WPF, GP, and male. The population density was changed to successfully induce WPF (Ji et al. 2019), and low temperature and short photoperiod were employed to produce GPe and males in cotton aphid (Additional file 1: Fig. S1). All three-wing morphs were induced from Fig. 1 Annual life cycle of cotton aphid. Cotton aphid has three representative winged morphs in its life cycle. In spring, the overwintering egg (a) hatches and develops into the fundatrix (b), which gives birth to fundatrigena (c) on primary host plants. The fundatrigena undergoes several generations and then produces wing parthenogenetic female (d), which migrates to second hosts, such as cotton, in late spring. During summer, cotton aphid alternatively chooses morphs between wing or wingless parthenogenetic females (e) over a dozen of generations to adapt to changing conditions, such as population density or host plant quality. In autumn, exclusively winged GP (f) and male (g) are produced successively and migrate to the primary hosts. Sexual female (h) is produced by GP and oviposits overwintering fertilized eggs after mating with the male (i). All these morphs above were captured in wild field. Scale bars, 0.2 mm nymphal WLPF through one or two generations. When newly-born WLPF nymphs (G1) were reared under SD condition, they exclusively produced GPe and males upon reaching adulthood (Fig. 2a). In other words, 100% of their offspring (G2) were winged. In addition, the experiments confirmed that GPe preceded males among the progeny of an individual G1 mother aphid in SD condition. The G1 mother had two reproduction cycles, each lasting for 7 days (Fig. 2a). The progeny sex ratio of the G1 mother was GP biased (Fig. 2b). Approximately 14 GPe were frequently produced at the first half of the reproductive time of the G1 mother, while 9 males were born during the second half period of the G1 aphid (Fig. 2b).

Morphology of the three-wing morphs of cotton aphid
Compared with WLPF, WPF, GP and males demonstrated similar body morphology, including full wings and melanization of the head and thorax (Fig. 3a). However, their reproductive system was obviously discriminated against each other. The numbers of embryos in WPF and GP were fewer than those in WLPF (Fig. 3a). In addition, the embryos in the ovaries of GP were green, while those in WPF and WLPF were yellow (Fig.  3a). By contrast, males had testes and accessory glands in the abdomen (Fig. 3a). As for wings, the shape and vein of the forewing and underwing were similar among the three-wing morphs (Fig. 3a). The length of the forewing and underwing of WPF is 1.60 and 0.92 mm, which is significantly shorter than that of GP (2.07 and 1.23 mm) and male (2.05 and 1.25 mm) (Fig. 3b). As for the body length, no significant differences were observed among WLPF, WPF, and GP, the exception ofexcept the males (Fig. 3c). In addition, WLPF gave birth to the largest number of offspring at 53, followed by WPF at 27 and GP at 8 (Fig. 3d). The internal genital characters of GP and male induced indoors were the same as those captured in the field (Additional file 2: Fig. S2), thereby indicating that the protocol for GP and male induction in this study was reliable. Even though these winged morphs were induced by different environmental factors with different reproductive systems and fecundities, they all came from WLPF (Additional file 1: Fig. S1). Regardless of the intricate wing-inducing stimulations, the winged morphs are probably regulated by several shared genes or signaling pathways in wing plasticity development.

Transcriptome assembly of the four morphs in cotton aphid
The cDNA libraries of the four morphs were sequenced using Illumina sequencing technology, and 984.87 M clean reads were obtained with Q30 of 90.91% and GC content of 42.20%. The clean reads were assembled into 46 501 unigenes with N50 length of 2 876 bp (Table 1). All unigenes matched previously described reads with more than 92.27% coverage, including no less than 31.25% of unique mapped reads (Table 1). Raw reads were submitted to NCBI SRA under the accession number PRJNA544300. Assembled unigenes were individually searched against public databases. In Nr, Nt, Swiss-Prot, KEGG, KOG, Pfam, and GO databases, a total of 33 998, 33 217, 26 857, 27 779, 25 938, 28 069, and 16 554 unigenes were annotated, respectively (Additional file 3: Fig. S3). Taken together, a total of 37 090 (79.76%) unigenes returned a substantial result in at least one of the searched databases by using this approach (Table 1).

DEGs in winged and wingless aphids
Based on the cutoff criteria (FC > 2, P < 0.001), 10 867 DEGs, including 7 270 upregulated and 3 597 downregulated, were identified in GP compared with those in WLPF (Fig. 4a). When the male was compared with WLPF, 11 137 upregulated and 8 197 downregulated DEGs were identified (Fig. 4a). For WPF versus WLPF, 5 067 DEGs, including 3 187 upregulated and 1 880 downregulated, were identified in a previous analysis (Ji et al. 2019) and showed here as well. The DEGs mentioned above were analyzed using a Venn diagram to identify the shared or separate genes potentially involved in the wing differentiation of the three-winged morphs (Bardou et al. 2014). A total of 769, 1 951, and 5 822 genes were upregulated exclusively in WPF, GP, and male compared with those in WLPF, respectively (Fig.  4b). Correspondingly, 701, 468 and 4 922 genes were downregulated exclusively in WPF, GP, and male compared with WLPF, respectively (Fig. 4b). The winged morphs also demonstrated 2 335 shared DEGs, including 1 658 upregulated and 677 downregulated, compared with WLPF (Fig. 4b). These shared DEGs probably contribute to wing differentiation in all winged morphs.
KEGG enrichment analysis of shared DEGs in the threewing morphs KEGG enrichment analysis was conducted separately on the shared upregulated and downregulated DEGs above to explore the pathways potentially involved in wing differentiation regulation in all three-wing morphs in the cotton aphid. Compared with WLPF, all the winged morphs showed several upregulated pathways including Fig. 3 Bodies, internal genitals, and fecundity of four morphs in cotton aphid. a Body, internal genital, forewing and underwing morphologies of WLPF, WPF, GP, and male obtained in laboratory condition. ov, ovariole; ag, accessory gland; t, testes. Scale bars, 0.2 mm. b Length of forewing and underwing of WPF, GP, and male. c Body length of WLPF, WPF, GP, and male with the exception of antennae and wings. d Fecundity of WLPF, WPF, and GP. The offspring of individual WLPF, WPF, and GP were counted. Asterisk indicates significant differences between two morphs (Student's t test, **, P < 0.01). Different letters indicate significant differences (P < 0.05) obtained using Tukey's multiple range tests. The length of wings and body of more than 25 individual adults for each morph were measured using a SteREO Discovery V8 microscope. More than 30 individual adults for each morph were observed daily for counting their offspring until death. The values in all panels represent mean ± SE. WLPF, wingless parthenogenetic female; WPF, winged parthenogenetic female; GP, gynopara insulin signaling, estrogen signaling, PPAR signaling, glycolysis/gluconeogenesis, pyruvate metabolism, starch and sucrose metabolism, fructose and mannose metabolism, glycerophospholipid metabolism, steroid biosynthesis, MAPK signaling, AMPK signaling, FoxO signaling, etc. (Fig. 5). These upregulated pathways were clustered into categories of the endocrine system, carbohydrate metabolism, lipid metabolism, and signal transduction (Fig. 5). By contrast, compared with WLPF, all winged morphs exhibited several downregulated pathways related to transcription and nucleotide metabolism, including RNA polymerase, pyrimidine metabolism, and purine metabolism (Fig.  5). The significantly enriched KEGG pathways on the shared DEGs are listed in Additional file 6: Table S2. Even though the upregulated and downregulated signaling pathways potentially involved in wing mode plasticity in the cotton aphid were identified, further studies are needed to determine their mechanism in wing dimorphism.

Validation of RNA-Seq data
Thirty-two upregulated and 12 downregulated genes that spanned the range of differential expression and were enriched in the identified pathway above were selected to validate the statistical analysis with RNA-Seq via RT-qPCR and to confirm the results of gene expression profiling. The FCs of all selected genes (either upregulated or downregulated) in RT-qPCR were consistent with the results in RNA-Seq (Additional file 7: Table S3).
Analysis of the Pearson correlation coefficient indicated a significant positive correlation (P < 0.000 1) between the results of RT-qPCR and RNA-Seq in all three-wing morphs compared with WLPF (Fig. 6). Therefore, the RNA-Seq data in this study were reliable.

DEGs related to the insulin signaling pathway
The insulin signaling pathway participates in wing differentiation and development in several insects (Ding et al. 2017;Xu et al. 2015). In the present study, it was highly upregulated in all three-wing morphs compared with WLPF ( Fig. 5), implying the importance of insulin signaling in wing dimorphism in A. gossypii. Fifteen insulin signaling-related genes, which were upregulated in all three-wing aphids in the RNA-Seq data, were analyzed using RT-qPCR technology. These genes were significantly highly expressed in all three-wing morphs, with FC ranges of 1.37∼17.46 in WPF versus WLPF, 1.96∼9.02 in GP versus WLPF, and 3.51∼87.20 in male versus WLPF (Fig. 7). Next, the candidate genes potentially involved in wing differentiation will be silenced using RNA interference (RNAi) to identify their function in the regulation of wing plasticity in all three-wing morphs in the cotton aphid.

Induction methods of GPe and males in host-alternating aphids
The distinction of wing dimorphism between hostalternating and non-host-alternating aphid species probably results from their diverse survival strategies. The   former always alternates between primary (usually woody) host plants and secondary (herbaceous) hosts, while the latter are usually monophagous (Hardie 2017). Thus, in host-alternating species, such as A. gossypii, the GPe and males are exclusively winged to migrate between the primary and secondary hosts (Hardie 2017).
To obtain an outline of all alate morphs in the cotton aphid, WPF, GP, male, and all apterous morphs were captured in the field from 2017 to 2018. A highresolution image of the life cycle of cotton aphid is shown in Fig. 1, and this image is beneficial to understand the wing plasticity in A. gossypii. Several attempts have been made to induce GP and male indoors and elucidate the mechanism underlying wing dimorphism. When fourth instar apterous nymphs, wingless adults, or alate virginopara adults (G0, initial generation) were transferred from LD condition to SD condition, GPe and males were produced qualitatively in their second generation (G2) (Campbell and Tregidga 2005;Hales et al. 1989;Hardie 1980;Hong and Boo 1998;Kwon and Kim 2017;Liu et al. 2014;Takada 1988), costing two generations (Additional file 8: Table  S4). However, the harvesting protocol in these studies above was unclear, and the corresponding accurate number of the two wing morphs was lacking and uncontrollable. In the study, an optimized protocol was established by rearing G1 (first generation) in SD condition (Additional file 1: Fig. S1), and they produced GPe and males in adulthood quantitatively (Fig. 2). Compared with the results of those the studies above, one generation was enough to produce the two-wing morphs. Moreover, the average number of GPe and males produced by a WLPF in SD condition can be precisely calculated and controlled (Fig. 2).
The morphological characters of the bodies and internal genitals of the two-wing morphs were described in detail, followed by a comparison between WPF and WLPF in the cotton aphid (Fig. 3). The internal genitals of GPe and males induced indoor agreed with those captured in the field (Additional file 2: Fig. S2). The wing vein patterns of the three-wing morphs were first compared with each other in the cotton aphid, and no significant difference in morphology was observed except for the length (Fig. 3). Taken all together, the effective induction method of GPe and males in this study is reliable and could facilitate the research of wing phenotypic plasticity in host-alternating aphid species.

Transcriptome comparison and DEG analysis
Transcriptomes were assembled using WLPF, WPF, GP, and male in the cotton aphid, and a total of 46 501 unigenes were obtained (Table 1), which were approximated to 44 310 in a previous study . After the repetitive genes and genes without annotation were filtered out, 37 090 unigenes were yielded in the present study (Table 1), while only 11 350 unigenes were obtained in Liu et al's study. This dissimilarity possibly resulted from the difference in winged morph selection strategies. WLPF, GP, and sexual female were sequenced in Liu et al's study, while WLPF, WPF, GP, and male were collected in the present study. Several previously unreported gene transcripts or isoforms probably existed in males and WPF in the cotton aphid. In addition, compared with WLPF, GP had 7 270 upregulated and 3 597 Fig. 7 Upregulated DEGs associated with insulin signaling pathway in the winged morphs. A total of 15 upregulated genes related to the insulin signaling pathway were observed with similar expression profiles between RT-qPCR and RNA-Seq results in the pairwise comparison between the winged and wingless morphs. Each value is the mean of four replicates, and error bars indicate SEs. InR1, insulin receptor 1; PEPCK, phosphoenolpyruvate carboxykinase, cytosolic [GTP]-like; glgP, glycogen phosphorylase-like; FASN, fatty acid synthase; GYS1, glycogen synthase; Phkb, phosphorylase b kinase regulatory subunit beta; Phkg1, phosphorylase b kinase gamma catalytic chain; HK2, hexokinase-2; PP1R3C, protein phosphatase 1 regulatory subunit 3C; Pten, phosphatidylinositol 3,4,5-trisphosphate 3-phosphatase and dual-specificity protein phosphatase; CBL, calcineurin B-like; Smad4, mothers against decapentaplegic homolog 4; Mdm2, E3 ubiquitin-protein ligase Mdm2; and FBXO32, F-box only protein 32. Fold change, the relative expression of a gene in one morph to another downregulated DEGs, while only 741 upregulated and 879 downregulated genes were demonstrated by a previous study . The difference could be attributed to the various strategies adopted to eliminate the potential influence of embryos in the mother's ovaries. In particular, embryos were manually removed from GPe in Liu et al's study, while all offspring were born before adult GPe were collected in the present study (Additional file 4: Fig. S4). Compared with other studies on the cotton aphid wing differentiation, the present study identified the shared and exclusively DEGs in WPF, GP, and males, and compared them with those in WLPFs, respectively (Fig. 4).

Pathways potentially involved in wing differentiation of the three-wing morphs
Several transcriptional studies have focused on aphids to identify signal pathways involved in wing differentiation. WPF exhibited 1 663 DEGs compared with WLPF in the cotton aphid, and these genes were significantly enriched in the ribosome, pyruvate metabolism, proteasome, lipid metabolism, protein synthesis and degradation, RNA transport, and antigen processing and presentation (Yang et al. 2014). In the present study, pyruvate metabolism, and antigen processing and presentation were also upregulated in all three-wing morphs compared with WLPF (Additional file 6: Table S2), indicating that these two pathways probably contribute to wing differentiation not only in WPF but also in GP and male. GP had a total of 1 620 DEGs compared with WLPF, in which six upregulated and seven downregulated signaling pathways were enriched, including starch and sucrose metabolism, phototransduction, dorsoventral axis formation, Wnt, and Notch . In this study, the pathway of starch and sucrose metabolism were also upregulated in the three-wing morphs compared with WLPF (Additional file 6: Table S2), which may be due to the indispensability of energy for flight apparatus and flight behavior.
The above studies have advanced the understanding of wing dimorphism in the cotton aphid. However, Liu et al's study mostly highlighted reproductive polyphenism, while Yang et al's study only focused on wing plasticity in WPF Yang et al. 2014). In the present study, three-winged morphs in the cotton aphid were compared with each other or to WLPF. A total of 2 335 shared DEGs, including 1 658 upregulated and 677 downregulated, were identified in all three-wing morphs compared with WLPF (Fig. 4), and they were significantly enriched in 49 upregulated and 7 downregulated KEGG pathways (Additional file 6: Table S2). The upregulated pathways were clustered into categories of signal transduction, lipid metabolism, carbohydrate metabolism, and endocrine system (Fig. 5). Energy allocation is important for the trade-off between winged morph and wingless morph in aphids (Yang et al. 2014). Thus, the upregulated lipid metabolism and carbohydrate metabolism in the three-winged morphs was consistent with the previous observation of significantly increased triglyceride content in the winged morph versus the wingless morph (Shi et al. 2010). Compared with WLPFs, all three alate morphs had upregulated insulin signaling pathway (Fig. 5), which was proven to regulate wing differentiation and development in several insects, including Nilaparvata lugens, Laodelphax striatellus, Sogatella furcifera, and Blattella germanica (Abrisqueta et al. 2014;Xu et al. 2015;Xu and Zhang 2017). Considering the importance of insulin signaling in wing determination, the relative expression levels of 15 related DEGs were validated in all three wing morphs and compared with those in WLPF (Fig. 7). The results showed that insulin signaling was potentially involved in the wing differentiation of the three wing morphs in cotton aphid.

DEGs associated with insulin, flight muscle, and energy
Insulin receptor 1 (InR1) leads to long-winged morph if active and short-winged morph if inactive in three planthoppers (Xu et al. 2015). Moreover, the silencing of InR1 disrupts the nymph-adult transition of alate viviparous females in A. (Toxoptera) citricidus (Ding et al. 2017). WPF, GP, and male showed 6.03-, 2.23-, and 6.70-fold increases in InR1 compared with WLPF, respectively ( Fig. 7; Additional file 7: Table S3). This finding demonstrated the potential importance of InR1 in wing regulation of the three wing morphs in the cotton aphid. Phosphoenolpyruvate carboxykinase [GTP]-like (PEPCK) and glycogen phosphorylase-like (glgP) were highly expressed in alate A. citricida adults, and individual silencing of these two genes could result in underdeveloped wings in WPFs at rates of 58%∼79% (Shang et al. 2016). Likewise, these two genes were significantly upregulated in all three wing morphs compared with WLPF in cotton aphid (Fig. 7). The expression of phosphatidylinositol 3,4,5-trisphosphate 3-phosphatase and dual-specificity protein phosphatase (Pten), the gene controlling junction lengthening and stability during cell rearrangement in epithelial tissue and contributing to cell packing in the wing development (Bardet et al. 2013), was upregulated the expression by 4.54-, 2.17-, and 8.55-fold in WPF, GP, and male compared with WLPF, respectively (Fig. 7), suggesting its importance in wing plasticity of the cotton aphid.
Flightin, a phosphorylated myofibrillar protein essential for thick filament assembly and sarcomere stability in insect flight muscles (Reedy et al. 2000), exhibits increased transcript accumulation in the winged parthenogenetic morphs and males in A. pisum (Brisson et al. 2007) and in WPFs of A. gossypii (Yang et al. 2014) and Rhopalosiphum padi (Zhang et al. 2019b). The disproportionately high levels of flightin transcript in the winged versus the wingless morphs likely result from the presence or absence of indirect flight muscles, which play a causal role in the morphological divergence of the winged morphs (Brisson et al. 2007). In the present study, flightin was increased by 2.32-, 6.50-, and 11.75fold in WPFs, GPe, and males compared with WLPF, respectively (Additional file 7: Table S3). This finding implied the potential importance of flightin in flight muscle formation or wing differentiation in the three alate morphs in the cotton aphid. Odorant receptor coreceptor (Orco), which mediates the winged morph differentiation of parthenogenetic female in Sitobion avenae (Fan et al. 2015), was also highly expressed in all three alate morphs compared with WLPF in the cotton aphid (Additional file 7: Table S3).
Taken together, these shared upregulated genes in the three-wing morphs underlined the importance of the signaling pathways of insulin, energy generation, Orco, and flightin in wing differentiation in the cotton aphid. The functions of these genes in wing dimorphism in A. gossypii will be confirmed using RNAi in a future study. To the knowledge of the authors, this study was the first to examine the transcriptome-wide patterns of differential transcript accumulation associated with threewinged morphs in the cotton aphid. This study could provide a baseline for future studies on the molecular basis of wing differentiation in A. gossypii.

Conclusion
This study provides a high-resolution image of the life cycle of cotton aphid, with clearly visualized body outlines, including all morphs captured in the field. GPe and males were successfully induced indoors, and the morphological characters of the bodies and reproductive organs of these two-wing morphs were compared with those of WPF and WLPF for the first time. A total of 2 335 shared DEGs were identified in all three-wing morphs via comparative transcriptomic analysis. The signaling pathways potentially involved in wing differentiation in these-winged morphs were obtained including insulin signaling, FoxO signaling, fatty acid biosynthesis and degradation and so on. The expression levels of DEGs associated with insulin signaling, flight muscle formation, and energy generation were validated via RT-qPCR. These findings could improve our understanding of wing dimorphism in the cotton aphid. Further studies on the functions of candidate genes via RNAi could provide a basis for developing genetic control strategies against this pest through the disruption of its migratory behavior.