Morphological characteristics, developmental dynamics, and gene temporal expressions across various development stages of Aphis gossypii sexual female

Abstract

Background

Aphis gossypii (Hemiptera: Aphididae) is a worldwide polyphagous phloem-feeding agricultural pest, and it can produce offspring by sexual or asexual reproduction. Compared with dozens of generations by parthenogenesis, sexual reproduction is performed in only one generation within one year, and little is known about the sexual reproduction of A. gossypii. In this study, sexual females of A. gossypii were successfully obtained through a previously established induction platform, and the morphological characteristics, developmental dynamics, and temporal gene expression were examined. Subsequently, signaling pathways potentially involved in regulating the growth, development, and reproduction of sexual females were investigated.

Results

The morphological observation showed that from the 1st instar nymph to adult, sexual females exhibited a gradually deepened body color, an enlarged body size, longer antennae with a blackened end, and obviously protruding cauda (in adulthood). The anatomy found that the ovaries of sexual females developed rapidly from the 2^nd instar nymph, and the embedded oocytes matured in adulthood. In addition, time-course transcriptome analysis revealed that gene expression profiles across the development of sexual females fell into 9 clusters with distinct patterns, in which gene expression levels in clusters 1, 5, and 8 peaked at the 2^nd instar nymphal stage with the largest number of up-regulated genes, suggesting that the 2^nd instar nymph was an important ovary development period. Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis revealed that a large number of genes in the sexual female adult were enriched in the TGF-beta signaling pathway and Forkhead box O (FoxO) signaling pathway, highlighting their important role in sexual female adult development and reproduction.

Conclusion

The morphological changes of the sexual female at each developmental stage were revealed for the first time. In addition, time-course transcriptomic analyses suggest genes enriched in the TGF-beta signaling pathway and FoxO signaling pathway probably contribute to regulating the development and oocyte maturation of sexual females. Overall, these findings will facilitate the regulating mechanism research in the growth and development of sexual females by providing candidate genes.

Introduction

Aphids are worldwide agricultural pests, harming multiple crops including cereal crops, Brassica crops, potatoes, cotton, vegetables, and fruits (van Emden et al., 2007). Most aphid species alternate parthenogenesis and sexual reproduction with varying seasons (Cuti et al., 2021). Parthenogenesis is characterized by its high reproductive efficiency and low consumption cost to produce considerable offsprings and can preserve stable genes. The harm of aphids to crops is mainly caused by parthenogenesis-induced rapid population expansion (Dixon, 1992).

Hence, considerable studies about growth and reproduction mechanisms in aphids mainly chose parthenogenetic females as the model other than sexual females. Kanturski et al. (2020) provided the first detailed morphological description of all body parts of parthenogenetic and sexual generations of adult pea aphids by scanning electron microscopy (SEM), which provides some reference points for quickly distinguishing between the two reproductive modes of pea aphids. In Acyrthosiphon pisum, silencing of the HAT gene increases the survival rate of offspring of parthenogenetic females, while silencing of the HDAC gene reduces the survival and delays the development of their offspring. The inactivation of the phenylalanine hydroxylase (PAH) enzyme affects fertility by reducing offspring amount and causes severe morphological defects in newborn nymphs (Simonet et al., 2016; Kirfel et al., 2020). For the brown citrus aphid (Aphis citricidus), the disruption of Vg and VgR expression in parthenogenetic females via RNAi delays the transition from nymphs to adults, prolongs the pre-reproductive period, but shortens the reproductive period, eventually decreases newborn nymph number (Shang et al., 2018). Knockout of SaEcR and SaUSP genes reduces the survival rate and fecundity of surviving aphids, demonstrating that SaEcR and SaUSP are important genes for the growth and development of parthenogenetic females in wheat aphid (Sitobion avenae), which can be used as RNAi targets for controlling this wheat aphid (Yan et al., 2016). In contrast, knowledge about sexual females especially in development and reproduction was few in aphids. No one has described the developmental process of sexual females fully so far. Liu et al. (2014) compared the morphology, ovary, and ovariole development of parthenogenetic females, gynopara, and sexual females in a comparative analysis of gene expression profiles in Aphis gossypii. Ji et al. (2023) only compared the antennae and body morphology of sexual females and males during the adult stage of Aphis gossypii, and did not refine the developmental process of sexual female aphids. Wieczorek et al. (2020) described the structure of the reproductive system of the sexual morphs form of Acyrthosiphon svalbardicum compared with Acyrthosiphon pisum and concluded that the unique characters include enormous fat body layer mainly in oviparous females, adhering to the structures of the internal reproductive system.

The life history of aphids mostly includes parthenogenesis and sexual reproduction, and sexual reproduction can enhance gene exchange among populations in different regions from different hosts so as to ensure the richness of population genetic structure (Tabata et al., 2016). At present, the research on the sexual female mainly focuses on morphological characterization, ovarian structure, and sex pheromone release (Miura et al., 2003; Kanturski et al., 2020; Koellner et al., 2022). The oviparous and viviparous adults of the pea aphid are almost morphologically identical (Liu et al., 2014). In A. pisum, the oviparous female (namely, sexual female) has a pair of ovaries, each of which contains seven ovarioles. In each ovary, there is usually only one oocyte developing at a time point, and the development of other oocytes is not synchronized (Wieczorek et al., 2019). A previous study has reported that round sex pheromone glands of A. pisum are located on the hind tibia of sexual females in adulthood, but not in parthenogenetic females. Moreover, genes encoding the downstream enzymes in the mevalonate pathway are highly expressed in the hind tibia of sexual females (Murano et al., 2018).

Few studies are available on the growth, development, and reproduction regulation of sexual females in aphid species. This is probably due to the fact that aphids only produce sexual females from late autumn to early winter with unknown birth sequences and uncertain mothers (Brisson et al., 2016; Lin et al., 2022). In Lachnus chosoni, the sexual generation shares the most features with the apterous vivipara and an aphid individual can only be confirmed as a sexual female when it produces eggs or be dissected (Kanturski et al., 2018). In addition, the sexual morphs induced indoors are still unstable and complex for aphid species. For example, in Schlechtendalia chinensis, the induction rate of sexual females substantially varies with temperature, with an induction rate of 89.6% at 7.5 °C and 5.3% at 18 °C, respectively, but at 22 °C, induction of sexual females failed (Liu et al., 2017). Campbell et al. (2005) found that transferring the winged Phorodon humuli adults reared under short-day conditions (12 h light/12 h dark) to long-day conditions (18 h light/6 h dark) results in the transformation from parthenogenesis to the generation of sexual female and male aphids, while the opposite conditions resulted in the production of a proportion of wingless parthenogenetic females, gynoparae, and males. Hong et al. (1998) raised Aphis spiraecola on Oenanthe javanica for three generations under the conditions of low temperature and short light to induce its sexual morphs, and found that males were produced in the second generation with an induction rate of 21%, and sexual females were produced in the third generation. In addition, Ji et al. (2023) have systematically studied the sexual morphs, feeding behaviors, and differentially expressed genes (DEG) between sexual females and males in adulthood. As mentioned above, though multiple sexual reproduction induction platforms have been established in various aphid species, the production and development processes of sexual females remain largely unknown, which results in the lack of a high-resolution picture of morphological characteristics and detailed description for sexual females, which prevents further study of sexual females in aphids. In contrast, we have established a stable and standardized sexual female induction platform and assembled a high-quality genome of cotton aphids at the chromosomal level, making A. gossypii an ideal model for studying the sexual female in aphids (Chen et al., 2022; Ji et al., 2021; Ji et al., 2023; Zhang et al., 2022).

In this study, we first retrieved and traced oviparous and sexual female aphid species observed in the field according to related records and descriptions in the book Aphids on the World’s Herbaceous Plants and Shrubs (Volume 2) (Blackman et al., 2006). Moreover, we investigated morphological changes, developmental dynamics, adult fecundity, and gene temporal expression of sexual females of A. gossypii in order to reveal the molecular regulation mechanism underlying their growth, development, and reproduction.

Methods and materials

Sexual female in aphids

The aphid species of egg-laying morphs reported were retrieved from the book Aphids on the World’s Herbaceous Plants and Shrubs (Volume 2) by keywords “oviparous female”, “ovipara”, “oviparae”, and “sexual female” (Blackman et al., 2006). Then the retrieval results are listed and categorized at the taxonomic levels.

Induction and feeding of sexual females of A. gossypii

A. gossypii was reared and purified for multiple generations in the laboratory of Institute of Cotton Research, Chinese Academy of Agricultural Sciences (Anyang City, Henan Province, China). The cotton seedlings (CCRI 49) were grown at (25 ± 1) ℃ with a relative humidity of 75% and a photoperiod of 14 h light/10 h dark. The sexual females were induced under low temperature and short-light conditions (18 ℃, 75% of relative humidity, and 8 h light/16 h dark photoperiod) according to the previously reported method (Chen et al., 2022; Ji et al., 2021, 2023).

Morphological characteristics and fecundity

The newborn sexual female within 12 h was transferred to a petri dish containing cotton leaves and 1.8% (mass fraction) agarose medium for the single culture. The instar of the sexual female was recorded according to molting times every 12 h, and a total of 30 sexual females were counted, while the body length, body width, and antennal length of those 30 sexual females were photographed and measured using a SteREO Discovery V8 microscope (Zeiss, Oberkochen, Germany) until adulthood. Meanwhile, the ovaries of sexual females at each development stage were isolated and observed in phosphate-buffered saline (PBS) buffer in parallel experiments. Eventually, the ovaries at the 2^nd, 3^rd, 4^th instar nymphs and adult sexual females were displayed except for the 1^st instar nymph due to the absence of obvious ovarian structure.

Preparation of RNA sequencing samples

The 1–4 instar nymphs (recorded as SF1, SF2, SF3, SF4, respectively) and adults (recorded as SFA) of sexual females of A. gossypii were collected with at least three replicates at each developmental stage with 50 individuals per replicate. Total RNA was extracted by Trizol^® reagent (Promega, Madison, WI, USA) according to the manufacturer’s instructions, and RNA quality was determined using spectrophotometer Nanodrop 2000 (Thermo, Wilmington, DE, USA) with 1% (mass fraction) agarose gel electrophoresis. The cDNA library construction and high-throughput sequencing were performed by Biomarker Technologies Company (Qingdao City, China).

Transcriptome assembly and gene annotation

After RNA-seq, the adaptor sequences and low-quality reads were removed to obtain clean reads. These clean reads were aligned against the published genome of A. gossypii (version ASM2018417v2) using HISAT2 (Kim et al., 2015; Zhang et al., 2022). The transcripts of each sample were reconstructed using String Tie and aligned with the original genome annotation information (Pertea et al., 2015). Then, the new genes were aligned against the NR (Non-Redundant Database), Swiss-Prot, COG (Cluster of Orthologous Groups), KOG (Clusters of orthologous groups for eukaryotic complete genomes), KEGG databases (Kyoto Encyclopedia of Genes and Genomes), and the functional annotation was performed using the Inter Pro integrated database. Using HMMER software, the predicted amino acid sequences of new genes were aligned against the Pfam database to obtain the annotation information (Eddy, 1998). All data were deposited in the NCBI SRA (sequence read archive) database with the accession number PRJNA1109341.

Identification and functional analysis of differentially expressed genes (DEGs)

The DEGs were identified in the pair-wise comparison between two neighboring development stages (SF1 vs. SF2, SF2 vs. SF3, SF3 vs. SF4, SF4 vs. SFA) using DESeq2 software with the thresholds of |log ₂ FC (fold change)|≥ 2 and false discovery rate (FDR) < 0.05 (Love et al., 2014). Using String Tie software, the gene expression level was calculated by the maximum flow algorithm and normalized as fragments per kilobase million (FPKM) (Trapnell et al., 2010; Pertea et al., 2015). In order to analyze the biological function of DEGs, KEGG pathway enrichment analysis was performed on DEGs with Q-value < 0.05 as the screening standard.

Gene temporal expression analysis

A fuzzy c-means algorithm was implemented to cluster the genes with similar expression patterns in various developmental stages of sexual females through Mfuzz software (Kumar et al., 2007). Then, the functions of genes categorized into those distinct clusters were predicted via KEGG pathways enrichment analysis as described above.

Validation of RNA-Seq data by real-time quantitative polymerase chain reaction (RT-qPCR)

In order to verify the accuracy of the transcriptome results, the relative expression of randomly selected 9 genes (PUF60, TYR, RNASEH2C, GRB2, G10, galE, ELOVL7, E1.10.3.3, and ABCC1) was quantified by RT-qPCR utilizing a Light Cycler 480 machine (Roche Diagnostics, Switzerland). The 20 μL RT-qPCR reaction system included 2 μL cDNA template, 7.2 μL nuclease-free water, 0.4 μL forward primer, 0.4 μL reverse primer, 10 μL 2 × TransStart® Top Green qPCR SuperMix (+ DyeI/ + DyeII) AQ131 (TransGen Biotech, Bejing, China). The RT-qPCR procedures were as follows: pre-denaturation 95 ℃ 5 min, followed by 40 cycles of 95 ℃ 5 s, 60 ℃ 15 s, and 72 ℃ 10 s. The melting curve was constructed to verify the specificity of amplification. Gene‐specific primers were designed using Primer Premier 6 and synthesized by Shanghai Sangon Biotech Co, LTD (Shanghai, China). In addition, the GAPDH gene was used as the internal control for gene expression level normalization (Ma et al., 2016). All primers used in this study are listed in Table 1. The relative expression level of each gene was calculated by the 2^−ΔΔC_t method (Livak et al., 2001).

Table 1 Transcriptome quantitative validation primers

Results

An overview of sex females in Aphididae

Sexual females were obtained in 1 096 Aphididae races, belonging to 25 subfamilies, 31 genera, and 113 races (Fig. 1, Table S1). At the subfamily level, sexual females were recorded mainly in Aphidinae (737), Calaphidinae (199), Lachninae (151), and Chaitophorinae (39) (Fig. 1 A). At the genus level, the most sexual females were found in Macrosiphini (585), followed by Aphidini (167), Eulachnini (126), and Calaphidini (49) (Fig. 1 B). At the race level, sexual females were observed primarily in Aphis (118), Cinara (117), Uroleucon (56), and Macrosiphum (55) (Fig. 1 C).

Fig. 1

Distribution of sexual females of A. gossypii at subfamily (A), genus (B), and race (C) levels. Only the top 20 taxa are shown

Morphology and fecundity of sexual females of A. gossypii

From the 1^st instar nymph to the adult, the body color of sexual females was gradually deepened from light green to dark green (Fig. 2 A), and the body size was gradually increased with the development, of which the mean body length at each development stage (1^st–4^th instar nymphs and adult) was 0.46 mm, 0.70 mm, 0.81 mm, 1.03 mm, and 1.26 mm, respectively, and the mean body width at each development stage was 0.23 mm, 0.35 mm, 0.41 mm, 0.56 mm, and 0.67 mm, respectively (Fig. 2 B, C). In addition, the mean antennal length of sexual females throughout the development period from the 1^st instar nymph to the adult was also gradually increased, which was 0.25 mm, 0.32 mm, 0.38 mm, 0.48 mm, and 0.59 mm, respectively (Fig. 2 D).

Fig. 2

Morphological characteristics dynamics of A. gossypii sexual females across the development. A Morphology of sexual females at five development stages. SF1 to SF4 represent the 1^st instar to 4^th instar nymphs of sexual females, respectively. SFA indicates sexual female adults. The scale bar is 200 μm. B–E Body length, body width, antenna length, and developmental duration of sexual females at 5 development stages were statistically analyzed by the one-way analysis of variance (ANOVA) to reveal the significance of differences between  each stage except for adulthood. The meanings of different letters (a, b, c, d, e) indicate significant differences. P-OP, pre-oviposition period; OP, oviposition period. F–I Ovaries of 2^nd, 3^rd, 4^th instar nymphs, and adult sexual females, respectively. The data were presented as mean ± standard errors of thirty biological replicates

Moreover, the mean duration time of each nymph stage (SF1–SF4) lasted for 2.37 d, 2.44 d, 2.13 d, and 2.68 d, respectively (Fig. 2 E), and the total nymph stage was 9.62 d. Interestingly, the mean adult life span was 21.21 d, including a longer pre-oviposition period of 14.95 d following a shorter oviposition period of 6.26 d (Table S2). Under laboratory conditions in this study, the average number of eggs laid by each sexual female without mating with male aphids was 3. Only trophoblast cells were observed in the whole ovaries of the 2^nd instar sexual female nymph (Fig. 2 F). In the ovaries of the 3^rd instar sexual female nymph, the oocytes became larger, but the trophoblast cells became smaller (Fig. 2 G). In the ovaries of the 4^th instar sexual female nymph, the oocytes took on a yellow color due to the deposition of vitellogenin (Fig. 2 H). Eventually, the ovaries of the sexual females were filled with mature oval oocytes in adulthood (Fig. 2 I).

Transcriptome assembly and gene annotation of sexual female

A total of 89.84 Gb clean reads were obtained from RNA sequencing of 15 samples, and the percentage of Q₃₀ bases was ≥ 92.78%. The clean reads of each sample were aligned to the specific reference genome (version ASM2018417v2), with an alignment rate ranging from 93.94% to 96.06% (Table S3). A total of 1 763 new genes were identified after blasting against the public databases (NR, Swiss-Prot, COG, KOG, and KEGG), of which 940 new genes were annotated eventually. Principal component analysis (PCA) and Pearson correlation analysis showed that the samples at different development stages exhibited obvious separation, but the samples from different replicates at the same development stage were close to each other (Fig. 3 A). Moreover, the Pearson correlation coefficient of gene expression levels at each development stage was greater than 0.9 (Fig. 3 B), indicating that our RNA-seq data were reliable and qualified for subsequent bioinformatics analysis.

Fig. 3

Principal component analysis (PCA) and Pearson correlation analysis heatmap of RNA-seq data. A PCA of samples at 5 development stages. SF1-1, SF1-2, and SF1-3, respectively, represent 3 biological replicates of SF1, and so on. Samples were clustered into 5 groups, respectively. SF1, SF2, SF3, SF4, and SFA represent the 1^st, 2^nd, 3^rd, 4^th instar nymphs and adult sexual females, respectively. B Pearson correlation analysis heatmap of 15 samples. In the heatmap, dark red represents a strong correlation, and blue denotes a weak correlation. Each column and row corresponds to the correlation between one sample and the other 14 samples (including itself)

Time-course transcriptome profiling across the development of sexual female

A total of 16 616 expressed genes across the development of sexual females fell into 9 different temporal pattern clusters, representing different expression dynamics, in which the values of zero expression were removed (Fig. 4 and Table S4). Cluster 4 contained the largest number of genes (2 473 genes), followed by cluster 6 (2 034 genes), cluster 3 (1 822 genes), cluster 8 (1 777 genes), cluster 2 (1 573 genes), cluster 9 (1 559 genes), cluster 1 (1 409 genes), cluster 5 (1 354 genes), and cluster 7 (826 genes). Moreover, the genes in cluster 3 and 4 were up-regulated during the development of sexual females, whereas the genes in clusters 6 and 9 were down-regulated significantly. In addition, most genes in clusters 1, 5, and 8 exhibited their highest expression level at the 2^nd instar nymph stage, while most genes in clusters 3, 4, and 5 displayed up-regulated expressions in the adult stage.

Fig. 4

Gene expression profiles across the development of sexual females. A total of 16 616 genes expressed in sexual females of cotton aphids fell into 9 distinct clusters of temporal expression patterns

KEGG pathways enrichment was performed to reveal the potential function of the genes in the above-mentioned 9 clusters (Fig. 5 and Table S5). The genes in cluster 1 were mostly enriched in pathways such as translation (ribosome) and sensory system (phototransduction-fly). The genes in cluster 2 were mostly enriched in the pathways such as energy metabolism (oxidative phosphorylation), lipid metabolism (fatty acid degradation, fatty acid elongation, glycerolipid metabolism, and others), carbohydrate metabolism (citrate cycle (TCA cycle), pentose and glucuronate interconversions, and ascorbate and aldarate metabolism), metabolism of cofactors and vitamins (retinol metabolism and folate biosynthesis), and amino acid metabolism (phenylalanine metabolism and tyrosine metabolism). The genes in cluster 3 were mainly enriched in pathways of RNA transport, replication, repair (DNA replication, mismatch repair, and homologous recombination), folding, sorting, and degradation (proteasome). The genes in cluster 4 included several vital signal transduction pathways (MAPK signaling pathway-fly, Wnt signaling pathway, Forkhead box O (FoxO) signaling pathway, Hippo signaling pathway-fly, and TGF-beta signaling pathway), transport and catabolism (endocytosis), environmental adaptation (circadian rhythm-fly), folding, sorting, and degradation (ubiquitin-mediated proteolysis), development and regeneration (dorso-ventral axis formation). The genes in cluster 5 were mainly involved in pathways of DNA replication, mismatch repair, base excision repair, and ribosome biogenesis in eukaryotes. The genes in clusters 6, 7, 8, and 9 were mainly enriched in signaling pathways related to neuroactive ligand-receptor interaction, replication and repair (nucleotide excision repair and Fanconi anemia pathway), and spliceosome and translation (aminoacyl-tRNA biosynthesis and ribosome). To verify the reliability of RNA-seq data, a total of 9 genes were randomly selected from 9 clusters with one gene per cluster, and their expression levels were validated by RT-qPCR. The results showed that the gene expression levels of these 9 genes were generally consistent with the results of RNA-seq analysis (Fig. 6), confirming the reliability of the RNA-seq data.

Fig. 5

KEGG pathway enrichment analysis of genes in 9 distinct gene expression clusters. The proportion of genes enriched in the target pathway indicates an enrichment degree. The color of the bubble from red to green indicates a gradually decreased P-value. The size of the bubble indicates the number of genes enriched in the pathway, of which the larger the bubble, the more the genes are enriched

Fig. 6

Verification of RNA seq results of 9 genes randomly selected from Cluster 1 to 9 by RT-qPCR. The gene expression level at the 1^st instar nymph stage was used as a control. Pearson correlation coefficients were calculated to compare results from RT-qPCR and RNA-seq data for each gene

Functional analysis of DEGs across the development of sexual female

A total of 2 511 DEGs were identified in the comparison of the first vs. the second instar sexual female nymphs, of which 1 287 and 1 224 genes were significantly up-regulated and down-regulated, respectively. A total of 1 643 DEGs were identified in the 2^nd vs. the 3^rd instar sexual female nymphs, of which 785 and 858 genes were significantly up-regulated and down-regulated, respectively. Out of 267 significant DEGs in the 3^rd vs. the 4^th instar sexual females, 99 genes were up-regulated, and 168 genes were down-regulated. Of 2 248 DEGs in the 4^th instars vs. adult sexual females, 663 and 1 585 genes were significantly up-regulated and down-regulated, respectively (Fig. 7 A). The Venn diagram showed that 14 DEGs were shared in all above 4 comparison groups, whereas 1 332, 672, 58, 1 001 DEGs were unique to the pairwise comparison of SF1 vs. SF2, SF2 vs. SF3, SF3 vs. SF4, SF4 vs. SFA, respectively (Fig. 7 B).

Fig. 7

Venn diagram of DEGs across the sexual female development process. A Number of DEGs in the pairwise comparison of SF1 vs. SF2, SF2 vs. SF3, SF3 vs. SF4, and SF4 vs. SFA. B Distribution of DEGs in the pairwise comparison of SF1 vs. SF2, SF2 vs. SF3, SF3 vs. SF4, and SF4 vs. SFA. SF1, SF2, SF3, SF4, and SFA represent the 1^st, 2^nd, 3^rd, 4^th instar nymphs and adult sexual females, respectively

Pathway enrichment analysis was performed to explore the roles of the up-regulated and down-regulated DEGs in the above four pairwise comparison groups (Fig. 8 and Table S6). In SF1 vs. SF2, upregulated DEGs were significantly enriched in replication and repair (DNA replication and homologous recombination), fatty acid degradation, and RNA transport pathways, whereas the downregulated DEGs were mostly enriched in lysosome, neuroactive ligand-receptor interaction, lipid metabolism (biosynthesis of unsaturated fatty acids and fatty acid elongation), various types of N-glycan biosynthesis, and carbohydrate metabolism (starch and sucrose metabolism, and galactose metabolism). In SF2 vs. SF3, the up-regulated DEGs were mainly significantly enriched in transport and catabolism (lysosome and autophagy-animal), lipid metabolism (fatty acid degradation and glycerolipid metabolism), retinol metabolism, and drug metabolism-cytochrome P450, whereas the down-regulated DEGs were mostly significantly enriched in neuroactive ligand-receptor interaction, Toll and Imd signaling pathway, and ABC transporters pathways. In SF3 vs. SF4, the majority of up-regulated DEGs were enriched in lipid metabolism (glycerolipid metabolism and fatty acid degradation), carbohydrate metabolism (pentose and glucuronate interconversions and ascorbate and aldarate metabolism), retinol metabolism pathways. However, there was no pathway significantly enriched with downregulated DEGs. In SF4 vs. SFA, the upregulated DEGs were mainly enriched in signal transduction pathways (MAPK signaling pathway-fly, TGF-beta, and FoxO signaling pathway), apoptosis-fly, and dorso-ventral axis formation, whereas the downregulated DEGs were primarily enriched in lysosome, lipid, amino acid, carbohydrate, and other synthetic and metabolic pathways.

Fig. 8

KEGG pathway enrichment analysis of significantly up-regulated and down-regulated DEGs in pairwise comparison of SF1 vs. SF2, SF2 vs. SF3, SF3 vs. SF4, and SF4 vs. SFA. SF1, SF2, SF3, SF4, and SFA represent the 1^st, 2^nd, 3^rd, 4^th instar nymphs and adult sexual females, respectively. Significantly enriched pathways were determined by the hypergeometric test at adjusted P < 0.05

Discussion

Dynamics of morphological characteristics across the development of A. gossypii sexual female

Although 1 094 aphid sexual females (also known as oviparous females or ovipara) have been observed and recorded, few studies have been conducted to characterize their morphology, growth, development, or reproduction. In this study, we investigated the morphological changes of sexual females of A. gossypii at 5 different developmental stages, and found that their body size enlarged gradually after each molting and the body color changed from light green to dark green. This conclusion is consistent with the results of Ji et al. (2023). The external morphological characteristics of sexual females were similar to those of parthenogenetic females except for the body color and size, which were also validated in this study (Liu et al., 2014; Ji et al., 2021). Büning (1985) found that in the ovaries of oviparous aphids, the nurse cells enlarged by endomitosis (n = 2⁸n–2¹⁰n), whereas in parthenogenetic aphids, the nurse cell nuclei remain small (n = 2²n–2⁴n) and that in parthenogenetic aphids, the previtellogenic growth of oocytes is highly inhibited, and vitellogenesis and chorionogenesis are blocked totally. However, in this study, we observed the ovarian development and oocyte maturation across the development of sexual females.

Signaling pathways related to ovarian development of sexual female

The DEG number of sexual females showed a “V” trend during the development process, namely, the number of up-regulated and down-regulated DEGs decreased with the development of nymphs, and both dropped to the lowest in the comparison of SF3 vs. SF4, then increased sharply in SF4 vs. SFA. In SF1 vs. SF2, the number of up-regulated DEGs was the largest (Fig. 7 A). These up-regulated DEGs were mainly enriched in fatty acid degradation, DNA replication, homologous recombination, base excision repair, mismatch repair, and RNA transport pathways (Fig. 8 and Table S6). In addition, gene expression profile analysis showed that the genes in clusters 1, 4, and 8 had high expression levels in SF2, and these genes were mainly involved in environmental information processing, genetic information processing, organismal systems, and cellular processes (Fig. 5 and Table S5). These results were supported by our dissection analysis results that the 2^nd instar nymph period was a critical period for the start of ovarian development of sexual females (Fig. 2). The above results indicated that these pathways probably jointly promoted the development of the ovaries of sexual females. However, related mechanisms remain to be unexplored.

It was worth noting that in SF3 vs. SF4, the up-regulated DEGs were mainly enriched in the pathways of terpenoid backbone biosynthesis, longevity regulating pathway-multiple species, and retinol metabolism. Sun et al. (2018) found that the RNAi of the FPPS gene responsible for the terpenoid backbone biosynthesis in aphids can significantly reduce the body size and fecundity of A. gossypii. Wang et al. (2022) have reported that the retinol metabolism pathway is significantly associated with the maturation of the Exopalaemon carinicauda ovaries. Combining these previous reports and our observation that these two pathways (terpenoid backbone biosynthesis and retinol metabolism) were enriched with up-regulated DEGs, we concluded that these two pathways might be related to the ovary development of sexual females in A. gossypii.

The TGF-beta signaling pathway has been reported to be associated with the embryonic development of Tribolium castaneum (Gao et al., 2023). FoxO signaling pathway is responsible for large-scale Vg synthesis required for simultaneous maturation of multiple eggs in oviparous insects (Wu et al., 2020). In this study, we also found that a large number of up-regulated DEGs in SF4 vs. SFA were mainly enriched in the TGF-beta signaling pathway and FoxO signaling pathway (Fig. 8). The above reports and our findings collectively indicated that these two pathways significantly enriched at the adult stage of the sexual female aphid might be related to the embryonic development and egg maturation of sexual female, which needs to be further confirmed by more experimental evidence in future.

Conclusion

In this study, the morphological characteristics, development dynamics, ovarian maturation process, and gene temporal expression of sexual females were investigated, with cotton aphid A. gossypii as a model insect. Through time-course transcriptomic analysis, it was found that the changes in differentially expressed genes were consistent with the changes in external morphology and ovary development of sexual females. Subsequent pathway enrichment analysis provided several vital signaling pathways potentially involved in regulating the development and ovary maturation such as the TGF-beta signaling pathway and FoxO signaling pathway in sexual females, of which roles should be confirmed in the future.