Skip to main content

Comparative transcriptional analysis and identification of hub genes associated with wing differentiation of male in Aphis gossypii

Abstract

Background

Aphis gossypii Glover (Hemiptera: Aphididae), a worldwide polyphagous phloem-feeding agricultural pest, has three wing morphs (winged parthenogenetic female, gynopara, and male) in the life cycle. The exclusive males could fly from summer hosts to winter hosts, which are essential for gene exchanges of cotton aphid populations from different hosts or regions. However, the molecular mechanism of wing differentiation of male in A. gossypii remains unclear.

Results

Morphological observation of male A. gossypii showed that there is no distinct difference in the external morphologies of the 1st and 2nd instar nymphs. The obvious differentiation of wing buds started in the 3rd instar nymph and was visible via naked eyes in the 4th instar nymphal stage, then adult male emerged with full wings. According to morphological dynamic changes, the development of wings in males were divided into four stages: preliminary stage (the 1st instar to 2nd instar), prophase (the 3rd instar), metaphase (the 4th instar), anaphase (the 5th instar). Results of feeding behavior monitoring via EPG (electrical penetration graph) technology indicated that although the male cotton aphids had strong desire to feed (longer duration of C 55.24%, F 5.05% and Pd waves 2.56%), its feeding efficiency to summer host cotton was low (shorter E1 3.56% and E2 waves 2.63%). Dynamic transcriptome analysis of male aphid at 5 different developmental periods showed that in the 3rd instar nymph, the number of up-regulated DEGs was significant increased, and time-course gene transcriptional pattern analyses results also showed that numerous genes categorized in clusters 3, 5, and 8 had the highest expressed levels, which were consistent with morphological changes of wing buds. These results indicate that the 3rd instar nymph is the critical stage of wing bud differentiation in males. Furthermore, through pathway enrichment analysis of DEGs and WGCNA, it revealed that the neuroactive ligand-receptor interaction, Ras signaling pathway, dopaminergic synapse, circadian entrainment and the corresponding hub genes of PLK1, BUB1, SMC2, TUBG, ASPM, the kinesin family members (KIF23, KIF20, KIF18-19) and the novel subfamily of serine/threonine (Aurora kinase A and Aurora kinase B) probably played an important role in the critical stage of wing bud differentiation.

Conclusion

This study explored morphological changes and genes transcriptional dynamics males in cotton aphid, revealed the phenomenon of low feeding efficiency of winged males on summer host cotton, and identified key signaling pathways and potential hub genes potentially involved in wing bud differentiation of male in A. gossypii.

Background

Aphis gossypii Glover (Hemiptera: Aphididae), a worldwide agricultural pest with extremely complex life history, distributed in more than 171 countries and feeding on 65 plant genera, causes huge losses to crop yield and quality every year (Willis 2017). In recent decades, excessive use of chemical pesticides, over large-area cultivation of Bt cotton and global warming have led to continual outbreak of piercing-sucking pests especially A. gossypii in China (Chen et al. 2013; Lu et al. 2010; Sun et al. 2015). The damage spread of cotton aphid was closely related to three winged morphs, in which winged male was indispensable for life cycle (Ji et al. 2021). In northern China, A. gossypii differentiates into winged males, then flies to winter hosts to mate with sexual females which lays fertilized overwintering eggs in autumn. The wing formation of male promotes gene exchanges of colonies from different geographical regions and hosts, which is crucial for the survival and reproduction of A. gossypii (Kwon et al. 2017). Considerable attempts have been made to successfully induce male A. gossypii under low temperature and short light conditions (Ji et al. 2021). However, the dynamic morphological characteristics, the development process, the feeding behavior, and particularly the molecular mechanism of wing differentiation in male A. gossypii remains enigmas.

The regulation mechanism of male wing differentiation was well studied in pea aphid Acyrthosiphum pisum, in which the sex-linked locus on the X-chromosome named aphicarus determines the wing production of male (Braendle et al. 2005; Ogawa et al. 2013). That is, male wing polymorphism exists in A. pisum and is under the control of exclusively genetics. However, the differentiation of wings in male A. gossypii was quite different: all males had wings (Lee et al. 2015). Both males of A. pisum and A. gossypii had only one X-chromosome, but their wing differentiation patterns were obviously distinct, which were related to the biological characteristics of these two species: A. pisum was monophagous and no-host-alternating, whereas A. gossypii was polyphagous and host-alternating. Long-distance flights were not fully essential for the former in autumn, hence wings in males of A. pisum were alternative. In contrast, males are exclusively winged to migrate between primary host (usually woody) and the second host (usually herbaceous) in the latter (Ardie et al. 2017). Moreover, only 12 species belong to the no-host-alternating category own male wing polymorphism, which account for less than 2% of the recorded aphid species (Blackman et al. 2007). Therefore, the differentiation of male wings of cotton aphid is more representative in Aphididae to some extent. Studies on male wing differentiation of cotton aphid will benefit for understanding male wing production in host-alternating aphids and developing genetic control strategies against this pest through the disruption of male migratory.

Studies on aphid wing differentiation also have focused on parthenogenetic aphid. The generation of winged morph is mainly influenced by environmental factors such as crowding, photoperiod, nutrition and natural enemies, which is called non-genetic polymorphism. Traditional histological methods were widely used to explore wing differentiation in parthenogenetic aphid (Shull 1938; Johnson et al. 1960; Ganassi et al. 2005a, b). Histological observation results showed that the development of wing is inherent, whether or not they eventually develop into winged form, whereas wingless morphology is achieved by altering the original developmental pathway in the embryonic or nymphal stage (Braendle et al. 2006). The endogenous factors that regulate the wing differentiation of parthenogenetic aphid include ecdysone, octopamine, and insulin. Activation or inhibition of these three signaling pathways can affect the percentage of winged aphid separately (Wang et al. 2016; Ding et al. 2017). In contrast, although the generation of male A. gossypii can be induced by low temperature, crowding and short photoperiod, the corresponding process of wing bud differentiation is still not clear.

In the study, we manually counted and analyzed the distribution of wingless or wing male in recorded aphids derived from literatures firstly. Then, the morphological changes of male A. gossypii at different developmental stages were systematically observed and feeding behavior of males were monitored via electrical penetration graph technology (EPG). Subsequently, gene expression patterns clustering, differential expressed genes (DEGs) identification, KEGG enrichment, and weighted correlation network analysis were integrated to explore pivotal genes and signal pathways potentially involved in the regulation of wing differentiation in male of cotton aphid. Taken together, this study forms a basis for deciphering the molecular mechanisms underlying wing production of males in the cotton aphid and also provides valuable resources for future studies on male wing formation in host-alternating aphid species.

Materials and methods

Investigation of male aphid wing in Aphididae

The aphid species in which males were observed in field were counted from the book of Aphids on the World's Herbaceous Plants and Shrubs through keywords searching (Blackman et al. 2007). The retrieval results were listed in MS Excel and categorized manually at subfamilies, races, and genera levels, respectively.

Male induction and morphological observation

Aphis gossypii colony, collected originally in Anyang, China, was reared on cotton seedlings under controlled laboratory conditions for more than 50 generations before the start of all experiments (25 ± 1℃, 75% relative humidity, 16L: 8D photoperiod). Males were induced by rearing newly-born viviparous nymphs (the first generation, G1) at short day (SD) conditions as our previous methods (Ji et al. 2021). Offspring (the second generation, G2) produced by G1 mothers about 8 days later will develop into males. The duration time of each nymph stage and adult longevity were calculated with more than 20 males. The body morphology of male at different developmental stages and adult male testes were visualized using a SteREO Discovery V8 microscope (Zeiss, Germany), in which the testes were isolated in phosphate-buffered saline (PBS) buffer. In addition, the body and antennae lengths as well as head widths in each of male developmental stage were measured with more than 30 males.

Males feeding monitoring via electrical penetration graph recording

Electrical penetration graph (EPG) apparatus (DC-EPG Giga-8, Wageningen, Netherlands), placed in Faraday metal shield, were used to monitor the feeding behavior of male A. gossypii on cotton (CCRI 49) seedlings. The experiments were according to a previously described protocol with some improvements as follows (Elisa et al. 2016). New molting adult males were connected to the flexible gold wires (10 μm diameter × 4 cm length) through their dorsal thorax with the water-soluble conductive silver glue. After staved for one hour, each individual was placed on a leaf of cotton seedling at four-leaf stage. The EPG signal of each male was monitored continuously for 12 h via software Style + d, and the waveforms were analyzed by using software style + a. During the 12-h recording, six waveforms described by Tjallingii (1988) and Prado and Tjallingii (1994), were identified as follows: (1) Np (non-penetration) waveform appears due to stylet withdrawal; (2) waveform C (the probing period waveform consists of aphid mouth needles probing between mesophyll cells or outside the sieve molecules in phloem); (3) Pd waveform (potential drops, the aphid's mouth needles probe inside the living cells); (4) F waveform (the aphid mouth needles encounter mechanical resistance outside the plant tissues and cells); (5) E1 waveform (the aphids' mouth needles secrete saliva into the sieve elements); (6) E2 waveform (the aphid mouth needles feed passively within the screen tube). EPG variables were processed and calculated using EPG-Excel data workbook (Sarria et al. 2009). A total of 12 replicates were set with each male.

Sample preparation for RNA-seq

Whole bodies of males at different developmental stages including the 1st to the 4th instar nymphs and adults were collected, respectively. Total RNA of males in each developmental stage was extracted separately using TRIzol® reagent (Promega, USA) following the manual. The RNA concentration was assessed by a spectrophotometer Nanodrop 2000 (Thermo, USA), and the RNA quality was assessed using Agilent 2100. Synthesis of cDNA and Illumina library generation were completed at Beijing Genomics Institute (BGI) (Shenzhen, China) using BGISEQ-500 sequencing. There were three biological replicates for each development stage of males, and each biological replicate contained 100 males, respectively.

Transcriptome assembly and gene annotation

After sequencing, the raw reads were obtained and subsequently the adapter sequences and low-quality reads were filtered to acquire clean reads. The RNA-seq clean reads were aligned to the A. gossypii genome (Quan et al. 2019) using HISAT2 (Kim et al. 2015). Transcripts of each sample were reconstructed using StringTie and the new transcripts were identified using Cufflinks (Pertea et al. 2015; Trapnell et al. 2012). Those new protein-coding transcripts were predicted by CPC and aligned to protein databases, including Nr, KEGG by blastx with the cutoff E-value of 10E-5 (Altschul et al. 1990).

Differentially expressed genes identification

Gene expression levels in each library were calculated by fragments per kilobase per million reads (FPKM) method using Bowtie2 and RESM (Langmead et al. 2012; Dewey and Li 2011). Fold change (FC), the relative expression level of a gene in one stage to another, and P-value adjusted in multiple tests by false discovery rate (FDR), were used to determine the differentially expressed genes between two stages by DEseq2 (Love et al. 2014; Reiner et al. 2003). Ultimately, FC > 2 and adjusted P < 0.001 were the thresholds to determine significant differences in gene expression. All raw RNA-Seq read data are deposited in the NCBI Short Read Archive under the BioProject accession code PRJNA807005. The expression profiles were assigned to 9 clusters by using the fuzzy c-means algorithm implemented in the R package Mfuzz (v2.34.0, -c 9 -m 1.25) (Lokesh et al. 2007). Significantly enriched KEGG pathways were identified through hypergeometric tests at the criteria of P < 0.05 by function phyper (https://search.r-project.org/CRAN/refmans/DPQ/html/phypers.html). Pearson correlation coefficients of all gene expression level between each two samples were calculated by function cor and reflected in the form of heat maps (https://cran.r-project.org/web/packages/corrplot/vignettes/corrplot-intro.html). Princomp was used for PCA analysis (https://search.r-project.org/CRAN/refmans/ggfortify/html/fortify.prcomp.html).

Gene co-expression network construction and visualization

WGCNA as a typical algorithm is widely used in gene co-expression network identification (Horvath et al. 2008). We constructed Weighted Gene Co-expression Network by the WGCNA package in R software (Langfelder et al. 2008). The R package and its source code and additional material is freely available at https://horvath.genetics.ucla.edu/html/CoexpressionNetwork/Rpackages/WGCNA/. The experimental procedures referred to the previous method (Zhang et al. 2021). The resulting modules were represented and visualised using Cytoscape_v3.8.0 (http://www.cytoscape.org/).

qRT-PCR for validation of RNA-seq data

The relative expression of genes was quantified by qRT-PCR utilizing a Light Cycler 480 machine (Roche Diagnostics, Switzerland). The reactions included 2 μL cDNA template, 7.2 μL nuclease-free water, and 0.4 μL primer F/R. qRT-PCR procedure was set as follows: 94 ℃ for 30 s, 45 cycles of 94 ℃ for 5 s, 55 ℃ for 15 s, and 72 ℃ for 10 s. Melting curve was conducted to verify the specificity of amplification. Gene‐specific primers were designed with Primer Premier 6 and GAPDH gene was used as the internal control to normalize expression level (Gao et al. 2017). All primers are listed in Additional file 1: Table S1. The relative expression level of each gene was calculated by the 2−ΔΔCt method (Livak et al. 2001). Three biological replicates were assayed for statistical analysis.

Results

Outline of male wing morph in Aphididae and cotton aphid

Male aphids were observed in 731 species of Aphididae, which belong to 12 subfamilies, 8 races, 153 genera (Additional file 2: Table S2). At the subfamily level, male aphids were recorded mainly in Aphidinae (672), Saltusaphidinae (16), and Calaphidinae (15) (Fig. 1A). At the race level, male aphids were observed primarily in Macrosiphini (490) and Aphidini (183) (Fig. 1B). At the genera level, male aphids were mostly found in Aphis genus (137), followed by genera Uroleucon (56) and Macrosiphum (48) (Fig. 1C). 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 (Fig. 1, Additional file 2: Table S2). At the subfamily level, aphids in Tamaliinae (3), Eriosomatinae (2), Greenideinae (2), Pterastheniinae (1), and Taiwanaphidinae (1) belonged to Class I, those in Chaitophorinae (10), Israelaphidinae (4), Lachninae (4), Anoeciinae (1) belonged to Class II. Emphatically, 376, 284, 12 species in Aphidinae belonged to Class I, II, III, respectively (Fig. 1A). At the race level, aphids in Eriosomatini (2) and Greenideini (2) belonged to Class I, those in Siphini (10), Lachnini (2), Tramini (2) belonged to Class II. Specifically, 69, 109, 4 species in Aphidini and 307, 175, 8 in Macrosiphini belonged to Class I, II, III, respectively (Fig. 1B). At the genus level, aphids in Therioaphis (12) and Diuraphis (6) 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 families of Taiwanaphidinae, Eriosomatinae, Greenideinae, Pterastheniinae, and Tamaliinae subfamilies were all winged, and those in Anoeciinae, Israelaphidinae, Lachninae, and Chaitophorinae subfamilies 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.

Fig. 1
figure 1

Distribution of aphid species with observed male at subfamily (A), race (B) and genus (C) level. Some races are nameless and therefore not counted in the corresponding race. The statistics of genus level only the top 10 genera are showed in the figure. The male aphids with winged (alatae), wingless (apterous), and alternate wing (alatae and apterous) were defined as Class I, II, III

Morphological changes especially wing bud differentiation and reproductive organ of males were showed in Fig. 2A, 2B. No distinct difference in the external morphologies between the first- and second-instar male nymphs were observed. The presumptive wing bud apophyses were visible on both sides of the second- and third-thoracic sections in the 3rd instar nymph and can be observed clearly by naked eyes in the 4th instar nymph. Adult male emerged with full wings, sclerotized and hardening thorax, as well as black streaks on the yellow abdomen. Taken together, wing differentiation in male of cotton aphid can be divided into four stages: preliminary stage (the 1st instar to 2nd instar), prophase (the 3rd instar), metaphase (the 4th instar), anaphase (the 5th instar). Besides, the duration time of the 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 d (Fig. 2C). The body length of the 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, 0.50, 0.65, 0.92, 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 the 1st to 4th instar nymphs until the significant differences were observed in the 3rd, 4th instar nymph and adult (Fig. 2 and Additional file 3: Table S3).

Fig. 2
figure 2

Morphological changes and biological parameters of male cotton aphid in different developmental stages. A Morphological changes of male cotton aphid. M1, M2, M3, M4 and MA represent the 1st, 2nd, 3rd, 4th instar nymphs and adult, respectively. B Reproduction system of male cotton aphid. C Developmental duration, D Body length, E Antenna length and F Head width of winged male aphid in different developmental stages. Data were showed as mean ± SE

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 aphids are showed in Table 1. In 12 h, the average frequency of probing, duration time and the average duration of each wave pattern of male cotton aphids on summer host cotton were significantly different. Pd waveform contains the most probing numbers (242.64), followed by waveform of 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 was much higher than waveform of 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 (1 100.54 s), C (488.34 s), F (443.80 s), E2 (190.58 s), E1 (166.65 s), Pd (112.65 s). Obviously, the mean duration times of probing waves (C waveform, F waveform, Pd waveform) was much longer than those of feeding waves (E1 waveform and E2 waveform).

Fig. 3
figure 3

Six basic EPG waveform maps of male cotton aphid feeding on cotton. Np (non-penetration) waveform appears due to stylet withdrawal; waveform C (the probing period waveform consists of aphid mouth needles probing between mesophyll cells or outside the sieve molecules in phloem); Pd (potential drops, the aphid's mouth needles probe inside the living cells); F (the aphid mouth needles encounter mechanical resistance outside the plant tissues and cells); E1 (the aphids' mouth needles secrete saliva into the sieve elements); 6) E2 (the aphid mouth needles feed passively within the screen tube)

Table 1 The statistics analysis of EPG parameters

Time-course transcriptome profiling of male A. gossypii

RNA sequencing was performed to monitor the transcriptomic dynamics across of male cotton aphid development process. A total of 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 (Additional file 4: Table S4). Principal components analysis (PCA) and Pearson correlation coefficient analysis showed that the samples during different developmental stages exhibited significant separation, whereas the 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 datas were reliable for further analysis.

Fig. 4
figure 4

Principal components analysis and sample relationship heatmap of RNA-seq data. Samples from M1 (M1_1, M1_2, M1_3), M2 (M2_1, M2_2, M2_3), M3 (M3_1, M3_2, M3_3), M4 (M4_1, M4_2, M4_3) and MA (MA_1, MA_2, MA_3) were clustered into 5 groups, respectively. M1, M2, M3, M4 and MA represent the 1st, 2nd, 3rd, 4th instar nymph and adult, respectively (A). In the heatmap, dark green represents strong correlation, and light green denotes weak correlation. Each column and row correspond to one sample’s relationships with the other 18 samples including itself (B)

Clustering gene time course expression of male was implemented by Fuzzy c-means algorithm through Mfuzz software. 11 936 genes involved in male cotton aphid development were divided into 9 distinct clusters of temporal patterns which representing different expression kinetics (Fig. 5 and Additional file 5: Table S5). Among these, cluster 8 includes the maximum gene number of 1 732, followed by cluster 3 (1 674), cluster 7 (1 566), cluster 4 (1 519), cluster 2 (1 301), cluster 5 (1 295), 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 those in cluster 1, 3, 5, 6, 8, and 9 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 the 3rd instar nymph stage,  those in cluster 1 and 9 exhibited their highest expression levels in the 4th instar nymph stage, which were the two vital periods for wing wud formation and fully wing emergence in male, respectively.

Fig. 5
figure 5

Gene expression profiles across the development of male. Fuzzy c-means algorithm was implemented to cluster gene expression patterns in all developmental stages of male through Mfuzz software. The 11 936 genes expressed in male cotton aphid were divided into 9 distinct clusters of temporal patterns representing different expression kinetics

KEGG signaling pathways enrichment was performed to reveal potential function of genes classified into the 9 clusters above (Fig. 6, Additional file 5: Table S5). 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 and 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.

Fig. 6
figure 6

Significantly enriched KEGG signaling pathways corresponding to 9 distinct gene expression clusters. Significantly enriched pathways were determined at adjusted P < 0.05 through the hypergeometric test. Rich factor is the proportion of genes enriched into the target pathway

DEGs identification and function analysis across the male development

Compared with the 1st instar male, 921 genes were significantly up-regulated and 803 genes were significantly down-regulated in the 2nd instar male. When male nymph transformed from the 2nd instar to the 3rd instar, up to 3 229 significantly differentially genes were identified, in which 2 796 and 433 were up- and down-regulated, respectively. Moreover, with the drastic change of wing bud morphology, the expression levels of 765 and 2 932 genes were significantly higher and lower in the 4th instar male nymph than those in the 3rd instar male nymph, respectively. After the last molting, male adults became fully winged, in which 1 032 and 3 586 genes were significantly up- down-regulated in comparison to the 4th instar nymph, respectively (Fig. 7A). Interestingly, the number of up-regulated genes in the development of males were increased and peaked at the pairwise between the 3rd instar and the 2nd instar (M3 vs M2), then decreased dramatically (M4 vs M3). Overlaps of DEGs among different pairwise comparison in five developmental stages of male were exhibited via Venn diagram (Fig. 7B). And 361, 593, 1 171, and 1 572 genes were exclusively significantly differentially expressed in comparison of the 2nd instar and the 1st instar (M2 vs. M1), the 3rd instar and the 2nd instar (M3 vs. M2), the 4th instar and the 3rd instar (M4 vs. M3), adult and the 4th instar (MA vs. M4), respectively.

Fig. 7
figure 7

Pairwise comparison of DEGs across the male development process. A The number of DEGs in the M2 vs. M1, M3 vs. M2, M4 vs. M3 and MA vs. M4. B Distribution of DEGs in the M2 vs. M1, M3 vs. M2, M4 vs. M3 and MA vs. M4. M1, M2, M3, M4 and MA represent the 1st, 2nd, 3rd, 4th instar nymphs and adult, respectively

The up- and down-regulated DEGs in the four pairwise comparison groups above were performed KEGG pathway enrichment separately (Fig. 8, Additional file 6: Table S6). 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 the 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 with the 3rd instar nymph, the 4th instar male has up-regulated 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.) (Additional file 6: Table S6). In comparison of the 4th instar nymph, the up-regulated DEGs adult males 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.

Fig. 8
figure 8

Significantly enriched KEGG signaling pathways of differentially expressed genes derived from pairwise comparison of M2 vs. M1, M3 vs. M2, M4 vs. M3, and MA vs. M4. Significantly enriched pathways were determined at adjusted P < 0.05 through the hypergeometric test

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 7 645 co-expressed DEGs were divided into 17 modules (Fig. 9A and Additional file 7: Table S7). Genes with highly similar expression patterns are clustered together into the same module. Meanwhile, the pattern of gene expressions was often associated with phenotypic changes. The magenta module (70 genes) and salmon module (42 genes) were highly related to the 1st and 2nd instar nymphs. 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 the 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) were obviously associated with the 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) were primarily related to the 5th instar nymph (Fig. 9B).

Fig. 9
figure 9

Network analysis dendrogram showing co-expression modules identified by weighted gene co-expression network analysis (WGCNA) of differentially expressed genes (DEGS). (A) Dendrogram plot with color annotation. The genes in the same branch could be assigned to different modules. (B) Heatmap plot of module-trait relationships. Red represents high adjacency (positive correlation) and blue represents low adjacency (negative correlation). The color scale on the right shows module-trait correlation from -1 (blue) to 1 (red)

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 showed that 31 genes were highly interactivity in the 1st and 2nd instar nymphs; 99 genes were highly interactivity in the 3rd instar nymph; 43 genes were highly interactivity in the 4th instar nymph. 53 hub genes were the most frequent interactions in the 5th instar nymph stage (Fig. 10). Top 10 hub genes as ranked by intra-module connectivity were showed in Table 2. In the 1st and 2nd instar nymphs, 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, 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 qRT-PCR were significantly in accord with those obtained by RNA-seq mostly with the R2 > 0.6 (Fig. 11).

Fig. 10
figure 10

Co-expression network of relevant modules and corresponding related phenotypes in specific developmental periods. Node size: larger indicates a more interactive, smaller indicates a less interactive. Node color: Red indicates up-regulated gene, Green indicates down-regulated gene

Table 2 The top 10 high interactive gene
Fig. 11
figure 11

Quantitative qRT-PCR of 9 DEGs in male cotton aphids. The genes expression level at M1 stage was used as a control for normalization. Pearson correlation coefficients were calculated by comparing results between qRT-PCR and RNA-seq data for each sample

Discussion

Morphology variation and wing bud differentiation of winged male

Males in Aphididae were usually neglected except for A. pisum, and only 731 species were recorded among more than 4 000 aphids (Blackman et al. 2007), of which 402, 317, 12 were exclusively winged (55%), wingless (43.36%), alternative wing (1.64%) (Fig. 1, Additional file 2: Table S2), which were defined as Class I, II, III here, respectively. Obviously, aphids of Class I such as A. gossypii were the representative species, which deserved to be explored further. In aphids, environmental conditions experienced during embryonic or nymphal stages determine the divergent adult phenotypes (Sutherland 1969). To better understand wing differentiation of male cotton aphid, we systematically observed their wing development dynamics (Fig. 1). Their morphological changes across different stages were basically consistent with those observed in A. pisum and Megoura viciae (Ganassi et al. 2005a, b; Benedetti et al. 1991). In A. gossypii and A. pisum, obvious morphological changes of wing bud were found in the 3th instar nymph, and the wing buds located on thoracic sections of the 4th instar male can be observed clearly by naked eyes (Ishikawa et al. 2008). However, wing buds of M. viciae in the 4th instar nymph were similar to those in the 3rd instar, in which only the cells were more regularly arranged in two layers with again an evident slot between them. Specially, although no significantly external morphological changes were observed in the 1st and 2nd instar nymphs of male cotton aphid, whereas a weeny wing bursa appeared on the prothorax of the 2nd instar male nymph, implying the development of wing primordia, which was consistent with previous reports that internal tissue structures already show differentiation from the 2nd instar onward (Ganassi et al. 2005a, b; Benedetti et al. 1991; Ishikawa et al. 2008). These results suggest that the divergent adult phenotypes were probably determined in embryonic or neonatal stage.

Feeding behavior of winged male A. gossypii

The results of feeding behavior analysis showed that male cotton aphid produced 6 kinds of waveforms on cotton seedlings (Fig. 3). The puncture waveform of male mouth needle outside the phloem (Pd, C, F waveform) contains the most probing numbers and longest duration time, which showed strong willingness to probe for feeding (Table 1). The puncture waveform of male needle in phloem (E1 and E2 waveform) contained the least number of puncture times and shortest duration time (Table 1), which indicated that although male has a strong desire to ingest, the feeding efficiency was lower. Consistent results have been reported in the comparison of the probing behaviors between alate and apterous of M. persicae on Solanum tuberosum, in which feeding activity of winged morph was lower (Boquel et al. 2011). Similar studies on aphid of A. fabae showed that the winged morph had higher feeding activity on the winter host (Powell et al. 2010). We hypothesized that lower feeding efficiency of male cotton aphid adults to the summer host (cotton) is probably due to enough energy that had been stored in nymphal stages or preference to winter host in adulthood. In fact, missions of wing male aphids are to migrate to winter host and mate with ovipara to complete life cycle in turn. Exactly, the passive feeding waveform E2 has the least number of waves and the shortest duration, which also hints that winged male A. gossypii adults have low adaptability to summer host cottons.

Signaling pathway in critical period of wing differentiation of male A. gossypii

Transcriptomic dynamics across male cotton aphid development were obtained via RNA-seq. The largest number of up-regulated DEGs occurred during the development from the 2nd instar to 3rd instars (Fig. 7A), and time-course gene transcriptional pattern analysis results also showed that numerous genes in clusters 3, 5, and 8 had the highest expressed level in the 3rd instar stage (Fig. 5). Actually, dramatically morphological changes, especially the wing bud apophyses, can be observed from the 2nd instar to the 3rd instar (Fig. 2A). Hence, we speculated these highly expressed genes in the 3rd instar above probably determined the wing bud differentiation in male. Compared with the 2nd instar nymph, the up-regulated DEGs in the 3rd instar nymph were enriched in neuroactive ligand-receptor interaction and Ras signaling pathway, dopaminergic synapse, and circadian entrainment (Fig. 8, Additional file 6: Table S6). Furthermore, highest expressed genes in the 3rd instar (cluster 3, 5, 8) were mostly involved in genetic information processing, followed by organismal systems and environmental information processing (Fig. 5, Additional file 5: Table S5). Obviously, those signaling pathway accelerate wing bud differentiation and development collaboratively and should get receive more attention in future.

From the 3rd, 4th instar to adult, the number of down-regulated DEGs were larger than up-regulated DEGs (Fig. 7A). Interestingly, the up-DEGs were mostly related to organismal systems, lipid metabolism, carbohydrate metabolism, energy metabolism (Additional file 6: Table S6). In addition, the down-regulated DEGs were mostly associated with genetic information processing including translation, replication and repair (Additional file 6: Table S6). These results indicated that the ability of cell and some tissue differentiation decreases with the maturation of aphid internal organs, and the entirely metabolism strengthen meanwhile. Specially, in MA vs. M4, the large number of up-regulated DEGs were mainly enriched in the fatty acid metabolism, fat digestion and absorption, PPAR signaling pathway, steroid hormone biosynthesis, glycerolipid metabolism and so on (Fig. 8). These signaling pathways were also reported significantly up-regulated in alate brown citrus aphid Toxoptera citricida and in long winged Sogatella furcifera (Shang et al. 2016; Liang et al. 2018). Moreover, silencing of lipid synthesis and degradation genes in T. citrida leading to underdeveloped wings revealed that both lipid metabolism and energy supply are essential for aphid wing bud differentiation (Shang et al. 2016; Liang et al. 2018).

Candidate hub genes of wing differentiation in males

WGCNA were used to identify the specifically expressed hub genes potentially involved in wing differentiation of male in A. gossypii among 7 645 DEGs, and the top 10 ones are shown as the order of interactivity in Table 2. Among the top 10 nodes in the younger larval stages (the 1st and 2nd instar nymphs), the CDC45, POLD1, RFA1, RFA2, RRM1, CTF4, ASPM, and the mini-chromosome maintenance (MCM) gene family members (Mcm4, Mcm5, Mcm6) were identified as candidate hub genes. CDC45 is an essential gene for the initiation of DNA replication and elongation (Luis and Ponting 2011). Cdc45 can bind to Mcm2-7 and GINS proteins to form a complex that is competent to unwind double stranded DNA and regulate DNA replication (Petojevic et al. 2015). The phenotype of lacking functional CDC45 gene was reported in the mice, and the loss of CDC45 gene function impaired proliferation of the inner cell mass (Yoshida et al. 2001). The MCM gene family play an important role in eukaryotes DNA replication, and silencing of Mcm4 results in inhibition of de novo DNA synthesis and polyploidization, and arrested oocyte maturation and ovarian growth in the migratory locust (Guo et al. 2014). The POLD1 gene (also known as CDC2) encodes the catalytic subunit of eukaryotic DNA polymerase-δ, and previous studies have shown that the POLD1 mutation was significantly impeded the growth and proliferative ability of the cells (Liu et al. 2020; Song et al. 2009). In addition, The RFA1, RFA2, RRM1, CTF4, and ASPM, are involved in eukaryotes DNA replication and cell proliferation (Longhese et al. 1994; Lei et al. 2012; Villa et al. 2016; Pan et al. 2008). In summary, all 10 candidate hub genes in younger larval stages (the 1st to 2nd instar stage) are involved in DNA replication and cell proliferation. Combined with the morphological study above, these candidate hub genes probably respond to the environment-induced wing plasticity in the embryonic stage, the 1st and 2nd instar nymphal stage. However, it remains to be determined how they respond to the environment-induced and regulate the early development of wing buds during the preliminary stage.

In the 3rd instar nymph, the wing buds begin to develop further and the obvious morphological changes of wing bud were observed. The genes related to cell division, sister chromatid cohesion, mitotic entry such as PLK1, BUB1, SMC2, TUBG, ASPM, the kinesin family members (KIF23, KIF20, KIF18-19) and the novel subfamily of serine/threonine (Aurora kinase A and Aurora kinase B) were identified as candidate hub genes potentially involved into wing development in the prophase (Table 2). However, these core genes involved in cell division and sister chromatid separation have been poorly studied in insects. Hence, how they regulate the wing differentiation need to be further explored. Studies in mammalian cells have shown that the loss of PLK1 function leads to a short delay entry into mitosis and subsequently often die by apoptosis (Sumara et al. 2004). In the Plk1-depleted cells, Aurora B inhibitor can silence the spindle checkpoint by stabilizing microtubule kinetochore interaction, and the cell exiting mitosis without chromosome segregation and cytokinesis (Sumara et al. 2004). Bub1 as a multi-task protein kinase plays an important role in mammalian chromosomal instability (Kawashima et al. 2010). The SMC2, TUBG, ASPM, the kinesin family members (KIF23, KIF20, KIF18-19) and the novel subfamily of serine/threonine (Aurora kinase A and Aurora kinase B) has been reported to be mainly associated with chromosome segregation, condensation and cell cycle in mammals (Pan et al. 2008; Strunnikov et al. 1995; Joshi et al. 1992; Zhu et al. 2020; Martin et al. 2013; Goldenson et al. 2014). To sum up, we hypothesized that these genes probably regulate the further development of wing buds by regulating cell division, sister chromatid cohesion and mitotic entry.

In the 4th instar nymph, the wing buds located on thoracic sections of male can be observed clearly by naked eyes. RSPH4A, KPNA2, DNALI, TOM70, TEKT1, DYX1C1, DNAH, and three kinase genes (ADK, HK, pfkA) showed high network interaction (Table 2). The RSPH4A, DNALI, and DNAH were reported to be closely related to the movement of cilia and sperm axoneme, and the loss of genes function led to male and female subfertility (Castleman et al. 2009; Rashid et al. 2006; Whitfield et al. 2019). Similarly, loss of TEKT1 gene function in Bactrocera dorsalis was also reported to lead to male infertility (Sohail et al. 2019). Interestingly, the DNALI1 and TEKT1 were highly expressed in male testis (Rashid et al. 2006; Sohail et al. 2019). Combining the functions of these genes and the high interactivity of these genes in the 4th instar nymphal stage suggests that the male cotton aphid not only undergo the further development of wings, but also undergo the development of the male reproductive system. In addition, these genes related to the glucose homeostasis and energy metabolism (HK and pfkA) also play a vital role in providing energy for flight muscle (Storey 1980; Siedler et al. 2012; Dixon et al. 1993). Transcriptome studies of Sogatella furcifera have shown consistent results that genes associated with energy metabolism and flight muscle is crucial in plasticity development of wing (Gao et al. 2019). Another studys in brown planthopper and S. furcifera have shown that insulin signaling genes also play an important role in the differentiation of wing type (Gao et al. 2019; Xu et al. 2015). In recent years, juvenile hormone and ecdysone signaling have also been shown to be involved in the development of wing plasticity (Zhang et al. 2019; Vellichirammal et al. 2017).

After last molting, adult male emerged with full wings. The genes associated with energy metabolism such as YWHAB_Q_Z, PCK, ACADM, and MDH2 showed high network interaction (Table 2). On the whole, the energy utilization of insects in flight can be divided into three types: sugar utilization type, fat utilization type and both of them. Locusts use carbohydrates and fat as the energy supply for initial flight and continuous flight, respectively (Rowan et al. 1979). However, when the flight energy is transformed from carbohydrate to fat, the carbohydrate is not simply exhausted, and a certain carbohydrate reserve can provide energy for the sudden flight behavior of insects (Michel et al. 1986; Liu et al. 2007; Mayer et al. 1969). In the study, genes associated with gluconeogenesis (YWHAB_Q_Z) (Ji et al. 2021) and fatty acid metabolism (PCK, ACADM, and MDH2) (Schmidt et al. 1981; Zhang et al. 2015) showed high network interaction implied that these genes played an important role in flight energy supply in winged male adult of cotton aphid. In particular, genes related to fatty acid biosynthesis and metabolism, namely ACADM and MDH2, are crucial in long periods of flight for finding primary hosts. TSSK1 and CHDH are two genes associated with male reproduction (Sohail et al. 2019; Johnson et al. 2010). The high interactivity of these reproductive-related genes also corresponded to the previous reports that the morphological plasticity of insect was accompanied by reproductive plasticity in response to environmental changes (Tufail et al. 2010).

Conclusion

The study explored morphological and transcriptional dynamics of male cotton aphids in the development, and revealed the phenomenon of low feeding efficiency of winged male cotton aphids on summer hosts. Dynamic transcriptome analysis of male aphid at 5 different developmental periods showed that the 3rd instar nymphs were the critical stage of wing bud differentiation in males. Pathway enrichment analysis of DEGs and weighted gene co-expression network analysis revealed key signaling pathway and candidate hub genes that may be involved in wing differentiation in males of cotton aphid.

References

Download references

Acknowledgements

We thank the editors and reviewers for all their constructive comments on our work. This research was funded by National Natural Science Foundation of China (No. 32102214) and Special Fund for Basic Public Welfare Research of Institute of Cotton Research of CAAS (No. 1610162020020604). Studies were carried out in laboratories at the Institute of Cotton Research of CAAS, Anyang, China.

Author information

Authors and Affiliations

Authors

Contributions

Ji JC, Ma XY, Cui JJ, and Luo JY conceived and designed the research; Huangfu NB and Ji JC collected the samples. Ji JC, Shi QY, and Chen LL conducted experiments. Ji JC and Huangfu NB analyzed the data. Huangfu NB and Ji JC wrote the paper. All authors have read and approved the manuscript.

Corresponding authors

Correspondence to MA Xiaoyan or JI Jichao.

Ethics declarations

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.

Competing interests

All authors declare that they have no conflict of interest.

Supplementary Information

Additional file 1: Table S1

The designed primers for qRT-PCR.

Additional file 2: Table S2

The taxonomic statistics of male aphid.

Additional file 3: Table S3

The statistics of developmental period and morphological data of male aphid.

Additional file 4: Table S4

Summary statistics of RNA-seq results in developing A. gossypii.

Additional file 5: Table S5

Gene clusters and enriched pathways.

Additional file 6: Table S6

Enriched KEGG pathways in comparison between different developmental stages.

Additional file 7: Table S7

Gene information statistics of different modules of WGCNA cluster.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

HUANGFU, N., SHI, Q., CHEN, L. et al. Comparative transcriptional analysis and identification of hub genes associated with wing differentiation of male in Aphis gossypii. J Cotton Res 5, 22 (2022). https://doi.org/10.1186/s42397-022-00130-x

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s42397-022-00130-x

Keywords

  • Aphis gossypii
  • EPG
  • Wing bud differentiation
  • Dynamic transcriptome
  • WGCNA