Transcriptome profiling of maternal stress‐induced wing dimorphism in pea aphids

Abstract Wing dimorphism, that is, wingless and winged forms, can be induced by maternal stress signals and is an adaptive response of aphids to environmental changes. Here, we investigated the ecological and molecular effects of three kinds of stress, namely crowding, predation, and aphid alarm pheromone, on wing dimorphism. These three stressors induced high proportion of up to 60% of winged morphs in offspring. Transcriptome analysis of stress‐treated female aphids revealed different changes in maternal gene expression induced by the three stressors. Crowding elicited widespread changes in the expression of genes involved in nutrient accumulation and energy mobilization. Distinct from crowding, predation caused dramatic expression changes in cuticle protein (CP) genes. Twenty‐three CP genes that belong to CP RR2 subfamily and are highly expressed in legs and embryos were greatly repressed by the presence of ladybird. By contrast, application of alarm pheromone, E‐β‐farnesene, caused slight changes in gene expression. The three factors shared a responsive gene, cuticle protein 43. This study reveals the adaptive response of aphids to environmental stresses and provides a rich resource on genome‐wide expression genes for exploring molecular mechanisms of ecological adaptation in aphids. OPEN RESEARCH BADGES This article has earned an Open Data Badge for making publicly available the digitally‐shareable data necessary to reproduce the reported results. The data is available at https://doi.org/10.5061/dryad.55b2b15.


| INTRODUC TI ON
Phenotype plasticity is a smart adaptive strategy universally acquired by many species in long-term evolutionary processes. In a changing environment, organisms possessing high plasticity in their behavior, morphology, and physiology would exhibit great fitness in survival and reproduction (Chevin, Lande, & Mace, 2010;Simpson, Sword, & Lo, 2011). Aphids (Hemiptera: Aphididae), such as the pea aphid (Acyrthosiphon pisum), are the representative species that survives an adverse and fluctuating environment through exhibiting such strong phenotypic plasticity. In warm spring and summer, the parthenogenetic wingless females can produce large numbers of wingless offspring to rapidly increase their population. Much of the aphids' energy is devoted to reproduction (Dixon & Howard, 1986;Zhang, Wu, Wyckhuys, & Heimpel, 2009). By contrast, wingless females produce winged offspring to migrate to new well-nourished plants when their host plant deteriorates and turn to sexual reproduction in autumn. The winged morphs input more energy on the development of wings and flight muscle rather than on reproduction (Brisson, Davis, & Stern, 2007). Therefore, the winged morphs present a smaller body size and less oviposition than the wingless morphs (Mackay & Wellington, 1975;Xu, Liu, Zhang, & Wu, 2011). The sexual female aphids produce eggs for overwintering through female and male copulation (Shingleton, Sisk, & Stern, 2003). Therefore, the interconversion between wing dimorphism and reproductive dimorphism can confer a selective advantage or disadvantage depending on environmental conditions. The various alternative phenotypes adopted by aphids to cope with the stressful environment make it an excellent model in studying adaptive dimorphism.
Abiotic and biotic environmental factors determine wing dimorphism in parthenogenetic female aphids. As a biotic factor, food source plays an important role in transgenerational wing dimorphism in aphids. Wingless females feeding on deteriorating plant sources promote production of winged offspring in Aphis craccivora and A. pisum (Johnson, 1966;Sutherland, 1969aSutherland, , 1969b. Supplying wingless aphids with a continuous flow of artificial diet with poor nutrition can increase the percentage of winged offspring in Myzus persicae (Harrewijn, 1976). In addition to nutrient substances, defensive and harmful substances may also influence transgenerational wing dimorphism in aphids. For example, precocene, a plant secondary metabolite, stimulates the production of winged offspring in A. pisum, which might be attributed to permanent juvenile hormone deficiency (Fridmancohen & Pener, 1980;Hardie, 1986).
Maternal A. fabae aphids tended by the mutualistic ant Formica fusca produce a high proportion of wingless offspring, and this phenomenon might be related to a juvenile hormone-related chemical from the ants (Kleinjan & Mittler, 1975). Moreover, the high-density signaling increases the proportion of producing winged individuals in Aphis glycines, Megoura crassicauda, and A. pisum (Ishikawa, Gotoh, Abe, & Miura, 2013;Martinez & Costamagna, 2018;Wilkinson et al., 2016). Such maternal stress-induced wing dimorphism has been a research hotspot for more than a century. However, whether or not these factors modulate the production of wing dimorphism through a common or distinct mechanism remains elusive.
The studies of wing morph differentiation and development have greatly laid a foundation of understanding the condition-dependent wing dimorphism (wing polyphenism) in viviparous female aphids (Brisson, 2010). Wing differentiation occurred at the early postembryonic instars, that is, at the first to second instar when the wing primordia and flight muscle appeared to develop, and subsequently, the flight muscles were degenerated in the destined wingless nymphs (Grantham & Brisson, 2018;Ishikawa, Hongo, & Miura, 2008). The wingless nymphs might transfer the nutrition and energy from the degraded flight muscle to postembryonic development, suggesting a tradeoff between dispersion and reproduction traits.
Recent breakthrough findings have unveiled some specific molecular mechanisms underlying the regulation of maternal stressinduced wing dimorphism. Transcriptome analysis of crowding stress-exposed pea aphids revealed prominent expression changes in genes associated with odorant binding, neurotransmitter transport, hormonal activity, and chromatin remodeling (Vellichirammal, Madayiputhiya, & Brisson, 2016). Specifically, crowding stress repressed the production of three monoamines, namely serotonin, dopamine, and octopamine in brain. The titers of these monoamines might signal developing embryos to be winged or wingless. Another study showed that the number of winged offspring was decreased by RNAi knockdown of the expression of a key octopamine synthesis enzyme, the tyramine β-hydroxylase (TβH; Wang, Zhang, Zhang, Tian, & Liu, 2016). The result partly confirmed the role of octopamine in the regulation of wing dimorphism in the pea aphids.
A subsequent study demonstrated that ecdysone signaling pathway is also involved in the transgenerational wing determination (Vellichirammal, Gupta, Hall, & Brisson, 2017). Females injected with ecdysone or its analog produced more winged offspring. On the contrary, ecdysone signaling suppressed by RNA interference targeting the ecdysone receptor (EcR) or by an EcR antagonist decreased the proportion of winged offspring. Insulin-related peptide 5 gene (Apirp5) regulated the alternation of wing morphs in pea aphids by affecting some physiological phenotypes such as body weight, embryo size, and carbohydrate and protein contents . The insulin signaling pathway could be also involved in the regulation of wing development in pea aphid. In addition, two laterally transferred viral genes, Apns-1 and Apns-2, could contribute to the production of winged offspring (Parker & Brisson, 2019). Despite these findings, the molecular mechanism for the regulation of maternal stress-induced wing dimorphism largely remains unclear.
This study aims to examine the effects of three representative external stimulations, that is, crowding, predation, and alarm pheromone, on the production of winged aphids and to further investigate the molecular mechanisms underlying maternal stress-induced wing dimorphism. Our findings may advance our understanding on the ecological adaptive mechanism of insect wing dimorphism under stressful environments.

| Insect
Green pea aphids (Acyrthosiphum pisim; Figure 1a) collected from Yunnan province, China in 2010, were reared on 2-week-old broad bean (Vicia faba) seedlings (Strain Linchan No. 7, Linxia Seed Company), in an incubator at 23 ± 0.5°C, 70% ± 5% relative humidity, and a 16L:8D photoperiod. The plants were watered daily and fertilized once a week. A parthenogenetic female was randomly selected to start a wingless clone at low density, approximately 20 aphids per seedling, to eliminate cross-generational effects (Sutherland, 1969a). All maternal aphids used in the experiment were from the same wingless clone and grown to the 8th day until maturity.
Adults of the two-spot ladybird, Adalia bipunctata, were collected from their natural habitats in Dongling Mountain, Beijing, China in July, 2015 and bred in groups of five in 8.5 cm diameter Petri dishes at 25 ± 0.5°C and a 16L:8D photoperiod. The beetles were fed with live pea aphids and readily laid clusters of eggs.

| Maternal stress I: Crowding
For this experiment, 1, 2, 5, 10, 30, or 50 adult individuals were placed on a leaf with the petiole wrapped with wet medical absorbent cotton and sealed by parafilm to maintain freshness in a Petri dish (8.5 cm) to manipulate pea aphid maternal density. A layer of filter paper dampened with a few drops of distilled water was placed on the inside bottom of the Petri dish to keep the moisture. Each density treatment was repeated 10 times. After 24 hr, the maternal adults in every treatment were transferred to a new Petri dish with a fresh leaf to continue the density treatment for another 24 hr. The newly born nymphs on the 2nd day were collected and transferred to a new seedling. After 6 days, the offspring reached the 4th instar stage when winged or wingless morphs are easily distinguishable by the presence of wing bud, and the winged offspring were then counted ( Figure 1b).

| Maternal stress II: Predation
Thirty adult aphids and a two-spot ladybird adult were placed together on a bean seedling covered with a transparent cylindrical bottle (diameter 7 cm, height 28 cm). The top of the bottle was covered with cotton gauze to prevent the ladybird from escaping. After 3 days, the ladybird and the surviving maternal aphids were removed from the plants (Weisser et al., 1999). In the control group, no presence of any predator was found. All offspring produced in 3 days F I G U R E 1 Experimental setup of maternal stress-induced transgenerational wing dimorphism in pea aphids. (a) Wingless and winged green pea aphids. (b) Experimental design for the induction of aphid wing dimorphism by three maternal stressors, that is, crowding, predation, and alarm pheromone, E-β-farnesene (EBF). To explore aphid density that caused the most significant crowding effect, we placed females (n = 1, 2, 5, 10, 30, or 50) in a Petri dish (8.5 cm) for 2 days, collected the offspring produced on the 2nd day, and counted the number of winged offspring after 6 days. To investigate the induction effect of predation, we placed one ladybird and 30 maternal aphids on a seedling for 3 days and collected all the offspring. No ladybird was placed in the control group. To test the effect of EBF, we placed five aphids in a Petri dish (8.5 cm) under EBF (5 µl, 100 ng/µl) or hexane stimulation five times per day (10:00, 12:00, 14:00, 16:00, and 18:00) for successive 2 days. The number of winged offspring on the 2nd day was counted 6 days later were left on the bean seedlings until they reached the 4th instar for the determination of the wing form and subsequent calculation of the percentage of winged offspring (Figure 1b). Each treatment was repeated 10 times.

| Maternal stress III: Alarm pheromone (EBF)
Five aphids were placed in a Petri dish (8.5 cm) with a fresh leaf that was treated as described above to keep its freshness. For EBF exposure, a filter paper strip (1 × 2 cm), which was added with 5 µl of 100 ng/µl EBF (Sigma-Aldrich, dissolved in HPLC-grade hexane) five times per day (10:00, 12:00, 14:00, 16:00, and 18:00) for 2 days, was placed in the Petri dish (Kunert, Otto, Röse, Gershenzon, & Weisser, 2010). The aphids in control groups received 5 µl of hexane at the same time points. The treatment and control Petri dishes were kept in different incubators to ensure that no alarm pheromone affect the control groups. The leaves were replaced daily. All offspring produced on the 2nd day were counted and transferred to fresh bean seedlings until they reached the 4th instar for the determination of the wing form ( Figure 1b). Each treatment was repeated 10 times.

| RNA extraction and preparation of Illumina sequencing libraries for RNA-Seq
After 2 days of stress exposure, the parent aphids were sampled and snap frozen in liquid nitrogen and stored at −80°C until RNA preparation. Fifty parent aphids were collected for each of four replicate samples for each treatment (i.e., stress exposure or control). Their offspring were reared on fresh bean seedling until they reached the 4th instar to confirm the efficacy of wing induction. The frozen adult aphids were disrupted and homogenized, and the mixture was oscillated by vortex blending for 2 min. RNA isolation and purification were then immediately performed on each sample using an RNeasy Mini Kit (Qiagen) according to the manufacturer's protocol. The purity of RNA was measured using Nanodrop 2000 (Thermo Scientific), and the integrity of RNA was determined by agarose gel electrophoresis and an Agilent 2100 TapeStation analysis (Aligent).
The mRNA was enriched and purified from the 20 µg total RNA of each sample by poly (T) magnetic beads, and then, the purified mRNA was cut into short clips using fragmentation buffer. The short clips were used as templates and reverse transcribed with the mixture of Superscript II (Invitrogen), Rnase H and DNA polymerase I (Illumina) to synthesize double-strands cDNA. After purification by QiaQuick PCR Purification Kit (Qiagen), cDNA was ligated with adenine at 3′-end and Illumina PE adapter. Adaptor-ligated cDNA was then subjected to PCR amplification. The cDNA library quality was determined by Agilent 2200 Bioanalyzer (Aligent). cDNA libraries were sequenced on the Illumina HiSeq 3000 by "paired-end" method at RiboBio Co. Ltd. A total of 24 libraries, including four repeats for each treatment, were constructed.
The original image data collected from the Illumina HiSeq 3000 platform were converted into raw sequencing reads. Raw reads were processed in order to obtain clean reads by filtering and removing reads with adaptors, short reads, and low-quality reads. Then, clean reads were produced and assembled into unigene using Trinity (Grabherr et al., 2011). Tophat2 was adopted to map clean reads to the pea aphid genome (AphidBase Official Gene Set v2, http://www. aphid base.com/) with parameter as follows: (-read-mismatches = 2,read-gap-length = 2; Trapnell, Pachter, & Salzberg, 2009). ANNOVAR was used to annotate the aligned transcripts in such an order of priority: exon, splicing region, intron, and intergenic region (Wang, Li, & Hakonarson, 2010). The gene model gff file was downloaded from the pea aphid genome.
The clean data were then transformed into RPKM (Expected number of Reads Per Kilobase of transcript sequence per Millions base pairs sequenced) in order to normalize the relative expression level of the matched unigenes (Mortazavi, Williams, McCue, Schaeffer, & Wold, 2008). Differentially expressed genes were determined by setting a fold change cutoff at least 1.5 and q value cutoff 0.05. Replicates with the low pairwise Pearson correlation coefficient (r < .7) between each others were excluded for further analysis (Table S1). Then, the average value of relative expression level of all the validated biological replicates in each treatment was used for transcriptome analysis.

| Bioinformatics analysis
Enrichment analysis for the supplied gene list was carried out based on an algorithm presented by GOstat (Beissbarth & Speed, 2004), with the whole annotated gene set as the background. The p-value was approximated using the chi-square test. Fisher's exact test was used when any expected value of count was below 5. Gene Ontology (GO) analysis including three ontologies, "Biological process," "Cellular component," and "Molecular function," provided a standardized gene functional classification method of all DEGs. The WEGO software was used for GO functional classification, and GO terms with a corrected p-value <.05 were defined as significantly enriched GO terms (Audic & Claverie, 1997). Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis was used to identify significantly enriched metabolic pathways or signal transduction pathways in DEGs based on the database at a criteria of p-value < .05 (Kanehisa et al., 2008).
Other graphs were plotted with GraphPad Prism5.

| Quantitative real-time PCR (qPCR)
Total RNA was extracted from the antenna, head, gut, cuticle, embryo, and leg using TRIzol reagent (Invitrogen) to test the tissue-specific distribution of cuticular protein (CP) genes of A. pisum.
The relative expression of mRNAs was quantified by qPCR utilizing a Light Cycler 480 instrument (Roche) in a 10 µl reaction volume consisting of 5 µl SYBR Green 1 Master Mix (Roche), 1 µl cDNA template, 3 µl H 2 O, and 0.5 µl primer F/R. qPCR procedure was set as follows: 95°C for 2 min, 45 cycles of 95°C for 30 s, 60°C for 30 s, and 72°C for 30 s. Gene-specific primers were designed with Primer Premier 5 according to the corresponding gene sequences, and unique amplification was verified by melting curve analysis (Table   S2). Three biological replicates were assayed for statistical analysis.
The 16S gene was selected as the internal control gene (Yang, Pan, Liu, & Zhou, 2014), and the relative expression level was analyzed by 2 −ΔΔCT method.

| Statistics
Data on the percentages of winged progenies were arcsine square root [X′ = arcsin(sqrt(x))] transformed and analyzed by SPSS 21.0 software (SPSS Inc). Significant differences between treatments were tested by one-way ANOVA or t test (p < .05).

| The three types of maternal stress-induced prominent wing dimorphism in pea aphids
We investigated whether or not three types of environmental stress, that is, crowding, predation, and alarm pheromone, induced wing dimorphism in pea aphid offspring. We first examined the effect of crowding stress by rearing wingless aphid mothers at different densities before tallying their winged offspring count. Figure 2a depicts the percentage of winged offspring produced by the females on the 2nd day of their crowding treatment. The solitary aphids did not produce any winged offspring. Low density (two aphids/leaf) had no remarkable effect on wing induction, leading to the winged progeny accounting for only 3%. However, the crowding density of five aphids/leaf significantly enhanced the effects (Student's t test, p < .0001) of this stressor and increased the percentage of winged offspring to 25%.
The percentage of winged offspring peaked at 60% for the group density of 10 aphids/leaf (Student's t test, p < .0001, compared with that of solitary aphids). Further increase in group density of up to 30 and 50 aphids/leaf also induced a high percentage of winged morphs, but the percentage was not higher than that with the density of 10 aphids/leaf. These results indicate that the production of winged aphids can be induced by a certain range of maternal crowding density, and excessive crowding does not enhance the induction effect.
Second, we investigated whether or not the presence of an aphid predator, that is, A. bipunctata, has a maternal effect on wing dimorphism. Pea aphid adults were subjected to predation stress by rearing 30 adult wingless aphids and one ladybird together on a seedling for 3 days. The percentage of winged offspring was significantly higher in mother aphids with predation exposure (25%) than that in the control without predation exposure (5%; p < .0001; Figure 2b).
Finally, the effect of the exposure of mother aphids to alarm pheromone (i.e., EBF) on the production of winged offspring was determined. The aphids reared at a density of five aphids/leaf responded strongly to EBF. The percentage of winged offspring in mother aphids with EBF exposure was significantly higher (29%) compared with that in the hexane control without EBF exposure (16%; p < .0001; Figure 2c). These results indicate that EBF stimulates winged morph production in crowded aphids.
F I G U R E 2 Effect of maternal stressors crowding, predation, and alarm pheromone on offspring wing dimorphism in pea aphids. (a) Effect of maternal density on the proportion of wing dimorphism. Female aphids were crowded at six densities for 2 days. Winged aphids were tallied on the 2nd day of crowding. Different letters indicate significant differences (p < .05) using Tukey's multiple range tests. (b) Effect of predation by a ladybird on offspring wing dimorphism. Thirty aphids were cohoused with one ladybird on a leaf for 3 days in each biological replicate. Winged aphids were tallied in the 3 days of treatment. (c) Effect of maternal exposure to alarm pheromone (EBF) on wing dimorphism. Five aphids on a leaf were exposed to 100 ng/µl EBF or hexane. Production of wing morphs was examined on the 2nd day after 2 consecutive days of treatment. Asterisk indicates significant differences between treatment and control (Student's t test, p < .05). Values in all panels represent mean ± SE. Ten biological replicates were used

| Transcriptome analysis reveals divergent molecular responses to the three maternal stress factors
To unveil the molecular mechanisms for the regulation of maternal stress-induced wing dimorphism, we sequenced the transcriptomes of the stress-exposed mother aphids that have produced winged morphs and of the nonstressed aphids. We excluded those biological replicates having low correlation of expression among replicates from the same treatment and included four replicates for crowding analysis, three replicates for predation analysis, and two replicates for EBF analysis.
A total of 12,417, 12,463, and 12,546 expressed transcripts were assembled in crowding, predation, and EBF treatments, respectively (Tables S3-S5 and Figure 3a). Different numbers of genes were significantly elicited by the three stressors, and most of the genes were induced by maternal crowding. A total of 489 genes were upregulated, and 290 genes were downregulated in response to crowding ( Figure 3b). In predation treatment, 183 genes were upregulated, and 118 genes were downregulated (Figure 3b).
Only a small number of DEGs (36), including 17 annotated DEGs, were detected in response to the EBF signal ( Figure 3b). The majority of these genes presented a moderate (i.e., 1.5 < fold change < 2) but significant (q < 0.05) change in their expression level (Table 1). Among the annotated DEGs, only cp43 and AF4/FMR2 family member 4 were upregulated. Among the downregulated DEGs, cp45 expression was highly and significantly induced by the EBF signal. Molecular chaperones, including three heat shock protein 68-like genes and one heat shock protein 70-like gene, were significantly downregulated. Other downregulated DEGs include titin, microtubule-associated protein futsch, nesprin-1, GPI-anchored adhesion, DNA ligase 1, aromatic-L-amino acid decarboxylase, longitudinal lacking protein, transcriptional regulator ATRX homolog, and lachesin (Table 1).
The majority of genes appeared a stressor-specific expression pattern. Venn diagram shows that 92%, 82%, and 50% of the DEGs in response to crowding, predation, and EBF, respectively, were stressspecifically expressed (Figure 3c). Only four common DEGs, including a CP gene cp43 and three unannotated genes, LOC100575566, LOC100574209, and LOC100574390 (Figure 3c  We further analyzed the functional involvement of the crowding-responsive genes and pathways by KEGG analysis. Three main functional categories, that is, lipid metabolism (7 pathways), amino acid metabolism (6 pathways), and carbohydrate metabolism (7 pathways), were significantly enriched (p < .05; Table S6). The major pathways involved in lipid metabolism are represented by glycerolipid/glycerphospholipid metabolism and cutin, suberine, and wax biosynthesis (Figure 4b). The subcategories with the largest number of DEGs in amino acid metabolism are responsible for the metabolism of starch, sucrose, ascorbate, and aldarate and the interconversions between pentose and glucuronate. Most of the DEGs involved in carbohydrate metabolism belong to the subcategories that control the metabolism of glycine, serine, threonine, arginine, proline, and tryptophan. These pathways are related to nutrient accumulation and energy mobilization. GO analysis revealed the involvement of transmembrane transport of anion, amino acid, organic anion, and intercellular signal transductions in the crowding-induced transgenerational wing production (Table S6 and Figure 4a).

| Cuticular protein genes were predominantly enriched and highly repressed in response to predation stress
We performed GO term analysis of the maternal transcriptomes induced by ladybird predation. In the molecular function domain, the structural molecule activity and structural constituent of cuticle were predominantly enriched (Figure 5a). In addition, genes functioning in neurotransmitter secretion, synaptic vesicle, and cell secretion were also significantly enriched and were closely associated with neuronal signaling and secretion (Table 2, Figure 5a and Table   S7). All the 28 annotated genes enriched in these terms belong to CP families (Table 2). KEGG analysis revealed significantly enriched four pathways, in which the annotated DEGs were all CP genes except for cprr 1-2 (Table S7 and Figure 5b).
We also performed a comprehensive analysis on these CP genes in pea aphid. There are 23 CPs among the 50 top DEGs based on their q value of expression difference in response to predation.
Expression level-q value correlation analysis showed that the majority of these CP genes presented high expression level relative to other DEGs (Figure 6a). RNA-seq data also indicated that the expression levels of the 23 CP genes in mother aphids were all significantly repressed by the presence of the ladybird (Figure 6b).
We further examined the tissue-specific expression of the 23 CP genes in pea aphids by quantitative PCR (see raw data in Table S8).
All the CP genes were highly expressed in the leg tissues. Twentythree CP genes, except gene cp23, were also highly expressed in embryos (although lower than that in legs; Figure 6c). These CP genes exhibit generally low expression in outer cuticles.
Finally, we also characterized the protein sequences of these CP genes. MUSCLE alignment of the predicted domain sequences indicated that the protein sequences of the 23 CP genes were highly similar, having 24 identical amino acid residues (Figure 6d). Protein family analysis based on the CuticleDB indicated that the 23 CP genes all belonged to RR2 subfamily and all presented a sequence motif that was homologous to the RR-2 motif reported in Anopheles sinensis (Figure 6e; Liu et al., 2017).

| D ISCUSS I ON
This study reported the prominent effects of different maternal stress factors, that is, crowding, predation, and alarm pheromone, on the induction of transgenerational wing dimorphism in aphids.
Increasing physical contact of pea aphids in a certain range of density caused many winged progenies (Lees, 1967;Martinez & Costamagna, 2018). Previous studies have revealed that contact between conspecific aphids or stimulation with a paintbrush induces the production of winged offspring (Johnson, 1965;Purandare, Tenhumberg, & Brisson, 2014). Predation and alarm pheromone are strong cues that can drive aphids to walk around and physically was more sensitive to ecological change and tended to produce more winged offspring under stresses. In addition, the content of carbohydrate and protein in the third-instar nymphs of the winged morph, whose wing primordia beginning to grow rapidly, is significantly higher than that in the wingless morph (Guo, Jiang, Yi, Liu, & Zhang, 2016). Thus, these two nutrients are vital for the development of primordia during this period. The nutrient accumulation and energy mobilization involved in the regulation of transgenerational wing dimorphism in A. pisum may represent a maternal adaptive strategy, in which resources are highly devoted to the development of winged embryos in a crowded population.
The results also imply an active involvement of amino acid catabolism in the regulation of wing dimorphism. In KEGG pathway for valine, leucine, and isoleucine biosynthesis, several genes were upregulated, including a gene encoding threonine dehydratase (LOC100165866), an allosterically controlled enzyme specific for Lisoleucine synthesis, and branched-chain-amino acid transaminase (LOC100167587) that catalyzes the final reaction of leucine, isoleucine, and valine biosynthesis (Umbarger, Neidhardt, & Curtiss, 1996).
The gene (LOC100161005) encoding pyrroline-5-carboxylate reductase that catalyzes the final step in proline synthesis was also upregulated. By contrast, the gene LOC100161119 encoding ornithine F I G U R E 5 Functional classification of GO and KEGG for DEGs in predation treatment. (a) GO classification of DEGs in predation treatment. The x-axis shows rich factor (Cluster frequency), and the y-axis shows significantly enriched GO terms. The color of symbol represents different −log 10 (Corrected p-value). The size of symbol represents the number of DEGs in each GO term. (b) KEGG classification of DEGs in predation treatment (p < .05).
These studies provide evidence for a close link between amino acid metabolism and transgenerational wing dimorphism in A. pisum.
In conclusion, our study provides a comprehensive understanding of the molecular genetic basis for the determination of winged descendants. The expression of CP genes related to transgenerational wing morphs can be induced by predation, crowding, and EBF stresses. The results imply a common regulatory mechanism underlying the transgenerational effects of these three stressful factors. However, stress factor-specific mechanisms are widespread in the regulation of wing form dimorphism. The responses to various stressors enable pea aphids to cope with environment changes and ensure their population persistence in growing seasons. Our comparative analyses have profound implications in understanding the divergent evolution of the regulation of phenotypic plasticity under distinct selection pressures.

ACK N OWLED G M ENTS
We thank Xieshuang Wang and Xiaolong Wang for photographing the winged and wingless aphids. This work was supported by the National Science Foundation of China (Grant No. 31572315 to L.C. and No. 31872304 to B.C.).

CO N FLI C T O F I NTE R E S T
None declared.

AUTH O R CO NTR I B UTI O N S
Li Chen and Bing Chen conceived the study, designed the experiments, analyzed the data, and revised the manuscript. Lin Hu performed the experiments on wing dimorphism induction and RNA sample collection. Wanying Gui analyzed the transcriptome data, performed qPCR, and wrote the manuscript.

DATA AVA I L A B I L I T Y S TAT E M E N T
The transcriptome data were deposited in NCBI SRA: PRJNA507326.

O PE N R E S E A RCH BA D G E S
This article has earned an Open Data Badge for making publicly available the digitally-shareable data necessary to reproduce the reported results. The data is available at https ://doi.org/10.5061/ dryad.55b2b15.