Identification and Integrated Analysis of MicroRNA and mRNA Expression Profiles During Agonistic Behavior in Chinese Mitten Crab (Eriocheir sinensis) Using a Deep Sequencing Approach

As a commercially important species, the Chinese mitten crab (Eriocheir sinensis) has been cultured for a long time in China. Agonistic behavior often causes limb disability and requires much energy, which is harmful to the growth and survival of crabs. In this paper, we divided crabs into a control group (control, no treatment) and an experimental group (fight, agonistic behavior after 1 h) and then collected the thoracic ganglia (TG) to extract RNA. Subsequently, we first used a deep sequencing approach to examine the transcripts of microRNAs (miRNAs) and messenger RNAs (mRNAs) in E. sinensis displaying agonistic behavior. According to the results, we found 29 significant differentially expressed miRNAs (DEMs) and 116 significant differentially expressed unigenes (DEGs). The DEMs esi-miR-199a-5p, esi-let-7d, esi-miR-200a, and esi-miR-200b might participate in the regulation of agonistic behavior by mediating neuroregulation and energy metabolism. Focusing on the transcripts of the mRNAs, the renin–angiotensin system (RAS) Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway might be involved in the regulation of agonistic behavior through glucose metabolism as this pathway was significantly enriched with DEGs. Besides, an integrated analysis of the miRNA and mRNA profiles revealed that the retinoid X receptor (RXR) was also involved in visual signal transduction, which was important for agonistic behavior. In addition, four vital agonistic behavior-related metabolic pathways, including the cAMP signaling, MAPK, protein digestion and absorption, and fatty acid metabolism pathways, were significantly enriched with the predicted target unigenes. In conclusion, the findings of this study might provide important insight enhancing our understanding of the underlying molecular mechanisms of agonistic behavior in E. sinensis.


INTRODUCTION
Chinese mitten crab (Eriocheir sinensis) has been loved by Chinese people and cultured in many areas of China for a long time (Cheng et al., 2018). Agonistic behavior is a common phenomenon in crustaceans and usually decreases the crabs' integrity, survival, and growth (Laranja et al., 2010;Harlıoglu et al., 2013;Romano and Zeng, 2016). In our previous studies, we reported the agonistic behavior of E. sinensis and further investigated the regulation of serotonin (5-HT) and dopamine (DA) in agonistic behavior through the cAMP signaling pathway (Pang et al., 2019a;Yang et al., 2019). First, we had divided the agonistic behavior in E. sinensis into three stages: approach, contact, and fight; and then, we found a sign in the concentration of 5-HT, DA, cAMP, and PKA after agonistic behavior (Yang et al., 2019). Subsequently, we injected the 5-HT, DA, or cAMP to crabs and observed the effect of agonistic behavior after injection. The results showed that 5-HT and DA could regulate the behavior through the cAMP-PKA signaling pathway (Pang et al., 2019a). However, the mechanism underlying agonistic behavior in E. sinensis is unclear, and research investigating this topic is lacking. An understanding of the regulatory mechanisms underlying the switching of agonistic behavior is essential for elucidating the behavioral plasticity of animals and ultimately decreasing the occurrence of agonistic behavior.
MicroRNAs (miRNAs) constitute a class of non-coding RNAs with a length of 18-22 nt that can regulate gene expression and usually play an important role in many biological processes (Bartel, 2004;. According to reports, miRNAs could mediate neurogenesis at many key steps in vertebrate and invertebrate (Li and Jin, 2010) and regulate animal behaviors through nervous systems (Lane et al., 2018). For example, miR-375 and miR-200b could regulate attraction/aversion behavior via the dopaminergic and GABAergic systems in amphibians (Dulcis et al., 2017). Another study showed that the overexpression of let-7d perfected the anxiolytic and depressant-like behaviors by targeting the DA D3 receptor (Bahi and Dreyer, 2018). Researchers have also found that miR-96, which could inhibit the 5-HT1B receptor, affected aggressive behavior in mice (Jensen et al., 2009). Although many reports described the effects of miRNAs on animal behaviors, studies investigating the effects of miRNAs on aggressive behavior in crustaceans were lacking. Combined with our previous researches, 5-HT and DA as two neurotransmitters played key roles in nervous systems and regulated the agonistic behavior in E. sinensis. So to research the relationship between miRNAs and nervous systems would undoubtedly be of great help to understand the regulatory mechanisms of agonistic behavior. In addition, agonistic behavior usually was accompanied by energy metabolism (Briffa and Elwood, 2002). We also found that the glucose level has a significant upregulation after the fight in E. sinensis (Yang et al., 2019). There were reports that miRNAs could regulate the energy metabolism in animals; for example, miR-199a-5p could participate in the regulation of glucose metabolism (Wu et al., 2015;Yan et al., 2017;Esteves et al., 2018). Therefore, the relationship between miRNAs and energy metabolism was valued in research on agonistic behavior. The central nervous system (CNS) plays a major role in regulating several behaviors in some invertebrates, including agonistic behavior (Momohara et al., 2018), reproduction behavior (Thongbuakaew et al., 2019), and phototactic behavior (Rivetti et al., 2018). Besides, the thoracic ganglia (TG) is an important nerve tissue that participates in the formation of the CNS in crustaceans (Mulloney and Smarandache-Wellmann, 2012). The TG is among the many places where many neurotransmitters transfer biological information, such as DA, 5-HT, octopamine (OA), and tyramine (TA) (Vomel and Wegener, 2008). In our previous studies, we had done a lot of research on TG, including in vitro culture and RT-qPCR (Pang et al., 2019a). Based on the measurements of the TG, we found that 5-HT and DA could regulate the agonistic behavior of E. sinensis through the TG (Pang et al., 2019a;Yang et al., 2019). Therefore, investigating the TG was suitable for analyzing the mechanism underlying agonistic behavior in E. sinensis.
In this study, we obtained the transcriptomes of miRNAs and mRNAs after agonistic behavior by using Illumina HiSeq TM 2000 technology. Also, we analyzed the functions of the significant differentially expressed miRNAs (DEMs) and mRNAs in the TG and performed an integrated analysis of their relationships. We aimed to reveal the underlying molecular mechanisms involved in the regulation of agonistic behavior in E. sinensis.

Animals and Sampling
Male Chinese mitten crabs (E. sinensis) with a body weight of 25.56 ± 4.73 g were obtained from the Shuxin crab base in Chongming Island, Shanghai (China). The crabs were maintained in separate opaque tanks of 29.0:18.0:19.5 cm (length:width:height) for at least 7 days under single-rearing conditions before the behavioral experiments, and the tanks were filled with thoroughly aerated freshwater to a depth of 12 cm by a circulating system. The intact crabs were reared individually to avoid social contact. A basal diet was used to feed the crabs once daily from 18:30 to 19:00.
The crabs were divided into a control group (control) and an experimental group (fight). In the experimental group, two crabs with a weight difference of 1-4% were paired in a new tank of 20.0:15.5:19.5 cm (length:width:height). The tank was filled with 10 cm of water and divided into equal halves by an opaque partition. Before each pairing, each pair of crabs was placed on either side of the partition. After 10 min, the partition board was removed, and the agonistic behavior of the crabs was observed for 1 h using a high-definition camera (H.264 DVR) based on our previous research (Yang et al., 2019). We ensured that agonistic behavior occurred within 1 h, and then, the TG were collected for deep sequencing. In the control group, the crabs were placed on either side of the partition for 70 min, and then, their TG were harvested.

RNA Extraction and Quality Check
The total RNA was extracted using a TRIzol reagent (Sangon Biotech Co., Ltd.) according to the manufacturer's instructions. After RNA extraction, every two crabs which were collected from one pair were mixed as one sample to subsequent detection. Qubit 2.0 (Invitrogen, United States) was used to determine the concentration and purity of the total RNA. The integrity of the RNA and genome contamination were examined by agarose gel electrophoresis.

Small RNA Library Construction and Illumina Sequencing
Based on the above results, three samples from the fight group (named F1, F2, and F3) and another three samples from the control group (named C1, C2, and C3) were used to construct six small RNA libraries. In each library, a 6 µl mixture contained approximately 1 µg or more total RNA, and 1 µL 3 adapter and DEPC H 2 O were used to prepare the RNA mix; then, the mixture was incubated at 72 • C for 2 min. Subsequently, 1.8 µl 10 × T4 RNA ligase buffer, 3.5 µl 25 mM MgCl 2 , 2.2 µl PEG8000 (50%), 0.5 µl RNase inhibitor, and 1 µl ligase 2 were added to the RNA mix. Then, 15 µl of the mixture was incubated at 22 • C for 1 h, 1 µl RT primer was added to the above mixture, and the reaction conditions used were as follows: 75 • C for 5 min; 37 • C for 30 min, and 25 • C for 15 min. To ligate with the 5 adapter, 16 µl of the above hybrid production was opened with a secondary structure at 70 • C for 2 min, and then, 1 µl ATP (10 mM), 1 µl 5 adapter, and 1 µl ligase 1 were added. The obtained reaction solution was incubated for 1 h at 20 • C. Then, the adaptor-ligated small RNAs were reverse transcribed to create the complementary DNA (cDNA) constructs. Subsequently, the cDNA constructs were amplified by PCR using an Illumina small RNA primer set and Phusion polymerase under the following conditions: 95 • C for 5 min; 15 cycles of 94 • C for 10 s, 55 • C for 10 s, and 72 • C for 15 s; 72 • C for 5 min; and 4 • C for ∞. The productions of PCR were detected by 12% polyacrylamide (PAGE). The reads were obtained with a single end sequencing pattern at Sangon Biotech Co., Ltd. (Shanghai, China) using Illumina HiSeq TM 2000.

MRNA Library Construction and Sequencing
Simultaneously, we also selected three samples from the fight group (named F4, F5, and F6) and another three samples from the control group (named C4, C5, and C6) to construct six mRNA libraries. The samples used for the mRNA transcriptome analysis were prepared according to a VAHTS TM mRNA-seq V2 Library Prep Kit for Illumina R (Vazyme Biotech Co., Ltd., Nanjing, China). Poly-T oligolinked magnetic beads were used to purify the PolyA mRNA from the total RNA, and the intact mRNA was broken into fragments with a bead washing buffer and a metal bath. The aforementioned mRNAs were used as templates to synthesize the first-strand cDNA. Then, the second-strand cDNA was synthesized using the second-strand buffer and second-strand enzyme mix. Subsequently, the cDNA was end-repaired, a base was added to the 3 end, and the cDNA was amplified by PCR. Finally, the constructed cDNA libraries were sequenced with Illumina HiSeq TM 2000.

Analysis of miRNA and mRNA Sequence Data
The adapter sequences (sequence: TGGAATTCTCGGGTGCC AAGGAACTC) were removed from the raw data using the cutadapt software. In addition, the low-quality reads (a base quality less than Q20, the average quality of four consecutive bases less than Q20 and reads shorter than 17 nucleotides) were filtered using trimmomatic (Bolger et al., 2014). Finally, the clean reads of 17-35 nt were used in the subsequent analysis.
For the animal miRNAs, miranda (Betel et al., 2008) was used to predict the target genes. The predicted miRNA targets were based on a threshold value for the parameters as follows: S ≥ 150, G ≤ −30 kcal/mol. The target genes were predicted based on the E. sinensis mRNA transcriptome sequence obtained in this study as a reference genome.
Regarding the mRNA sequencing, raw reads from Illumina HiSeq TM 2000 may contain sequencing primers and low-quality sequences, which can affect the analytical quality. Therefore, the raw reads were cleaned through the following three steps: (a) linger sequences were discarded; (b) Q < 20 (Q = −10 log 10 E) bases were removed; and (c) read lengths shorter than 25 bp were discarded. Then, the clean reads were used for the de novo assembly using Trinity. Regarding the annotation, the assembled final unigenes with screening condition were annotated using the NCBI nucleotide sequence (NT) and NCBI protein nonredundant (Nr) databases. The unigenes were also classified using the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases.

Differential Expression Analysis of miRNAs and mRNAs
The expression levels of the miRNAs were determined by using RPM (normalized to determine the expression in transcripts per million). To identify the significant DEMs between two groups, the read counts of the samples were analyzed as input data using edgeR (Xiaobei et al., 2014). The screening conditions were as follows: p-value < 0.05 and | logFoldChange| > 1. The expression levels of the genes were determined by using the TPM (transcripts per million) method. An analysis of the differential expression levels across the samples was performed using DEGseq, and the significant differentially expressed genes (DEGs) were identified (q-value < 0.05, | FoldChange| > 2). Then, the significantly enriched terms were obtained by mapping the DEGs to the corresponding KEGG pathways.

Quantitative Real-Time PCR-Verified miRNAs
Five DEMs were randomly chosen for a quantitative realtime PCR (RT-qPCR) analysis to validate the accuracy of the RNA-seq results. The total RNA (500 ng), which was returned after sequencing (three samples per group), was used as a template to synthesize first-strand cDNA with a miRNA 1st Strand cDNA Synthesis Kit (by stem loop: 5 -GTCGTATCCAGTGCAGGGTCCGAGGTA TTCGCACTGGA TACGAC-3 , Vazyme Biotech) according to the manufacturer's instructions . Specific primers were designed based on the DEMs with a miRNA design software (Vazyme designed) ( Table 1).
The RT-qPCR reaction system contained 5 µl of 2 × miRNA Universal SYBR qPCR Master (Vazyme Biotech), 1 µl of diluted first-strand cDNA, 3.6 µl of PCR-grade water, 0.2 µl of specific primers, and 0.2 µl of mQ Primer R (5 -AGTGCAGGGTCCGAGGTA TT-3 ) . The mixtures were run under the following conditions: 95 • C for 5 min; 40 cycles of 95 • C for 10 s and 60 • C for 30 s; and 95 • C for 15 s, 60 • C for 1 min, and 95 • C for 15 s on an ABI 7500 Real-Time PCR System (Life Technology, United States). We used snRNA as an internal control . The DEM expression levels were calculated using the 2 − Ct method (Pang et al., 2019b).

STATISTICAL ANALYSIS
The data are expressed as the mean ± SD, and a T-test analysis of variance was used to compare two groups. A P-value < 0.05 indicated a statistically significant difference.

Overview of Small RNA Sequences in E. sinensis
To explore the role of miRNAs in the regulation of agonistic behavior in E. sinensis, we obtained 34,367,222 clean reads from the control group and 28,150,606 clean reads from the fight group using Illumina HiSeq TM 2000 (Supplementary Table S1). The small RNA sequences with lengths of 17-35 nt were the most enriched in 21-23 nt (Figure 1). As shown in Supplementary Table S2, the clean reads in the six small RNA libraries were composed of rRNAs, tRNAs, snRNAs, snoRNAs, and miRNAs.

Identification of Mature miRNAs
The miRNA clean reads were blasted against the mature sequences of known miRNAs included in the miRbase (Kozomara and Griffiths-Jones, 2014), which revealed a total of  893 mature miRNAs in E. sinensis, and the number of mature miRNAs in each library is shown in Supplementary Table S2.

DEM Analysis and RT-qPCR Verification
The levels of the DEMs between the control and fight groups were compared according to RPM values ≥ 5 as determined by a linkage hierarchical cluster analysis. We obtained 29 DEMs, namely, 8 upregulated and 21 downregulated DEMs, compared to the control group, indicating that these miRNAs might be involved in agonistic behavior-specific regulation (Figure 3).
Five DEMs were selected to verify the RNA-Seq results. The results of the RT-qPCR analysis confirmed that esi-miR-194-5p and esi-miR-192 were significantly downregulated in the fight group (Figure 2). Additionally, esi-let-7b-3p and esi-miR-750 were significantly upregulated in the fight group. We also found that esi-miR-146b-5p does not exhibit a significant difference in expression between the control and fight groups. Our results confirmed the data from the RNA-Seq.

De novo Assembly and Unigene Functional Annotation
From RNA-Seq, we obtained 145,309,464 raw reads from the fight group and 142,210,006 raw reads from the control group. The clean reads were selected by excluding reads that did not conform to the requirements. Thus, there were 136,737,066 and 139,889,948 clean reads in the fight and control groups, respectively (Supplementary Table S1), and these reads were used for de novo assembly. Then, 217,281 unigenes were obtained from six mRNA libraries. In total, 18,950 unigenes had a length greater than 1,000 bp, and 53,333 unigenes were greater than 500 bp in length (Figure 4).
Besides, the unigenes were used for further functional analysis. In total, 30,458 unigenes were annotated by the GO database and assigned to the following three groups: molecular function, biological process, and cellular component; these groups were further divided into 57 categories. Most unigenes were concentrated in the binding, cellular process, and cell categories. The signaling pathways of 8,447 unigenes were analyzed using the KEGG KAAS online pathway analysis tool 1 . In addition, all unigenes were mapped to the following five processes: organismal systems (25.64%), metabolism (25.13%), genetic information (16.16%), cellular processes (15.90%), and environmental information processing (15.17%). Most unigenes were associated with signal transduction, translation, transport, and catabolism; carbohydrate metabolism; and endocrine system pathways.

DEG Analysis
We performed an analysis to detect significant differences in expression levels across the samples and found 101 upregulated significant DEGs and 15 significant downregulated DEGs between the fight and control groups. Furthermore, 1 http://www.kegg.jp/blastkoala/ significant DEGs were assigned to KEGG pathways. The significantly enriched pathways (q-value < 0.05) included the renin-angiotensin system (RAS) (Figure 5), NOD-like receptor signaling pathway, cytosolic DNA-sensing pathway, renin secretion, and ribosome biogenesis in eukaryotes.

MiRNA Target Gene Prediction and Functional Annotation
Miranda was used to predict the target genes based on the mRNA transcriptome. In our work, in total, 12,819 unigenes were targeted by 647 miRNAs. In addition, 43 miRNAs were predicted to target more than 100 genes; among these, miR-6914-5p was predicted to target 2,258 genes, which is the largest number observed in this study. Also, we found that 3,688 unigenes were targeted by two or more miRNAs. Furthermore, 15 DEGs targeted by 23 miRNAs were predicted ( Table 2).
In addition, the predicted target unigenes of 29 DEMs were shown in the Supplementary Material (S1).
To analyze the functions of the target genes among the different groups, most genes were classified into 44 different GO terms. The terms male meiosis, apical cortex, and ARF guanyl-nucleotide exchange factor activity were the most highly enriched in the biological process, cellular component, and molecular function categories, respectively, between the control and fight groups. In addition, 259 target genes were assigned to 169 KEGG pathways. The top 20 enriched pathways are shown in Figure 5 and include the cAMP signaling pathway, MAPK signaling pathway-fly, protein digestion and absorption, and fatty acid biosynthesis.

DISCUSSION
In this study investigating agonistic behavior in E. sinensis, we reveal the underlying molecular mechanisms through an analysis of miRNA and mRNA expression profiles using a deep sequencing approach. We found 29 miRNAs and 116 unigenes with significant differential expression in the fight group. In the 29 DEMs, we did not find any that was characterized in crustaceans, but the mse-miR-281 (Manduca sexta) (Zhang et al., 2012), lgi-miR-750 (Lottia gigantea) , sme-miR-315-5p (Schmidtea mediterranea) (da Silva et al., 2015), cin-miR-4171-5p (Ciona intestinalis), and nve-miR-2025-3p (Nematostella vectensis) were characterized in invertebrates. And then, we analyzed the functions of the DEMs and found that they play some key roles in the regulation of animal behaviors. For example, esi-miR-199a-5p as a DEM could induce neuronal apoptosis by targeting the HIF-1α gene and then affect learning and memory behaviors in rats (Yan et al., 2017). Besides, miR-199a-5p was believed to participate in the regulation of glucose metabolism (Wu et al., 2015;Yan et al., 2017;Esteves et al., 2018). Agonistic behavior usually was accompanied by energy consumption, especially glucose metabolism (Briffa and Elwood, 2002). Due to the high conservation of miRNAs (Lewis et al., 2005), we assumed that miR-199a-5p may FIGURE 4 | KEGG pathway of the renin-angiotensin system. In addition, the red color indicates upregulated genes, and the yellow color indicates genes with no significant differences in expression. ACE is angiotensin converting enzyme. regulate crabs' agonistic behavior by mediating neurogenesis and energy metabolism.
Based on the analysis of the mRNA transcripts, we also found some DEGs which could regulate energy metabolism. For example, we observed that the KEGG pathway of the RAS was significantly enriched with DEGs ( Figure 5). The RAS, which consists of angiotensinogen, angiotensin-converting enzyme (ACE) isoforms, and angiotensin peptides and their receptors, performs many important functions, including glucose homeostasis (Luther and Brown, 2011), electrolyte balance (Jackson et al., 2018), and features of a neurotransmitter (Phillips et al., 1979) in mammals. In the RAS, angiotensinogen is transformed to form angiotensin (Ang) I under renin action and is then cleaved by ACE to form Ang II. Ang II can induce metabolic effects by acting on cell surface G protein-coupled receptors (White et al., 2019). In fact, the RAS has been shown to mediate glucose metabolism (Gamba et al., 2019). As previously discussed, agonistic behavior is usually accompanied by glucose metabolism in crustaceans. Therefore, there is new evidence suggesting that the RAS may participate in the regulation of agonistic behavior in crabs. Based on the above studies, ACE, which was significantly upregulated in our results, plays a key role in the RAS pathway. A series of studies show that ACE is related to glucose and lipid metabolism (Vaughan et al., 2016;Glover et al., 2019). In addition, ACE modulates behavioral and metabolic responses to diet in Drosophila melanogaster (Glover et al., 2019). Moreover, ACE may alter anxiety-like behavior (Wang L.A. et al., 2018) and cognitive ability (Gamba et al., 2019) by mediating the RAS. Therefore, we propose that the ACE gene plays a key role in the RAS pathway to alter energy metabolism, thereby affecting agonistic behavior in E. sinensis. Focusing on the neuroregulation in behavior, research show that let-7d played an important role in regulating the attention deficit hyperactivity disorder; Wu et al. (2010) found that an increased level of let-7d could decrease tyrosine hydroxylase activity by targeting galectin-3 and affect DA metabolism in rats (Wu et al., 2010). A recent study also showed that through a D3R target-mediated mechanism, let-7d improved anxiety and depression disorders (Bahi and Dreyer, 2018). The above evidence suggested that let-7d is involved in the upstream and downstream regulation of DA. According to previous studies, DA is a key neurotransmitter that may affect agonistic behavior in crustaceans. For example, Briffa and Elwood (2002) confirmed that the circulating levels of DA significantly decrease after a fight in Pagurus bernhardus, and the same result was found in shore crab (Carcinus maenas) (Sneddon et al., 2000). In addition, an injection of DA could improve catechol-O-methyltransferase (COMT) gene expression, and COMT might be related to aggressiveness in the swimming crab (Portunus trituberculatus) . In our previous two studies, we also found that the DA level is significantly decreased during agonistic behavior and suggested that an injection of DA reduced agonistic behavior in E. sinensis (Pang et al., 2019a;Yang et al., 2019). In contrast to let-7d, two DEMs, that is, esi-miR-200a and esi-miR-200b, exhibited a significant downregulation in the fight group. However, these DEMs also played important roles in regulating the dopaminergic system (Dulcis et al., 2017;Wu et al., 2018). Wu et al. (2018) proved that miR-200a could influence Parkinson's disease by mediating the DA2 receptor, which regulates the cAMP-PKA signaling pathway. In addition, miR-200b is associated with methamphetamine (METH) addiction, and METH could induce behavioral changes by changing the DA level in rats (Sim et al., 2017). Thus, the significant change in esi-let-7d, esi-miR-200a, and esi-miR-200b indicated that these factors may influence agonistic behavior in E. sinensis by mediating the DA regulation mechanism. However, the different expression levels among esi-miR-200a, esi-miR-200b, and esilet-7d suggest that these factors may play different roles in the regulation of agonistic behavior.
In our results, we predict that the E. sinensis retinoid X receptor (RXR) as a DEG was targeted by five miRNAs ( Table 2). In the eyes, retinoic acid (RA) is responsible for phototransduction through binding the G-protein receptors (Uray et al., 2016). The activity of RA can be regulated by RXR, which is important for the proper function of the retina (Janssen et al., 1999). RXR can regulate the development of cone photoreceptors in many animals, and the first step of the phototransduction of vision requires cone photoreceptors (Hennig et al., 2008). To the best of our knowledge, the visual systems are widely used for animal communication, including color discrimination (Chiou et al., 2011), individual recognition , and other animal behaviors . Undoubtedly, differences in visual signals can also induce changes in the behavioral responses of crustaceans. For example, mantis shrimp (Squilla empusa) displayed varying degrees of agonistic behavior according to vision differences (Wortham-Neal, 2002). In summary, we propose that visual signals were produced when the two crabs met, and then, the signals were conducted into the animals' CNS by receptors that may contain RXR. Subsequently, the command for agonistic behavior is issued by the CNS. Of course, the behavior may be regulated by complex physiological processes, such as neurotransmission, energy metabolism, and immunoregulation. Based on our results, genes other than those discussed above exhibited significant expression differences. Due to the length of this paper, we only analyzed the DEGs considered to perform considerable functions in the regulation of agonistic behavior. The other DEGs will be further explored in future studies.
Analyzing the predicted target unigene enrichment KEGG pathways (Figure 5), we found that the cAMP signaling pathway is the most enriched. The cAMP is a second messenger that can be activated by neuromodulators and then affect protein kinase A (PKA) activity (Buckley et al., 2016;Momohara et al., 2016). According to reports, the cAMP-PKA signaling pathway was detected to be essential for mediating loser and winner effects in agonistic behavior (Momohara et al., 2016). Our previous studies also showed that an injection of CPT-cAMP into crabs inhibited agonistic behavior in E. sinensis (Pang et al., 2019a). Therefore, the results of this study further demonstrate that the cAMP signaling pathway participates in the regulation of agonistic behavior. And the cAMP pathway played an important role in the energy metabolism and the conduction of nerve signals in crustaceans (Pang et al., 2019a;Yang et al., 2019). Also, we found that the MAPK pathway was significantly enriched with unigenes. The MAPK pathway is a three-tiered cascade (MAP3Ks-MAP2Ks-MAPKs), and the MAPK family includes extracellular signal−regulated kinase (ERK), P38, and c−Jun NH2−terminal kinase (JNK) (Malki et al., 2014). Several studies have reported that the MAPK pathway is linked to highly aggressive behavior (Malki et al., 2014;Malki et al., 2016a,b). In MAPK signaling cascades, ERK seems to be mainly associated with aggression (Malki et al., 2014). ERK is a key factor in the MAPK pathway that has been revealed to be involved in energy metabolism (Kuo et al., 2019). In fact, our results further show that the protein digestion and absorption pathway and fatty metabolism pathway are significantly enriched with DEGs. As previously mentioned, energy metabolism may be regulated by ERK or other key factors involved in agonistic behavior. To sum up, this study suggested that the regulation between miRNAs and genes could mediate the energy metabolism and the neurogenesis to regulate the agonistic behavior in E. sinensis.

CONCLUSION
In this paper, we performed an integrated analysis of miRNA and mRNA expression profiles and found that the miRNAs and mRNAs were involved in the behavior by mediating the progress of neurogenesis and energy metabolism in E. sinensis. We propose that DEMs and DEGs, mainly including esi-miR-199a-5p, esi-let-7d, the esi-miR-200 family, RXR, and ACE genes, may regulate agonistic behavior. Additionally, the miRNAs and mRNAs mainly focus on the RAS, cAMP, MAPK, and energy metabolism pathways after a fight. These results reveal the underlying molecular mechanisms of agonistic behavior and could guide researchers in investigating behavioral plasticity in the Chinese mitten crab.

DATA AVAILABILITY STATEMENT
The datasets generated for this study can be found in the https://dataview.ncbi.nlm.nih.gov/object/PRJNA565886?reviewer =ugafl4uq46l71e8ro649glsvhh.

ETHICS STATEMENT
The animal study was reviewed and approved by the Animal Experiments Ethics Committee of Shanghai Ocean University.

AUTHOR CONTRIBUTIONS
The experiments of this study were designed by YP and XY. YS, XS, and JL assisted YP and LH to complete several animal experiments. Results analysis and manuscript writing were made by YP, and then were reviewed and edited by XY. The funding resources came from YC. All authors read and approved the final manuscript.