Transcriptional and Translational Relationship in Environmental Stress: RNAseq and ITRAQ Proteomic Analysis Between Sexually Reproducing and Parthenogenetic Females in Moina micrura

Moina micrura is a kind of small-bodied water flea within the family Moinidae. Similar to Daphnia, M. micrura could also switch its reproduction mode from parthenogenetic female (PF) to sexual female (SF) to adapt to the external environment. To uncover the mechanisms of reproductive switching in M. micrura, we used both RNA-Seq and iTRAQ analyses to investigate the differentially expressed genes (DEGs) and their protein products between SF and PF in M. micrura. A total of 1665 DEGs (702 up-regulated, 963 down-regulated) and 600 differentially expressed proteins (DEPs) (102 up-regulated, 498 down-regulated) were detected in SF. Correlation analyses indicated that 31 genes were expressed significantly differentially at both transcriptomic and proteomic levels, including 15 up-regulated genes and 16 down-regulated genes in SF. Meanwhile, our data also showed that 528 DEPs have discordant expression at transcript level, implying post-transcriptional (including translational) regulation. These top up-regulated genes and their protein products in SF were mainly grouped into the globin-related family, vitellogenin-related family, cuticle-related family, Hsp-related family and methyltransferases-related family, which were all involved in the reproductive switching in Daphnia. In contrast, a cluster of orthologous groups revealed that up-regulated genes and their protein products in PF were strongly associated with the metabolic process, which may be responsible for rapid population proliferation in M. micrura.


INTRODUCTION
Moina micrura is a kind of small-bodied crustaceans within the family Moinidae. It is referred to as water flea but smaller than its more well-known cousins: the larger Daphnia magna and the medium-sized Daphnia pulex. M. micrura is the very ideal live bait for fish larvae and economically important crustaceans, and could serve as an indicator for evaluating water quality in the aquatic ecosystem (Arauzo and Valladolid, 2003). Similar to other water fleas, M. micrura utilizes a special reproductive strategy known as cyclical parthenogenesis, in which parthenogenesis (asexual reproduction) and sexual reproduction can be switched when the environmental characteristics change (Hobaek and Larsson, 1990;Rojas et al., 2001). When the environmental conditions are ideal and appropriate, parthenogenesis is occurred to produce a large number of parthenogenetic embryos. When facing deteriorated environmental conditions (e.g., shortened day length, low or high temperature, lack of nutrients and overpopulation), M. micrura could undergo an unusual transition from asexual to sexual reproduction and produce resting eggs to live through the difficult days. This reproductive strategy is an adaptation to changes in external environment; parthenogenesis contributes to the rapid proliferation in the favorable seasons, while sexual reproduction results in the increase of the genetic diversity and survival through poor environmental conditions (Barton and Charlesworth, 1998;Toyota et al., 2016). Therefore, M. micrura is also a great model organism for elucidating the molecular mechanism of reproductive switching in cladocerans. However, little is known about the precise mechanism of reproductive switching in M. micrura.
Previous studies have demonstrated that the reproductive switching of cladocerans was closely related with the changes of environmental factors including the photoperiod, temperature, food, and population density. For example, excessive temperature could induce the emergence of male and produce a large number of resting eggs (Green, 1965;Gordo et al., 1994). Short light period could significantly increase the proportion of SF and males (Baer and Owens, 1999;Toyota et al., 2015). Food is another effective factor in the sex determination of cladocerans, and its quantity and quality have a direct effect on the reproductive switching of cladocerans (Arnold, 1971). Population density also plays an important role in the reproductive switching of cladocerans. The amount of resting eggs usually shows a positive correlation with population density (Banta and Brown, 1929). In natural waters, male and SF are usually found when the population density of cladocerans is high (Carvalho and Hughes, 1983). Besides these ecological factors, several key genes have also been detected to play a pivotal role in the regulation of reproductive switching of cladocerans. For example, Kato et al. isolated five important genes involved in reproductive switching in D. magna, such as DM-domain (Kato et al., 2008), transformer (Kato et al., 2010), and doublesex (Kato et al., 2011). In addition, basing on the EST sequence analysis, other reproductive-specific genes were also detected, such as doublesex/mab-3 domain, ubulin alpha chain, induction gene (CSP) and heat shock protein (HSP) gene (Xu et al., 2009). These genes, which are closely related to perceiving the changes of environment, may play a key regulatory role in parthenogenesis of cladocerans. However, the reproductive switching is a complicated process in cladocerans. Little is known about which genes regulate reproductive switching, and what the interactions between transcriptional and translational control of gene expression are.
High-throughput sequencing technologies empower us to analyze alterations of an organism's global molecular profile in response to environmental stresses. RNA sequencing (RNA-Seq) is a transcript quantification technique and has been widely and successfully applied to gene annotation, expression analysis and transcript profiling (Yang et al., 2011). As well, proteomics is a high-throughput approach to address gene functions that cannot be offered by genome sequences and is the most direct way to confirm the function of a gene (Ahsan et al., 2008). Isobaric tags for relative and absolute quantitation (iTRAQ) is the second generation of gel-free proteomic analysis that allows for a more accurate quantitation of proteins than low-throughput proteomic methods, such as traditional two-dimensional gel-based approaches. Recently, transcriptome analysis has been used to determine the DEGs between PF and SF in D. pulex (Toyota et al., 2015) and D. similoides (Zhang et al., 2016), respectively. However, as far as we know, there is no study to conduct proteomic profiling to examine the proteins involved in the reproductive switching under environmental stresses.
The aim of this study is to further investigate the key genes and their protein products involved in reproductive switching. To achieve this, using M. micrura as model, both RNA-Seq and iTRAQ technology are used to identify the differentially expressed mRNAs and proteins between SF and PF. These findings will provide important information on the molecular mechanism of the reproductive switching in M. micrura.

Sample Preparation
Moina micrura was isolated from the South Lake (Wuhan city, China), and the pure lines were established by monoclonal culture in our laboratory. Healthy parthenogenetic organisms were cultivated under a 12 h light/12 h dark photoperiod at 25 • C, and fed chlorella for 1 year. According to the biological characteristics of M. micrura, reproduction switching occurs when the population density reaches a certain level. Healthy SF and PF were collected and their reproductive state was confirmed, respectively, using OLYMPUS BA200 microscope (Figure 1), and put into cryogenic vials (100 per vial). Then samples in the cryogenic vials were washed three times with sterilized water. After that, water was drained and the samples were immediately frozen in liquid nitrogen and stored at −80 • C until further processing. In this study, three tubes of SF and PF were collected for transcriptomic and proteomic sequencing, respectively.

Total RNA Isolation and cDNA Library Construction
Total RNA was extracted from SF and PF using TRIzol reagent (Invitrogen, United States) according to the manufacturer's instructions, and genomic DNA was removed using DNase I (TaKara). The Nanodrop 2000 spectrophotometer was used to assess sample purity and RNA concentration, and the quality of RNA was analyzed on an Agilent 2100 Bioanalyzer using the RNA 6000 Nano kit (Agilent Technologies, Santa Clara, CA, United States). All samples displayed a 260/280 ratio greater than 2.0 and RNA integrity numbers (RIN) ≥7.5. Before the construction of library, ribosomal and viral RNA were removed and poly(A)+ mRNA was isolated with magnetic Oligo-dT beads (Invitrogen, United States), and then 10 µg of total RNA for each sample was used for the construction of libraries by using a TruSeq TM RNA sample preparation Kit from Illumina (San Diego, CA, United States), all performed according to the Illumina protocol. Thereafter, libraries were sequenced using Illumina HiSeq TM 2000 system at Majorbio Biotech Co., Ltd. (Shanghai, China). We chose a read depth of 600 million 150bp single end reads. All raw-sequence read data were deposited in NCBI Sequence Read Archive (SRA 1 ) with accession number SRP136759.

Transcriptome de Novo Assembly and Bioinformatics Analysis
The raw paired end reads were trimmed and quality controlled by SeqPrep 2 and Sickle 3 with default parameters. Then clean data from the samples (SF and PF) were used to do de novo assembly with Trinity 4 and all the assembled full-length sequences were named unigenes (Grabherr et al., 2011). All the unigenes were predicted and used for Blastx search and annotation against the NR, Swiss-Prot, COG, and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases (Camacho et al., 2009). Functional annotation by gene ontology (GO) terms was analyzed by using Blast2GO software (Conesa et al., 2005), and GO functional classification for unigenes was analyzed by using WEGO software (Ye et al., 2006). The read counts were further normalized into FPKM values. The fold changes were calculated by using RSEM software v 1.2.7 (Li and Dewey, 2011), and the DEGs were analyzed using the R Bioconductor package, edgeR (Robinson et al., 2010). The P-value set the threshold for the differential gene expression test. The threshold of the P-value in multiple tests was determined by the value for the false discovery rate (FDR) (Benjamini and Yekutieli, 2001). DEGs were screened with a cut-off conditions of fold change (FC) > 2.0 and FDR ≤ 0.001.

Protein Preparation and iTRAQ Labeling
The M. micrura samples were sonicated in ice with appropriate amount of lysis buffer (8 M urea, 0.3% SDS) and protease inhibitors (Thermo, United States). The proteins were reduced with 10 mM DTT (final concentration) at 56 • C for 1 h and alkylated by 55 mM IAM (final concentration) in the darkroom for 1 h. The reduced and alkylated protein mixtures were precipitated by adding 4× volume of chilled acetone at −20 • C overnight. After centrifugation at 4 • C, 30,000 g, the pellet was dissolved in 0.5 M TEAB (Applied Biosystems, Milan, Italy) and sonicated in ice. After centrifugation at 30,000 g at 4 • C, an aliquot of the supernatant was taken to determine the protein concentration with a 2-D Quant Kit (GE Healthcare). The proteins in the supernatant were kept at −80 • C for further analysis. More details for iTRAQ labeling are provided in the Supplementary Supporting Information.

Analysis of the Identified Proteins
The raw data obtained from LC/LC-MS/MS analysis were processed using Proteome Discoverer TM Software 2.1 (Thermo, United States). The search parameters included uniprot-organism6669-D. pulex as protein database, Iodoacetamide as Cys alkylation, oxidation (M), acetyl (Protein N-Terminus) and iTRAQ8plex (Y), as the potential dynamic modifications, and iTRAQ8plex (K), iTRAQ8plex (N-Terminus), as the static modifications. Functional annotations of the proteins were conducted with the Blast2GO (Gene Ontology) program against the non-redundant protein database (NR; NCBI), and the differentially annotated proteins were further assigned to KEGG database. The significant p-value of the difference between the samples was calculated by using the t-test function in the R language. For identifying significantly up-or down-regulated proteins, the threshold values of SF/PF or PF/SF ratios were ≥1.50 or ≤0.67 (≥1.5-fold) and the p-value < 0.05 in both two iTRAQ analyses.

Correlation Analysis Between Transcriptome and Proteome Data
To assess the potential relevance of quantitative information between mRNA and proteins, the cut off values (DEGs: |Fold change (FC)| > 2.0 and FDR ≤ 0.001; DEPs: p-value < 0.05 and |FC| > 1.50) were used to screen the subsets of mRNA and proteins with apparent expression. Then, we used RNA-seq data as a searchable database, all identified protein sequences were analyzed and queried with the RNA-seq data. To further dig out the information from proteomes, the correlation analysis was conducted between DEPs and the whole transcript database.

Real-Time Quantitative PCR (RT-qPCR) Validation
Total RNA of M. micrura was extracted from SF and PF samples with TRIzol reagent (Invitrogen, United States), respectively. The Nanodrop 2000 spectrophotometer was used to assess sample purity and RNA concentration, and 10 µg total RNA from PF and SF were reversely transcribed by PrimeScript RT reagent kit (Takara, Dalian, China), respectively. The RT samples obtained were subjected to qPCR using a ABI 7500 real-time PCR system (Applied Biosystems, United States) with specific primers for DEGs, respectively (see Supplementary Table S1 for primer sequences and PCR condition). In these experiments, the comparative Ct method (2 − Ct ) of relative quantification (Livak and Schmittgen, 2001) was used to analyze the real-time quantitative PCR data, using QuantStudio TM Real-Time PCR Software, Version 1.2. Parallel real-time PCR for GAPDH was also conducted to serve as the internal control. The sizes of the amplified products were confirmed through gel electrophoresis.

Data Analysis
All the data were subjected to statistical analysis with SPSS 19.0 software (IBM Corp, Armonk, NY, United States). Physiological data were analyzed using Student's t-test and the Tukey method for one-way ANOVA at 95% confidence. Real-time RT-PCR results were evaluated with the Student's t-test. The data were expressed as the means ± standard errors of the mean (SEM). Differences between experimental groups were considered as significant at P < 0.05.

RNA-Seq Data and Identification of DEGs
To survey the key genes in reproductive switching, a high-throughput transcriptome was used to compare mRNA expression profiles between SF and PF. In the present study, the cDNA libraries for SF and PF were established and sequenced. The numbers of clean reads obtained from SF and PF of M. micrura were 77,923,354 and 68,509,476, respectively, and the read length was 145 bp and all Q30 were higher than 96%. After quality filtering, 69,857 unigenes were acquired, in which the average length was 995 bp. For these unigenes, 20,131 (28.8%), 24,282 (34.8%), 8,807 (12.6%), 24,543 (35.1%), and 13,213 (18.9%) were annotated to NR, Swiss-prot, KEGG, COG, and GO database, respectively. Furthermore, fragments per kilobase of transcript per million fragments (FPKM) analysis showed that total 1665 unigenes (FC > 2.0 and FDR < 0.001) were differently expressed between SF and PF. Compared to PF, 963 genes were up-regulated and 702 genes were down-regulated in SF (Figure 2A).

GO and KEGG Pathway Enrichment Analysis of DEGs
To understand the distribution of DEGs at a macro level, the set of DEGs between the SF and PF were mapped in accordance with the GO terms and the KEGG pathway, respectively. GO enrichment analysis showed that the up-regulated genes in SF were highly related to cuticle (e.g., carbohydrate binding, polysaccharide binding, chitin binding, chitin metabolic process, and structural constituent of cuticle), hemoglobin (e.g., iron ion binding, heme binding and tetrapyrrole binding), oxidative stress (e.g., oxidoreductase activity and oxidation-reduction process), and Vg (e.g., lipid transporter activity) ( Figure 3A). Yet, most down-regulated genes in SF were grouped in the metabolic process (e.g., macromolecule/heterocycle/aminoglycan metabolic process), primary metabolic process (e.g., protein metabolic process, proteolysis and carbohydrate/polysaccharide metabolic process), catalytic activity (e.g., hydrolase activity and transferase activity), and peptidase activity (e.g., endopeptidase activity and serine-type peptidase/endopeptidase activity) ( Figure 3A). KEGG pathway enrichment analysis demonstrated that most up-regulated genes in SF were widely related with organismal systems (e.g., longevity regulating pathway, PPAR signaling pathway etc.) and metabolism (e.g., citrate cycle, pentose and glucuronate interconversions etc.). In addition, organismal systems (e.g., protein digestion and absorption, lysosome etc.) and environmental information processing (e.g., neuroactive ligand-receptor interaction and sphingolipid signaling pathway) were found to be the most-enriched pathways in down-regulated genes ( Figure 4A).

Protein Identification and Quantitation
Given the advantage of high resolution and high-quality accuracy, the Q-Exactive mass spectrometer was used to detect the DEPs between SF and PF. The iTRAQ analysis of M. micrura global proteome revealed 2352 protein hits in Proteome Discoverer. The error in the matching of the peptide segment found in the database was below 0.05Da ( Figure 5A). The molecular weight of most proteins was shown in Figure 5B. As shown in Figure 5C, the peptide length was around in 9-20 aa, among which 9-15 aa interval was the peak area. The numbers of peptides identified into the proteins was exhibited in Figure 5D, which showed that most of protein sequence coverage was at 1-40%.

Identification and Enrichment Analysis of DEPs
A total of 600 DEPs were identified between SF and PF (FC ≥ 1.5 and p-value < 0.05), among which included 102 up-regulated proteins and 498 down-regulated proteins in SF ( Figure 2B). GO analysis showed that most DEPs were grouped into cellular component (e.g., cell part, macromolecular complex, organelle etc.), binding, cellular process, and metabolic process ( Figure 3B). Additionally, most up-regulated DEPs in SF were preferentially related to  developmental process and transcription regulator activity ( Figure 4A). In addition to the GO analysis, KEGG pathway enrichment analysis was also used to further elucidate the putative functions of these DEPs. The results demonstrated that up-regulated proteins in SF were mainly involved in the genetic information processing, especially translation process relating to RNA transport, ribosome and ribosome biogenesis in eukaryotes, while only ribosome was found to be the most-enriched pathways in down-regulated genes in SF ( Figure 4B).

Correlation Analysis Between Transcriptome and Proteome Data in SF and PF
As shown in Figure 6A, a total of 2352 genes were identified at both protein and transcript levels. However, only 72 DEPs were also detected in DEGs (Figure 6B). The Pearson correlation coefficient was −0.1542 (p < 0.01, Figure 6C), which indicated that the relative expression quantity of the genes at transcript level correlated with protein level positively but poorly. In addition, most of the genes were concentrated in the center of the coordinate, indicating that these genes were consistent at transcriptome and proteome levels. The GO analyses for both DEGs and DEPs demonstrated that the major GO terms were focused on "binding, " "metabolic process, " "catalytic activity, " "cellular process, " and "single-organism process, " respectively ( Figure 6D). Furthermore, within the most significantly enriched pathways, most DEGs were amassed in "metabolic pathways, " "ribosome, " and "protein digestion and absorption, " while the major-enriched pathways in DEPs were "ribosome, " "ribosome biogenesis, " "RNA transport, " and "spliceosome" (Figure 6E).
The comparison between transcriptome and proteome displayed that the expression patterns of proteins showed a general correlation with their mRNA expression at transcript level. Furthermore, our data also revealed significantly discordant regulation between transcriptional and translational levels. We found 59 up-regulated (FC > 1.5) (Supplementary Table S11, with function annotation) and 469 down-regulated proteins in SF (Supplementary Table S12, FC > 1.8) at proteomic level but in absence of changes in mRNA expression at transcript level ( Figure 6F). In addition, 4 up-regulated (Supplementary  Table S13) and 11 down-regulated proteins (Supplementary Table  S14) showed abundance changes in the opposite direction at transcript level ( Figure 6F).
Besides transcriptional regulation, it is clear that gene expression is also modulated at post-transcriptional levels.

DISCUSSION
Similar to Daphnia, M. micrura could also switch their reproduction mode from parthenogenesis to sexual reproduction to adapt to the external environment. However, little is known about the molecular mechanisms of reproductive switching. In the present study, based on the transcriptomic and proteomic analyses, we tried to investigate DEGs and DEPs in M. micrura (SF vs. PF) and identify the critical genes involved in reproductive switching.
Gene ontology and KEGG analyses demonstrated that up-regulated and down-regulated genes in M. micrura (SF vs. PF) had different enrichment terms. The most frequently enriched terms in SF were cuticle (e.g., carbohydrate binding, polysaccharide binding, chitin binding, chitin metabolic process, and structural constituent of cuticle), hemoglobin (e.g., iron ion binding, heme binding and tetrapyrrole binding), and oxidative stress (e.g., oxidoreductase activity and oxidation-reduction process). Chitin is involved in cuticle formation in insects and is the most important barrier tissues against any environmental stresses (Pesch et al., 2015). Hemoglobin could allow the organisms to survive in low oxygen stress (Pirow et al., 2001). These results, taken together, suggest that the up-regulated genes in SF play an important role in producing a protective cuticle structure to resist unfavorable conditions, and increasing oxygen transport and storage in respond to hypoxia. In contrast, most basic metabolic processes (e.g., "metabolic process, " "primary metabolic process, " "macromolecule metabolic process, " and "protein metabolic process") were the most frequently enriched groups in PF. As well known, PF can produce a large number of parthenogenetic embryos which could develop directly into PF without fertilization (Arbačiauskas, 2004), these results imply that up-regulation of various metabolic processes in SP might play an important role in nutrient uptake and further contribute to the growth of neonates.
Correlation analysis between transcriptomics and proteomics exhibit that the enrichment of proteins in different GO terms are consistent with the transcriptome results. However, the pathways enriched in proteins are especially grouped in genetic information processing, while various metabolic processes are dominated in transcriptome level. In addition, our present study demonstrates that 31 genes are significantly differentially expressed at both transcriptomic and proteomic levels, including 16 up-regulated genes and 15 down-regulated genes in SF. Although transcriptional and translational regulation are generally coordinated, we also identify a lot of genes that oppose this trend. A total of 543 DEPs in SF are showed differential expression pattern in transcript level, implying the involvement of post-transcriptional regulation. Previous study has also demonstrated that post-transcriptional processing determined steady-state protein levels (Chen and Harmon, 2006). To further explore the effects of post-transcriptional regulation in reproductive switching, we also dig out some proteins that went against the overall trend of concordant regulation at mRNA and protein levels, including Dmagvtg1, Ferh, Cat-1, Ndufaf5, Pxt, cuticle protein 7 (N/A c ), Cygb2, endocuticle structural glycoprotein SgAbd-2 (N/A l ), Cyp301a1, Pcp20, Ecu02_0100, ferritin heavy chain (N/A m ), and Apod. These post-transcriptionally regulated genes may play a more important role in reproductive switching than those transcriptionally regulated genes in M. micrura. Finally, the critical genes and their protein products involved in reproductive switching are analyzed as follows.

Oxygen-Transport Genes
The globins are a superfamily of heme-containing globular proteins, involved in binding and/or transporting oxygen, specifically for iron ion binding, heme binding, and oxygen transport (Baldi et al., 1994;Culbertson, 2010). In the present study, the globin members in SF were found to contain one significantly up-regulated gene (hemoglobin) at transcriptomic and proteomic levels and three significantly . The x-axis shows the mean value of the expressed protein in the group, and the y-axis indicates the mean value of the FPKM of the transcript in the same group. The expression of the transcript and the protein is taken as the logarithmic value. Each point represents a protein and its associated transcripts, rho in the legend represents the pearson correlation coefficient, p is the correlation test p-value. when rho > 0, it is called negative correlation; rho < 0, said positive correlation; rho = 0, called zero correlation, that is, no correlation; | rho| larger, the greater the correlation between the two groups. up-regulated genes (Cygb2, Ferh and ferritin heavy chain) at post-transcriptional level. The hemoglobin is one of the most important globin protein, which plays a key role in oxygen transport. In invertebrates, hemoglobin concentrations are low in most cases. However, hypoxia can cause a significantly increase in the hemoglobin concentrations of freshwater zooplankton   organisms, which indicates that the induction of hemoglobin synthesis could improve the tolerance of the daphnids to low environmental dissolved oxygen levels, and allows the organisms to survive under this stress (Pirow et al., 2001). Given that environmental oxygen levels are relatively low in the occurrence of sexual reproduction in cladoceran, it is further demonstrated that the dissolved oxygen levels may be an important environmental factor for the reproduction switching of cladoceran. Similar to hemoglobin, the Cygb2 existing in proteomic data also belongs to the globin family, which is significantly up-regulated in SF. Since previous studies have demonstrated that Cygb2 could be up-regulated under hypoxic conditions and played a key role in oxygen-requiring reactions (Fago et al., 2004;Schmidt et al., 2004), we speculate that Cygb2 may be also involved in intracellular oxygen storage or transfer in M. micrura. In addition to hemoglobin and Cygb2, Ferh and ferritin heavy chain also played an important role in ferric iron binding and transport. Larade and Storey (2004) reported that the concentration of ferritin could be rapidly increased in response to stresses such as anoxia, which implies that ferritin is an acute phase protein for hypoxia (Beck et al., 2002). As we know, environmental signals could regulate a variety of key physiology processes in animals, the globin-related genes involved in oxygen transport were all up-regulated in SF in response to hypoxia. These results, taken together, indicate that low oxygen content may lead to reproductive switching of M. micruras. On the other  hand, the high expression levels of the globin-related genes in SF may also indicate that SF needs more oxygen to finish sexual reproduction and to produce resting eggs.

Vitellogenin-Related Genes
Vitellogenin (Vg/Vit-2) and vitellogenin-2 (Vg2) are both the members of the Vg family that play a key role in the formation of yolk proteins in cladocerans (Kato et al., 2004). Vg is the major precursor of the egg-yolk proteins, vitellins (Vn), which provide energy for embryonic development in oviparous organisms (Matozzo et al., 2008). In the present study, the expression levels of Vg gene and vitellogenin-2 gene were both significantly higher in SF compared to PF. Similar results were also reported in D. similoides (Zhang et al., 2016). In addition, we also found that Dmagvtg1 was significantly up-regulated in SF at post-transcriptional level. Previous study has reported that Dmagvtg1 encoded Vg, and the precursor to Vg was the most abundant polypeptide in D. magna eggs (Kato et al., 2004).
These results indicate that the SF could produce more vitellins to package into the resting eggs for embryonic development. In addition, similar to Vg, the F40A3.3 was also highly expressed in SF, which was annotated as lipid transport or lipid transporter activity, and played an important role in lipid binding. Hence, it may be also involved in the synthesis of egg-yolk proteins in SF.

Cuticle-Related Genes
In the present study, seven genes and their protein products (cuticle protein 6, cuticle protein 7, endocuticle structural glycoprotein SgAbd-2, endocuticle structural glycoprotein Sgabd-8, Pcp20, Pcp27, and Cyp301a1), related to the synthesis of cuticle or chitin and the formation of cuticle, were significantly up-regulated in SF compared to PF. Of these genes, the expression levels of cuticle protein 7, endocuticle structural glycoprotein SgAbd-2, Pcp20 and Cyp301a1 were found significantly higher in SF compared to PF at post-transcriptional level. The cuticle proteins and chitin are both involved in the formation of cuticle, which is the most important barrier against the environmental stresses in insects (Neville et al., 1976;Liu et al., 2014). In cladocerans, the cuticle also consists of cuticle protein and chitin, which can withstand adverse conditions of external environment (Repka et al., 1995;Liu et al., 2014). A recent study demonstrated that D. magna could develop an array of morphological changes in the presence of Triops cancriformis, including changes of carapace morphology and cuticle hardening (Otte et al., 2014). Besides, another study suggested that Cyp301a1 was involved in the formation of the adult cuticle in Drosophila melanogaster (Sztal et al., 2012). Taken together, the high expression of cuticle-related genes and proteins indicates that these genes could help SF to undergo a series of corresponding changes in cuticle structure as a response to the adverse external environmental conditions. Furthermore, the SF needs to produce resting eggs to survive in unfavorable environmental conditions, and those resting eggs are packed by a black ephippium with a thick cuticle. These findings indicate that the up-regulated cuticle-related genes in SF might also play a key role in the formation of ephippium.

Hsp-Related Genes
In this study, three members of heat shock protein (Hsp)-related family (Hsp-16.2, ECU02_0100, and Cryaa) were also highly detected in SF. Hsp-16.2 encodes a small heat shock protein (sHsp) which plays an important role in protecting proteins against stress-induced aggregation (Pozsgai et al., 2007). Previous study has demonstrated that the anti-apoptotic effect of Hsp16.2 was mediated by the activation of Hsp90, with which Hsp16.2 binds (Bellyei et al., 2007a). Furthermore, overexpression of Hsp16.2 can increase lipid raft formation, thus helping to stabilize the plasma membrane (Bellyei et al., 2007a,b). ECU02_0100 is an Hsp70-related protein, which is a member of Hsp70s. Hsp70s play a vital role in stress tolerance and survival under adverse conditions (Mikulski et al., 2011). They can also promote the correct folding of newly synthesized proteins (Mayer and Bukau, 2005). Cryaa can encode alpha-crystallin A chain (αAC), which belongs to the alpha-crystallin (αC). Similar to Hsp16.2, αC has obvious sequence homology with the amino acid sequence of 90-100 residues at the C terminus of the sHSP family (Klemenz et al., 1991). It is also demonstrated that the expression of αC increase significantly under heat stress and oxidative stress, which suggests that the αC is also a molecular chaperone (similar to sHSP) in response to multiple environmental stresses. These findings, as a whole, indicate that the high expression of three Hsp genes in SF can protect the M. micrura against various environmental stress stimuli (e.g., low or high temperature, hypoxia) by maintaining cellular homeostasis. In addition, these findings also indicate that Hsp proteins may be the key proteins that deal with environmental stresses and maintain the homeostasis, which play an important role in the success of reproductive switching.

Methyltransferase-Related Genes
Methyltransferases are a large group of enzymes including protein methyltransferases, DNA/RNA methyltransferase, natural product methyltransferases, and non-S-Adenosyl methionine (SAM) dependent methyltransferase. For the protein methyltransferase, previous studies have demonstrated that farnesoic acid O-methyltransferase was a major contributing enzyme in the synthesis of methyl farnesoate in decapod and branchiopod crustaceans (Ruddell et al., 2003;Li et al., 2010). Recent study has reported that methyl farnesoate played a key role in regulating the reproductive switching in water flea (LeBlanc and Medlock, 2015). These findings suggest that the up-regulated genes encoding methyltransferases (Mettl10 and Ndufaf5) in SF may also be related to the reproductive switching in M. micrura. On the other hand, the DNA methyltransferase can catalyze the transfer of a methyl group to DNA. DNA methylation, a key component of genetic regulation, occurs primarily at the 5-carbon of the base cytosine, forming 5 methylcytosine, and does not entail a change in DNA sequences (Lan et al., 2010). The presence of DNA methylation has been detected in D. magna, which suggested that potentially epigenetic effect might occur in this species (Vandegehuchte et al., 2009). In this study, the Mettl10 and Ndufaf5 gene products were both highly detected in SF, which suggested that the DNA methylation might also be related to the reproductive switching in M. micrura.

Genes Encoding Some Peptidases, Peroxidase and Lipocalin
Our iTRAQ results showed that the genes encoding some peptidases, peroxidase and lipocalin were all highly detected in SF compared to PF. Cathepsin L-like proteinase was up-regulated in Helicoverpa armigera during larval molting (Wang et al., 2010). Chorion peroxidase was involved in the formation of a rigid and insoluble egg chorion in Aedes aegypti (Li et al., 1996). We speculate that both cathepsin L-like proteinase and chorion peroxidase may be involved in the development of embryos and juveniles under adverse condition in M. micrura. Apolipoprotein D was an acute response protein with a protective and was involved in the mechanisms regulating protection from oxidative stress (Ganfornina et al., 2008). In the present study, apolipoprotein D was also significantly up-regulated in SF, which may contribute to protecting M. micrura from oxidative stress. These findings, taken together, suggest that the up-regulated cathepsin L-like proteinase, chorion peroxidase and apolipoprotein D in SF may play a role in reproductive switching in M. micrura.
In summary, we identified several critical DEGs and DEPs involved in reproductive switching in M. micruras basing on the analysis of transcriptome and proteome. We found post-transcriptional regulation may play an important role in reproductive switching in M. micrura. Those top up-regulated genes and their protein products in SF mainly belonged to globin-related family, Vg-related family, cuticle-related family, Hsp-related family and methyltransferase-related family, which indicates that these genes may play a key role in the reproductive switching in M. micrura. In contrast, the clusters of orthologous groups reveal that up-regulated genes in PF are strongly associated with the metabolic process, which may be responsible for rapid population proliferation of PF. Our comprehensive transcriptome and proteome data provide a valuable resource and offer novel insights into understanding post-transcriptional regulation in reproductive switching in M. micrura. Furthermore, this study also provides many key candidate genes and their protein products related to reproductive switching for further functional studies in M. micrura.