Insights into wing dimorphism in the worldwide agricultural pest Aphis gossypii , the host-alternating aphid CURRENT

Background: Three wing morphs exists in the life cycle of the worldwide pest Aphis gossypii, i.e., wing parthenogenetic female (WPF), gynopara (GP) and male, which were produced mostly by crowding and host quality, photoperiod, loss of X chromosome, respectively. However, the shared molecular mechanism underlying their wing differentiation remains an enigma. Here we firstly induced gynoparae and males indoors and compared the characters of these wing morphs in body, internal genitals and fecundity. Then we identified the shared and separate differentially expressed genes (DEGs) and signaling pathways potentially involved in the wing morphs regulation in WPF, GP and male compared to wingless parthenogenetic female (WLPF). Results: Newly-born nymphs reared in short photoperiod condition exclusively produce gynoparae and males in adulthood successively, in which the sex ratio is gynoparae biased. Compared with WLPF, three wing morphs have similar morphology in bodies but is obviously discriminated in the reproductive system and fecundity. Built upon our previous study, 37 090 annotated unigenes were obtained from libraries constructed by the four morphs above through RNA-sequencing, in which 10 867, 19 334 DEGs were identified in pairwise comparison of GP vs. WLPF, Male vs. WLPF, respectively. Furthermore, 2 335 shared DEGs including 1 658 up- and 677 downregulated were obtained in these wing morphs compared to WLPF. The 1 658 shared up-regulated DEGs were enriched in multiple signaling pathways including insulin, FoxO, MAPK, strarch and sucrose metabolism, fatty acid biosynthesis and degradation which hint their key roles in the regulation of wing plasticity in cotton aphid. Gene expression levels were validated by using Pearson’s correlation (r) and potential roles of 15 DEGs related to the insulin signaling pathway in cotton wing dimorphism were discussed. Conclusions: The results of this study establish a solid foundation for deciphering molecular mechanisms underlying the switch between wingless and wing morphs in the cotton aphid and provide valuable resources for future research on the host-alternating aphid species. the of molecular mechanisms wing thick filament and sarcomere stability in Drosophila flight Cell Biol

alternative wing dimorphism, including two categories, long or short wing, winged or wingless morphologies [1,2]. The former is commonly observed in plant hoppers or crickets, in which host plant nutrition affects the wing developmental plasticity and insulin signal pathway participates in the regulation of wing differentiation [3][4][5]. In contrast, the later mainly exist in aphid species, which consists of fully winged and flight capable versus fully wingless and flight-incapable morphs [6,7].
However, the molecular mechanism of wing morphs selection in aphids is still not well understood.
Considerable attempts have been made to uncover the molecular mechanism underlying wing morphs switch in Acyrthosiphon pisum, a non-host-alternating aphid species. Two types of wing dimorphisms exist in this species, including environmental wing polyphenism and genetic wing polymorphism [6,8]. Environmental induced wing polyphenism occurs among parthenogenetic females and was regulated by ecdysone or octopamine [9,10]. By contrast, genetically determined wing polymorphism has been found only in males, in which the sex-linked locus on the X-chromosome named aphicarus controls the wing differentiation [11,12]. These excellent studies in pea aphid have greatly advanced our understanding of the molecular mechanism underlying wing phenotypic plasticity in the non-host-alternating aphid species. However, not all aphid species own these two types of wing dimorphisms, especially dimorphic males. For example, only about 10% of European species have been recorded with winged and wingless males [13]. Furthermore, hostalternating/heteroecious aphids are often more economically important than no-host-alternating species, as the secondary hosts usually herbaceous species, and frequently a crop plant [14]. More distinctive is, males in all host-alternating aphids are exclusively winged [14] .
In addition to the fully winged male, another exclusively winged morph gynopara, capable of giving birth to sexual female offspring, exists in host-alternating aphid species [2,[15][16][17][18]. In other words, three wing morphs exists in these species, namely wing parthenogenetic female (WPF), gynopara (GP) and male, which is extremely different from the non-host-alternating aphid species. Several shared or separate mechanism probably participate in wing differentiation of these three wing morphs. To date, it is well known that WPF, GP and male was induced by crowding, poor plant nutrition, shortened photoperiod, lower temperature [18][19][20]. However, less well known is the molecular mechanism underlying wing differentiation regulation of these three wing morphs in host-alternating aphids.
The representative of host-alternating aphids, Aphis gossypii Glover, also known as cotton aphid or melon aphid, caught our attention. This aphid species is the top 10 agricultural pests which distributes in 171 countries and damage various crops including those in Cucurbitaceae, Malvaceae, Solanaceae, and Rutaceae [21]. In this species, WPFs appearing in spring and summer facilitate population expansion among various host plants and different regions, while gynoparae and males produced in autumn can fly from second hosts to seek winter hosts and promote gene communication of colonies from different secondary hosts or regions [17,18]. As an economically important destructive pest in agriculture, horticulture and forestry, several induction strategies for WPF, gynopara and male in this species have been obtained so far [18,20,22,23]. In a word, cotton aphid is a compelling laboratory model to study the molecular mechanism of wing dimorphism phenomenon in host-alternating aphids.
Recently, transcriptome sequencing has offered powerful tools for genes expression profiles analysis, differentially expressed genes (DEGs) selection, function genes mining [24]. In our previous study, the effects of postnatal crowding on wing morphs induction was identified and several signaling pathways potentially involved the wing differentiation of WPF were obtained in cotton aphid [25]. Here, in order to identify the shared important genes and signaling pathways probably involved in wing differentiation of the three wing morphs in cotton aphids, we induced gynoparae and males in the laboratory and employed a comparative transcriptome analysis approach to identify differentially expressed genes (DEGs) in gynoapra, male and wingless parthenogenetic female. Moreover, through a conjoint analysis with DEGs identified in the comparison of WPF vs. WLPF [25], thousands of shared DEGs between these three wing morphs and WLPF were obtained. On the basis, the putative roles of these shared DEGs and significantly enriched signal pathways potentially involved in the wing morphs switch of cotton aphid were discussed. Taken together, our results provide the most valuable resources to date for gene expression profiles in three wing morphs of the host-alternating aphid, A.
gossypii. Our finding also could benefit the understanding of molecular mechanisms in wing mode switch of aphid species.

Results
The production of gynoparae and males The annual life cycle of cotton aphid is characterized by three representative wing morphs including wing parthenogenetic female, gynopara and male (Fig. 1). Wing morphs of aphids were mostly induced by density, host plant quality, temperature and photoperiod. Here we used low temperature and short day as wing-inducing stimulation to produce gynoaprae and males in cotton aphid. When newly-born nymphs (G1) were reared under the SD conditon (18 ± 1 ℃, 75% relative humidity, 8L:16D photoperiod), they exclusively produced two wing morphs, gynoparae and males, in adulthood ( Fig. 2A; Fig. 3A). In other words, 100% of their offspring (G2) were winged. Besides, our experiments confirm that gynoparae precede males among the progeny of an individual G1 mother aphid ( Fig. 2A) and show that the sex ratio in the offspring is gynoparae biased (Fig. 2B). Specifically, about 14 gynoparae were frequently produced at the first half of the reproductive time of G1 mother, while 9 males were born subsequently during the second half period of G1 aphid (Fig. 2B).

Morphology Of Three Wing Morphs Of Cotton Aphid
Compared with WLPF, the WPF, GP and male have similar morphology in bodies, including fully wings, heavier melanization of head and thorax (Fig. 3A). However, the reproductive system of these three wing morphs is obviously discriminated. In comparison to WLPF, number of embryos in WPF and gynopara were both fewer (Fig. 3A). Besides, embryos in the ovaries of gynopara are green while those in WPF and WLPF are yellow (Fig. 3A). In contrast, testes and accessory gland exist in the abdomen of male (Fig. 3A). There are no significant differences between WLPF, WPF and gynopara in body length with the exception of male (Fig. 3B). In addition, WLPF give birth to most offspring with the number of 53, followed by WPF with the number of 27 and gynopara with the number of 8 ( Fig. 3C). Internal genitals characters of gynopara and male induced in controlled indoor were the same as those captured in the wild field ( Fig. S1: A, gynoapra; B, male). Anyway, though these wing morphs were induced by different environmental factors with different reproductive system and fecundity, they probably were regulated by several shared genes or signaling pathways in wing differentiation process. 069, 16 554 unigenes were annotated in above database, respectively (Fig. S1). Taken together, a total of 37 090 (79.76%) unigenes returned a significant result at least in one of the searched databases using this approach (Table 1).

DEGs In Winged And Wingless Morphs
Based on the cutoff criteria (FC > 2, P < 0.001), 10 867 DEGs including 7 270 upregulated and 3 597 downregulated genes were identified in GP compared to WLPF (Fig. 4A). When male was compared to WLPF, 11 137 upregulated and 8 197 downregulated DEGs were identified (Fig. 4A). For WPF vs. WLPF, 5 067 DEGs including 3 187 upregulated and 1 880 downregulated genes were identified in our recently study [25], which was showed here as well. In order to identify the shared or separate genes potentially involved in wing differentiation of the three wing morphs in cotton aphids, DEGs above were analyzed by Venn diagram. Interestingly, 769, 1951, 5 822 genes were up-regulated exclusively in WPF, GP and male compared to WLPF, respectively (Fig. 4B). 701, 468, 4 922 genes were downregulated exclusively in WPF, GP and male compared to WLPF, respectively (Fig. 4B). Above all, 2 335 shared DEGs including 1 658 up-and 677 downregulated were identified in all these three wing morphs compared to WLPF (Fig. 4B). These shared DEGs probably contribute to the wing dimorphism regulation in cotton aphid.

KEGG enrichment analysis of shared DGEs in three wing morphs
To explore the pathways potentially involved in the wing differentiation regulation in all three wing morphs in cotton aphid, the shared up-and downregulated DEGs above were performed KEGG enrichment analysis separately. Compared to WLPF, several pathways were up-regulated in all wing morphs, including insulin, FoxO, fatty acid biosynthesis and degradation, metabolism of glycerophospholipid, fructose and mannose, starch and sucrose, steroid biosynthesis and so on (Fig. 5). These up-regulated pathways were clustered into categories of signal transduction, lipid metabolism, carbohydrate metabolism, endocrine system (Fig. 5). By contrast, several pathways related to nucleotide metabolism and transcription were downregulated in these wing morphs compared to WLPF, including RNA polymerase, purine and pyrimidine metabolism (Fig. 5). Though the up-and downregulated signaling pathways potentially involved in wing mode plasticity in cotton aphid were identified, how they work in wing dimorphism is needed to be elucidated in further studies yet.
All significantly enriched KEGG pathways performed on shared DEGs were listed in Table S2.

RT-qPCR Validation
To confirm the results of gene expression profiling, we chose 32 up-and 12 downregulated genes that spanned the range of differential expression to validate our statistical analysis with RNA-Seq by RT-qPCR method (Table S3). Fold changes of all these selected genes (either up-or down-regulated) in RT-qPCR were consistent with the results in RNA-seq (Table S3). Pearson correlation coefficient analysis by SPSS 19.0 software indicated a significant positive correlation (P < 0.0001) between the results of RT-qPCR and RNA-seq in all these three wing morphs compared to WLPF (Fig. 6). Therefore, the RNA-Seq data in our study were reliable.

DEGs Related To Insulin Signaling Pathway
Insulin signaling pathway participates in wing differentiation and development in several insects [4,26]. Interestingly, it was highly expressed as well in all these three wing morphs compared to WLPF in cotton aphid (Fig. 5), which implies the importance of insulin signaling in wing dimorphism in A.
gossypii. Thus, 15 insulin signaling related genes which were upregulated in all these three wing aphids in the RNA-seq data were analyzed through using the method of RT-qPCR. These genes were significantly high expressed in all these three wing morphs with the fold change ranges of 1.37 ~ 17.46 in WPF vs. WLPF, 1.96 ~ 9.02 in GP vs. WLPF, 3.51 ~ 87.20 in male vs. WLPF (Fig. 7). Next, silencing of these candidate genes potentially involved in wing differentiation by using RNA interference (RNAi) will be carried out to identify their function in the regulation of wing plasticity in all these three wing morphs in cotton aphid.

Discussion
Exclusively winged gynopara and male in A. gossypii Insects produce alternative wing morphs in response to environment changing [1], in which winged or wingless morphs were mostly observed in aphids species. Considerable studies on wing phenotypic plasticity mostly have been carried out on pea aphid, the non-host-alternating species, in which two wing morphs environmentally induced or genetically determined, namely WLPF and male, existed [8,27]. However, three wing morphs existed in cotton aphid, the host-alternating species, in which males were exclusively winged and an extra special wing morph, gynopara, were reported in this species [18]. Distinction of wing dimorphism between host-alternating and no-host-alternating aphid species probably result from their diverse survival strategies. The former always alternates between primary (usually woody) host plants and secondary (herbaceous) hosts, while the latter were usually monophagous [28]. Hence, in host-alternating species A. gossypii, gynoparae and males are exclusively winged to migrate between the primary and secondary hosts. By contrast, not all males were winged in Acyrthosiphon pisum, possibly because this species was not entirely necessary to migrate in winter [8]. Certainly, winged parthenogenetic females were produced for population expansion in both host-alternating and non-host-alternating species in summer [28].
Besides, male wing differentiation in pea aphid were determined by a single locus on the X chromosome called aphicarus [8]. However, molecular mechanism about wing differentiation of gynopara and male in cotton aphid were still not well known. Species exhibiting three alternative wing morphs provide valuable models for studying wing differentiation in aphids and some shared physiological or molecular baselines probably exist. Here, methods for gynopara and male induction in laboratory were established by shortening photoperiod in cotton aphid (Fig. 2). Compare to previous studies [18,29], one generation induction was enough to produce these two wing morphs with a high effectiveness. Furthermore, morphological characters in body and internal genitals were firstly described and discriminated in three wing morphs (Fig. 3). All these results broaden our understanding about wing phenotypic plasticity in aphid species.

Transcriptomes Comparison And DEGs Analysis
We firstly assembled transcriptomes by WLPF, WPF, gynopara and male in cotton aphid and obtained a total of 46 501 unigenes (Table 1), which was approximated to the number of 44 310 in previous study [20]. After filtering out the repetitive genes and genes without annotation, 37 090 unigenes yielded in our study (Table 1), while only 11 350 unigenes were obtained in Liu's study. This may result from the difference in wing morphs selection strategies. Only WLPF, gynopara and sexual female were sequenced in Liu's study, while extra male and WPF were collected in our study meanwhile. Several previously unreported gene transcripts or isoforms probably existed in male and WPF in cotton aphid. In addition, compared to WLPF, gynopara owned 7 270 up-regulated and 3 597 down-regulated DEGs while there were only 741 up-and 879 down-regulated genes in previous study [20]. Different strategies were adopted to eliminate the potential influence of embryos in the mother's ovaries: embryos were manually removed in Liu's study while all offspring were born before adult cotton aphid were collected in our study. Compared to these excellent studies performed on cotton aphid wing differentiation, we firstly identified the shared and exclusively DEGs in WPFs, GP, male compared to WLPFs, respectively (Fig. 4).
Pathways potentially involved in wing differentiation of three wing morphs Several transcriptional studies have focused on aphids to seek signal pathways involved in the wing differentiation. 1 663 DEGs were identified in WPF compared to WLPF in cotton aphid, which were significantly enriched in ribosome, pyruvate metabolism, proteasome, lipid metabolism, protein synthesis and degradation, RNA transport, antigen processing and presentation [29]. Our results were consistent with those finding, in which pyruvate metabolism, antigen processing and presentation were up-regulated in all these three wing morphs compared to WLPF as well (Table S2). 1 620 DEGs were identified in gynoparae compared to WLPF, in which 6 up-and 7 downregulated signaling pathways were enriched, including starch and sucrose metabolism, phototransduction, dorso-ventral axis formation, Wnt, Notch [20]. Likewise, pathway of starch and sucrose metabolism were upregulated in these three wing morphs (Table S2). This probably results from the indispensability of energy for flight apparatus and flight behavior.
Those studies advanced our understanding of wing dimorphic in cotton aphid. However, Liu's study mostly was highlighted on reproductive polyphenism while Yang's study only performed on WPF [20,29]. In this study, three wing morphs in cotton aphid were firstly compared with each other or to wingless parthenogenetic female. 2 335 shared DEGs including 1 658 up-and 677 downregulated were identified in all three wing morphs compared to WLPF (Fig. 4), which were significantly enriched in 49 up-and 7 downregulated KEGG pathways (Table S2). The upregulated pathways were clustered into categories of signal transduction, lipid metabolism, carbohydrate metabolism, endocrine system ( Fig. 5). Energy allocation is important for trade-off between wing morph and wingless morph in aphids [29]. Thus, the upregulated lipid and carbohydrate metabolism in three wing morphs is consistent with the previous observation of the significantly higher triglyceride content in winged morph versus the wingless morph [23]. Besides, compared to WLPFs, all these three alate morphs had up-regulated insulin signaling pathway (Fig. 5), which had been proved regulating wing differentiation and development in several insects including Nilaparvata lugens, Laodelphax striatellus, Sogatella furcifera, Blattella germanica [4,30,31]. Considering the importance of insulin in wing determination, relative expression levels of 15 related DEGs in RNA-seq data were validated in all three wing morphs compared to WLPF (Fig. 7), which implied insulin is potentially involved in the wing differentiation of 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 active in three planthoppers [4]. Besides, silencing of InR1 disrupts nymph-adult transition of alate viviparous females in A. (Toxoptera) citricidus [26]. InR1 were increased significantly to 6.03, 2.23, 6.70 folds in WPF, gynopara, male compared to WLPF, respectively (Fig. 7, Table S3). This hints the potential which plays a causal role in the morphological divergence of the wing morphs [33]. Particularly, flightin were increased to 2.32, 6.50, 11.75 folds in WPF, gynopara, male compared to WLPF in our results likewise, respectively (Table S3). This imply the potential importance of flightin in flight muscle formation or wing differentiation in three wing morphs in cotton aphid. Besides, phosphoenolpyruvate carboxykinase [GTP]-like (PEPCK) and glycogen phosphorylase-like (glgP) were expressed highly in alate A. citricida adults, and silencing of these two genes individually can all resulted in underdeveloped wings in WPFs at the rates of 58-79% [35]. Likewise, these two genes were significantly upregulated in all three wing morphs compared to WLPF in cotton aphid (Fig. 7). In addition, odorant receptor co-receptor (Orco) which mediates winged morph differentiation of parthenogenetic female in Sitobion avenae [36], was also higher expressed in all three alate morphs compared to WLPF in cotton aphid (Table. S3).
Taken together, these shared upregulated genes in three wing morphs underline the importance of signaling pathways of insulin, energy generation, orco, flightin in wing differentiation in cotton aphid.
The functions of these genes in wing dimorphism in A. gossypii will be confirmed by RNAi strategy in our future study. To our knowledge, this is the first study examining transcriptome-wide patterns of differential transcript accumulation associated with three wing morphs in cotton aphid, which will provide a baseline for future studies on molecular basis of wing differentiation in A. gossypii.

Conclusion
A. gossypii is the top 10 worldwide destructive sap-sucking pest which distributes in 171 countries with various hosts and vector of more than 50 plant viruses [21]. The ecological success of this pest was mostly due to three wing morphs in the life cycle, namely WPF, gynopara and male. In this study, we induced gynoparae and male successfully indoors and compared the morphological characters of these three wing morphs in body and reproduction organ for the first time. Then comparative transcriptomic analysis were employed and the shared 2 335 DEGs were identified in all these three wing morphs compared to parthenogenetic wingless female. At last, potentially signaling pathways involved in the regulation of wing differentiation in all these wing morphs of cotton aphid were achieved by KEGG enrichment analysis. Besides, the annual life cycle of cotton aphid was described and morphological characters of different morphs captured in the wild field were visualized clearly.
Hence discoveries above in this study will advance our knowledge about wing dimorphism in cotton aphid and provide basis for developing genetic control strategies against this pest through the disruption of its migratory behavior.

Induction of gynoparae and males
A. gossypii colony, collected originally in Anyang, Henan, China, was reared on cotton seedlings at Institute of Cotton Research of China Academy Agricultural Sciences under controlled laboratory conditions for more than 50 generations before the start of all experiments (25 ± 1 ℃, 75% relative humidity, 16L:8D photoperiod). Gynoparae and males were induced by rearing newly-born viviparous nymphs (the first generation, G1) at short day (SD) conditions (18 ± 1 ℃, 75% relative humidity, 8L:16D photoperiod). Offspring (the second generation, G2) produced by G1 aphids were translated to a separate cotton seedling daily and their morphs were identified in adulthood. Progeny numbers of G1 mothers were counted every day.

Morphological Characters And Fecundity
WLPFs were obtained by rearing aphids solitarily, while WPFs were generated by rearing aphids in crowding condition [25]. Then, morphological characters of the WLPF, WPF, GP and male in body and internal genitals at their adult stages were visualized using a SteREO Discovery V8 microscope (Zeiss, Germany). In addition, the body length (from head to abdomen, antennae and wings were not included) of each morph was measured and the offspring numbers of WLPF, WPF, GP were counted. Furthermore, internal genitals characters of gynopara and male captured in the wild field were visualized as well to use as a reference.

Sample Preparation For RNA-seq
Gynoparae and males used for RNA sequencing were selected in adulthood, in which adults of gynoparae were gathered individually after they produced all offspring to eliminate the potential influence of embryos embedded in their ovaries. There were four biological replicates for each morph.
Each biological replicate contained 50 gynoparae or male adults, respectively. The samples were immediately frozen in liquid nitrogen and stored at -80℃.

Transcriptome Assembly And Gene Annotation
Total RNA was isolated from the pooled whole body of aphid samples using TRIzol reagent (Promega, USA) according to the manufacturer's instructions. The purity and integrity of RNA were assessed using an Agilent 2100 Bioannlyzer (Agilent, America). cDNA libraries were sequenced in paired-end modes on BGISEQ-500 system at Beijing Genomics Institute (BGI) (Shenzhen, China). 8 libraries were constructed with four biological replicates per morph. After sequencing, the raw reads were obtained and subsequently the adapter sequences and low-quality reads were filtered to acquire clean reads. Ultimately, FC > 2 and adjusted P < 0.001 were the thresholds to determine significant differences in gene expression. Gene ontology (GO) enrichment and significantly enriched KEGG pathways were identified through hypergeometric tests at the criteria of P < 0.05 by function phyper in R software.

RT-qPCR For Validation Of RNA-seq Data
Single-stranded cDNA was reverse-transcribed synthesized using 1 µg of RNA from each sample with (Thermo, America) using GoTaq® qPCR Master Mix (Promega, USA), initiated at 95 °C for 5 min, followed by 40 cycles of 95 °C for 5 s, 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-ΔΔCt method [43] and normalized by the housekeeping gene GAPDH [44]. Primers used for RT-qPCR validation were showed in Table S1. Statistical analyses were performed by Student's t-test using SPSS 19.0 software.
Significant difference was considered at P < 0.05. Values were reported as mean ± SE.

Ethics approval and consent to participate
We declare that our study was conducted in strict compliance with ethical standards. Collection of A. gossypii in this study did not require permits because it is a common cotton pest in China.

Consent for publication
Not applicable.

Availability of data and materials
All data generated or analyzed during this study are included in this published article (and its Supplementary Information files). Transcriptome datasets in this study are available at NCBI project PRJNA544300 with accession number SRP199362, and SRA with accession number SRS4812382.       Values in all panels represent mean ± SE. WLPF, wingless parthenogenetic female; WPF, wing parthenogenetic female; GP, gynopara.       Table S2.xlsx Table S1.docx Table S3.docx Figure. S2.tif