Outline of male wing morph in Aphididae
Male aphids were observed in 731 species of Aphididae, which belong to 12 subfamilies, 8 races, 153 genera (Table S2). At the subfamily level, male aphids were recorded mainly in Aphidinae (672), Saltusaphidinae (16) and Calaphidinae (15). At the races level, male aphids were observed primarily in Macrosiphini (490) and Aphidini (183). At the genera level, male aphids were mostly found in Aphis genus (137), followed by genera Uroleucon (56) and Macrosiphum (48). Furthermore, males were exclusively alate in 402 species, exclusively apterous in 317 species, and alternative wing (alate or apterous) in 12 species. We defined these three aphid species as Class I, II, III, in which males were winged, wingless, alternate wing, respectively. At the subfamily level, aphids in Tamaliinae (3), Eriosomatinae (2), Greenideinae (2), Pterastheniinae (1), Taiwanaphidinae (1) were belonged to Class I, in Chaitophorinae (10), Israelaphidinae (4), Lachninae (4), Anoeciinae (1) were belonged to Class II, in Calaphidinae (Class I, 14; Class II, 1) and Saltusaphidinae (Class I, 3; Class II, 13) were belonged to Class I or II. Emphatically, 376, 284, 12 species in Aphidinae belong to Class I, II, III, respectively (Fig. 1A). At the race level, aphids in Eriosomatini (2) and Greenideini (2) were belonged to Class I, in Siphini (10), Lachnini (2), Tramini (2) were blonged to Class Ⅱ, in Calaphidini (Class I, 14; Class II, 1) were belonged to Class I or II. Specifically, 69, 109, 4 species in Aphidini and 307, 175, 8 in Macrosiphini were belonged to Class I, II, III, respectively (Fig. 1B). At the gunus level, aphids in Therioaphis (12) and Diuraphis (6) were belonged to Class I and II, respectively; 136 and 72 species scattered in 11genera were categorized into Classed I, II, respectively; 114, 121, 12 species in 6 genera were classified into Class I, II, III, respectively (Fig. 1C). Obviously, males in Taiwanaphidinae, Eriosomatinae, Greenideinae, Pterastheniinae, Tamaliinae subfamilies were all winged, and males in Anoeciinae, Israelaphidinae, Lachninae, Chaitophorinae subfamilies were all wingless. At the races level, males in Eriosomatini, Greenideini races were all winged, and males in Siphini, Lachnini, Tramini races were all wingless. At the genera level, all males are wingless mainly in Therioaphis and Rhopalomyzus, and all males are winged mainly in Diuraphis, Lipaphis and Thripsaphis. Taken together, males in 402, 317, 12 observed aphid species were exclusively alate, apterous, and alternative wing, respectively (Fig. 1, Table S2)
Dynamic changes of male cotton aphid morphological characteristics
Morphological changes especially wing bud differentiation of males at different developmental stages were examined via stereomicroscope and shown in Fig. 2A. No distinct difference in the external morphologies between the first- and second-instar male nymphs were observed, except for a weeny wing bursa appeared on the prothorax of 2nd instar male nymph, implying the development beginning of wing primordia. At 3rd instar stage, the presumptive wing bud apophyses were visible on both sides of the second- and third-thoracic sections. Subsequently, the wing buds located on thoracic sections of 4th instar male can be observed clearly by naked eyes. After the last molting, adult male emerged with fully wings, sclerotized and hardening thorax, as well as black streaks on the yellow abdomen. According to these morphological characteristics, the development of wings of male cotton aphid can be divided into four developmental stages: preliminary stage (1st instar to 2nd instar), prophase (3rd instar), metaphase (4th instar), anaphase (5th instar).
Reproductive organ of male was presented in Fig. 2B: two accessory glands exist under the testes, and the length of male reproductive system is 0.83–1.07 mm. The duration time of 1st, 2nd, 3rd, 4th instar nymphs of male is 2.59 d, 2.65 d, 3.26 d, 4.5 d, respectively, and the longevity of male adult is 11.33 day (Fig. 2C). The body length of first-, second-, third-, fourth-instar nymphs and adult were 0.69, 0.83, 0.94, 1.21, 1.67 mm, respectively (Fig. 2D), the corresponding antenna length were 0.41 mm, 0.50 mm, 0.65 mm, 0.92 mm, 1.43 mm, respectively (Fig. 2E), and the head width were 0.20, 0.22, 0.24, 0.25, 0.37 mm, respectively (Fig. 2F). Generally, the duration times, body length, antenna length and head width of male increases slowly from 1st to 4th instar nymphs until the significant differences were observed between 4th instar nymph and adult (Fig. 2).
Male A. gossypii feeding behavior monitoring via EPG
Males feeding on cotton seedings showed typical feeding behavior patterns that have been previously documented with 6 clear probe events followed successively by Np waveform, C waveform, Pd waveform, E1 waveform, E2 waveform, F waveform, respectively (Fig. 3). Feeding behavior parameters of male cotton aphid are shown in Table 1. In 12 h, the average frequency of probing, duration time and the average duration of each wave pattern of male cotton aphid on secondary host plants were significantly different. Pd waveform contains the most probing numbers (242.64), followed by C, Np, E1, F, E2, with the numbers of 53.71, 33.07, 14.21, 13.21, 6.29, respectively. The longest total duration time was C waveform (397.76 min, 55.24%), which were much higher than Np (231.93 min, 32.21%), F (36.37 min, 5.05%), E1 (25.66 min, 3.56%), E2 (18.96 min, 2.63%) and Pd (18.41 min, 2.56%). Furthermore, the mean duration times are in order as follows: Np (1100.54 s), C (488.34 s), F (443.80 s), E2 (190.58 s), E1 (166.65 s), Pd (112.65 s). On the whole, the mean duration times of probing waves (C waveform, F waveform, Pd waveform) was much longer than that of feeding waves (E1 waveform and E2 waveform).
Time-course transcriptome profiling of male A. gossypii
RNA sequencing was performed to monitor the transcriptomic dynamics across of male cotton aphid development process. 918. 37 M clean reads were obtained with an average Q30 percentage of 88.41% from 12 nymphal and 3 adult samples, in which the mean genome and gene set mapping percentage of clean reads were 87.74% and 78.63%, respectively (Table S3). PCA and Pearson correlation coefficient analysis showed that the samples during different developmental stages exhibited significant separation, whereas variability was low between different replicates in the same stage (Fig. 4A and 4B). Specifically, the Pearson correlation coefficient of intra-group gene expression levels in each stage was all greater than 0.9, which indicated that our RNA-seq data were reliable for further analysis.
Clustering gene time course expression of male was implemented by Fuzzy c-means algorithm through Mfuzz software (Fig. 5 and Table S4). 11936 genes involved in male cotton aphid development were divided into 9 distinct clusters of temporal patterns which representing different expression kinetics (Fig. 5 and Table S4). Among these, cluster 8 includes the maximum gene number of 1732, followed by cluster 3 (1674), cluster 7 (1566), cluster 4 (1519), cluster 2 (1301), cluster 5 (1295), cluster 1 (994), cluster 9 (961), cluster 6 (894). Moreover, genes in cluster 2 were upregulated as a whole in male development, and genes in cluster 4 and 7 were downregulated mostly, whereas cluster 1, 3, 5, 6, 8 and 9 represent genes displaying a bi-modal expression pattern in male development. Furthermore, in the last molting of male (eclosion), genes in cluster 6 and 2 were abruptly upregulated, while genes in cluster 1 and 9 are downregulated dramatically. In addition, genes in cluster 3, 5 and 8 mostly expressed the highest levels in 3rd instar nymph stage, in cluster 1 and 9 exhibited their highest expressed levels in 4th instar nymph stage, which were the two vital periods for wing wud formation and fully wing emergence in male.
KEGG signaling pathways enrichment was performed to reveal potential function of genes classified into the 9 clusters above (Fig. 6, Table S4). Genes in cluster 2 were enriched in pathways of membranes transport (ABC transporters), lipid metabolism (fatty acid degradation, steroid biosynthesis), energy metabolism (oxidative phosphorylation), endocrine system (adipocytokine signaling pathway, aldosterone synthesis and secretion, PPAR signaling pathway, renin secretion), digestive system (cholesterol metabolism, fat digestion and absorption, gastric acid secretion), and environmental adaptation (thermogenesis). Cluster 3 and 4 genes were primarily associated with genetic information processing. Specifically, genes in cluster 3 involved in signaling pathways such as transcription (basal transcription factors, spliceosome), translation (RNA transport), “folding, sorting and degradation” (ubiquitin mediated proteolysis, protein processing in endoplasmic reticulum, proteasome), and genes in cluster 4 contributed to several vital signaling pathways such as replication and repair (base excision repair, DNA replication, mismatch repair), transcription (RNA polymerase), translation (ribosome), cell growth and death (cell cycle). Genes in cluster 5 were related to pathways of the neuroactive ligand-receptor interaction and circadian entrainment. Cluster 6 and 7 genes were primarily taken part in signal transduction, in which calcium signaling pathway, phospholipase D signaling pathway, and insulin secretion were enriched in the former, while mTOR signaling pathway, Jak-STAT signaling pathway were enriched in the latter. At last, genes in cluster 9 participated in pathway of citrate cycle (TCA cycle) and carbon metabolism.
DEGs identification and function analysis across the male development
Compared to 1st instar male, 921 genes were significantly up-regulated and 803 genes were significantly down-regulated in 2nd instar male. Then, when male nymph transformed from 2nd instar to 3rd instar, up to 3229 significantly differentially genes were identified, in which 2796 and 433 were up- and down-regulated, respectively. Moreover, with the drastic change of wing bud morphology, the expression level of 765 and 2932 genes were significantly higher and lower in 4th instar male nymph than in 3rd instar male nymph, respectively. After the last molting, male adult became fully winged, in which 1032 and 3586 genes were significantly up- down-regulated in comparison to 4th instar nymph, respectively (Fig. 7A). Interestingly, the number of up-regulated genes in the development of male were increased and peaked at the pairwise between 3rd instar and 2nd instar (M3 vs M2), then decreased. Overlaps of DEGs among different pairwise comparison in five developmental stages of male were exhibited via Venn diagram (Fig. 7B). 361, 593, 1171 and 1572 genes were exclusively significantly differentially expressed in comparing of 2nd instar and 1st instar (M2 vs. M1), 3rd instar and 2nd instar (M3 vs. M2), 4th instar and 3rd instar (M4 vs. M3), adult and 4th instar (MA vs. M4), respectively.
The up- and down-regulated DEGs in the four pairwise comparison groups above were performed KEGG pathway enrichment separately (Fig. 8, Table S5). In M2 vs. M1, mostly up-regulated DEGs were significantly enriched in organismal system (protein digestion and absorption, longevity regulating pathway, relaxin signaling pathway, ovarian steroidogenesis, dopaminergic synapse), cellular processes (Oocyte meiosis), while the down-regulated DEGs were involved in immune system (Toll and Imd signaling pathway). For M3 vs. M2, the mainly up-regulated DEGs were significantly enriched in environmental information processing (neuroactive ligand-receptor interaction, Ras signaling pathway), followed by organismal systems (dopaminergic synapse, circadian entrainment), while the down-regulated DEGs were significantly enriched in metabolism of phenylalanine and tyrosine, and steroid hormone biosynthesis. Compared to 3rd instar nymph, 4th instar male has up-regulated citrate cycle, oxidative phosphorylation, thermogenesis, citrate cycle and carbon metabolism, whereas the regulated DEGs were primarily associated with genetic information processing of transcription, translation, replication and repair (i.e. DNA replication, mismatch repair, homologous recombination, RNA transport, etc.). In comparison to 4th instar nymph, the up-regulated DEGs adult male were mainly enriched in mentalism of lipid and carbohydrate (such as fatty acid metabolism, steroid hormone biosynthesis, pentose and glucuronate interconversions), digestive system (i.e. salivary secretion, fat digestion and absorption), endocrine system (i.e. signaling pathway of PPAR and GnRH), the down-regulated DEGs were significantly enriched in genetic information processing signaling pathway, such as Ribosome and RNA polymerase.
Weighted gene co-expression network analysis
To investigate the key modules of genes associated with wings formation in male, weighted gene co-expression network analysis was preformed, in which 7645 co-expressed DEGs were divided into 17 modules (Fig. 9A and Table S6). Genes with highly similar expression patterns are clustered together into the same module. Meanwhile, the pattern of gene expression was often associated with phenotypic changes. The magenta module (70 genes) and salmon module (42 genes) was highly related to 1st instar nymph and 2nd instar nymph. The blue module (165 genes), brown module (90 genes), pink module (71genes), black module (76 genes), purple module (59 genes) and red module (81 genes) was positively correlated with 3rd instar nymph. The grey 60 module (29 genes), turquoise module (167 genes), green module (85 genes) and lightcyan module (30 genes), brown module (90 genes), pink module (71genes) was obviously associated with 4th instar nymph. The grey 60 module (29 genes), turquoise module (167 genes), the midnightblue module (30 genes), yellow module (87 genes), tan module (46 genes), cyan module (32 genes) and greenyellow (52 genes) was primarily related to 5th instar nymph (Fig. 9B).
To further identify the hub genes associated with phenotypic, a co-expression network was constructed and the connectivity patterns and hub genes of the relevant phenotypes were exhibited (Fig. 10). The Cytoscape software visualization result shown that 31 genes were highly interactivity in 1st instar nymph and 2nd instar nymph; 99 genes were highly interactivity in 3rd instar nymph; 43 genes were highly interactivity in 4th instar nymph. 53 hub genes were the most frequent interactions in 5th instar nymph stage (Fig. 10). Top 10 hub genes as ranked by intra-module connectivity were showed in Table 2. In the 1st instar nymph and 2nd instar nymph, the most highly connected gene was MCM4 (DNA replication licensing factor), followed by MCM5 (DNA replication licensing factor), CDC45 (cell division control protein 45), POLD1 (DNA polymerase delta subunit 1). In the 3rd instar nymph, the PLK1 (polo-like kinase 1) and AURKA were the most highly interactivity gene. In the 4th instar nymph, among the 10 top hub genes, such as, ADK (adenosine kinase), RSPH4A (radial spoke head protein 4A), HK (hexokinase) and KPNA2 (importin subunit alpha-2) were identified as the highly interactive gene. In the 5th instar nymph, 4CL (4-coumarate--CoA ligase) was the most highly interactivity gene, followed by YWHAB_Q_Z (14-3-3 protein beta/theta/zeta), PCK (phosphoenolpyruvate carboxykinase (GTP)) and SRPK3 (serine/threonine-protein kinase SRPK3). The gene expression patterns of 9 hub genes from the co-expression network analysis validated by RT-qPCR were significantly in accord with those obtained by RNA-seq mostly with the R2>0.6 ( Fig. 11).