Genomic Characterization and Expression of Juvenile Hormone Esterase-Like Carboxylesterase Genes in Pacific White Shrimp, Litopenaeus vannamei

The sesquiterpenoid methyl farnesoate (MF), a juvenile hormone (JH) analog, plays important roles in many physiological processes of crustaceans, such as morphogenesis, molting and reproduction. Juvenile hormone esterase-like (JHE-like) carboxylesterase (CXE) is a key enzyme in MF degradation, playing a significant role in regulating MF titer. However, its function is barely known in shrimp. In this study, a total of 21 JHE-like CXEs (LvCXEs) were characterized in Pacific white shrimp Litopenaeus vannamei, based on the full genome and multi-transcriptomic data. LvCXE has a conserved triplet catalytic site (Ser-Glu-His) and a characteristic GxSxG motif. Most LvCXEs were highly expressed in the hepatopancreas, which was the main site for MF degradation. LvCXEs containing a GESAG motif showed a specific expansion in the L. vannamei genome. Those GESAG-containing LvCXEs presented differential expressions at different larvae stages and different molting stages of L. vannamei, which suggested their potential functions in development and molting. Additionally, when the transcription level of CXEs was inhibited, it could lead to failed molt and death of L. vannamei. When we further detected the expression levels of the key ecdysone responsive transcription factors including LvE75, LvBr-C, LvHr3 and LvFtz-f1 after the CXE inhibitor was injected into L. vannamei, they all showed apparent down-regulation. These results suggested that the expansion of LvCXEs in the L. vannamei genome should contribute to the regulation of metamorphosis at larvae stages and frequent molting during the growth of L. vannamei.


Introduction
Carboxylesterase is a superfamily of multifunctional enzyme which ubiquitously exists in animals, plants and microbes [1,2]. Based on sequence similarity and substrate specificity, carboxylesterases can be divided into eight subfamilies-α-esterases (Ae), β-esterases (Be), juvenile hormone esterases (JHE), gliotactins (Gli), acetylcholinesterases (AChE), neurotactins (Nrt), neuroligins (Nrl) and glutactins (Glt) [3]. Among them, JHE is considered to be the predominant enzyme involved in degradation of juvenile hormone (JH), which is a group of sesquiterpenoid compounds with pleiotropic functions in development, metamorphosis, molting, growth, reproduction and pheromone biosynthesis of insects [4,5]. JHE converts JH to JH acid and then regulates JH titer at appropriate levels, and thereby directs development and metamorphosis in insects [6,7]. The increase of JHE activity at late developmental stage is one of the important biochemical events that leads to pupation [8,9].

Characterization of CXE Genes in L.vannamei genome
Among 51 CXE genes, a total of 21 LvCXEs held a complete open reading frame (ORF), which were named as LvCXE1-LvCXE21, respectively. The predicted proteins of these LvCXEs varied from 488 to 669 amino acid residues. The estimated Mw ranged from 54.3 to 74.5 kDa, with pI of 4.53 to 8.98 (Table 1). The exon numbers of all LvCXEs varied from 8 to 13 (Figure 1). LvCXE2 and LvCXE3 were distributed in scaffold883, LvCXE4 and LvCXE5 located in scaffold962, LvCXE9 and LvCXE10 were distributed in scaffold1839. Other LvCXEs scattered in L. vannamei genome (Table 1).

Structure Analysis of LvCXEs
A total of 15 of 21 LvCXEs contained ten motifs (Figure 2), indicating that their protein sequences showed high similarity. Among ten motifs, motif 5, 4, 1, 3, 6 contained the RF, GG, DQ, GxSxG, E/D motif, respectively, which were mentioned in Table 2. Except for the incomplete LvCXE12, all LvCXEs were expected to contain a signal peptide sequence (17-26 aa) at the N-terminal ( Table 2). The multiple alignment revealed that all LvCXEs had a typical domain organization, including a conserved catalytic triad (Ser-Glu-His) and a carboxylesterase-specific GxSxG motif. A total of 14 LvCXEs contained the GESAG motif, and the other LvCXEs were different. In addition, LvCXEs had RF (arginine, phenylalanine) and DQ (aspartic acid, glutamine) residues at the same position as insect JHEs (Figure 3 and Table 2).

Structure Analysis of LvCXEs
A total of 15 of 21 LvCXEs contained ten motifs (Figure 2), indicating that their protein sequences showed high similarity. Among ten motifs, motif 5, 4, 1, 3, 6 contained the RF, GG, DQ, GxSxG, E/D motif, respectively, which were mentioned in Table 2. Except for the incomplete LvCXE12, all LvCXEs were expected to contain a signal peptide sequence (17-26 aa) at the N-terminal ( Table 2). The multiple alignment revealed that all LvCXEs had a typical domain organization, including a conserved catalytic triad (Ser-Glu-His) and a carboxylesterase-specific GxSxG motif. A total of 14 LvCXEs contained the GESAG motif, and the other LvCXEs were different. In addition, LvCXEs had RF (arginine, phenylalanine) and DQ (aspartic acid, glutamine) residues at the same position as insect JHEs ( Figure 3 and Table 2).

Phylogenetic Analysis
The phylogenetic analysis was conducted to investigate the evolutionary relationship between LvCXEs, insect JHEs and other crustacean CXEs. The results showed that all LvCXEs were clustered within the monophyletic clade of crustacean CXE ( Figure 4). Meanwhile, JHE genes of insect were

Phylogenetic Analysis
The phylogenetic analysis was conducted to investigate the evolutionary relationship between LvCXEs, insect JHEs and other crustacean CXEs. The results showed that all LvCXEs were clustered within the monophyletic clade of crustacean CXE ( Figure 4). Meanwhile, JHE genes of insect were clustered in another monophyletic clade. Additionally, most LvCXEs containing the GESAG motif were closely related to each other in the phylogenetic tree, suggesting that GESAG-containing CXE genes expanded specifically in L. vannamei. Furthermore, most of JHE genes of insect contained the GQSAG motif rather than GESAG.

Spatial and Temporal Expression of LvCXEs
RNA-seq data was used to analyze the spatial expression patterns of LvCXEs in 16 different tissues of the adult L. vannamei. The hierarchical clustering analysis showed that the expression patterns of LvCXEs can be divided into 4 modules ( Figure 5). Module 1, including LvCXE2, LvCXE9, LvCXE10, LvCXE15 and LvCXE21, was found to be expressed in almost all tissues detected, with exception of ovary, testis and muscle. Module 2, including LvCXE12, LvCXE17, LvCXE19, was also expressed in all tissues detected, but the expressions was lower than group 1. Module 3, including LvCXE3, LvCXE4, LvCXE5, LvCXE8, LvCXE14, LvCXE16 and LvCXE18, showed expression in specific tissues. For example, LvCXE4, LvCXE5 and LvCXE18 were almost specifically expressed in stomach, gill and antennary gland, respectively. Module 4, including LvCXE1, LvCXE6, LvCXE7, LvCXE11, LvCXE13 and LvCXE20, was highly expressed in hepatopancreas.
We also assessed the expression profiles of LvCXEs at early developmental stages of L. vannamei. Obviously, the expression levels of all LvCXEs could be clustered into two modules. Module 1 included LvCXE4, LvCXE6, LvCXE9, LvCXE10, LvCXE11, LvCXE13, LvCXE15, LvCXE16, LvCXE19 and LvCXE21, in which contained seven GESAG-containing genes. They showed high expressions from zoea I (Z1) to post-larva (P1) stage. In detail, their expressions increased until mysis II (M2) stage and decreased in mysis III (M3) stage, which was the end of metamorphosis. Other LvCXEs were clustered in module 2, whose expression levels were relatively low during the whole early developmental stages ( Figure 6). were closely related to each other in the phylogenetic tree, suggesting that GESAG-containing CXE genes expanded specifically in L. vannamei. Furthermore, most of JHE genes of insect contained the GQSAG motif rather than GESAG.

Spatial and Temporal Expression of LvCXEs
RNA-seq data was used to analyze the spatial expression patterns of LvCXEs in 16 different tissues of the adult L. vannamei. The hierarchical clustering analysis showed that the expression patterns of LvCXEs can be divided into 4 modules ( Figure 5). Module 1, including LvCXE2, LvCXE9, LvCXE10, LvCXE15 and LvCXE21, was found to be expressed in almost all tissues detected, with exception of ovary, testis and muscle. Module 2, including LvCXE12, LvCXE17, LvCXE19, was also expressed in all tissues detected, but the expressions was lower than group 1. Module 3, including LvCXE3, LvCXE4, LvCXE5, LvCXE8, LvCXE14, LvCXE16 and LvCXE18, showed expression in specific tissues. For example, LvCXE4, LvCXE5 and LvCXE18 were almost specifically expressed in stomach, gill and antennary gland, respectively. Module 4, including LvCXE1, LvCXE6, LvCXE7, LvCXE11, LvCXE13 and LvCXE20, was highly expressed in hepatopancreas. We also assessed the expression profiles of LvCXEs at early developmental stages of L. vannamei. Obviously, the expression levels of all LvCXEs could be clustered into two modules. Module 1 included LvCXE4, LvCXE6, LvCXE9, LvCXE10, LvCXE11, LvCXE13, LvCXE15, LvCXE16, LvCXE19 and LvCXE21, in which contained seven GESAG-containing genes. They showed high expressions from zoea I (Z1) to post-larva (P1) stage. In detail, their expressions increased until mysis II (M2) stage and decreased in mysis III (M3) stage, which was the end of metamorphosis. Other LvCXEs were clustered in module 2, whose expression levels were relatively low during the whole early

Expression Patterns of LvCXEs During Molting
All LvCXEs were detectable with dynamic expression patterns during the molting cycle. Module 1, including 16 LvCXEs, LvCXE1, LvCXE2, LvCXE4, LvCXE5, LvCXE6, LvCXE7, LvCXE9, LvCXE10, LvCXE11, LvCXE13, LvCXE16, LvCXE17, LvCXE18, LvCXE19, LvCXE20 and LvCXE21, showed similar expression patterns that the expression levels were high at C and D0 stage, decreased in D3 stage, and increased rapidly in D4 stage during the molting period (Figure 7, Supplementary Figure  S1). Furthermore, 11 of 14 GESAG-containing LvCXEs exhibited molt-dependent expressions, indicating that these expanded LvCXEs might be crucial for the molt cycle of L. vannamei. However, the results of hierarchical clustering analysis showed that the expression patterns of module 2, including LvCXE3, LvCXE8, LvCXE12, LvCXE14 and LvCXE15, were distinct from module 1. They exhibited low expressions in the C-D3 stage and rapid high expressions in D4 stage.

Expression Patterns of LvCXEs during Molting
All LvCXEs were detectable with dynamic expression patterns during the molting cycle. Module 1, including 16 LvCXEs, LvCXE1, LvCXE2, LvCXE4, LvCXE5, LvCXE6, LvCXE7, LvCXE9, LvCXE10, LvCXE11, LvCXE13, LvCXE16, LvCXE17, LvCXE18, LvCXE19, LvCXE20 and LvCXE21, showed similar expression patterns that the expression levels were high at C and D0 stage, decreased in D3 stage, and increased rapidly in D4 stage during the molting period (Figure 7, Supplementary Figure S1). Furthermore, 11 of 14 GESAG-containing LvCXEs exhibited molt-dependent expressions, indicating that these expanded LvCXEs might be crucial for the molt cycle of L. vannamei. However, the results of hierarchical clustering analysis showed that the expression patterns of module 2, including LvCXE3, LvCXE8, LvCXE12, LvCXE14 and LvCXE15, were distinct from module 1. They exhibited low expressions in the C-D3 stage and rapid high expressions in D4 stage.

OTFP Inhibition
To further verify the relationship between LvCXEs and the molting process of L. vannamei, we conducted OTFP inhibition experiments. The real-time PCR results revealed that the expression levels of 13 LvCXEs were significantly down-regulated in hepatopancreas after OTFP injection ( Figure 8A). The average inhibition efficiency was about 55%, and the highest was 88%. The expression levels of 6 LvCXEs showed no significant difference or were not detected. The molting times of each group were recorded and the cumulative molting rate was calculated every 24 h after injection. The results showed that a total of 14 shrimps of control group, accounting for 24.1%, molted in 12 days.
However, none of them molted successfully in the OTFP inhibition group. Moreover, eight shrimps died due to their exoskeletons being separated incompletely from the body, which resulted in unsuccessful molting process ( Figure 8B). Additionally, many key ecdysone response genes were also affected by OTFP inhibition. The expressions of E75 (ecdysone-induced protein 75), Br-C (Broad-Complex), Hr3 (hormone receptor 3) and Ftz-f1 (Fushi tarazu factor-1) gene were all significantly down-regulated after OTFP inhibition ( Figure 8C). Compared with the control group, the relative expressions of E75, Br-C, Hr3 and Ftz-f1 gene in OTFP inhibition group reduced by about 91%, 97%, 93% and 79%, respectively.

OTFP Inhibition
To further verify the relationship between LvCXEs and the molting process of L. vannamei, we conducted OTFP inhibition experiments. The real-time PCR results revealed that the expression levels of 13 LvCXEs were significantly down-regulated in hepatopancreas after OTFP injection ( Figure 8A). The average inhibition efficiency was about 55%, and the highest was 88%. The expression levels of 6 LvCXEs showed no significant difference or were not detected. The molting times of each group were recorded and the cumulative molting rate was calculated every 24 h after injection. The results showed that a total of 14 shrimps of control group, accounting for 24.1%, molted in 12 days. However, none of them molted successfully in the OTFP inhibition group. Moreover, eight shrimps died due to their exoskeletons being separated incompletely from the body, which resulted in unsuccessful molting process ( Figure 8B). Additionally, many key ecdysone response genes were also affected by OTFP inhibition. The expressions of E75 (ecdysone-induced protein 75), Br-C (Broad-Complex), Hr3 (hormone receptor 3) and Ftz-f1 (Fushi tarazu factor-1) gene were all significantly down-regulated after OTFP inhibition ( Figure 8C). Compared with the control group, the relative expressions of E75, Br-C, Hr3 and Ftz-f1 gene in OTFP inhibition group reduced by about 91%, 97%, 93% and 79%, respectively. and LvHr3 after OTFP inhibition. Significant differences in gene expression levels between two treatments were shown with a star (*) at p < 0.05 or two stars (**) at p < 0.01. 'ns' represents no significant difference.

Discussion
Previous studies have showed that a specific expansion of the CXE gene family was found in the LvHr3 after OTFP inhibition. Significant differences in gene expression levels between two treatments were shown with a star (*) at p < 0.05 or two stars (**) at p < 0.01. 'ns' represents no significant difference.

Discussion
Previous studies have showed that a specific expansion of the CXE gene family was found in the L. vannamei genome through comparative genomic analysis [36]. Since we found seven, two and zero CXE genes in the genome of Procambarus virginalis, Parhyale hawaiensis and Daphnia duplex, L. vannamei has the largest number of CXE genes. The expansion of the gene family is supposed to strengthen the specific biological process or contribute to the adaptive evolution [36,37]. Hence, comprehensive analysis of the CXE gene family would not only study their function systematically, but also illustrate the relationship between its expansion and specific biological characteristics.
The features of sequences determine their biological functions. All LvCXEs identified in this study had conserved motifs and sites similar to those of insect JHEs or other crustacean CXEs. Furthermore, all CXE genes of crustaceans were clustered into a clade of the phylogenetics tree. These results indicated that all LvCXEs belong to the CXE family. Low bootstrap values of some branches may be caused by low sequence similarity between CXEs and JHEs, since there are no large conserved domains. In general, JHE containing a GQSAG motif can degrade JH in insects [38,39]. In this study, only LvCXE5 contained the GQSAG motif, but it was not clustered with insect JHEs in the phylogenetic tree and not expressed in the hepatopancreas or gonads. It indicated that LvCXE5 was not an ortholog of insect JHE. Meanwhile, previous studies found that CXE genes containing the GESAG motif were proposed to have esterase activity for MF in crustaceans [30,32]. Consistently, 14 LvCXEs contained GESAG in our study. Compared with JHE genes containing GQSAG and esterase activity for JH, the difference may be caused by the divergence between crustaceans and insects. Hence, we speculated that LvCXEs containing a GESAG motif may degrade MF in L. vannamei. However, since the study on CXE and MF in crustaceans is limited, it is arbitrary to conclude that all LvCXEs are involved in the degradation of MF just according to their sequence characteristics.
To explore the preliminary function of LvCXEs, we first analyzed their expressions profiles in various tissues. MF esterase activity has been detected in the hepatopancreas and gonads of many crustaceans [40,41], which were considered to be the main sites for MF inactivation. In the present study, seven LvCXEs exhibited high expression in the hepatopancreas, but only LvCXE17 was expressed in the ovary. This is probably because the shrimp samples of different tissue transcriptomes are not in the stage of gonadal development. Interestingly, MF esterase activity has not been detected in the hemolymph of crustaceans so far [30], but we firstly found that three LvCXEs showed a high expression in hemocyte. Additionally, we firstly found that about half of LvCXEs were expressed in the antennal gland, which provided clues for the significant role of MF and the antennal gland in osmoregulation [12,42]. Meanwhile, our study also provided clues for better understanding the diverse functions of CXE genes and MF in crustaceans.
Similar to the functions of JH in insects, MF inhibits the metamorphosis of crustaceans [16,17]. Hence, the expression levels of LvCXEs were analyzed to further investigate their preliminary function during the developmental stage. Consistent with the increase of JHE activity in the late developmental stage of insects [8,9], 10 LvCXEs, including seven GESAG-motif containing genes, began expression in the nauplius, and were highly expressed in the zoeae and mysis stage gradually in L. vannamei. These may decrease MF titer and enable shrimp larvae metamorphosis to occur successfully during the early developmental stage. Furthermore, the decreased expression levels of LvCXEs at the post-larvae stage may restore MF titer. The above results reveal that LvCXEs possess similar function to the JHE gene of insects during early development. However, L. vannamei has evolved a pattern that multiple LvCXEs control and compensate functionally to maintain normal metamorphosis simultaneously.
The study of the function of LvCXEs during the molting cycle is of great significance since MF plays an important role in the molting process. In our study, most LvCXEs presented molt-related expressions, suggesting the potential relationship between LvCXEs and the molting process of L. vannamei. Under normal conditions, the expression levels of LvCXEs in module 1 increased significantly at D4 stage, which might result in a rapid decrease of MF titer and lead to the decrease of ecdysone in hemolymph. Similarly, the levels of MF, which fluctuate consistently with ecdysone [20][21][22], rise gradually at the early pre-molt stage but rapidly decline at the late pre-molt stage in the hemolymph of M. rosenbergii [43] and P. clarkia [15]. However, the injection of OTFP in D3 stage inhibited the transcription of LvCXEs and then might lead to the increase of MF titer rapidly, which elevated ecdysone titer in hemolymph abnormally. Furthermore, the synthesis of ecdysone in Y organ was then inhibited owing to the negative feedback, which was supported by the result that the expression levels of transcription factors LvE75, LvBr-C, LvHr3 and LvFtz-f1 were significantly down-regulated. Hence, the above results confirm that the high expression of LvCXEs in D4 stage is necessary for the molting cycle, and the CXE gene family plays an important role in the regulation of frequent molting in L. vannamei. However, MF exposure or MF injection accelerated the molting process of crustaceans according to previous studies [15,18]. We guess that they did not perform the MF treatment near molting, because they did not check the molting stage of the samples before the experiment. Hence, we speculate that MF can accelerate the molting process at the inter-molt or early pre-molt stage, but hinder the molting process near molting. Furthermore, it is also of concern that the expression patterns of module 2, including LvCXE3, LvCXE8, LvCXE12 LvCXE14 and LvCXE15, were different from module 1. Besides, the expression levels of LvCXE3, LvCXE8, LvCXE12 and LvCXE14 were very low not only in the hepatopancreas, but also at the early development stages. The above results showed that the function of these five LvCXEs should be different or divergent after the gene family expansion. However, there is no significant difference in the amino acid composition of their domains and key sites. It is speculated that the differences in their expression patterns are caused by different promoter regions. Of course, it needs more evidence to support this speculation.

Phylogenetic Analysis
The full-length amino acid sequences of CXE from crustaceans and JHE from insects were collected from GenBank databases. All collected sequences were aligned using MUSCLE program (v3.8.31) [52] with default parameters. Then the best-fit model for constructing the phylogenetic tree was selected by ProtTest 3.4.2 (parameter: -all-distributions -F -AIC -BIC -tc 0.5) [53]. A maximum likelihood (ML) phylogenetic tree of LvCXEs was generated by RAxML 8.0.26 [54] with a WAG+I+G+F model. Meanwhile, the ML phylogenetic tree of arthropod JHEs/CXEs was generated with a WAG+G+F model. All trees were tested with 100 bootstrap replicates and finally visualized using iTOL [55].

Expression Profiles
In previous studies, the RPKM (reads per kilobases per million reads) values of all transcripts of 20 early development stages, 8 molting stages and 16 adult tissues were calculated [44][45][46]. The RPKM values of all LvCXEs were extracted and then normalized with log2 conversion. Heat maps were then created and clustered with hierarchical clustering using TBtools software [51].

OTFP Inhibition
To investigate the relationships between LvCXEs and the molting process of L. vannamei, a total of 120 shrimps (about 5 cm for body length) in D0-D4 stage were selected and maintained in filtered seawater at a temperature of 25 ± 1 • C with continuous aeration. Before the experiment, all shrimps were kept in the aquarium for 3 days to acclimate them to the laboratory conditions. They were then randomly selected into a control group and experimental group. For each group, three replicates were conducted. The control group was injected with 10 µL mixture of ethanol and PBS (volume ratio 3:7) and the experimental group was injected with 10 µL 4.19 × 10 −3 mol/L OTFP solution, which was diluted by the mixture of ethanol and PBS (volume ratio 3:7). The molting times of each group were recorded and the cumulative molting rate was calculated every 24 h after injection. To further quantify the expression levels of some molt-related genes, 20 shrimps in D3-D4 stage were randomly selected into a control group and experimental group. The control group was injected with 10 µL mixture of ethanol and PBS (volume ratio 3:7) and the experimental group was injected with 10 µL 4.19 × 10 −3 mol/L OTFP solution, which was diluted by the mixture of ethanol and PBS (volume ratio 3:7). At 48 h postinjection, the hepatopancreas of four shrimps from each group was collected separately for RNA extraction.

RNA-Isolation and qRT-PCR
Total RNA was extracted using RNAiso Plus reagent (TaKaRa, Kyoto, Japan) according to the manufacturer's instructions. Agarose electrophoresis and NanoDrop 2000 (Thermo Fisher Scientific, Waltham, MA, USA) were used to detect the quality and the concentration of RNA, respectively. First-strand cDNA was synthesized using 1 µg of total RNA with a PrimeScript RT reagent kit (TaKaRa). According to the instruction, the mixed primer included Oligo d(T) primer and random hexamer primer was used and the genomic DNA had been removed by DNase treatment for cDNA synthesis. cDNA was synthesized according to the following procedure-37 • C for 1 h, and 85 • C for 5 s. The cDNA samples were stored at −80 • C for further use.
SYBR Green-based quantitative real-time PCR (qRT-PCR) was performed to detect the expression levels of LvCXEs and four ecdysone response genes in the hepatopancreas using THUNDERBIRD TM SYBR ® (TOYOBO, Osaka, Japan). Four technical replicates were conducted for each gene to eliminate the systematic errors. The 18S rRNA was used as an internal standard. Primer sequences were listed in Supplementary Table S1. qPCR was performed with an Eppendorf Mastercycler ep realplex (Eppendorf, Hamburg, Germany) using the following program-denaturation at 95 • C for 15 min; 40 cycles of 95 • C for 15 s, annealing temperature for 15 s, and 72 • C for 20 s. The PCR melting-curve was used to check the specificity of the PCR product. In addition, relative expression levels were calculated with the formula 2 −∆∆C T using the comparative C T method [56]. An unpaired two-tailed t test was used for statistical analysis by R (version 3.3.3). p value < 0.05 was considered statistically significant.

Conclusions
In summary, a total of 21 complete LvCXEs were identified in the L. vannamei genome. Their genome localizations, gene structures and conserved domains were comprehensively analyzed. The CXE genes with a GESAG motif expanded specifically in the L. vannamei genome, which increased the diversity of gene structures and functions. At the end of the early developmental stage, the high expression levels of LvCXEs may enable shrimp larvae to develop successfully. The high expression levels of LvCXEs in D4 stage are essential for the molting cycle. In contrast, the abnormal function of LvCXEs will lead to the failure of the molting process. In conclusion, the expansion of the CXE gene family is related to precise regulation for early development and frequent molting in L. vannamei.