DOI:10.32604/biocell.2022.016563
| BIOCELL DOI:10.32604/biocell.2022.016563 | |
| Article |
Transcriptome profile analysis of the accompanying migratory parasitic wasp Aenasius bambawalei (= Aenasius arizonensis girault) (Hymenoptera: Encyrtidae): Genes related to fertilization involved at different stage of ovary development
1Guangxi Key Laboratory of Biology for Crop Diseases and Insect Pests, Guangxi Academy of Agricultural Sciences, Nanning, 530007, China
2Zhejiang Institute of Landscape Plants and Flowers, Zhejiang Xiaoshan Institute of Cotton & Bast Fiber Crops, Hangzhou, 311251, China
3State Key Laboratory for Managing Biotic and Chemical Threats to the Quality and Safety of Agro-Products, Institute of Plant Protection and Microbiology, Zhejiang Academy of Agricultural Sciences, Hangzhou, 310021, China
*Address correspondence to: Jun Huang, junhuang1981@aliyun.com
Received: 23 March 2021; Accepted: 06 May 2021
Abstract: Age-related declines in fertilization success have been reported for a wide range of species. The fertilization of parasitic wasps is closely related to egg production and sperm storage. Aenasius bambawalei (Hymenoptera: Encyridae) is a key parasitic wasp of the important invasive mealybug Phenacoccus solenopsis (Hemiptera: Pseudococcidae). The female offspring ratio of this parasitic wasp was declined with parental age in mass rearing under laboratory conditions. To investigate the molecular mechanisms involved in the reproduction of A. bambawalei, an extensive analysis of the impact of age on transcriptome profile of mated ovaries of this wasp was performed by comparing the gene expression profiles of various maternal ages: the early stage (ES), the intermediate stage (IS), and the advanced stage (AS). In total, 358 differentially expressed genes were identified, with 17.60% (63 genes) of the changes associated with greater expression in fertilization. Moreover, the expression of serine protease 47 precursor, serine protease inhibitor 3/4, glucose dehydrogenase, fatty acyl-CoA reductase 1-like, major royal jelly, and acyl-CoA delta (11) desaturase-like was detected by real-time fluorescence quantitative (RT-qPCR). The results showed that fertilization related genes exhibited a stage-specific pattern. Egg production and sperm storage genes in A. bambawalei were significantly modified in the transcriptome, providing a starting point for the genetic dissection of fertilization.
Keywords: Phenacoccus solenopsis; Egg production; Sperm storage; Transcriptomics; Fertilization; Aenasius bambawalei
Age-related declines in fertilization success have been reported for a wide range of species, such as flies (Bonduriansky and Brassil, 2002), beetles (Jehan et al., 2021) and humans (Kidd et al., 2001; Lane et al., 2014). In Hymenoptera, males are produced from haploid oocytes and only female offspring emerge from diploid eggs (Damiens et al., 2003; Bressac et al., 2008; Heimpel and de Boer, 2008; Nguyen et al., 2013). Consequently, the fertilization of oocytes is a key regulator of the sex ratio (Hardy, 1994; Nguyen et al., 2013), which is very important for mass rearing and utilization of parasitic wasps (Hardy, 1994). However, few mechanisms responsible for these phenomena have yet been established (Jones and Elgar, 2004).
The fertilization of oocytes mainly depends on two physiological parameters: the pattern of egg production (Jervis and Kidd, 1986; Rivero-Lynch and Godfray, 1997) and the dynamics of sperm storage (Cardozo et al., 2020; Matsuzaki and Sasanami, 2017; Chevrier and Bressac, 2002). For parasitic wasps, there are two strategies of egg production: species that have all their egg complement mature upon emergence are termed proovigenic, while those that mature eggs during the adult stage are termed synovigenic (Flanders, 1950). Therefore, for synovigenic eggs, the number of mature eggs varies according to adult age (Espinosa and Virla, 2018). The dynamics of sperm storage significantly influence the fitness of parasitic wasps (Cardozo et al., 2020). However, the effects of maternal age on sperm storage and eventually on fertilization are ambiguous (Damiens et al., 2003).
The mealybug Phenacoccus solenopsis Tinsley (Hemiptera: Pseudococcidae) is a highly invasive and polyphagous pest (Arya et al., 2020; Fand and Suroshe, 2015; Fand et al., 2011). Aenasius bambawalei Hayat (= Aenasius arizonensis Girault) (Hymenoptera: Encyrtidae) is an accompanying parasitoid wasp of this mealybug (Ahmed et al., 2015; Zhang et al., 2016). This parasitic wasp is an ideal biocontrol tool because of its host-specificity, fast reproduction compared to the host, and female-biased sex ratio (Fand et al., 2011). Most importantly, it plays a significant role in keeping mealybug populations under control (Ashfaq et al., 2010; Hayat, 2009; Wang et al., 2010). We observed that Aenasius females were synovigenic in our pilot study. In addition, the female offspring sex ratio of this parasitic wasp declined with parental age in mass rearing. This phenomenon is not propitious for efficient use of this parasitic wasp. A previous study demonstrated that the use of deep RNA-sequencing technology on dissected ovaries is a powerful method to investigate the molecular mechanisms involved in hymenopteran reproduction (Niu et al., 2014). However, an extensive analysis of the impact of age on the transcriptome profile of mated ovaries of A. bambawalei was needed.
Here, we aimed to identify the genes that control egg production and sperm storage. First, we explored the ovary gene expression changes in A. bambawalei females of three ages: early (ES), intermediate (IS) and advanced (AS) stages. Then, we compared the gene expression profiles of different maternal ages and identified the genes that might be related to egg production and sperm storage. The results will help researchers to understand the underlying mechanisms associated with fertilization of this parasitic wasp, design better biological control strategies, and develop monitoring tools for parasitism rates in the field (Calla et al., 2015).
The ovary transcriptomes of mated A. bambawalei at three female ages were analyzed to measure the gene expression changes associated with egg production and sperm storage.
Insect rearing and sample preparation
We conducted studies in areas where P. solenopsis and A. bambawalei occur, and no specific permits were required. The invasive mealybug P. solenopsis (Hemiptera: Pseudococcidae) has become a highly polyphagous pest and has been reported to damage >200 plant species from approximately 24 countries in tropical and subtropical regions of the world (Fand and Suroshe, 2015). The initial P. solenopsis used in this experiment was from a common green plant, Torenia fournieri (Scrophulariales: Scrophulariaceae), in the suburbs of Hangzhou (30°09’14”N; 120°13’30”E), Zhejiang Province, China in May 2011. Mealybugs have been maintained on their original host, Gossypium hirsutum (Malvales: Malvaceae) under laboratory conditions ever since. A. bambawalei (= Aenasius arizonensis Girault) (Hymenoptera: Encyrtidae) is an accompanying parasitic wasp of P. solenopsis (Ahmed et al., 2015; Zhang et al., 2016). The wasp colony was acquired from the above P. solenopsis colony and reared on colonies of P. solenopsis in a laboratory rearing system at the Zhejiang Academy of Agricultural Sciences. Female adults of P. solenopsis (the most suitable host stage for A. bambawalei) and A. bambawalei adults of different ages were used in this study.
The rearing approaches of A. bambawalei were similar to those described by Huang et al. (2017) and Zhang et al. (2016). Briefly, a leaf containing 30 female adults of P. solenopsis was first placed in the rearing system. Then, six pairs of newly emerged (1–2 days) and mated adults of A. bambawalei (Zhang et al., 2016) were introduced into the rearing system. Finally, the rearing system was reversed on a cup of water to prevent leaves from becoming desiccated. To ensure that each studied female was inseminated and by a similar number of males, we set up the following steps: in a 10-mL plastic finger tube, one female wasp and one male wasp were placed without mating, and at the same time, an absorbent cotton was dipped in 10% honey water and provided to supplement the nutrition of parasitic wasps. The success of mating can be judged by observing the mating behavior of parasitoids (King and Kuban, 2012). Given that nutrition supplementation can prolong the longevity of female wasps (Harvey et al., 2012), a 10% honeydew solution was supplied daily for A. bambawalei. The laboratory cultures of A. bambawalei and P. solenopsis were maintained in a growth chamber (27 ± 1°C, 14L:10D, and 75 ± 5% relative humidity (RH)).
Because the signal for male or female development appears to be already present in the egg upon oviposition, ovaries of three different female ages were chosen as the experimental materials. The ovaries of A. bambawalei were classified according to Jeanne’s classification method (Jeanne, 1972) with some modifications. The ovarian development process was divided into three stages: an early stage (ES, without or only partially visible developed oocytes, Fig. 1A), an intermediate stage (IS, containing two or more mature oocytes, Fig. 1B), and an advanced stage (AS, containing mature oocytes and apoptotic cells, Fig. 1C).
Figure 1: Schematic plot of Aenasius bambawalei ovaries at different ages. (A) Ovary from the 1st d, ES (early stage, without or only partially visible developed oocytes), (B) Ovary from the 3rd d, IS (intermediate stage, containing two or more mature oocytes), (C) Ovary from the 7th d, AS (advanced stage, containing mature oocytes and apoptotic cells). ov-ovary, do-developing oocytes, s-spermatheca, mo-mature oocyte, ac-apoptotic cell.
Then, the female adults at ES (the 1st d), IS (the 3rd d), and AS (the 7th d) were dissected under a microscope to obtain the ovaries. For the transcriptomics experiment, we used a total of 3 treatment combinations with n = 300 individuals per treatment combination (total n = 900 experimental females). These ovary combinations were placed into 1.5-mL microcentrifuge tubes, flash-frozen in liquid nitrogen, and stored at −80°C before RNA extraction.
RNA extraction, library preparation and sequencing
Total RNA from the A. bambawalei ovary was extracted using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer’s standard protocol. RNA integrity and quality were assessed using a bioanalyzer (Agilent 2100 Technologies). The samples for transcriptome analysis were prepared with an Illumina kit according to the manufacturer’s recommendations. Briefly, mRNA was purified from 3 µg of total RNA with oligo (dT) magnetic beads. After purification, the mRNA was fragmented, and the cleaved RNA fragmentswere used to synthesize first-strand cDNA with random hexamers and SuperScript II reverse transcriptase (Invitrogen). Then, second-strand cDNA was synthesized with DNA polymerase I and RNase H. These cDNA fragments were subjected to an end repair process and adapter ligation. Fragments 300–350 bp in length were excised and then enriched for 12 PCR cycles to generate the final cDNA library on the HiSeq 2000 Platform (Illumina) of CapitalBio Corporation, Beijing, China.
Transcriptome data processing and assembly
The raw images generated by sequences were converted into nucleotide sequences (called raw reads) using the Illumina GA Pipeline 1.6. Clean reads were obtained by removing adaptor reads, empty reads, and low-quality sequences (reads with unknown sequences ‘N’), and the threshold of length cutoff was more than 30 bp. Then, the high-quality reads were de novo assembled using Velvet 1.2.103 software to construct unique consensus sequences (Zerbino and Birney, 2008), and the kmer length was 57 bp. Finally, the completeness of the transcriptome was evaluated by Benchmarking Universal Single-Copy Orthologs (BUSCO) (Simão et al., 2015).
Gene expression profile establishment
The trimmed Solexa transcriptome reads were mapped onto the unique consensus sequences using Bowtie (Langmead et al., 2009). A Perl script (Toolbox for Biologists, TBtools v0.6673) was used to process the mapping results and to generate the gene expression profile.
Functional annotation and classification
Gene annotation was performed through BLASTn and BLASTx (E values ≤ E−5) (Altschul et al., 1997) searches against the Swiss-Prot (Altschul and Gish, 1996) and GenBank databases (E values ≤ E−10) and the best matched sequence was selected based on the gene annotation. GetORF from EMBOSS software was used to analyze reliable coding sequences (CDSs). Unigenes were assigned functional annotations by sequence similarity comparisons against the Clusters of Orthologous Groups of proteins database (Tatusov et al., 1997; Tatusov et al., 2003) with BLAST (Behrouzian and Buist, 2002) (E values ≤ E−10). The metabolic pathway was constructed based on the KEGG database using BLASTX (E values ≤ E−10) (Kanehisa et al., 2006).
Interpro domains (Mulder et al., 2003) were annotated by InterProScan Release 4.8 (Zdobnov and Apweiler, 2001), and functional assignments were mapped onto Gene Ontology (GO) (Harris et al., 2004). WEGO (Ye et al., 2006) was used to perform GO classification and to draw a GO tree. The datasets are available at the NCBI Short Read Archive (SRA) with the accession number: PRJNA273425.
Detection of differentially expressed unigenes and qPCR analysis
Tests for differential expression were carried out using DEGseq. The value of reads per kilobase per million reads (RPKM) was calculated for gene expression analysis. DEG-seq (Altschul and Gish, 1996) was implemented to identify differentially regulated genes between two samples using the two unpaired MA-plot-based methods to detect and visualize the gene expression differences where significant P-values were less than 0.01 or the expression ratio was more than 2 pairwise comparisons. In addition, the genes exhibiting differential expression that have a role in fertilization were detected by referring to previous studies (Pannebakker et al., 2013; Cook et al., 2015).
RNA extraction and integrity and quality assessment were performed according to the transcriptome sequencing procedure. Then, 1 µg of RNA was reverse transcribed to cDNA by SuperScript IV (Invitrogen, USA). Finally, qPCR was conducted in a 25-µL reaction mixture (12.5 µL of 2 × SYBR qPCR mixture, 1 µL of DNA Template, 0.5 µL of forward and reverse primers and ddH2O to the final volume) following the protocol for the SYBR Green Supermix Kit (Takara, Japan) on a BIO-RAD CFX96TM Real-Time System (Bustin et al., 2009). The reference gene was elongation factor 1α (Ef1α). Primers for the serine protease 47 precursor, glucose dehydrogenase, fatty acyl-CoA reductase 1-like, serine protease inhabitor 3/4, major royal jelly protein-like 7 precursor, and acyl-CoA delta (11) desaturase-like were designed for this study using primer-BLAST (Tab. 1). The qPCRs were set up as follows: 95°C for 3 min, 40 cycles of 95°C for 15 s and 61°C for 45 s. At the end of each run, melt curves with a gradual increase in temperature were performed (72°C for 10 s, 61°C for 5 s and 95°C for 5 s) to check for non-specific amplification. Dissociation curves were checked following the PCR to ensure the efficiency of PCR is 90% to 110%. We calculated the relative expression levels following the 2-ΔΔCt method. Error bars represent the means ± standard deviations from 3 biological replicates. The expression profiles of various genes at different ages were determined by one-way analysis of variance (ANOVA).
Differentially expressed genes were used to generate cluster diagrams with Cluster 3.0 (http://bioservices.capitalbio.com/xzzq/rj/3885.shtml) using the hierarchical method. This was followed by uncentered correlation, and average linkage parameters were chosen to calculate the distances in genes and samples.
Illumina sequencing and reads assembly
To obtain an overview of the ovary transcriptomes of A. bambawalei of different ages, RNA was extracted from female adults at the ES, IS, and AS stages and was then sequenced using the Illumina sequencing platform. Sequencing of normalized cDNA libraries resulted in approximately 33 Gbp of sequence data (GenBank accession No. PRJNA273425, Suppl. Tab. S1-1). Approximately, 17,011,204 to 19,407,760 paired-end (PE) raw reads were obtained from the ovaries of ES, IS, and AS stages, and the total length were 3,436,263,208 bp to 3,920,367,520 bp (Tab. 2, Suppl. Tab. S1-2). Using Velvet software for de novo assembly, all reliable reads from the samples of different ages were first assembled into contigs. Subsequently, the contigs were pooled together and then further assembled into scaffolds. Then, the scaffold sequences were assembled into clusters for consensus analysis. Finally, a total of 51,448 nonredundant (nr) sequences (transcripts), ranging from 100 to 24,805 bp were generated, and the average length was 1,915.36 bp (Suppl. Tab. S1-3, Fig. 2A). The number of high-quality reads was 32,626,310 to 37,128,892 in ES, IS, and AS stages, and the proportion of thelength of high-quality reads to the total length of the raw reads was 92.95% to 93.49% (Tab. 2, Suppl. Tab. S1-2). The value of N50 was 3,661 bp (Suppl. Tab. S1-4).
Figure 2: Unigene traits. (A)The distribution of the unigene length. The proportion of the sequences with matches (with a cutoff E-value of 1.0E−5) in the NCBI nr database was greater among the longer assembled sequences. (B) Characteristics of the homology search of query sequences aligned by BLASTx to the nr database. (B–1) E-value distribution of unigene BLASTx hits in the nr database with an E-value threshold of 1.0E−5. (B–2) Similarity distribution of the top BLAST hits for each sequence. (B–3) Species distribution of the first BLAST hits for each sequence with a cutoff E-value of 1.0E−5.
We analyzed the length distribution of all nr sequences after assembly. The proportions of the shorter and longer nr sequences were both greater than those of intermediate nr sequences. Specifically, the sequences of 100–500 bp and the sequences longer than 2000 bp showed matching ratios of 33.1% and 34.4%, respectively, whereas the matching ratios of the sequences of 500–1000 bp and 1000–2000 bp decreased to 14.3% and 8.2%, respectively (Fig. 2A).
Comparison with the GenBank nr database and the Swiss-Prot database revealed that most unigenes (31,269, 60.78%) were annotated (E-value < 10−5) (Suppl. Tab. S1-4). The rate of homologous sequences with E-values ranging from 1.0E−5 to 1.0E−50 was 52.78%, while the rate of homologous sequences with E-values less than 1.0E−50 was 47.22% (Fig. 2B-1). The percentage of sequences with a similarity higher than 80% distribution was 19.10%, while the majority of the sequences showed a similarity range from 20% to 80% (Fig. 2B-2). The species with the best BLASTx matching result (first hit) was Nasonia vitripennis (31.58%), followed by Tachypleus tridentatus (7.27%), Millerozyma farinose (3.95%), and Lactobacillus jensenii (3.95%) (Fig. 2B-3).
Annotation of predicted proteins
GetORF from EMBOSS software was used to analyze the reliable coding sequences (CDSs) of the 51,448 nr sequences and to further translate them into functional proteins. Among the 51,448 nr sequences, 42,790 unigenes (83.20%) contained ORFs, and the average ORF length was 1,181.10 bp (Suppl. Tab. S1-5). Comparisons with the GenBank nr and Swiss-Prot databases revealed that 31,269 nr sequences showed higher comparability with known gene sequences in existing species (with a cutoff E-value of 1.0E−5) (Suppl. Tab. S1-4). The proportion of annotated protein was large among the shorter sequences (100–500 bp) (Fig. 3A). Within the known genes, the average size of the encoded proteins was 1,134 amino acids. In addition, 8,658 nr sequences (16.80%) showed no significant similarity to any other protein sequences and were not annotated as unknown protein or putative proteins in the public database (Suppl. Tab. S1-5).
Gene Ontology (GO) and Clusters of Orthologous Groups (COG) classification
GO assignments were used to classify the functions of the annotated genes. Based on sequence homology, the sequences could be categorized into three different GO trees. Then, the three main categories were further classified into 46 functional groups (Suppl. Tab. S2-1). In the three main categories (cellular component, molecular function, and biological process) of the GO classification, the dominant groups were the “Cell/cell part”, “Binding”, and “Cellular process” groups. The genes with high expression levels were from the categories of “the Metabolic Process”, and “the Catalytic Activity”, whereas among the genes in the categories of “Viral Reproduction”, “Protein Tag”, “Viron Part”, and “Multi/-Organism Process”, only a few genes were expressed (Fig. 3B). We further performed phylogenetic classification with clusters of the COG database. In total, 11,058 genes were matched and then grouped into 25 functional categories (Suppl. Tab. S2-2, Fig. 3C). Among these categories, the cluster for “General function prediction only” represented the largest group (2,282, 26.06%), followed by “Posttranslational modification, protein turnover, and chaperones” (1,100, 9.95%), and “Translation, ribosomal structure and biogenesis” (835, 7.55%). Categories such as nuclear structure (9, 0.08%), cell motility (13, 0.12%), and RNA processing and modification (102, 0.92%) represented the smallest groups (Suppl. Tab. S2-2, Fig. 3C).
Figure 3: The traits of annotated proteins. (A) Distribution of the protein length. (B) Histogram presentation of Gene Ontology classifications. The results are summarized in three main categories: cellular location, molecular function, and biological process. The right y-axis indicates the number of genes in a category. The left y-axis indicates the percentage of a specific category of genes in that main category. (C) Histogram presentation of clusters of orthologous groups (COG) classification. A total of 11,508 sequences had COG classifications among the 25 categories.
To identify the biological pathways active in A. bambawalei, we mapped the annotated sequences into the reference canonical pathways in Kyoto Encyclopedia of Genes and Genomes (KEGG). Among the 51,448 nr sequences, approximately 4,225 transcripts could be grouped into six categories of 303 known metabolic or signaling pathways (Suppl. Tab. S2-3). The pathways with the most highly expressed genes in different stages included amoebiasis (12/157), human diseases, neurodegenerative diseases, Parkinson’s disease (14/235), and metabolism, energy metabolism, oxidative phosphorylation (14/276) (P < 0.01) (Tab. 3). These annotations provided a valuable resource for investigating the specific processes, functions, and pathways of A. bambawalei.
Analysis of differentially expressed genes (DEGs)
The age of A. bambawalei, i.e., the different stages after eclosion, led to significant changes in gene expression. From the 51,448 genes tested (Suppl. Tab. S1), 358 genes were differentially expressed (Suppl. Tab. S4-1). Notably, referring to previous literature, 63 genes exhibited changes possibly associated with the egg production and sperm storage of A. bambawalei, or 17.60% of the total expression changes (Suppl. Tab. S4). Pairwise analyses revealed that 128 genes showed significant expression differences (P < 0.001) between ES and IS, 124 genes showed significant expression differences (P < 0.001) between ES and AS, and 274 genes showed significant expression differences (P < 0.001) between IS and AS (Suppl. Tab. S4-2). In ES, IS, and AS, the number of upregulated genes were 59, 68, and 80, respectively, while the numbers of downregulated genes were 69, 56, and 194, respectively (Suppl. Tab. S4-2).
In ES vs. IS, some genes possibly related to egg production and sperm storage, such as, heat shock protein 90 (1), heat shock cognate protein 70 (1), serine protease 47 precursor (2) and serine protease inhibitor 3/4 (7) were upregulated and serine protease 43 precursor (1) was downregulated (Tab. S4-1). In ES vs. AS, heat shock cognate protein 70 (1) and serine protease 47 precursor (1) were still upregulated and fatty acyl-CoA reductase 1-like (8), acyl-CoA desaturase 1-like (1), major royal jelly protein-like 7 precursor (1), acyl-CoA Delta (11) desaturase-like isoform (1) and acyl-CoA Delta (11) desaturase (1) were downregulated (Suppl. Tab. S4-1). In IS vs. AS, serine protease 52 precursor (1) and serine protease 43 precursor (1) were upregulated and fatty acyl-CoA reductase 1-like (8), acyl-CoA desaturase 1-like (2), acyl-CoA delta (11) desaturase-like isoform (1), acyl-CoA delta (11) desaturase (1), serine protease inhibitor 3/4 (8), serine protease 47 precursor (2), major royal jelly protein 1 (3), major royal jelly protein-like 7 precursor (6), and glucose dehydrogenase [acceptor]-like (3) were downregulated (Suppl. Tab. S4-1). Moreover, the qPCR results confirmed that these genes exhibited stage-specific expression patterns (Fig. 4). Serine protease 47 precursor and serine protease inhibitor 3/4 were highly expressed in the ES stage, while glucose dehydrogenase, fatty acyl-CoA reductase 1-like, major royal jelly, and acyl-CoA delta (11) desaturase-like were upregulated in the AS stage. In addition to these differentially expressed genes, some other genes potentially related to sex allocation were highly expressed in all stages, such as purity-of-essence (poe), neurocalcin, scavenger receptor class B member 1, bcl-2-related ovarian killer gene, and uniquinone biosynthesis (Suppl. Tab. S5).
To understand the biological functions of these genes, we annotated the most highly upregulated pathways on the metabolic pathway map. Notably, specific enrichment was observed in pathways likely related to fertilization, such as the “genetic information processing”, “human disease”, and “organismal systems” (Suppl. Tab. S4).
Figure 4: Relative expressions levels of serine protease 47 precursor (A), glucose dehydrogenase (B), fatty acyl-CoA reductase 1-like (C), serine protease inhibitor 3/4 (D), major royal jelly protein-like 7 precursor (E), and acyl-CoA delta (11) desaturase-like (F) at different ages. Error bars represent the means ± standard deviations from 3 biological replicates. One-way ANOVA was used to determine the significant differences with different lowercase letters (a–c) (P < 0.05).
In the present study, we characterized the transcriptome data of A. bambawalei ovaries using the Illumina HiSeq sequencing platform to detect the genes likely related to age and fertilization. Our results indicated that from a total of 51,448 genes tested, 358 genes were differentially expressed in this parasitic wasp. Notably, referring to previous literature, 63 genes exhibited changes possibly associated with the egg production and sperm storage of A. bambawalei, or 17.60% of the total differentially expressed changes. Egg production is one of the important aspects of parasitic wasp fitness (Jervis and Kidd, 1986; Rivero-Lynch and Godfray, 1997), and the egg load strategy has a key role in age-related egg production regulation (Flanders, 1942). Our pilot study indicated that the egg load increases with female age but later declines. A likely explanation is that the increase in the proportion of females with resorbing oocytes occurs due to host deprivation in synovigenic species (Flanders, 1942). This is consistent with our observation that the Aenasius females are synovigenic. Additionally, the mated female Aenasius of different ages in our study expressed highly conserved genes critical for reproduction across several insect species. Among them, some genes exhibited stage-specific expression patterns, consistent with the physiological functions of these genes (Belancio et al., 2015). For example, in the ES stage, putative proteolysis regulators (serine proteases 47 precursor and serine protease inhibitor 3/4) and heat shock proteins (HSPs) were highly expressed. These genes could cleave inactive molecules into their active forms (Fenker and Stanfield, 2015), participate in oogenesis (Ambrosio and Schedl, 1984; Neuer et al., 2000) and ovulation (LaFlamme et al., 2012), and enhance oogenesis (Pilpel et al., 2008) in a number of species. Therefore, these results could predict that the high expression levels of these genes in the ES stage might be relate to the oogenesis of A. bambawalei.
In addition to egg production, the dynamics of sperm storage are another factor influencing the fitness of parasitic wasps (Chevrier and Bressac, 2002). Glucose dehydrogenase (GLD), expressed in the spermathecal duct, is associated with the storage and utilization of sperm in Drosophila melanogaster (Iida and Cavener, 2004). We confirmed that GLD is indeed upregulated in the AS stage (Tab. S2, Fig. 4) and discovered that a total of three GLD genes were differentially expressed in this stage. Given that sperm storage and release are very important for parasitic female wasps (Cook et al., 2015), we suggest that this gene family may have a role in the sperm storage of Aenasius females. At the same time, some genes that could change the energy utilization and metabolism associated with reproductive behavior (Schmitzová et al., 1998; Imjongjirak et al., 2005) were upregulated in older Aenasius specimens, such as the major royal jelly protein-like 7 precursor (MRJPs) (Thompson et al., 2006), fatty acyl-CoA reductase 1-like, and acyl-CoA delta (11) desaturase (Swanson et al., 2001; Pannebakker et al., 2013). The widespread presence of these genes across distant species supports the hypothesis of an important role of such molecules in the female Aenasius. However, some other genes potentially associated with egg production and sperm storage did not vary in all stages. These genes include the positive and the negative genes associated with parasitic wasp fitness. For example, the scavenger receptor class B member 1 gene associated in the process of future oogenesis to replace the eggs being deposited in N. vitripennis (Iida and Cavener, 2004; Cook et al., 2015) and spermatogenesis (Casado et al., 2012). The neurocalcin gene is present in spermatogenic cells (Jankowska et al., 2014) and may be a part of the fertilization machinery from humans to mice (Jankowska et al., 2010). The purity-of-essence (poe) gene is potentially associated with sterility in the oviposition of N. vitripennis (Pannebakker et al., 2013) and is a late block to spermatogenesis in D. melanogaster (Fabrizio et al., 1998). This is concordant with the fact that hormone-regulated cyclic cell turnover occurred in all stages.
The genes highly expressed in all stages can be interpreted from two aspects. First, because many genes are expressed regardless of an organism’s developmental stage (Pannebakker et al., 2013), most of the highly expressed genes in A. bambawalei showed no significant expression differences between ES, IS, and AS. Second, one of the crucial limitations of behavioral transcriptomics at the moment is that the annotations of genes are associated with molecular functions at a cellular or tissue level that may seem a long way from the regulation or control of a focal behavior (unless a key developmental gene or signal pathway is identified as a hit) (Pannebakker et al., 2013). Therefore, to understand the biological functions of these genes, we annotated the most highly upregulated pathways on the metabolic pathway map. Specific enrichment was observed in pathways such as the “genetic information processing”, “human disease”, and “organismal systems”. These pathways will be further investigated in future functional research on sex allocation related genes. Many candidate novel genes that may be involved in egg production and sperm storage mechanisms are worthy of further investigation. By studying and regulating these genes, we could artificially control the offspring sex ratio of A. bambawalei in the future.
We found differences in highly expressed genes at three stages of A. bambawalei. The highly expressed genes were related to genetic information processing, human disease, and organismal systems. The differences may be attributed to the aging of organisms. Other genes related to egg production and sperm storage were highly expressed in all three stages of A. bambawalei. Although there were no statistically significant differences, this information will be valuable for analyzing the reproductive tactics of A. bambawalei. The results of this study suggest that further analyses of the transcriptomes of A. bambawalei ovaries at different stages could clarify the roles of the many important genes related to egg production and sperm utilization.
In conclusion, the different stages of A. bambawalei females led to significant gene expression in ovaries. Among the differentially expressed genes, 294 genes (approximately 69.55% of the total differentially expressed genes) exhibited changes likely associated with the egg production and sperm storage of A. bambawalei. In the future, by studying and regulating these genes, we might be able to artificially control the offspring sex ratio of A. bambawalei to some extent.
Acknowledgement: The authors are deeply grateful to Prof. Wang Xiaowei of the Zhejiang University, for his valuable critical comments and suggestions.
Availability of Data and Materials: All data generated or analyzed during this study are included in this published article (and its supplementary information files).
Supplementary Material: The supplementary material is available online at DOI: https://10.32604/biocell.2021.016563.
Authors’ Contributions: The authors confirms contribution to the paper as follows: conceived and designed the experiments: JZ, JH, YT, XL; contributed reagents, materials, analysis tools and data: YT, XL; performed the experiments: JZ; analyzed and interpreted the data: JZ, JH; draft manuscript preparation: JZ, JH. All authors reviewed the results and approved the final version of the manuscript.
Ethics Approval: This research does not involve human data and we conducted studies in areas where P. solenopsis and A. bambawalei occur, and no specific permits were required.
Funding Statement: J Zhang was supported by National Natural Science Foundation of China (31801801), the Zhejiang Provincial Public Welfare Technology Research Program (LGN20C140004) and Foundation of Guangxi Key Laboratory of Biology for Crop Diseases and Insect Pests (2015-KF-4), and J Huang was supported by National Natural Science Foundation of China (31772234).
Conflicts of Interest: The authors declare that they have no conflict of interest to report regarding the present study.
SUPPLEMENTARY TABLE S1. The total tested genes.
SUPPLEMENTARY TABLE S2. GO, COG, KEGG.
SUPPLEMENTARY TABLE S3. Signaling pathways.
SUPPLEMENTARY TABLE S4. Differentially expressed genes in three stages.
SUPPLEMENTARY TABLE S5. Genes likely related to fertilization but not differentially expressed.
 | This work is licensed under a Creative Commons Attribution 4.0 International License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. |