Understanding molt control switches: Transcriptomic and expression analysis of the genes involved in ecdysteroidogenesis and cholesterol uptake pathways in the Y-organ of the blue crab, Callinectes sapidus

The crustacean molting process is regulated by an interplay of hormones produced by the eyestalk ganglia and Y-organs (YO). Molt-inhibiting hormone and crustacean hyperglycemic hormone released by the sinus gland of the eyestalk ganglia (EG) inhibit the synthesis and secretion of ecdysteroid by the YO, hence regulating hemolymph levels during the molt cycle. The purpose of this study is to investigate the ecdysteroidogenesis pathway, specifically genes linked to changes in ecdysteroid levels occurring at early premolt (ePM). To this end, a reference transcriptome based on YO, EG, and hepatopancreas was de novo assembled. Two genes (cholesterol 7-desaturase Neverland and cytochrome p450 307a1-like Spook) involved in ecdysteroidogenesis were identified from the YO transcriptome using sequence comparisons and transcript abundance. Two other candidates, Hormone receptor 4 and probable cytochrome p450 49a1 potentially involved in ecdysteroidogenesis were also identified. Since cholesterol is the ecdysteroid precursor, a putative cholesterol carrier (Apolipoprotein D-like) was also examined to understand if cholesterol uptake coincided with the increase in the ecdysteroid levels at the ePM stage. The expression level changes of the five candidate genes in the YO were compared between intermolt (IM) and induced ePM (iePM) stages using transcriptomic analysis. Expression analysis using qPCR were carried out at IM, iePM, and normal ePM. The increase in Spook and Neverland expression in the YO at the ePM was accompanied by a concomitant rise in ecdysteroid levels. The data obtained from iePM stage were congruent with those obtained from the normal ePM stage of intact control animals. The present findings support the role of Halloween genes in the ecdysteroidogenesis and molt cycle in the blue crab, Callinectes sapidus.

Introduction Crustacean molting, required for somatic growth, is regulated by hormones produced by the eyestalk ganglia (EG) and Y-organs (YO). The ganglia, located within the eyestalk of decapod crustaceans, contains the primary endocrine tissue, the X-organ medulla terminalis-sinus gland complex. This tissue produces and releases the four crustacean hyperglycemic hormone family members: molt-inhibiting hormone (MIH), crustacean hyperglycemic hormone (CHH), mandibular-organ inhibiting hormone, and gonad/vitellogenesis-inhibiting hormone [1][2][3]. Together with a multifunctional CHH, MIH suppresses YO activity, keeping ecdysteroid (Ecd) in the hemolymph low at the intermolt (IM) stage [3][4][5][6][7]. The injection of MIH and CHH into eyestalk-ablated animals can substitute for eyestalk tissue itself, suppressing hemolymph Ecd [4,5,7]. Hemolymph Ecd levels in turn influence MIH transcription in EG of the blue crab, Callinectes sapidus [8,9], possibly through tissue-specific isoforms of ecdysone/retinoid X receptors [10]. An interaction between the two tissues for molt control warrants further study on what changes occur in the YO from IM to premolt (PM), specifically leading to increased hemolymph Ecd.
With cholesterol as the precursor molecule, the YO produces two inactive forms of Ecd: αecdysone and 25-deoxyecdysone (25dE), which are converted in the peripheral tissues into the active forms of 20-hydroxyecdysone (20E) and 25-deoxy-20 hydroxyecdysone or Ponasterone A (PonA). The hydroxylation reactions are catalyzed by cytochrome P450 mono-oxygenases encoded by shade. The conversion process is partial, indicated by the measurable levels of PonA, 3-dehydro-20-hydroxyecdysone (3d20E), 20-hydroxyecdysone (20E) and 25-deoxyecdysone (25dE) in the hemolymph [4]. In C. sapidus, PonA is the principal Ecd in the hemolymph of PM animals, followed by 20E. As alluded to earlier, eyestalk neuropeptide MIH and/ or CHH in Ecd synthesis is a key regulator because eyestalk-ablated animals at the PM stage carry 20E as the major Ecd in the hemolymph, not PonA [4].
Hemolymph Ecd concentrations being cyclic repeatedly rise and fall during the molt cycle of several crustacean species including the blue crab, C. sapidus [4], the American lobster, Homarus americanus [11], and the tanner crab, Chionoecetes bairdi [12]. Although the hemolymph Ecd concentrations differ by species, the described pattern is that animals at the IM stage have low levels of Ecd ranging from 0 ng/mL-100 ng/mL, while those at PM have much higher levels with 200 ng/mL to 2500 ng/mL [11]. What exactly triggers the changes in Ecd concentrations remains to be unraveled in decapod crustaceans.
Since Ecd synthesis starts with cholesterol as the precursor, cyclic fluctuation of Ecd concentrations during the molt cycle implies changes in cholesterol levels [13][14][15]. Hence, cholesterol uptake and transport are considered the critical, rate-limiting step for Ecd synthesis in arthropods, as they do not synthesize cholesterol de novo [16]. In vertebrates, cholesterol uptake and transport are facilitated by the lipoproteins such as Low-Density-Lipoprotein or apolipoprotein (Apo), through a receptor-mediated endocytic pathway, Low-Density-Lipoprotein-receptor [17]. Insect lipoprotein is involved in cholesterol uptake by the insect prothoracic gland [18], suggesting that cholesterol-binding protein and its receptor expression may be positively associated with cholesterol transport and Ecd synthesis during the molt cycle. A hormone that is an equivalent of insect prothoracic gland hormone has not yet been found in crustaceans.
In this study, we aimed to determine the increase in the YO activity, especially changes in the expression level of genes involved in crustacean Ecd synthesis during the transition from IM to PM, specifically when the C. sapidus physiologically commit to the next molt cycle. To this end, we first compared the transcriptomic data of the YO of intact animals with those of eyestalk-ablated ones (induced early premolt iePM stage) and identified the genes that are potentially involved in early ecdysteroidogenesis and cholesterol transport/ uptake. We then measured the distribution of these genes across all crab tissues. And finally, we measured the gene expression in the YO of animals at different molt stages and in eyestalk-ablated animals, together correlating with the Ecd and cholesterol concentrations in the hemolymph.

Animal culture
The blue crab, C. sapidus, was produced in the Blue Crab Hatchery, Aquaculture Research Center (ARC, Institute of Marine and Environmental Technology, Baltimore, MD, USA). The animals were reared in individual compartments (15 cm × 15 cm) in a closed recirculating aquaculture system with 25 ppt artificial seawater at 22-23˚C [8]. They were fed daily with a piece of squid (10% to 20% of body weight) between 4-6 P.M. EST until they reached the target experimental size (>70 mm). All animals were sacrificed in the afternoon (2-6 P.M.). Water quality was examined daily by the ZooQuatic Laboratory LLC (Baltimore, MD).

Eyestalk ablation experiment
Adult C. sapidus females and males at IM (n = 5-7) were used for bilateral eyestalk ablation. In brief, hemolymph (50 μL) was directly withdrawn into a 28.5 gauge insulin syringe containing 50 μL of a modified crustacean anticoagulant [28] (pH 7.4) at a 1:1 ratio (volume fraction), for basal Ecd levels at t = 0 [29]. The animals were chilled by placing them on ice for 5 min, and both eyestalks were then removed. The bleeding through the cut was monitored, ensuring the survival of the animal. The ablated animals were returned to their compartment and maintained for 7 days, with hemolymph sampling on days 3 and 7. After sampling at t = 7 days post-ablation, the animals were sacrificed for tissue dissection. The YOs were dissected in icecold diethylpyrocarbonate-treated crustacean saline [30] under a stereomicroscope (Leica), frozen immediately on dry ice and stored at -80˚C until further analysis. EG and hepatopancreas were also collected in male C. sapidus following the same procedure. manufacturer's recommendations. RNA concentrations were measured using a Nanodrop. RNA samples from 3 to 5 individuals were pooled and submitted to Macrogen (www. macrogenusa.com) for RNA sequencing, where the quantity and quality of RNA for each sample was checked using a Bioanalyzer (Agilent).

De novo assembly and quantification
Four RNAseq read sets were generated: EG, hepatopancreas, and the YO at IM and iePM. The quality of the raw reads was evaluated with FastQC. The raw reads were pre-processed using Trimmomatic to remove reads containing adapter, reads containing poly-N, and low-quality reads [32]. Clean reads were assembled as a reference transcriptome using Trinity (v2.6.6) with default parameters [33]. The N50 parameter and the % mapped reads to the assembly using bowtie2 were assessed. A clustered transcriptome was built using CD-HIT-EST v4.6.6 with 0.95% identity [34]. The completeness of the clustered and unclustered assemblies were evaluated using Benchmarking Universal Single-Copy Orthologs (BUSCO) analysis [35]. BUSCO v3.0 was used with the dependencies augustus-3.2.1 and hmmer-3.2.1 to compare the transcriptomes against the arthoropoda_odb9 lineage dataset with an e-value threshold of 1e-05. Trimmed reads obtained from each tissue were mapped back to the reference transcriptome using Bowtie2 v2.1.0. The abundance of transcripts and genes were evaluated using RSEM v1.2.28 and calculated as Transcripts Per Million (TPM) values [36]. A transcript was considered expressed using the threshold of TPM � 1. Unique and shared transcripts were identified using an UpSet plot [37] generated by the web-based tool (https://gehlenborglab. shinyapps.io/upsetr/).

Functional annotation and candidate genes
Transcript abundance was used to explore the YO transcriptomes (IM and iePM). The most highly expressed transcripts of the YO at iePM (TPM � 100) were aligned against the subset taxonomy group "Arthropods" of the NCBI non-redundant (nr) protein database of using blastx-fast with an e-value cut off of 1.0e-5 using Blast2go suite in OmicsBox [38]. Contigs were functionally annotated following Blast2go workflow using Gene Ontology and completed by an eggnog annotation in OmicsBox [39]. The genes of interest were further curated manually using BLAST. The genes involved in the ecdysteroid synthesis and cholesterol uptake and transport that were found based on GO terms and sequence descriptions were selected for further analysis.

Sequence homology and phylogenetic tree analysis
The mRNA sequences were translated using the ExPasy (https://web.expasy.org). The predicted open reading frames were aligned against invertebrate sequences for homology and phylogeny. Phylogenetic tree was generated using Phylogeny.fr [40]. Sequences were aligned using MUSCLE, GBLOCKS site selection, and presented fasta format as supplemental data (S1-S5 Files). For each alignment, the sequences from Callinectes sapidus, were compares using blastp to find the best hit to: P. trituberculatus, S. paramamosain, H. americanus, P. monodon, P. vannamei, Hyalella azteca, Daphnia magna, Apis mellifera, and Drosophila melanogaster, followed by alignment and phylogenetic analysis as described above. The annotations were confirmed by comparisons with the canonical genes described in D. melanogaster. In one case, for the Apolipoprotein D-like gene a single, canonical sequence from D. melanogaster was not available, so metazoan outgroups were used.

Tissue distribution in an adult female C. sapidus: Primers and sequence verification
The expression and sequence of these genes were confirmed using endpoint PCR after candidates were identified in the YO transcriptome. A pair of primers for each gene (Table 1) were generated using NCBI Primer-BLAST based on the transcriptomic sequences. PCR templates were cDNA containing 12.5 ng total RNA equivalent from 17 different tissues in a single adult female C. sapidus at IM and ovarian stage 2. The PCR conditions were as follows: initial denaturation 94˚C for 2.5 min followed by 35 cycles at 94˚C for 30 s; 58˚C for 30 s and 72˚C for 1.5 min. The PCR products were analyzed on 1.5% agarose gel and stained with ethidium bromide. The PCR amplicons with the expected size were cloned and sequenced for confirmation as described [41,42].

RNA extraction, cDNA synthesis and quantitative RT-PCR assay
Total RNAs of YO from the prepuberty females at different molt stages (IM n = 8, ePM n = 6, PM n = 8) were extracted as described above and were quantified with a Nanodrop Lite Spectrophotometer (Thermo Scientific). Reverse transcription with 1 μg to 1.5 μg total RNA was carried out using a Takara-PrimeScript RT reagent Kit with gDNA Eraser (Clontech). After diluted to a final concentration of 12.5 ng/μL, the quality of each cDNA sample was examined by amplifying the AK gene [8,10]. The standards for the genes were generated as described [43][44][45]. Briefly, a plasmid DNA containing each gene PCR product (generated using primers in Table 1) was generated and diluted to range from 200 to 2 million copies/well. The average amplification efficiency of the standards and samples were -3.521 ± 0.152 (n = 12) and R 2 value of 0.995 ± 0.005 (n = 12, S1 Table). For expression levels of each gene (the primers listed in Table 2), the cDNA samples were assayed in duplicate using Fast SYBR Green Master Mix (Applied Biosystems) on an Applied Biosystems 7500 Real-Time PCR instrument. The melting curve analysis was included at the end of each assay, with expected melting temperature (Tm) ranged from 78 (ApoD), 80 (Spook and cyp49a1) to 81 (Neverland and HR4).

PLOS ONE
Understanding molt control switches: Transcriptomic and expression analyses of the Y-organ in the blue crab

Cholesterol assay
Free cholesterol and cholesteryl ester were measured using a fluorometric Amplex Red Cholesterol Assay Kit by following the manufacturer's instructions (Invitrogen). Briefly, 10 μL of hemolymph (a 1:1 ratio of hemolymph: anticoagulant [28]) from each animal at the following molt stages (IM, n = 8; ePM n = 15; PM, n = 11; IM n = 4; and iePM, n = 5) was diluted in 40 μL of reaction mix previously prepared, with the cholesterol standard ranging from 1 μg/ mL to 8 μg/mL. The Amplex Red reagent was added to the wells, and the plates were incubated in the dark for 30 min at 37˚C. Fluorescence was measured at 560 nm with an excitation at 590 nm immediately after incubation (SpectraMax M5). All samples, blanks, and standards were measured in duplicate.

Statistical analysis
Statistical analyses were performed using the R software package (v3.4.3). When two groups were compared, Student's t-test was calculated and significance was noted as follows:" � " P<0.05 and " �� " P<0.01. For multiple group comparisons, a non-parametric Kruskal-Wallis test followed by Dunn post-hoc test with Holm adjustment was performed. The significance at P<0.05 was noted as" � " or letters.

Sequencing and read assembly
A reference transcriptome was de novo assembled using Trinity, including the YO transcriptome at IM and iePM stages, eyestalk ganglia and hepatopancreas (total 161M reads). A total of 377,335,498 bases were assembled in 260,172 Trinity genes and 362,308 Trinity transcripts (Table 3). BUSCO analysis of the transcriptome against the arthropod database of 1,066 genes (Fig 1) revealed 93% complete with 1% missing. Among the complete genes, the majority (68%) were duplicated. The redundancy of the Trinity assembly was then filtered using CD-HIT-EST, resulting in 260,148 Trinity genes (99% retained) and 313,761 transcripts (87% retained) ( Table 3). The BUSCO analysis of the resulting Trinity_95 transcriptome improved the completeness to 96%, decreased duplicates from 67.8% to 26.5%, increased single-copy genes from 32.2% to 73.5% and reduced missing BUSCO to 0.5% (Fig 1). The Trinity_95 assembly was then used as the reference transcriptome for further analysis. Individual read sets of each tissue were mapped back to the Trinity_95 reference transcriptome: 33,758,864 and 52,953,709 reads for the YO at IM and iePM, respectively Table 2. List of primers that were used for the qPCR assays.

Functional annotation and genes involved in the molting process
The data from the YO transcriptomes were presented in the context of the molting process. The transcriptome at the iePM stage was investigated by focusing on the 608 highly expressed sequences with TPM>100 (max TPM of 269,633.82 for IM and 174,208.76 for iePM). Contigs involved in the ecdysteroidogenesis, together with cholesterol transport and uptake pathways, were searched using the following GO terms: "Process: ecdysteroid biosynthetic process," " ecdysone biosynthetic process", "Function: cholesterol binding," and "Process: regulation of cholesterol transport." Among the ten genes retrieved (Table 4), four belonged to the Halloween genes: Neverland (cholesterol 7-desaturase; accession No MW556747), Spook (cyp307a1-like; accession No MW556746), Disembodied (cyp302a1), and Shadow (cyp315a1) based on reciprocal best hits when using Drosophila and the blue crab transcriptome.
At the iePM, Spook was the most highly expressed Halloween gene as it was the 3rd most abundant transcript in the iePM transcriptome. Three of four Halloween genes (Spook, Neverland, and Disembodied) showed two-fold increase in the YO from IM to iePM (Table 4). Shadow presented two isoforms with different expression patterns. Three genes involved in cholesterol transport and uptake (Apolipoprotein D-like, accession No MW556748; D-βhydroxybutyrate dehydrogenase, and NPC1) were identified among the most expressed genes in the transcriptome (Table 4). NPC1 presented the highest increase at iePM with a mean fold

PLOS ONE
Understanding molt control switches: Transcriptomic and expression analyses of the Y-organ in the blue crab change of 4.2 followed by Apolipoprotein D-like and D-β-hydroxybutyrate dehydrogenase, a fold change of 3.7 and 1.7, respectively. Genes involved in ecdysteroidogenesis or cholesterol transport and uptake (Hormone receptor 4, accession No MW556749 and Adrenodoxin, mitochondrial) also increased from IM to iePM (1.9 and 5.3, respectively). Additionally, another cytochrome P450: probable cytochrome P450 49a1 (accession No MW556750) was also identified as highly expressed in the transcriptome with a mean fold change of 2.1.

Phylogenetic tree and sequence homology of candidate genes
Five of the ten genes described above were chosen for further investigation in C. sapidus YO based on their function and expression in the YO transcriptomes. These genes were involved in the two first steps of ecdysteroidogenesis (Neverland and Spook), an additional cytochrome p450 (probable cytochrome p450 49a1), one gene involved in the cholesterol transport (Apolipoprotein D-like), and Hormone receptor 4.
Neverland's Rieske iron-sulfur domains (cd03531, pfam00355) were identified in the C. sapidus' sequence (Table 5). Both Spook and Probable cytochrome p450 49a1 displayed conserved domains form the p450 superfamily, including the shared conserved domains pfam00067 and COG2124. Apolipoprotein D-like function of transport of small molecules was Table 4

. Abundance of genes in transcript per million (TPM) and their corresponding fragments per kilo base per million mapped reads (FPKM) involved in ecdysteroid synthesis, cholesterol uptake and transport in the YO transcriptomes of C. sapidus at intermolt (IM) and induced early premolt (iePM).
When a gene presented several transcripts, the abundance of both transcripts is shown. confirmed by the conserved domain lipocalin-like domain (pfam08212). The Hormone receptor 4 was characterized by the ligand binding domain of hormone receptors (smart00430 and pfam00104), including the ligand binding domain of orphan nuclear receptor Ecdysoneinduced receptor DHR4-like (cd06953) and a zinc finger conserved domain (cd00202, pfam00320, COG5641). The putative amino acid sequence of Neverland, Spook, Hormone receptor 4, and probable cytochrome p450 49a1 showed the highest similarity of the gazami crab, P. trituberculatus sequences (Fig 3A-3E). Such high identity was reflected in the phylogenetic trees where C. sapidus grouped with the other crab sequences in all five trees. The clade of H. americanus with the two crab species was found in four trees, albeit without strong support, followed by the two shrimp species. However, in the Hormone Receptor 4 phylogeny, the H. americanus sequence was placed outside the crab shrimp clade, albeit without bootstrap support. Callinectes sapidus, the other decapods, and an amphipod species, Hyalella Azteca, were grouped closely in well-supported clades (96 to 100% support).

PLOS ONE
Understanding molt control switches: Transcriptomic and expression analyses of the Y-organ in the blue crab The putative amino acid sequence of Neverland, Spook, Hormone receptor 4, and probable cytochrome p450 49a1 showed the highest similarity of the gazami crab, P. trituberculatus sequences (Fig 3A-3E). This high identity was reflected in the phylogenetic trees where C. sapidus grouped with the other crab sequences in all five trees. The clade of H. americanus with the two crab species was found in four of the trees, albeit without strong support, followed by the two shrimp species. However, in the Hormone Receptor 4 phylogeny, the H. americanus sequence was placed outside the crab shrimp clade, albeit without bootstrap support. Callinectes sapidus, the other decapods and an amphipod species, Hyalella azteca were grouped closely in well-supported clades (96 to 100% support). Interestingly, Daphnia magna sequences were separated from crustaceans and were located close to insects (Fig 3A-3E). Drosophila melanogaster sequences were located as an outgroup in four of five genes (Fig 3A, 3B, 3D and  3E). For the Apolipoprotein D-like sequence, D. melanogaster and Apis mellifera sequences were not conserved in alignments, so the outgroup of Mus musculus and Homo sapiens were used.
Interestingly, Daphnia magna sequences were separated from crustaceans and were located close to insects (Fig 3A-3E). Drosophila melanogaster sequences were located as an outgroup in four of five genes (Fig 3A, 3B, 3D and 3E). For the Apolipoprotein D-like sequence, D. melanogaster and Apis mellifera sequences were not conserved in alignments, so the outgroup of Mus musculus and Homo sapiens were used.

Candidate genes tissue distribution in the blue crab C. sapidus
The expression of five candidate genes were confirmed using endpoint PCR across 17 different tissues (Fig 4) with Spook, Neverland, Hormone receptor 4, probable cytochrome p450 49a1, and Apolipoprotein D-like showing expression in the YO. Interestingly, the expression of Spook and Neverland was also noted in the EG, thoracic ganglia complex, and brain. Hormone receptor 4, Spook, and probable cytochrome p450 49a1 were expressed in the ovaries. Apolipoprotein D-like had a wide distribution among the tissues tested.

Changes in circulating ecdysteroid and cholesterol in the hemolymph of intact and ablated animals at different molt stages
The Ecd levels in the hemolymph of intact normal animals were determined at three stages: IM, ePM, and mid-PM. The hemolymph had significantly low levels of Ecd at IM with 4.8 ng/

Candidate gene expression in the YO after eyestalk-ablation
An eyestalk-ablation experiment was carried out to examine if the elevated Ecd was related to gene expression using a qPCR assay with the primers listed in Table 2. Spook presented the highest transcripts/μg total RNA in the YO (Fig 6A).

PLOS ONE
Understanding molt control switches: Transcriptomic and expression analyses of the Y-organ in the blue crab

Candidate gene expression in intact animals at different molt stages
The validation of upregulated gene expression at iePM was carried out with animals undergoing a natural molt cycle at IM (n = 8), ePM (n = 6) and PM (n = 8). In this experiment, two genes, Spook and probable cytochrome p450 49a1 (cyp49a1), had the highest expression in the YO at PM (Fig 7A).

Summary of the ecdysteroidogenesis pathway
Based on the results from this study, we proposed the following summary of the ecdysteroidogenesis pathway in C. sapidus. First, the cholesterol concentration increased at ePM (turquoise arrow), resulting in the initial rise in ecdysone synthesis and release, via the upregulation of Neverland (Nvd) and Spook (Fig 8) at ePM. Then this elevated ecdysone concentration in hemolymph upregulated Apolipoprotein D-like (ApoD) at PM (dark blue arrow) to further accelerate the transport of cholesterol for ecdysteroid synthesis. The increased ecdysteroid concentration also up-regulated other ecdysone responsive genes, including probable cytochrome p450 49a1 (cyp49a1) and Hormone receptor 4 (HR4) at PM.

Discussion
The present study describes the identification and changes in the expression of candidate genes involved in ecdysteroidogenesis and cholesterol uptake in YO of the blue crab, C. sapidus. The candidate genes were identified using RNAseq analysis of the YO transcriptomes at two different molt stages: IM and iePM. Their expression in the YO is discussed in the context of the increase of hemolymph Ecd and cholesterol at iePM and ePM.
The transcriptomic analysis allows the functional genomics study of the critical physiological processes of non-model organisms, including decapod crustaceans. Specifically, transcriptomic analysis has been adopted to unravel the molecular mechanisms governing molting in the whiteleg shrimp, Litopenaeus vannamei [46], the Chinese mitten crab, Eriocheir sinensis [47], the land crab, Gecarcinus lateralis [48][49][50] and the Danube crayfish, Pontastacus leptodactylus [51]. A first YO transcriptome of intact animals has been sequenced in the blue crab C. sapidus [52] but further investigation is required to better understand the mechanism of the molting process, specifically by which the levels of Ecd are regulated during the molt cycle.
In the current study, the sequencing data for gene discovery were generated using the libraries obtained from pooled YO RNA samples from multiple individuals per molt stage. It is not unusual to use pooled RNA samples for sequencing. A similar approach has been used for gene discovery and identifying specific candidate genes in C. sapidus [53], the greentail prawn, Metapenaeus bennettae [54], or the giant river prawn, Macrobrachium rosenbergii [55]. A nonreplicated study can be used as a pilot study to test a proposed hypothesis and be validated by targeted gene-expression or further replicated study [56]. The high quality of YO transcriptome is supported with~96% of completed BUSCOs for the clustered assembly that validates the completeness and accuracy of our transcriptomic data. The reduction of duplicated BUSCOs and the increase of single-copy BUSCOs confirm the use of CD-HIT-EST to reduce redundancy. These results agree with other BUSCO analyses carried out in crustaceans such as the greentail prawn, M. bennettae [54], and the black tiger shrimp, Penaeus monodon [57]. Moreover, the Trinity assembly presents a high N50 (2,854), which implies the presence of isoforms. The clustered Trinity_95 showing a reduced N50 (1,510) after CD-HIT-EST confirms the redundancy improvement, which is similar to the N50 (1537.79) of the YO filtered transcriptome of G. lateralis [48].
As expected, when using a TPM threshold of 1 more contigs are found in the YO at iePM than IM, corroborating the elevated YO transcriptional activity. Overall, the number of contigs (40,605 at IM and 46,240 at iePM, respectively) obtained in this study is similar to the filtered YO transcriptome in G. lateralis (48,590) [48] and to the 31,661 predicted protein coding region in the intact C. sapidus [52]. As a first exploration, our work was focused on highly abundant transcripts (TPM) in the YO at iePM. All the sequences of interest belonged to the 33,609 transcripts shared by the two transcriptomes and presented a higher abundance (TPM) at iePM than IM.
The 608 most abundant sequences (TPM>100) were screened to identify candidate genes. Through functional annotation and based on insect research, the genes possibly regulating ecdysteroid levels in YO were identified, specifically those that are involved in the ecdysteroidogenesis (Neverland, spook), ecdysone responsive genes (cyp49a1, Hr4), and involved in the cholesterol transport and uptake (ApoD).
The involvement of Halloween genes (cytochrome p450 family) and a Rieske domain oxygenase in ecdysteroidogenesis have been well-described in insects [20,21,[58][59][60][61]. Orthologous genes identified in Daphnia pulex [25] underscore the functional importance and conservation of this pathway in the arthropod phylum. Homology for these five genes is determined based on phylogenetic analysis, which shows in most cases a clear decapod and crustacean group. These sequencing trees are generally consistent with taxonomic relationships in most cases. The first step has been described similarly in insect and crustaceans. First, the 7,8-dehydrogenase, Neverland, converts cholesterol into 7dC. The presence of Neverland has been identified in the YO of the black tiger shrimp, P. monodon [62], and the land crab, G. lateralis [49]. The next step, from 7dC to Δ4-Diketol, remains to be unraveled in crustaceans and is referred as the 'Black box', but the Halloween gene Spook is involved in this step in insects [23].
The molt cycle is regulated by an interplay between ecdysteroid and neuropeptide hormones (such as MIH and CHH [8]). YO is the tissue for ecdysteroidogenesis in decapod crustaceans and displays the exclusive MIH binding site [63]. However, Neverland and Spook expression in C. sapidus is not exclusive to the YO. Expression of these genes was observed in other tissues such as the EG, thoracic ganglia complex, brain, and ovary. A similar pattern has been observed in P. trituberculatus [26], but further study of the localization and translation of mRNAs in these tissues is needed. Neverland and Spook have been identified in ovarian cells in other species such as the salmon louse, Lepoeophtheirus salmonis [64], the water flea, D. magna [65], the acari, Ornithodoros moubata [66], the mosquito, Anopheles gambiae [67], and the fruit fly, Drosophila melanogaster [68]. In D. magna, two tissue-specific paralogs of Neverland are found in the gut epithelium (Neverland1) and in the ovary (Neverland2) [65]. The expression of Neverland and Spook in ovarian tissue in C. sapidus could be maternal RNAs for Ecd synthesis translated during embryogenesis, as an earlier report states that crustacean ovaries contain maternal Ecd [69]. Alternatively, the ovary may be another site for ecdysteroidogenesis in C. sapidus for ovarian development that takes place multiple times during adulthood. Adult females experience puberty-terminal molt after which the high levels of MIH in the hemolymph inhibits the YO activity to maintain low hemolymph Ecd at the permanent IM [70].
Expression of the candidate genes possibly responsible for the changes in Ecd levels was further confirmed using qPCR assays. A significant upregulation of Neverland expression in the YO at iePM after eyestalk-ablation is also seen in the intact animals at the same molt stage. Interestingly, Neverland presents the highest expression at ePM followed by a decrease at PM. A similar expression pattern is observed during the molt cycle in the YO of the black tiger shrimp, P. monodon [62]. Since the conversion of cholesterol to 7dC is required for further Ecd production at PM, it is proposed that the translated Neverland may have a long half-life. The presence of MIH keeps the Ecd low during IM by suppressing YO activity. Hence, it is plausible to suggest that the elevated MIH levels at ePM may negatively affect Neverland expression [8]. The increase of cholesterol concentration observed in the hemolymph at ePM in C. sapidus supports its recruitment for ecdysteroidogenesis. Like Neverland, upregulation of Spook expression is seen in the YO at iePM. This induction after eyestalk-ablation was not observed in the YO of P. trituberculatus [26]. Despite the lack of statistical significance, Spook presents a higher expression in C. sapidus at ePM and PM in a natural molt cycle. A similar increase in ePM has been demonstrated in YO of G. lateralis [49]. The increased Ecd in the hemolymph at ePM (or iePM) and PM, confirms the production of the molting hormone by the YO. Interestingly, Spook exhibits the highest expression in both qPCR assays and transcriptome analysis in C. sapidus. This result implies that Spook or the 'Black box' step (from 7dC to Δ4-Diketol) may be a key step in ecdysteroidogenesis.
The function of cytochrome P450 includes not only steroid synthesis and breakdown but also clearance of fatty acids and xenobiotics. While the function of cyp49a1 remains to be elucidated in C. sapidus, cyp49a1 shows high expression in both YO and ovary together along with a slight expression in the spermatheca and heart. Similar to the expression of Neverland and Spook, the up-regulation of cyp49a1 expression is observed at iePM after eyestalk-ablation and at PM during a natural molt cycle. In A. gambiae, cyp49a1 expression was higher in the pupae stage rather than the adult mosquito [71] suggesting a role of cyp49a1 in arthropod development at a life-stage dependent manner. It is important to note in the present study that the tissues associated with the adult stage, i.e., ovarian and spermathecal tissues, were developing in the animals undergoing pubertal terminal molt. This suggests that cyp49a1 could be involved in the development of C. sapidus.
Ecdysteroid regulates arthropod (insects and crustaceans) growth by signaling through a heterodimer receptor of EcR and RXR/Ultraspiracle (USP). The binding of ecdysone to the EcR-RXR/USP triggers a cascade of Ecd-responsive genes including Hormone receptor 4 [59]. The involvement of Hormone receptor 4 seems to be limited to metamorphosis in Drosophila [72,73] and ecdysis in the German cockroach, Blatella germanica [74]. In the Chinese mitten crab, Eriocheir sinensis, Hormone receptor 4 expression differs between the zoeal and megalopal stages, highlighting its role in crustacean metamorphosis [75]. In C. sapidus, Hormone receptor 4 was mostly expressed in the YO and ovary with a low detection in the mandibular organ and brain. Its expression was induced in the YO after eyestalk-ablation. The change in expression, though non-significant during the prepurberty molt cycle, shows a slight increase at PM. In order to have a better understanding of the role played by Hormone receptor 4 in C. sapidus, it would be interesting to investigate expression immediately after eyestalk ablation.
Unlike vertebrates, arthropods do not de novo synthesize the cholesterol needed to produce Ecd. Transport and cholesterol uptake may be limiting steps for ecdysteroidogenesis. Cholesterol is transported in the hemolymph into lipoproteins [18,76], which are in turn associated with apolipoproteins and involved in the lipoprotein membrane integrity and in binding to receptors [76]. We identified Apolipoprotein D-like expression in different tissues in the blue crab including the YO. An increase in Apolipoprotein D-like expression by qPCR after eyestalk-ablation was confirmed by an up-regulation at PM in a natural molt cycle, implicating Apolipoprotein D-like in cholesterol uptake by the YO in C. sapidus.
Increased cholesterol concentrations in the hemolymph of C. sapidus at iePM and at ePM in a normal molt cycle provide the YO with the precursor needed to increase the production of Ecd. While Ecd production by the YO reaches its maximum at PM [4], the cholesterol levels remain similar in ePM and PM. This suggests that there may be some differences in the rate of cholesterol transport in the hemolymph and uptake by the YO. The YO at PM may increase cholesterol uptake while maintaining the necessary amounts of hemolymph cholesterol at ePM levels. This notion is supported by studies that demonstrated radiolabeled cholesterol accumulated the most in the YO while a small increase was observed in the hemolymph suggesting an active uptake by the YO [14,15].

Conclusions
This study presents an exploration of the blue crab C. sapidus YO transcriptome during normal and induced molt cycles to identify candidate genes in the ecdysteroid synthesis, cholesterol uptake and transport, and the characterization of their expression by qPCR. Ecdysteroidogenesis genes (Neverland) and members of the Halloween genes (Spook, Disembodied, and Shadow) have been identified in the YO transcriptomes, and the expression of the early ecdysteroid synthesis steps (Neverland and Spook) was characterized during a natural molt cycle. The expression of two ecdysone-responsive genes (probable cytochrome p450 49a1 and Hormone receptor 4) increased accordingly with the Ecd production. A similar pattern was observed with Apolipoprotein D-like and the cholesterol concentration in the hemolymph, suggesting its role in the uptake of the latter by YO. Overall, key genes with increased transcripts at ePM, along with elevated cholesterol in the hemolymph, corroborate the elevated Ecd synthesis by the YO as summarized in Fig 8, which overall supports the earlier finding of the presence of positive feedback of Ecd on YO [8]. The transcriptome and expression data are described in this paper. The function of each gene at the translation and protein level remains to be studied in the future.
Supporting information S1 Table. Quantitative RT-PCR standard curve parameters.