Abstract

   While trade-offs between flight capability and reproduction is a common
   phenomenon in wing dimorphic insects, the molecular basis is largely
   unknown. In this study, we examined the transcriptomic differences
   between winged and wingless morphs of cotton aphids, Aphis gossypii,
   using a tag-based digital gene expression (DGE) approach. Ultra
   high-throughput Illumina sequencing generated 5.30 and 5.39 million raw
   tags, respectively, from winged and wingless A. gossypii DGE libraries.
   We identified 1,663 differentially expressed transcripts, among which
   58 were highly expressed in the winged A. gossypii, whereas 1,605
   expressed significantly higher in the wingless morphs. Bioinformatics
   tools, including Gene Ontology, Cluster of Orthologous Groups,
   euKaryotic Orthologous Groups and Kyoto Encyclopedia of Genes and
   Genomes pathways, were used to functionally annotate these transcripts.
   In addition, 20 differentially expressed transcripts detected by DGE
   were validated by the quantitative real-time PCR. Comparative
   transcriptomic analysis of sedentary (wingless) and migratory (winged)
   A. gossyii not only advances our understanding of the trade-offs in
   wing dimorphic insects, but also provides a candidate molecular target
   for the genetic control of this agricultural insect pest.

   Keywords: Aphis gossypii, trade-off, migration, digital gene
   expression, wing polyphenism.

Introduction

   Trade-offs between life-history traits have been studied in a variety
   of organisms [41]^1^-[42]^3. Individuals capable of fly benefit from
   the ability, but it is also costly on energy for flight mechanism and
   flight behavior. This energy, alternatively, can be used for other
   aspects, such as body size and reproduction [43]^4. In the past forty
   years, trade-offs between flight and reproduction has been widely
   studied [44]^5^-[45]^9. Wing dimorphic insects have been studied
   extensively and have emerged as a great model system for the study of
   trade-offs in the past several decades [46]^8^-[47]^10. Wing
   polyphenism has been observed in a wide range of insect species,
   including Hemiptera, Orthoptera, and Coleoptera, to facilitate their
   range expansion [48]^8^, [49]^10^-[50]^12. Wing dimorphic insects
   display commonly dispersing and non-dispersing morphs [51]^13^,
   [52]^14. Flight morphs, with long wings, flight muscles and flight
   capability, need additional energy and time for the migratory flight.
   Flightless morphs, with reduced or absent wings and underdeveloped
   flight muscles, show an earlier onset of oviposition and enhanced
   reproductive output [53]^8^, [54]^15^-[55]^17. Trade-offs have been
   studied at behavioral, morphological, physiological and biochemical
   levels in wing dimorphic insects, such as crickets, migratory locusts,
   and aphids [56]^16^, [57]^18^, [58]^19. For instance, biochemical
   aspects and intermediary metabolisms are investigated in a wing
   dimorphic cricket, Gryllus firmus [59]^20^-[60]^22. However, the
   genetic basis of biochemistry and metabolism related to trade-offs
   remains unclear.

   There are trade-offs related to dispersal and reproduction in aphids
   between two wing morphs. The wing polyphenism in aphids is well-studied
   [61]^23^, [62]^24. Winged and wingless aphids exist in both sexual and
   parthenogenetic stages of its life cycle. The migrants in spring,
   summer and early autumn are parthenogenetic generations, and the
   migration normally occurs in spring and summer [63]^25^, [64]^26. Wing
   polyphenism occurs primarily among parthenogenetic females, whereas
   wing polymorphism, which is genetically determined, has been found only
   in males [65]^14^, [66]^24. Parthenogenetic females with different wing
   morphs are descended from the same ancestry. A single genotype can
   produce various phenotypes as a result of exposure to various
   environmental conditions [67]^27. Such phenotypic plasticity is one of
   the main reasons contributing to the ecological success of aphids and
   their increasing pest status [68]^26^, [69]^28. Plasticity is not only
   an adaptive behavioral strategy that has evolved in most insect orders
   for coping with complex and uncertain environments, but also a key
   factor that leads to population expansion and infestation.
   Wing-polyphenic aphids exhibited distinct morphological differences
   between two wing types of aphids, and showed differences in their life
   history, particularly in their behavioral and physiological traits
   [70]^14^, [71]^29^, [72]^30. The biochemical basis of wing polyphenism
   related to energy resources, including glycogen, trehalose, and lipids,
   was investigated in cotton aphid, Aphis gossypii [73]^30 and grain
   aphid, Sitobion avenae [74]^31. Total lipid and triglyceride were
   significantly higher in winged morphs, whereas glycogen, soluble
   protein and water were higher in wingless morphs in both species.
   Recent genetic studies [75]^32^-[76]^34, especially the release of a
   draft genome of pea aphid Acyrthosiphon pisum [77]^35, lays the
   foundation for a better understanding of trade-offs between flight and
   reproduction at the molecular level. Genome-wide transcriptome analyses
   provide a new avenue for the future functional characterization study.

   Digital gene expression (DGE) is a tag-based RNA sequencing approach
   that allows millions of RNA reads, including differentially expressed
   transcripts, to be identified in a sample with or without any
   annotation information [78]^36. For aphid wing polyphenism study, DGE
   has several advantages - a greater sequencing depth, detection of
   unknown transcripts, and sensitivity [79]^37. Although external factors
   such as environmental cues, population density and nutrition status
   have been studied for years, the genetic basis of trade-offs in the
   cotton aphids, A. gossypii, a major cotton pest worldwide, has yet been
   examined in details. The objectives of this study are 1) to reveal
   differential gene expression profiles between winged and wingless A.
   gossypii adults, and 2) to establish a link between transcriptomic
   profiles and biochemical differences observed in previous studies. Our
   studies provide molecular basis that underpins trade-offs between
   flight capability and reproduction in A. gossypii and will facilitate
   the sustainable cotton aphid management through the RNA
   interference-mediated genetic disruption of its migratory behavior.

Materials and Methods

Aphis gossypii maintenance and sample preparation

   A single A. gossypii female adult was collected from a cotton field
   (40°01′N, 116°17′E) at the China Agricultural University, Beijing,
   China, in July of 2009. Aphid colonies generated from this individual
   A. gossypii had been established in the laboratory for more than one
   year before being sacrificed in the subsequent experiments.
   Specifically, parthenogenetic adults were provisioned with cotton
   seedlings (three-leaf stage, cultivar: Zhong 35) and maintained in a
   growth chamber at 60% relative humidity, and 16: 8 L: D photoperiod
   with 24±1˚C in the light and 22±1˚C in the dark [80]^30^, [81]^34^,
   [82]^38. Winged forms were induced by crowding in cotton leaves two
   weeks after A. gossypii was transferred to new host plants. The fourth
   instars were placed into individual dishes, and observed daily until
   molting. A total of 50 winged or wingless A. gossypii adults were
   collected separately within two days after molting. Aphids were frozen
   in liquid nitrogen immediately after the collection. RNA was extracted
   using the RNeasy Mini kit (Qiagen). Genomic DNA was digested by
   RNase-Free DNase Set (Qiagen). The quality and quantities of RNAs were
   determined using a Nanodrop2000 and Agilent bioanalyzer (Agilent).

DGE library preparation and sequencing

   High-throughput DGE was performed by Solexa/Illumina sequencing
   technology at the Beijing Genomics Institute, Shenzhen, China. DGE
   libraries (winged and wingless) were prepared using the Illumina Gene
   Expression Sample Prep Kit. Briefly, poly(A)^+ RNA was enriched from 6
   µg of total RNA using oligomagnetic beads. Double-stranded cDNAs were
   directly synthesized on the poly(A)^+ RNA-bound beads and then digested
   with a 4 base recognition enzyme NlaIII. The fragmentized cDNAs were
   purified from the magnetic beads, containing 3' ends and a sticky end.
   Then the Illumina adaptor 1 (sense:
   5'ACACTCTTTCCCTACACGACGCTCTTCCGATC3') was ligated to the 5' ends of
   these cDNA fragments. Digestion of MmeI and ligation of the Illumina
   adaptor 2 (sense: 5'GATCGGAAGAGCGGTTCAGCAGGAATGCCGAG3') produced short
   DNA fragments, with two adaptors and 21bp tags started with the CATG
   sites. Later on, Primer GX1 and Primer GX2 were added for PCR. After 15
   cycles of linear PCR amplification, 95bp fragments were purified using
   6% TBE PAGE Gel electrophoresis. After denaturation, the
   single-stranded molecules were attached to the sequencing chip for
   sequencing via Illumina HiSeq™ 2000 System. During this process,
   adaptor 1 was used as sequencing primer. Each tunnel of chip (flow
   cell) generated millions of raw tags with sequencing length of 35 bp
   [83]^39^, [84]^40.

Mapping of DGE tags

   Clean tags were generated by removing 3' adaptor sequences, empty reads
   (with only 3' adaptor sequences), low quality tags (tags with unknown
   nucleotide ''N''), tags with erratic lengths, and tags with copy number
   of 1. A total of 88,851 A. gossypii ESTs were selected at NCBI on Oct,
   28, 2010. After eliminating low quality sequences and Phrap assembly
   (minmatch: 42; minscore: 40, repeat_stringency:0.96), 23,630 unigenes
   were obtained and used as reference dataset ([85]Supplementary
   Material: Table S1). Clean tags were aligned to the reference dataset
   with only one nucleotide mismatch tolerance. The number of unambiguous
   clean tag corresponding to each transcript was counted and normalized
   to TPM (number of transcripts per million clean tags) [86]^36^,
   [87]^41.

Identification of differentially expressed transcripts

   The differentially expressed transcripts between winged and wingless
   patterns were identified using an algorithm developed by Audic and
   Claverie [88]^42. P value corresponds to differential gene expression
   test. FDR (False Discovery Rate) is a method of determining the
   threshold of P Value in multiple tests and analysis through
   manipulating FDR value [89]^43. In this study, we considered a
   transcript as differentially expressed when there was at least a 2-fold
   difference in expression between the two samples with a FDR value≤
   0.001.

   To functionally annotate these unigenes, distinct sequences were
   searched by BLASTx against the NCBI non-redundant protein database
   [90]^44 with a cut-off E-value of 10^-5, and assigned Cluster of
   Orthologous Groups (COG) and euKaryotic Orthologous Groups (KOG)
   protein classes with rpstblastn algorithm against COG and KOG databases
   in NCBI (E-value ≤10^-5) [91]^45^, [92]^46. The BLASTx result was
   imported into Blast2GO software (version 2.6.4) to map the Gene
   Ontology (GO) categories [93]^47, and the Kyoto Encyclopedia of Genes
   and Genomes (KEGG) pathways [94]^48. The enriched GO terms were
   developed by Blast2Go using Fisher's Exact Test with Multiple Testing
   Correction (FDR<0.05)[95]^47. The enrichment analysis of pathway was
   performed to further understand genes biological functions.

qRT-PCR analysis

   PCR reactions were performed in a 20 μl final volume containing 2×
   Power SYBR Green PCR Master Mix buffer (Applied Biosystems), 1 μM of
   each specific primer and 1 μl of diluted cDNA. PCR amplification was
   performed on an ABI Prism 7300 Sequence Detection System (Applied
   Biosystems) under the following conditions: 10 min at 95 °C for
   polymerase activation, and 40 cycles of 95°C for 15s, 60°C for 30s, and
   72°C for 30s for amplification signal collection. A default melting
   curve was produced to confirm single gene-specific peaks without primer
   dimerization present. Housekeeping genes, β-actin and 18S rRNA were
   subjected to the reference gene selection experiment. Eventually,
   β-actin was used as a constitutive control to normalize the expressions
   of target genes using a -ΔΔCt method.

   To validate DGE results, a total of 20 transcripts were chosen from the
   top 10 differentially expressed transcripts of winged and wingless
   libraries, respectively, and three significantly enriched pathways. The
   gene-specific primers designed for reference and target genes were
   listed in [96]Supplementary Material: Table S2. Three independent
   replications were carried out for the qRT-PCR validation study, with
   approximately 30 winged or wingless females for each replication.

Results

Summary of digital gene expression sequencing

   Reference sequences were obtained from 88,851 ESTs, of which 36 low
   quality (less than 100bp after remove vector) sequences were discarded.
   23,630 unigenes, containing 12,507 contigs and 11,123 singletons, with
   a mean length of 715 bp, were assembled by Phrap. The generated
   unigenes were subjected to BLAST search and functional annotation using
   various databases.

   In total, 5.30 and 5.39 million raw tags were obtained from winged and
   wingless adults, respectively. Distribution details of tags were shown
   in [97]Supplementary Material: Figure S1. Clean tags were mapped to
   9,011 and 10161 unigenes, respectively (Table [98]1). Sequencing
   saturation analysis shows the number of detected genes almost ceases to
   increase when sequencing amount (total tag number) reaches 2 million
   ([99]Supplementary Material: Figure S2), consequently, 5 million tags
   are sufficient for the expression profiling of cotton aphid. The total
   number of tags sequenced in the two samples was 10,479,972. Tags with
   the same sequence were identified as a distinct tag, and their
   distributions in different wing morphs were documented. Tags in the two
   wing forms were compared in total and in distinct (Figure [100]1).
   Surprisingly, morph-specific tags (tags detected in one morph only)
   accounted for a very small proportion (0.85% in winged morph, 2.12% in
   wingless morph) of the total clean tags, while the percentage of
   morph-specific tags increased obviously among the distinct tags.

Table 1.

   Sequencing summary.
   Sequence Winged Morph Wingless Morph
   Total % Distinct % Total % Distinct %
   Clean Tag 5221975 - 44414 - 5257997 - 61655 -
   Mapped to gene 4781916 91.57 21675 48.80 4345336 82.64 27998 45.41
   Unambiguous tag 1329628 25.4 20748 46.71 2322979 44.18 26749 43.38
   Unknown tag 440059 8.4 22739 51.20 912661 17.36 33657 54.59
   [101]Open in a new tab

Figure 1.

   [102]Figure 1
   [103]Open in a new tab

   Morph-specific sequencing tags in the winged and wingless A. gossypii.
   Number and percentage of morph-specific tags were documented in total
   clean tags (A) and distinct clean tags (B), respectively.

Differentially expressed transcripts in winged and wingless aphids

   We obtained 1,663 differentially expressed transcripts from a total of
   8,308 tag-mapped unigenes. Among them, 58 transcripts were expressed
   significantly higher in winged aphids, while 1,605 transcripts had a
   higher expression in wingless morphs (|log[2] ratio|≥1; FDR value≤
   0.001). A total of 334 (20.69%) differentially expressed transcripts
   had a differential ratio higher than 2 (|log[2] ratio|≥2; Figure
   [104]2), in which 3 transcripts were detected in winged forms only and
   34 transcripts were expressed in wingless forms exclusively
   ([105]Supplementary Material: Table S3). Top ten differentially
   expressed transcripts and their ratios are summarized in Table [106]2.

Figure 2.

   [107]Figure 2
   [108]Open in a new tab

   Expression dynamics of differentially expressed transcripts between
   winged and wingless A. gossypii. Fold Change represents the log[2 ]of
   the ratio of gene expression level in winged in relation to wingless
   adult females. Transcript expressed at a constant level (ratio of 1)
   has a log2 (ratio) of zero, which can be interpreted as "equally
   expressed among winged and wingless aphids". Positive number represents
   highly expressed transcripts in winged adult females. Negative number
   represents highly expressed transcripts in wingless adult females.

Table 2.

   Top ten differentially expressed transcripts between winged and
   wingless A. gossypii.
   Gene ID log[2]FC* Description
   AG10278 12.5 Flightin-like
   gi|289097293 11 3-ketoacyl-CoA thiolase, mitochondrial-like
   AG4984 10.1 MPA13 allergen-like
   gi|289151497 8.3 N/A
   AG12533 8.1 Phosphate carrier protein, mitochondrial-like
   AG13078 5.7 Adenine nucleotide translocator-like
   AG6703 5.6 PXMP2/4 family protein 4-like
   AG11738 5.3 Tectonin beta-propeller repeat-containing protein-like
   AG12565 3.8 Hypothetical protein LOC100165045 isoform 1
   AG11798 3.5 Very long-chain specific acyl-CoA dehydrogenase
   mitochondrial-like isoform 1
   AG10938 -9.7 Seryl-tRNA synthetase, cytoplasmic-like
   AG2471 -9.3 Hypothetical protein LOC100571294
   gi|289129407 -9.2 N/A
   AG9276 -9.1 RNA-binding protein PNO1-like
   AG11819 -9 T-complex protein 1 subunit alpha-like
   gi|289076469 -9 N/A
   AG7395 -8.9 Exosome complex exonuclease MTR3-like
   gi|289138983 -8.7 N/A
   AG4684 -8.6 N/A
   AG3568 -8.6 N/A
   [109]Open in a new tab

   “*”: FC is the fold change of gene expression level between winged and
   wingless morphs. Positive and negative values, respectively, denote
   highly and lowly expressed gene expression in winged adult females
   relative to wingless individuals. Assembly details are provided in
   [110]Supplementary Material: Table S3.

Functional analysis of differentially expressed transcripts

   Following assignment of BLASTx top hits, 1,306 differentially expressed
   transcripts were functionally annotated with four additional databases.
   Among 8,308 sequenced transcripts, 3,184 transcripts were annotated
   with GO terms, including 2,318 involved in biological process, 2,547 in
   molecular function and 1,954 in cellular component. As to 1,663
   differentially expressed transcripts, 763 transcripts were assigned
   with 3,168 GOs and can be classified into 42 secondary level categories
   (Figure [111]3). The category of biological process consists of 573
   transcripts at 22 secondary levels. Molecular function category has 619
   transcripts at 9 secondary levels. Cellular component includes
   491transcripts in 11 secondary categories.

Figure 3.

   [112]Figure 3
   [113]Open in a new tab

   Gene ontology classification of transcripts differentially expressed
   between winged and wingless A. gossypii. A total of 763 genes (45.88%
   of the differentially expressed transcripts) were categorized into
   three Gene Ontology classes, cellular components, molecular functions,
   and biological processes, respectively.

   All differentially expressed transcripts were enriched in 16 GO group
   when compared to all the A. gossypii unigenes ([114]Supplementary
   Material: Table S4). 7 of them were very deep level (specific) GO terms
   (Figure [115]4). Only one term (membrane) was underrepresented. In
   addition, transcripts were annotated by COG and KOG databases, by which
   676 transcripts were categorized into 23 COG functional groups (Figure
   [116]5) and 1,144 transcripts into 25 KOG classes (Figure [117]5). With
   KOG classification, in addition to the category of general function
   prediction only, the number of transcripts in category A (RNA
   processing and modification), O (Posttranslational modification,
   protein turnover and chaperones), and U (Intracellular trafficking,
   secretion, and vesicular transport) were 802, 748 and 759,
   respectively, higher than the number in the other categories.

Figure 4.

   [118]Figure 4
   [119]Open in a new tab

   Differentially expressed transcripts enriched specific Gene Ontology
   terms. All 763 differentially expressed transcripts are enriched in 16
   GO groups. 7 deep level (specific) GO terms were gained. 6 terms were
   overrepresented, and one term was underrepresented.

Figure 5.

   [120]Figure 5
   [121]Open in a new tab

   Functional classification of differentially expressed transcripts
   between winged and wingless A. gossypii. Based on Clusters of
   Orthologous Groups (COG, Figure_5A), 676 genes (40.65% of the
   Differentially expressed transcripts ) were annotated and categorized
   into 23 functional classes. In comparison, 1144 genes (68.79% of the
   Differentially expressed transcripts ) were annotated and categorized
   into 25 specific classes within the euKaryotic Orthologous Groups (KOG,
   Figure_5B). A: RNA processing and modification, B: Chromatin structure
   and dynamics, C: Energy production and conversion, D: Cell cycle
   control, cell division, chromosome partitioning, E: Amino acid
   transport and metabolism, F: Nucleotide transport and metabolism, G:
   Carbohydrate transport and metabolism, H: Coenzyme transport and
   metabolism, I: Lipid transport and metabolism, J: Translation,
   ribosomal structure and biogenesis, K: Transcription, L: Replication,
   recombination and repair, M: Cell wall/membrane/envelope biogenesis, N:
   Cell motility, O: Posttranslational modification, protein turnover,
   chaperones, P: Inorganic ion transport and metabolism, Q: Secondary
   metabolites biosynthesis, transport and catabolism, R: General function
   prediction only, S: Function unknown, T: Signal transduction
   mechanisms, U: Intracellular trafficking, secretion, and vesicular
   transport, V: Defence mechanisms, W: Extracellular structures, Y:
   Nuclear structure, Z: Cytoskeleton.

   Pathway analysis of differentially expressed transcripts was performed
   using the Kyoto Encyclopedia of Genes and Genomes (KEGG) database.
   Among 1,663 differentially expressed transcripts, 878 transcripts could
   be annotated in KO based on sequence homologies and mapped to 211 KEGG
   pathways (Table [122]3). Pathway-oriented analysis showed that the top
   three significantly enriched metabolic pathways in DGEs were ribosome,
   pyruvate metabolism and proteasome. Differentially expressed
   transcripts were also enriched in pathways relating to protein
   synthesis and degradation, lipid metabolism, immunity, RNA transport,
   and some signaling pathways.

Table 3.

   Top ten enriched metabolic pathways among differentially expressed
   transcripts.
   Pathway ID Description DET^a Reference^b P-value^c Q-value
   ko03010 Ribosome 49 144 1.35E-14 2.85E-12
   ko00620 Pyruvate metabolism 21 83 8.96E-05 9.45E-03
   ko03050 Proteasome 15 53 0.000239 1.68E-02
   ko03040 Spliceosome 44 268 0.001539 7.15E-02
   ko04141 Protein processing in endoplasmic reticulum 42 254 0.001695
   7.15E-02
   ko00280 Valine, leucine and isoleucine degradation 15 65 0.002402
   8.45E-02
   ko00062 Fatty acid elongation in mitochondria 6 16 0.004098 1.24E-01
   ko04612 Antigen processing and presentation 15 75 0.009839 2.60E-01
   ko04146 Peroxisome 23 134 0.011384 2.67E-01
   ko01040 Biosynthesis of unsaturated fatty acids 9 38 0.014408 2.67E-01
   [123]Open in a new tab

   “^a”: Differentially expressed transcripts

   “^b”: Reference sequences were extracted from NCBI (see Materials and
   Methods for details).

   “^c”: P-value is calculated based on Hypergeometric Test. Pathways with
   Q value≤0.05 are significantly enriched in differentially expressed
   transcripts. Pathway ID can be found in KEGG database.

qRT-PCR validation study

   qRT-PCR analysis was carried out to validate differentially expressed
   transcripts identified by DGE. All 20 genes showed the similar
   expression pattern as indicated by the DGE analysis (Table [124]4),
   suggesting that DGE approach can be a powerful tool for the discovery
   of candidate genes in a high-throughput manner.

Table 4.

   Differentially expressed transcripts validated by qRT-PCR analysis.
   Gene ID Description Fold Change* SE
   DGE qRT-PCR
   AG10278 Flightin 5961.13 1631.77 159.36
   gi|289097293 Kat 2068.00 11.06 1.93
   AG4984 Fabp 1092.00 12.25 2.18
   AG12533 phosphate carrier protein, mitochondrial-like 270.69 129.31
   12.41
   AG13078 Ant 52.62 55.20 0.78
   AG11798 very long-chain specific acyl-CoA dehydrogenase,
   mitochondrial-like 11.63 12.14 2.33
   AG10938 SerRS -837.00 -3.32 0.21
   AG13124 ribosomal protein S28e-like 11.45 26.07 5.96
   AG13172 ribosomal protein LP0 -4.53 -2.03 0.54
   AG12143 ribosomal protein S27A -4.37 -3.42 1.07
   AG13161 40S ribosomal protein S4-like -4.15 -1.18 0.17
   AG11336 40S ribosomal protein S21-like -4.03 -2.46 0.18
   AG13187 ribosomal protein L15 -4.03 -3.76 1.16
   AG8620 acylphosphatase 11.18 26.64 6.96
   gi|289143076 acylphosphatase-2-like -9.64 -2.06 0.46
   AG8620 pyruvate kinase-like -7.07 -1.87 0.69
   AG10855 malate dehydrogenase, cytoplasmic-like -6.52 -2.70 0.99
   AG12539 pyruvate dehydrogenase E1 component subunit alpha,
   mitochondrial-like -4.16 -2.51 1.21
   AG9726 proteasome activator complex subunit 3-like -4.88 -1.69 0.18
   AG12585 proteasome inhibitor-like -4.32 -1.10 0.05
   [125]Open in a new tab

   “*”: Fold change of gene expression level between winged and wingless
   morphs. SE is the standard error of qRT-PCT. Positive and negative
   values, respectively, denote highly and lowly expressed gene expression
   in winged adult females relative to wingless individuals.
   Morph-specific genes are highlighted in shade. Assembly details are
   provided in [126]Supplementary Material: Table S3.

Discussion

DGE analysis

   We sequenced wing-polyphenic aphids using the DGE approach and obtained
   10.69 million raw base sequences, 10.48 million of which were clean
   tags. Clean tags were compared in different libraries. Interestingly,
   compared to the clean tags in non-specific forms, clean tags of
   morph-specific forms account for a small fraction of the total clean
   tags (Figure [127]1A), however, they comprise a higher level of
   percentage in distinct tags (Figure [128]1B), indicating that
   morph-specific tags have a high diversity. This result may suggest that
   expression level of common genes which maintain the basic growth and
   development is higher relative to specific genes that play crucial and
   multiple roles in trait variation. Typically, the reference sequences
   used by DGE analysis are selected from genome sequence information of
   evolutionarily related species to the target species, especially when
   there is no genome information available for the target species. In
   current DGE analysis, we chose EST sequence of A. gossypii as reference
   sequences, 46,766 tags were mapped to 18,141 transcripts (Table
   [129]1).

Differentially expressed transcripts

   Through DGE analysis, we obtained a total of 1,663 differentially
   expressed transcripts, in which 58 transcripts were expressed
   significantly higher in winged aphids and 1,605 transcripts had a
   higher expression in wingless A. gossypii (Figure [130]2). This
   discrepancy is a reflection of the flight-reproduction trade-offs, in
   which genes associated with reproduction is likely over-numbered the
   flight capability. Among these differentially expressed transcripts,
   Flightin and seryl-tRNA synthetase (SerRS) was highly expressed in
   winged and wingless A. gossypii, respectively. Transcripts related to
   lipid metabolism and energy production were found at higher expression
   levels in migratory morphs (Table [131]2, [132]Supplementary Material:
   Table S2). Differential energy allocation is an important part of
   trade-offs between flight capability and reproduction. This result is
   consistent with the previous observation of the significantly higher
   triglyceride content in winged aphid versus the wingless aphid
   [133]^30^, [134]^31.

   Associated with lipid synthesis and degradation, transcripts of
   3-ketoacyl-CoA thiolase (KAT), very long-chain specific acyl-CoA
   dehydrogenase, medium-chain specific acyl-CoA dehydrogenase and
   long-chain-fatty-acid-CoA ligase 1 were shown significantly higher
   expression levels in winged morphs. Furthermore, 15 highly expressed
   “mitochondrial like” transcripts in winged morphs may contribute to
   energy production and thought to be involved in aspects of the citric
   acid cycle, respiratory chain, ATP synthesis and transport. These
   results suggest that lipid is the fuel for aphid migration. Similarly,
   transcripts involved in energy production showed elevated
   transcriptions in winged morphs of A. pisum [135]^34 and the brown
   planthopper, Nilaparvata lugens [136]^40.

   Differentially expressed transcripts were significantly enriched in 7
   GO terms, and 5 of them, structural molecule activity,
   ribonucleoprotein complex, translation, translation factor
   activity/nucleic acid binding, mitotic spindle organization, were
   directly or indirectly associated with protein translation and
   assembly. This is consistent with the pathway enrichment analysis. The
   remaining 2 term, membrane and intracellular non-membrane-bounded
   organelle, were associated with organelle, in which a series of
   biological process occurs, including partial protein synthesis
   processes in endoplasmic reticulum, fatty acid elongate in
   mitochondria, peroxisome break down very long chain fatty acids.

   qRT-PCR validation study confirmed the results from DGE analysis (Table
   [137]4). Fully replicated qRT-PCR analysis, to certain extent,
   complements the lack of DGE replications. All 20 genes subjected to the
   qRT-PCR analysis showed the similar trends with the DGE results. In
   comparison to DGE analysis, however, the fold changes in gene
   expression of Flightin, Fabp, Kat and SerRS were less dramatic in
   qRT-PCR analysis. Such discrepancy is not without precedent [138]^40.
   Due to the computation needs (the lowest expression level of any gene
   in the DGE analysis, by default, is 0.01 even if it is not expressed),
   the differential gene expression level is typically overestimated in
   the DGE analysis. In this study, the expression of Fabp, Kat and SerRS
   were morph-specific, and consequently, the result from the DGE analysis
   is likely overestimated.

Pathway enrichment analysis

   Through the Kyoto Encyclopedia of Genes and Genomes (KEGG) annotation
   system, differentially expressed transcripts were significantly
   enriched (Q<0.05) in three pathways, ribosome, pyruvate metabolism and
   proteasome (Table [139]3). Of the differentially expressed transcripts
   5.58% were in the ribosome pathway. While in the whole reference
   background, only 1.72% of transcripts were annotated in this pathway.
   The ribosome is required for protein synthesis, which is critical for
   all aspects of cell growth and division [140]^49. They are usually
   recognized as stable expression genes. However, in the current study,
   48 highly expressed transcripts in wingless morph and 1 in winged morph
   were coded for the ribosomal proteins. Ribosomal genes are typically
   highly expressed genes. Previously reported research results showed
   gene products that are components of ribosomes were over-representation
   in A. pisum [141]^34. One 40S ribosomal protein S20 displayed a higher
   expression level in brachypterous morph of N. lugens than in
   macropterous adult females [142]^40. This evidence convinces us of
   difference in ribosomal proteins between two morphs of A. gossypii. Our
   data indicated that more protein synthesis was detected in the
   non-dispersing morphs, which results in faster developmental rates and
   higher fecundity [143]^14^, [144]^50.

   It is expected that pyruvate metabolism was a differentially expressed
   transcripts enriched pathway. 21 transcripts in this pathway were
   screened as differentially expressed transcripts. Pyruvate is a key
   intersection in several metabolic pathways, and known as the “hub” of
   carbohydrate, fats and proteins [145]^16^, [146]^46. Pyruvate can be
   generated from glucose through glycolysis [147]^16. In recent study,
   wingless morphs showed significantly higher soluble sugar and glycogen
   content than winged morphs during adulthood in A. gossypii [148]^30 and
   in S. avenae [149]^31. The result of gene expression was consistent
   with our previous research in biochemical composition. Most transcripts
   in this pathway were highly expressed in wingless morphs. Our results
   suggest that sugar or glycogen metabolism is faster in wingless morphs
   and those carbohydrates were used for development and production. This
   is consistent with the contention that sugar and glycogen is not the
   fuel of aphids for migration [150]^51.

   26S proteasome, working with ubiquitin, operates the energy-dependent
   regulated proteolysis process in eukaryotic cells. The immunotype
   proteasomes, whose catalytic subunits are replaced by the related
   counterparts, were discovered recently in vertebrates [151]^52. The
   protein degradation process including immunity may be diverse in
   different morphs. In previous study, soluble protein content were
   higher in dispersal morphs than in reproduction morphs in A. gossypii
   [152]^30 and S. avenae [153]^31. As such, proteasome for protein
   degradation and protein turnover increased.

Genes of interests

   Genes which we are interested in include Flightin, fatty acid binding
   protein gene (Fabp), adenine nucleotide translocase gene (Ant),
   abnormal wing discs1 (Awd1) and abnormal wing discs2 (Awd2)
   ([154]Supplementary Material: Table S2).

   The primary difference between winged and wingless adults is the
   existence of flight apparatus, including wing and indirect flight
   muscle. For winged adults, genes expressed in muscle were identified in
   this study, Flightin and Fabp. Flightin is a multiply phosphorylated,
   myofibrillar protein in Drosophila indirect flight muscles (IFM), and
   an essential part of the flight muscle contractile mechanism for thick
   filament assembly and sarcomere stability in Drosophila
   [155]^53^-[156]^56. Particularly, Flightin exhibits the greatest
   differential transcript accumulation in winged morphs of A. pisum
   [157]^34 and N. lugens [158]^40. Another gene coding FABP in muscle was
   one of the top ten most differentially expressed genes in winged
   morphs. FABP will increase the solubility of fatty acids and thus lead
   to a rapid transport through the cytosol, since the lipid is the
   preferred fuel for insect migrating [159]^57 In recent research, free
   fatty acid content was significantly higher in winged than wingless
   adults before flight behavior in A. gossypii [160]^30 and A. pisum
   [161]^31.

   ANT belongs to the mitochondria carrier protein which is embedded in
   the inner membrane of mitochondria and required for importing ADP into
   the mitochondrial matrix throughout the ATP synthesis [162]^33. The
   gene expression of Ant gene is significantly higher in winged morphs
   than wingless morphs in A. gossypii (Table [163]2). The expression
   pattern of Ag Ant is similar as that of Mp Ant [164]^33, suggesting the
   need of energy for winged morphs during the radiation to diverse food
   resources and habitats [165]^58.

   The awd1 and awd2 showed a significantly higher expression level in
   wingless morphs than winged morphs in A. gossypii. The awd of
   Drosophila encodes a nucleoside diphosphate kinase (NDPK). Complete
   loss of function of the awd gene causes lethality in Drosophila after
   the larval stage. Mutations in awd gene cause abnormal tissue
   morphology and aberrant differentiation. [166]^59^, [167]^60. The awd
   plays an important role in development and differentiation of insects.
   They are not involved in the presumptive wing-patterning gene network
   [168]^61. Little attention was given to awd genes in the previous
   studies related to wing polyphenism.

Perspectives and future research

   This study sheds lights on the genetic basis of trade-offs between
   flight capability and reproduction. Genes involved in lipid metabolism
   show higher expression levels in winged A. gossypii than wingless
   morphs. In contrast, genes involved in sugar metabolism are elevated in
   wingless morphs. Our combined results suggest aphids using lipid rather
   than sugar or glycogen as fuel during migration. Flight ability of
   alate aphids depends on the energy supply of lipids. KAT, FABP, and
   very long-chain specific acyl-CoA dehydrogenase seems to be more
   important for winged morphs. Silencing of these genes using RNA
   interference (RNAi) may disrupt their energy supply, and compromise
   their ability to fly. RNAi through injection [169]^62^, [170]^63 and
   feeding [171]^64 methods have been well established in A. pisum. It is
   germane to investigate the pest control potential of these genes,
   including Flightin, through genetically interrupting the dispersal of
   A. gossypii.

Supplementary Material

   Table S1-S4, and Fig.S1 - S2.
   [172]Click here for additional data file.^ (6.7MB, xlsx)

Acknowledgments