Mating-Induced Trade-Offs upon Egg Production versus Fertilization and Offspring’s Survival in a Sawfly with Facultative Parthenogenesis

Simple Summary Study of mating-induced trade-offs between reproduction and survival is conducive to provide evolutionary insights into reproductive strategies and aging. Using RNA sequencing and bioinformatics, we found that mating induced changes of genes and pathways related to reproduction and survival in females of a pine sawfly. Mating induced substantial downregulation on genes associated to immunity, stress response, and longevity. However, mating induced divergent reproductive response, with downregulation on genes related to egg production while upregulation on genes related to egg fertilization. Considering the nature of limited resources in adults, low fecundity and egg protection behavior in this sawfly, we suggest that mating triggers trade-offs between reproduction and survival in this insect and females of this species have evolved specific strategies to adapt to the living conditions, e.g., restrict whole fecundity to ensure higher fertilization and offspring’s survival. Abstract Investigation of mating-induced trade-offs between reproduction and survival is conducive to provide evolutionary insights into reproductive strategies and aging. Here, we used RNAseq and bioinformatics to reveal mating-induced changes of genes and pathways related to reproduction and survival in female Cephalcia chuxiongica, a pine defoliator with facultative parthenogenesis and long larval dormancy. Results showed that mating induced substantial downregulation on genes and pathways associated to immunity, stress response, and longevity. However, mating induced divergent reproductive response, with downregulation on genes and pathways related to egg production while upregulation on genes and pathways related to egg fertilization. Considering the nature of limited resources in adults, low fecundity, and egg protection behavior in C. chuxiongica, we suggest that mating triggers trade-offs between reproduction and survival in this insect and females of this species may have evolved specific strategies to adapt to the environmental and hosts’ conditions, e.g., restrict whole fecundity to ensure higher fertilization and offspring’s survival. Moreover, mating induced significant responses on genes and pathways that play important roles in vertebrate reproduction while their function in insects are unclear, such as the progesterone-mediated oocyte maturation pathway; the significant regulation after mating suggests that their function may be evolutionarily conserved in animal kingdom.


Introduction
Trade-offs between reproduction and survivorship is the central theme in evolutionary biology of senescence. Negative relationships between reproduction and survivorship are Insects 2021, 12, 693 2 of 18 often observed in organisms, in which increased fecundity may cause costs on somatic maintenance (such as immunity and stress responses) and longevity [1][2][3][4]. The resources allocation model suggests that the trade-offs between reproduction and survival is likely due to alternative allocation of limiting energetic resources because both reproduction and soma maintenance are energetically costly [1,5]. Genetic model suggests that there are antagonistic alleles that may promote negative genetic correlation between reproduction and survival, i.e., alleles may promote fecundity at the expense of somatic maintenance [2,6]. These two models are not mutually exclusive [2,7].
In insects, the endocrine network plays vital roles in reproduction and maintenance regulation, which mainly involves insulin-like/IGF-1 signaling (IIS), juvenile hormone (JH), 20-Hydroxyecdysone (20E), yolk precursor vitellogenin (Vg), and yolk proteins (YPs) [2]. JH and IIS are important mediators of egg maturation; improved JH and IIS levels promote the expression of Vg and YPs and the uptaking of Vg and YPs into oocytes [1]. However, elevated JH and IIS may inhibit immune responses, such as by reducing phenoloxidase activity and the production of antimicrobial peptides [1,8]. In contrast, 20E may promote immune responses, such as by upregulating the expression of antimicrobial peptide genes in Drosophila melanogaster [8,9]. Therefore, although the mechanism remains poorly understood, such opposite effects of different hormones seem to be intimately involved in the trade-offs between reproduction and survivorship [1,2].
Female insects are likely to face a limited protein supply during reproductive period because many of them do not feed on a protein source as adults while elevated egg production process requires higher protein supplies at this stage [10][11][12]. Female insects thus may have evolved behavioral and physiological strategies to allocate limited resources optimally between reproduction and survival [1,13]. Mating is a key switch for sexual reproduction in insects, which incurs major changes in the physiology and behavior in females [14]. Previous studies have found that mating did positively affect female reproduction and negatively affect female immune activities in some insects [1,5,15]. However, some studies also found that there is no difference in immune responses between mated and virgin insect females and in some insects mating even upregulates female immune responses [5]. These diverse findings may be due to the differences in measurement methods and immune indicators and the differences in species and mating systems. For example, males may directly suppress female immunity, thereby promoting sperm storage and egg fertilization by their sperms in female reproductive tracts [16]. However, mated females may need to have higher post-mating immune activity to provide defense against mating transferred infections, particularly in polygamous mating systems [5,17]. The heat shock response is essential for proteostasis and cellular health [4]. Similarly, studies on reproduction-induced heat shock response also found inconsistent results. In C. elegans, the heat shock response was inhibited at the onset of reproduction, probably due to competing requirements of the germline and soma [4]. In honey bees, however, increased reproduction did not cause loss of heat shock response in the reproductive queen [3]. Therefore, honey bee queens may possess an atypical uncoupling of the reproduction-maintenance trade-off, in which stress response may be maintained during reproduction [3]. Therefore, study on the trade-offs between reproduction and survival in different mating systems by using methods that can more comprehensively examine various indicators (such as high-throughput sequencing and bioinformatics) will provide deeper insights in this field.
Recent progress on high-throughput sequencing and bioinformatics has revolutionized our understanding of life and organism. Previous transcriptome analysis in a number of insect species have shown that mating can induce expression changes in genes related to reproduction, immunity, stress response, and longevity [18][19][20]. A resent transcriptomebased study in the sweet potato whitefly females found 434 deferentially expressed genes (DEGs) between mated and unmated groups, with many of them encoding binding-related proteins and genes associated with lifespan [20]. Another study in Spodoptera litura female moths showed a divergent response in DEGs in relation to reproduction and immunity [15]. These results suggest that trade-offs on reproduction versus survival were induced by Insects 2021, 12, 693 3 of 18 mating in some insects. In the present study, we studied genes and pathways related to reproduction and survival (immunity, stress response, and longevity) and mating-induced regulation in females of a pine sawfly, Cephalcia chuxiongica. This insect has a facultative parthenogenesis reproductive system and a long (19 months) larval dormancy stage [21]. Study of mating-induced regulation on reproduction and survivorship in such a special insect species is expected to achieve some new enlightenment on the evolution of reproductive strategies and aging.
The sawfly Cephalcia chuxiongica Xiao (Hymenoptera: Pamphiliidae) is a serious pest of pines in China [22,23]. Almost fifty Cephalcia species have been found all over the world, with most of them being important forest pests [24,25]. Cephalcia chuxiongica larvae feed on pine needles and pupate in the soil under the trees. The life cycle of this species is long (about 22 months) due to a 19-month larval diapause in the soil [22]. Similar to many other species from Hymenoptera, C. chuxiongica females perform sexual reproduction under normal conditions, while in the absence of males, they can perform parthenogenesis [22]. The male and female adults of C. chuxiongica can mate a few hours after eclosion; unmated females (virgin) may start to lay unfertilized eggs through parthenogenesis a few days (>3 d) after eclosion [21]. Sawflies do most of their feeding during the larval stage, and adults of some species may feed on pollen and nectar [26]. So far, there is no report on adult feeding habits in C. chuxiongica; adults even did not feed on the provided honey solution, suggesting that this species do not feed during the adult stage [21].

Sample Collection
Cephalcia chuxiongica adults were collected in August 2019 under infested pine trees in the Pinus yunnanensis forest, located on a small mountain (25.600166N, 103.439364E) near Xundian town in Yunnan Province, China. Cephalcia chuxiongica pupate and eclose underground during July to Sept. Eclosed adults crawl out from the soil for reproduction in the forest. To confirm adults virginity and age, newly eclosed adults were directly dug out from the soil below the pine trees in the morning [27]. Males and females were sexed according to morphological traits [23] and reared separately on the pine trees (on the pine needles covered by nylon mesh bag [30 cm diameter, 40 cm long] to avoid insect escape). As mentioned above, C. chuxiongica adults can mate soon after emergence and mating mostly happen at noon. Therefore, mating was allowed in the same day after collection during 12:00-14:00 by pairing males and females on the pine needles covered by bags with one pair per bag. Mating events (two insects engaged at the tip of the abdomen; mating duration is about 20 min) [22] were recorded. Males were removed after mating and the mated females were then individually reared in the same bag and their abdomens were sampled at 1 h, 6 h, and 24 h after mating. Abdomens from fifteen females were used as a replicate, and three replicates were used for each sampling time point. Virgin females at the same age as mated individuals were used as controls. All samples were placed in liquid nitrogen for immediate freezing after sampling and stored at −80 • C. One replicate from Virgin-1 h and one replicate from Virgin-6 h were discarded due to lower sequencing quality (see below) and thus finally, two replicates for Virgin-1 h and Virgin-6 h, and three replicates for other treatments (Virgin-24 h, Mated-1 h, Mated-6 h and Mated-24 h) were used for analysis.

Library Preparation and Sequencing
Total RNA was extracted from samples using Trizol reagent (Invitrogen Inc., Calsbad, CA, USA) and then was treated with the RNase-free DNase I to eliminate genomic DNA. The purity and concentration of RNA was assessed by using Qubit RNA Assay Kit (Life Technologies Inc., Grand Island, NY, USA) and the NanoPhotemeter spectrophotometer (Implen Inc., Westlake Village, CA, USA). The RNA integrity was checked by the Agilent Bioanalyzer 2100 system (Agilent Technologies Inc., Palo Alto, CA, USA). One microgram RNA per sample was used for the preparation of the sequencing libraries by using NEBnext Ultra RNA Library Prep Kit for Illumina (New England BioLabs Inc., San Diego, MA, USA) following the manufacturer's instructions and index codes were added to attribute sequences to each sample. The quality of each sample library was assessed using the Agilent Bioanalyzer 2100 system. Ultimately the library per sample was sequenced by Illumina Hiseq4000 platform (Majorbio Biotechnology Inc., Shanghai, China) and the paired-end reads were generated.

Quality Control and Assembly
The raw reads obtained were initially processed by trimming the adapter and low quality reads to produce clean reads. The clean reads were assembled using the Trinity software (version 2.5.1) to generate transcripts [28]. Then transcript analysis was performed to remove redundancies with TGICL software (version 2.1) and acquire unigenes without redundancy [29].

Differential Expression Analysis
Gene expression levels were determined as transcripts per million (TPM). The differential expression analysis between samples was performed using the edgeR R package (3.0.8). p-value was adjusted using q-value [30]. q < 0.05 and |log2(foldchange)| > 1 was set as the threshold for significantly differential expression.

Functional Annotation and Enrichment Analysis of Differentially Expressed Genes (DEGs)
Using the BLAST [31] software, the DEGs were compared with NCBI non-redundant protein (NR), Swiss-Prot protein database (Swiss-Prot), Gene Ontology (GO), KEGG Orthology (KO), Cluster of Orthologous Groups (COG), and Pfam databases to obtain annotation information about the DEGs. GO enrichment analysis of DEGs was implemented by using the GOSeq program and KEGG enrichment was performed using the KOBAS software. GO terms and KEGG pathways with q < 0.05 were significantly enriched in DEGs.

Validation by qRT-PCR
To verify the accuracy of the RNAseq data, 24 DEGs were used for qRT-PCR and the β-Tubulin (GeneBank ID: MZ028603) was used as a reference gene. Total RNA was extracted from above samples by using RNAiso plus (TaKaRa Inc., Dalian, China) and cDNA was synthesized using PrimeScript RT reagent Kit (Takara Inc., Dalian, China). The PCR was performed by QuantStudio 7 Flex (Thermo Fisher Scientific Inc., Waltham, MA, USA) with gene specific primers (Table S1) and the following program: 95 • C for 30 s, followed by 40 cycles of 95 • C for 5 s, 60 • C for 30 s, and dissociation. The 2 −∆∆CT method [32] was used to calculate the relative expression. Differences of gene expression levels between treatments were analyzed by ANOVA using SPSS 25.0. The rejection level was set at α < 0.05.

Sequencing and Assembly
By RNAseq using Illumina HiSeq4000 platform,~55,000,000 clean reads were obtained from each of the 16 sequenced libraries (Table S2). The percentages of Q20 and Q30 of all samples' clean reads ranged from 97.76% to 98.21% and from 94.10% to 95.21%, respectively. The biological replicates were highly correlated (Table S3), which affirmed the reproducibility of the RNASeq technology and biological replicates. The transcriptome raw reads were deposited to the NCBI SRA database (accession no.: SRR13991392~SRR13991407).  (Figure 2; Table S7). . Genes with significant differential expression were indicated by red dots (upregulated) and green dots (downregulated). Genes with no significant differential expression were represented by grey dots.  . Genes with significant differential expression were indicated by red dots (upregulated) and green dots (downregulated). Genes with no significant differential expression were represented by grey dots. SRR13991392~SRR13991407).

Outline of Mating-Induced Transcriptional Changes
There are 149, 320, and 1068 DEGs within Mated-1 h vs.  (Figure 2; Table S7). . Genes with significant differential expression were indicated by red dots (upregulated) and green dots (downregulated). Genes with no significant differential expression were represented by grey dots.  To better understand their functions, all these DEGs were annotated based on NR, Swiss-Prot, GO, KO, COG, and Pfam databases (Tables S4-S6) and then were submitted for GO (Table S8) and KEGG (Table S9) enrichment analysis. Most of these DEGs (1235/1537 = 80.4%) were successfully mapped to at least one of these databases (Tables S4-S6). A total of 309 GO terms were significantly enriched in these DEGs, with most of them (59.5%) being assigned to biological process (BP) terms, and then 17.2% and 23.3% being assigned to cellular component (CC) terms and molecular function (MF) terms, respectively (Figures 3a and 4; Table S8). Moreover, there are 51 KEGG pathways that were significantly enriched in these DEGs, which were assigned to six of Level I Categories (Figures 3b and 5;  Table S9). (1235/1537 = 80.4%) were successfully mapped to at least one of these databases (Tables S4-S6). A total of 309 GO terms were significantly enriched in these DEGs, with most of them (59.5%) being assigned to biological process (BP) terms, and then 17.2% and 23.3% being assigned to cellular component (CC) terms and molecular function (MF) terms, respectively (Figures 3a and 4; Table S8). Moreover, there are 51 KEGG pathways that were significantly enriched in these DEGs, which were assigned to six of Level I Categories (Figures 3b and 5; Table S9).    Mating-induced regulation on genes and pathways related to reproduction and survivorship were then studied in detail as follows based on the above annotation and enrichment analysis.

Transcriptional Changes at 1 h Post-Mating
At 1 h after mating, 68 genes were downregulated and 81 genes were upregulated in mated females compared to virgin females (Figure 1a). The log2FoldChange (LC) value of DEGs varied from −13.04 to 12.74 (Table S4).
Two reproduction-related genes were found within these DEGs, with one encoding follicle cell protein (downregulated) and one encoding cytochrome P450 (upregulated) ( Table 1). KEGG enrichment analysis of the 149 DEGs found one reproductive related pathway, the steroid hormone biosynthesis pathway, which enriched to one upregulated DEG (Table 2).
Two immunity related genes were found within these DEGs ( Table 3). Both of them encode antimicrobial peptides and both were downregulated in mated females compared to virgin ones. Six immunity related pathways were enriched based on the 149 DEGs, with four pathways enriched to downregulated DEGs and two pathways enriched to upregulated DEGs (Table 4).
No longevity and heat shock-related genes and pathways were found within these DEGs of M1 h-vs.-V1 h. Mating-induced regulation on genes and pathways related to reproduction and survivorship were then studied in detail as follows based on the above annotation and enrichment analysis.

Transcriptional Changes at 1 h Post-Mating
At 1 h after mating, 68 genes were downregulated and 81 genes were upregulated in mated females compared to virgin females (Figure 1a). The log2FoldChange (LC) value of DEGs varied from −13.04 to 12.74 (Table S4).
Two reproduction-related genes were found within these DEGs, with one encoding follicle cell protein (downregulated) and one encoding cytochrome P450 (upregulated) ( Table 1). KEGG enrichment analysis of the 149 DEGs found one reproductive related pathway, the steroid hormone biosynthesis pathway, which enriched to one upregulated DEG (Table 2).
Two immunity related genes were found within these DEGs ( Table 3). Both of them encode antimicrobial peptides and both were downregulated in mated females compared to virgin ones. Six immunity related pathways were enriched based on the 149 DEGs, with four pathways enriched to downregulated DEGs and two pathways enriched to upregulated DEGs (Table 4).
No longevity and heat shock-related genes and pathways were found within these DEGs of M1 h-vs.-V1 h. Plays a role in the storage of replication-dependent histone mRNAs and proteins during oogenesis. [53]

Transcriptional Changes at 6 h Post-Mating
Within the 320 DEGs of M6 h-vs.-V6 h, 206 were downregulated and 114 were upregulated in mated females compared to virgin ones (Figure 1b). The LC value of DEGs differed from -10.34 to 12.29 (Table S5).
Sixteen reproductive-related genes (all downregulated) were found within these DEGs (Table 1) One immunity-related gene was found within these DEGs (Table 3), which is the phenoloxidase-encoding gene (LC: −3.53). Three immunity-related KEGG pathways were significantly enriched to 12 downregulated DEGs (Table 4).
Still no longevity and heat shock related genes and pathways were found within these DEGs.

Transcriptional Changes at 24 h Post-Mating
Within the 1068 DEGs of M24 h-vs.-V24 h, 104 genes were downregulated and 964 genes were upregulated in mated females compared to virgin ones (Figure 1c). The LC value of DEGs changed from −10.50 to 14.64 (Table S6).
A total of 51 transcripts were annotated as reproductive related genes within these DEGs (Table 1), including seven Vitellogenin and Vitellogenin-like (all downregulated), one Vitellogenin-receptor (upregulated), 26 Zona pellucida protein-like (all upregulated), seven Cyclin-like (all upregulated), and ten other reproductive-related genes. Three reproductive related KEGG pathways were significantly enriched (  (3) progesterone-mediated oocyte maturation pathway, enriched to 13 upregulated DEGs. To get some evolutional clues on oocyte maturation, unigenes annotated to the progesterone-mediated oocyte maturation pathway (Table S10) were mapped to the same pathway of Xenopus (KEGG id: map04914) ( Figure 6). Results show about half of the genes on map04914 were found in the unigenes of C. chuxiongica. This study also revealed an additional 25 unmapped progesterone-mediated oocyte maturationrelated unigenes (Table S10).

Transcriptional Changes at 24 h Post-Mating
Within the 1068 DEGs of M24h-vs.-V24h, 104 genes were downregulated and 964 genes were upregulated in mated females compared to virgin ones (Figure 1c). The LC value of DEGs changed from −10.50 to 14.64 (Table S6).
A total of 51 transcripts were annotated as reproductive related genes within these DEGs (Table 1), including seven Vitellogenin and Vitellogenin-like (all downregulated), one Vitellogenin-receptor (upregulated), 26 Zona pellucida protein-like (all upregulated), seven Cyclin-like (all upregulated), and ten other reproductive-related genes. Three reproductive related KEGG pathways were significantly enriched (  (3) progesterone-mediated oocyte maturation pathway, enriched to 13 upregulated DEGs. To get some evolutional clues on oocyte maturation, unigenes annotated to the progesterone-mediated oocyte maturation pathway (Table S10) were mapped to the same pathway of Xenopus (KEGG id: map04914) ( Figure 6). Results show about half of the genes on map04914 were found in the unigenes of C. chuxiongica. This study also revealed an additional 25 unmapped progesterone-mediated oocyte maturation-related unigenes (Table S10). Three immunity related genes were found within these DEGs with all of them being downregulated (Table 3). Thirteen immunity-related KEGG pathways were significantly enriched, with ten of them being enriched to downregulated DEGs and three being enriched to upregulated DEGs (Table 4).
One longevity-related pathway was significantly enriched, which was related to four downregulated DEGs (Table 4).
Three immunity related genes were found within these DEGs with all of them being downregulated (Table 3). Thirteen immunity-related KEGG pathways were significantly enriched, with ten of them being enriched to downregulated DEGs and three being enriched to upregulated DEGs (Table 4).
One longevity-related pathway was significantly enriched, which was related to four downregulated DEGs (Table 4).

Validation of RNAseq by qRT-PCR
Eight DEGs from each group (M1 h-vs.-V1 h, M6 h-vs.-V6 h and M24 h-vs.-V24 h, respectively) were used to verify the accuracy of RNAseq, which include six reproductionrelated genes (DN64818_c0_g6, DN52313_c1_g5, DN65191_c2_g1, DN66978_c2_g1, DN78694_ c0_g1, and DN71364_c1_g1), three immunity-related genes (DN70573_c2_g1, DN60391_c5_g2 and DN63597_c1_g1), and one heat shock response gene (DN70765_c7_g9) (Figure 7). The expression levels of these genes measured by qRT-PCR were similar to the results of RNAseq analysis, which suggested that the RNAseq data were reliable.  Figure 7). The expression levels of these genes measured by qRT-PCR were similar to the results of RNAseq analysis, which suggested that the RNAseq data were reliable.

Discussion
Previous studies generally found that mating induces the upregulation of egg production (fecundity)-related genes, such as YPs and Vg [15,18,19]. In the present study, however, we found an opposite result, where all Vg encoding DEGs (five in M6h-vs.-V6h and seven in M24h-vs.-V24h) were found downregulated largely (Table 1). Moreover, other positive modulator of fecundity, including Follicle cell protein, Haemolymph juvenile hormone binding protein, UDP-glucosyltransferase, Apolipophorin III, and Matrix metalloproteinase, were also downregulated significantly (Table 1). In addition, the upregulation of Insulinase (which destroys or inactivates insulin) in mated females may also negatively affect female fecundity. These results suggested that mating downregulated reproduction in terms of fecundity (number of eggs) in C. chuxiongica.
On the contrary, we found a significant upregulation on genes related to oocyte maturation and embryogenesis (Table 1): (1) Seven Cyclins, which play roles in the process of oocyte maturation and the onset of embryogenesis [47][48][49]; (2) two Zygote arrest protein, which are ovary-specific maternal factors that play essential roles during the oocyte-toembryo transition [51]. Moreover, two oocyte maturation-related pathways, the oocyte meiosis pathway and the progesterone-mediated oocyte maturation pathway, have been significantly enriched upregulated DEGs (Table 2). Fully grown oocytes are arrested in

Discussion
Previous studies generally found that mating induces the upregulation of egg production (fecundity)-related genes, such as YPs and Vg [15,18,19]. In the present study, however, we found an opposite result, where all Vg encoding DEGs (five in M6 h-vs.-V6 h and seven in M24 h-vs.-V24 h) were found downregulated largely (Table 1). Moreover, other positive modulator of fecundity, including Follicle cell protein, Haemolymph juvenile hormone binding protein, UDP-glucosyltransferase, Apolipophorin III, and Matrix metalloproteinase, were also downregulated significantly (Table 1). In addition, the upregulation of Insulinase (which destroys or inactivates insulin) in mated females may also negatively affect female fecundity. These results suggested that mating downregulated reproduction in terms of fecundity (number of eggs) in C. chuxiongica.
On the contrary, we found a significant upregulation on genes related to oocyte maturation and embryogenesis (Table 1): (1) Seven Cyclins, which play roles in the process of oocyte maturation and the onset of embryogenesis [47][48][49]; (2) two Zygote arrest protein, which are ovary-specific maternal factors that play essential roles during the oocyte-toembryo transition [51]. Moreover, two oocyte maturation-related pathways, the oocyte meiosis pathway and the progesterone-mediated oocyte maturation pathway, have been significantly enriched upregulated DEGs (Table 2). Fully grown oocytes are arrested in the first meiotic prophase (Figure 6), which is a state of low metabolic activity without detrimental effects on subsequent embryogenesis [57]. At the end step of oocyte maturation, a process termed meiotic reinitiation or resumption, takes place before or during the time when the oocyte moves from the follicle to the oviduct (ovulation). Meiotic maturation involves the activation of various signal transduction pathways, which converge to activate maturation-promoting factors ( Figure 6). In vertebrate, the marks of meiotic maturation include [58]: (1) Resumption of meiosis I, including germinal vesicle breakdown, chromosome condensation, and spindle formation; (2) the transition from meiosis I to meiosis II; and (3) arrest in metaphase II. Meiosis II will complete after egg fertilization. In many invertebrates, however, oocytes maturation proceeds only to metaphase of meiosis I, which is when they are fertilized [58]. Therefore, the upregulation of genes and pathways in relation to oocyte meiosis and maturation in C. chuxiongica (Tables 1 and 2) are likely to function in meiotic resumption and maturation in fully grown oocytes, from which the egg is ready for the following fertilization.
The progesterone-mediated oocyte maturation pathway of the African clawed frog Xenopus laevis is the most intensively studied model system for meiotic maturation [58]. Therefore, we mapped unigenes that annotated to the progesterone-mediated oocyte maturation pathway to the same pathway of Xenopus (KEGG id: map04914) and found that about half of the genes on map04914 were present in the unigenes of C. chuxiongica ( Figure 6; Table S10). Further searching found an additional 25 unmapped progesteronemediated oocyte maturation-related unigenes in C. chuxiongica (Table S10). These results suggest that this pathway may also play roles in egg maturation of C. chuxiongica but factors and signaling pathways may be species-specific. In the present study, we also found the presence of Cyclin-B1 (KO id: K05868) and Xkid (KO id: K10403; kinesin-like DNA-binding protein, also known as Kinesin-like protein, KIF22) [57] in the enriched DEGs ( Figure 6, Table S10). These proteins are all dispensable for meiosis I entry but play essential roles in the progression from meiosis I to meiosis II [58]. As mentioned above, oocytes maturation in many invertebrates proceeds only to metaphase of meiosis I until they are fertilized. If it is the case in C. chuxiongica, then the upregulation of Cyclin-B1 and Xkid may be triggered by mating factors (such as male accessary protein and sperms) or fertilization, which then will promote the transition of eggs from meiosis I to meiosis II. This warrants further studies.
More interestingly, this study also found significant expression changes in a number of genes that may relate to fertilization and egg hatching, including a number of Zona pellucida protein and Ovastacin/Astacin-like (Table 1). The Zona pellucida (ZP) is a glycoprotein layer surrounding the plasma membrane of mammalian oocytes, which is a vital constitutive part of the oocyte and functions in primary binding and induction of the sperm acrosome reaction [42,44]. A family of 18 ZP-like protein encoding genes has been identified in D. melanogaster [59]. These genes are specifically expressed during embryogenesis and during differentiation of epithelial tissues. However, whether these ZP-like proteins also play roles in fertilization is still unknown in D. melanogaster or other insect species because little is known about the molecular mechanism of fertilization in insects [59]. Ovastacin plays a role in the post-fertilization block to sperm binding by cleaving ZP2 in the ZP to prevent polyspermy in mouse [52]. Some astacins may function as hatching enzyme in insects, which degrades the egg envelope to release the embryo [60]. The reproductiverelated regulation on these genes may suggest they are crucial in the reproduction process of insects.
Insect cytochrome P450 families contain a class of enzymes, which play diverse functions in detoxification and the biosynthesis of hormones [34]. In Drosophila, mating induced downregulation in six P450 genes and upregulation in 22 P450 genes. In B. tabaci, mating induced upregulation in two P450 genes [20]. In the present study, one P450 gene (CYP2C) was significantly upregulated in mated females shortly after mating (1 h postmating) and four P450 genes were downregulated in mated females at 6 h post-mating (Table 1). Upregulation of P450 shortly after mating may be a detoxification response as males may transfer slightly toxic seminal fluid during mating [61] and downregulation of P450 sometime after mating may be due to trade-offs between reproduction and survival.
Reproduction associated heat shock responses are inconsistent in different species, which is downregulated after the onset of reproduction in C. elegans [4] while in honey bees, increased reproduction did not cause loss of heat shock response in the reproductive queen [3]. In B. tabaci, mating led to the upregulation of Hsp68 and Hsp70 genes at different time points in mated females [20]. In the present study, we also found significant expression changes in thirteen Hsps-encoding transcripts 24 h post-mating, with some of them (four sHsps and four Hsp70s) being downregulated and others (three Hsp70s and two Hsp90s) being upregulated substantially (Table 3). Hsps act as molecular chaperones to improve organisms' survival, development, and reproduction under different stresses. Studies have shown that Hsps have multiple functions during reproductive process, including gamete protection [58], oocyte maturation ( Figure 6) [62], reproduction vs. survival trade-offs [63], etc. In addition, one longevity-associated pathway was significantly enriched, which related to four downregulated Hsp70s coding DEGs (Table 4). Therefore, mating induced different response in different Hsps may be due to their different functions in reproduction and survival.
Insect immune defenses are carried out through humoral and hemocyte responses, with the former playing defensive role through the synthesis of antimicrobial peptides and the latter through encapsulation and phagocytosis [64]. In the present study, four antimicrobial peptides encoding gene (Defensin, Hymenoptaecin, Megourin, and Fungal protease inhibitor) were found in C. chuxiongica and all of them were downregulated considerably in mated females (Table 3). In addition, two hemocyte responses related enzyme-encoding genes, Phenoloxidase and Lysozyme, were also found downregulated in mated females (Table 3). Lysozyme is involved in immune defense by hydrolyzing the bacterial cell walls whereas phenoloxidase participates in the synthesis of oxidative free radicals and the defensive melanization [1]. A total of 21 immunity related KEGG pathways were enriched, with sixteen of them being enriched to downregulated DEGs and five being enriched to upregulated DEGs. These results have suggested that the immunity activity has been downregulated substantially in mated females.
Above results and discussion have suggested that mating is also an essential switch in C. chuxiongica, a species with facultative parthenogenesis, which induced significant trade-offs between reproduction and survivorship. As mentioned above, the male and female adults of C. chuxiongica can mate a few hours after eclosion and mated females start to lay fertilized eggs in the subsequent day after mating, as well as start to show egg protection behavior; unmated females (virgin) may start to lay unfertilized eggs through parthenogenesis a few days (>3 d) after eclosion [21]. Therefore, the pattern of gene expression changes after mating is consistent with the post-mating behavioral and physiological changes. Moreover, the present study also suggested that mating may also induce trade-offs within reproduction on fecundity vs. fertility in C. chuxiongica females, in which females may downregulate oogenesis and egg production to restrict fecundity but upregulate meiotic maturation in fully grown oocytes and fertilization-related pathways to favor the following egg fertilization. C. chuxiongica mostly occurs in pine forests at high altitude and barren soil, its larvae are obligate pine defoliator, which mainly feed on Pinus yunnanensis [21]. C. chuxiongica has a long (19 months) larval dormancy stage and adults do not feed [21], which may force adults to use limited resources optimally on reproduction and survival. In addition, C. chuxiongica females lay eggs on the surface of pine needles [23], which will facilitate larval feeding after hatching but also will likely to incur attack by natural enemies. C. chuxiongica females may thus have evolved egg protection behavior. The fecundity of C. chuxiongica is about 50 eggs per female [22]. Therefore, compared to other high fecundity insect species, such as S. litura that lay more than one thousand eggs per female [65] and show upregulation of fecundity related genes after mating [15], C. chuxiongica may have evolved a different reproductive strategy, i.e., restrict whole fecundity while ensuring higher egg fertilization and offspring survival. Future studies to clarify the hypotheses established in this study by using other techniques, such as RNAi and 2D electrophoresis/MS, will help to provide deeper insights in this field.
The present study and previous studies [21][22][23] on C. chuxiongica suggest that this species has evolved multiple physiological and behavioral strategies (such as facultative parthenogenesis, long larval diapause, and trade-offs on reproduction and survival) to adapt to its living environment.  Institutional Review Board Statement: Not applicable.

Data Availability Statement:
The transcriptome raw reads have been deposited to the NCBI SRA database (accession no.: SRR13991392~SRR13991407). Other data generated or analyzed during this study were included in this article and its Supplementary Information.