t i l e Candidate Chemosensory Genes Identified in the Procambarus clarkii (Decapoda: Cambaridae) by Antennules Transcriptome Analysis

a vital the reproduction of However, information is Procambarus clarkii , important aquaculture we conducted the antennules transcriptome of Procambarus clarkii and identified 26 olfactory-related genes, including nine ionotropic receptors, two ionotropic glutamate receptors, ten variant ionotropic glutamate receptors, four cytochrome P450s, and one carboxylesterase. Our study suggested that ionotropic receptors were the main odorant receptors in Procambarus clarkii . Additionally, the key genes that are responsible for olfactory transduction in Procambarus clarkii were identified in the cAMP-mediated olfactory transduction pathway, including Cyclic nucleotide-gated cation channel beta-1, Adenylyl cyclase III, Calcium/calmodulin-dependent protein kinase II and Na + /Ca 2+ exchanger. As a whole, our study laid a solid foundation for further functional elucidations of olfactory molecular mechanism in Procambarus clarkii , and provided further insight for a better understanding of olfaction molecular mechanism in crustaceans.


INTRODUCTION
R ed swamp crayfish, Procambarus clarkii (Crustacea: Decapoda), has high production and economic values in China . In 2019, the total annual production of P. clarkii reached nearly 2,089,600 tons based on the Ministry of Agriculture data in China (www. yyj.moa.gov.cn). P. clarkii is an omnivorous species in freshwater, and adult crayfish in nature have a wide range of diets: Fish, shrimp, clam, algae, batata, and corn (Gherardi, 2006;Yue et al., 2010). Moreover, olfaction is the main factor affecting the ingestion of crayfish, and the antennules of P. clarkii have previously been reported as the main olfactory organ used to distinguish stimuli in the ingestion process (Kruangkum et al., 2019;Charles et al., 2016). However, the olfaction molecular mechanisms associated with the antennules of P. clarkii remain largely unknown.

O n l i n e F i r s t A r t i c l e
Olfaction plays a critical role throughout the life history of most animals (Du et al., 2018). Odor molecules can stimulate the olfactory organ to form an electric signal in the olfactory sensory neurons via various chemoreceptors and proteins, leading to diverse behaviors through a series of responses. The processes of olfactory detection are mediated by several functionally interrelated gene families: Odorant receptors (ORs) and ionotropic receptors (IRs) during chemosensory stimulus reception; odorant-degrading enzymes (ODEs) including antennalexpressed carboxylesterases (CXEs) and cytochrome P450s (CYPs) during degradation of redundant odorant molecules. ORs belonging to the G-protein-coupled receptors (GPCRs) family can form binding sites for odorant to transform binding chemical signals into neural signals (Benton et al., 2006). IRs, which are evolved from ionotropic glutamate receptors (IGluRs) in ancient protostomes, respond to both environmental and cellular signals and function similarly to ORs (Benton et al., 2009;Croset et al., 2010). A kind of special IRs, called IR coreceptors (e.g. IR25a, IR8a, and IR93a), is necessary and highly conserved in insects, and IR co-receptors can bind with IRs to take part in detecting odor (Abuin et al., 2011;Rytz et al., 2013). Furthermore, the extra odor molecules might affect the sensitivity of olfactory organs, ODEs can quickly degrade odor molecules and terminate short odor signals for continuous stimulation (Prestwich, 1987) to maintain the sensitivity of olfactory organs.
The major olfactory receptor of crustaceans was presented as the IRs identified in antennules (Daniel, 1997). Studies on the olfactory systems of crustaceans have been carried out on Decapoda such as Panulirus argus, Pagurus bernhardus, and Coenobita clypeatus (Charles et al., 2016). As confirmed by in situ hybridization, PargIR25a and PargIR93a were expressed in P. argus (Corey et al., 2013). In hermit crab, 18 IRs including IR25a and IR93a were identified in P. bernhardus, and IR25a, IR93a, and 27 divergent IRs including 22 variant ionotropic glutamate receptors (VIGluRs) were identified in C. clypeatus (Groh et al., 2013;Groh-Lunow et al., 2014). Despite the aforementioned studies, the researches on the olfactory receptors of Decapoda are still insufficient.
To get a better understanding of the molecular mechanisms underlying olfactory chemoreception in Decapoda, we selected P. clarkii as a representative economic shrimp and performed transcriptomics on antennules via Illumina sequencing technology. This study could help better understand the mechanisms of olfactory recognition in P. clarkii, and provide a solid foundation for further study on other Decapoda species.

Animals and sample preparation
The intact and healthy P. clarkii were obtained from Yancheng Sunshine Supermarket, Jiangsu, China. All P. clarkii were cultured in freshwater tanks (150 × 150 × 100 cm) with adequate aeration at 20±1℃. The experimental water quality parameters were monitored daily to maintain conditions of pH with 7.0 -7.5, salinity at 20 parts per thousand (ppt) and dissolved oxygen > 6 mg/L. In order to achieve the better effect of food induction, we chose a variety of plant foods and meat foods and observed shrimp feeding during the temporary rearing period. In a sevenday temporary rearing, P. clarkii were fed with algae, corn, duckweed, clam, shrimp, and fish at 19:00 daily. By comparing the amount of food left, it's found that P. clarkii preferred to eat corn (plant food) and clam (meat food). Therefore, corn and clam were selected as the inducer in the formal experiment. After 7 days of rearing in the aquarium, every 3 shrimps were randomly selected as a group, with a total of three groups of different inducers: water group (WG); herbivory group (HG); creophagy group (CG). The inducers (corn filtrate in HG; clam filtrate in CG and pure water in WG) were released 2 cm in front of the forehead of P. clarkii, and the antennules were immediately dissected from the individuals after 2 min. All samples were immediately frozen in liquid nitrogen, and stored at -80℃ until used.

Illumina sequencing and de novo assembly
Total RNA of antennules samples were extracted using Trizol reagent (Invitrogen, CA, USA). The purity and concentration of the RNA were detected by Nanodrop 2000, and the RNA integrity was detected by agarose gel electrophoresis. The libraries of tissues were constructed by TruseqTM RNA sample prep Kit, sequenced on an Illumina HiSeq 2000 platform (Illumina, San Diego, CA, USA). To obtain high-quality clean reads, trimmed adaptor sequences and low-quality reads were filtered from the RNA-seq raw reads. The transcriptome de novo assembly was carried out by Trinity with parameters set as default, only the longest transcript of each cluster was considered as a unique gene (unigene) for subsequent annotation.

Gene functional annotation
The unigenes were annotated by searching against the NCBI Non-Redundant protein database (NR), String, Swiss-Prot database, and EuKaryotic Orthologous Groups (KOG) databases, with a threshold of E-value ≤ 10 -5 . Based on the results of NR database annotation, the Blast 2GO program was employed to perform Gene Ontology (GO) annotation. The relationships among the unigenes O n l i n e

F i r s t A r t i c l e
and construct pathways were predicted by the Kyoto Encyclopedia of Genes and Genomes (KEGG) database.

Differential expression analysis and enrichment analysis
The expression level of each unigene was measured by using the fragments per kilobase of exon model per million mapped reads (FPKM) method using RSEM (Mortazavi et al., 2008). Significantly differentially expressed genes (DEGs) were selected with multiple check corrections based on the standard FDR < 0.05, |logFC| >= 1. To ensure the high quality of differentially expressed gene data, we eliminated the unigenes with obsessively low levels of expression (RPKM < 1). GO term enrichment and KEGG pathway functional enrichment analysis of DEGs was performed using Goatools (https://github.com/tanghaibao/ GOatools) and KOBAS (http://kobas.cbi.pku.edu.cn/ home.do), respectively. We defined the corrected p-value <= 0.05 to indicate significantly enriched (Minoru et al., 2008).

Construction of phylogenetic trees
To explore the relationships about olfactory receptors in Decapoda, olfactory receptors from different species (Litopenaeus vannamei, Eriocheir japonica sinensis, P. argus, C. clypeatus, P. bernhardus, and Palaemon serratus) were used for a comparison with our data. Amino acid sequence alignments were performed using MEGA5, with phylogenetic trees constructed using the neighborjoining (NJ) method with 10,000 bootstrap replications. Additionally, to ensure the accuracy of phylogenetic analysis, the phylogenetic trees of P. clarkii chemosensory genes were also constructed in IQ-TREE using the bestfitting substitution-model automatically with maximumlikelihood (ML) with 1,000 bootstrap replications. FigTree v 1.4.4, Adobe Illustrator CS4 and Photoshop C36 were used to visualize the phylogenetic trees.

Sequencing and assembly
The transcriptomic characteristics of P. clarkii were analyzed with high throughput Illumina sequencing technology. A total of 56,902,942, 51,010,266 and 49,809,984 raw reads were, respectively obtained from WG, HG, CG. After quality checks, 56,875,944, 50,939,572 and 49,744,810 clean reads were severally found for WG, HG, CG ( Table Ⅰ). Further assembly analysis showed that 755,807 transcripts contributed to 582,398 unigenes with a mean length of 479.89 bp ( Supplementary Fig. S1), the N50 length of 512 bp, and a GC content of 45.02% (Supplementary Table SⅠ).

Functional annotation
According to the results of annotation, 45,459 (25.61%) and 36,436 (20.52%) unigenes showed significant matches in the NR and Swiss-Prot protein databases, while 35,628 (20.07%) and 24,640 (13.88%) unigenes could be classified by GO and KEGG databases, respectively (Table II). GO annotation analysis suggested that the unigenes could be classified into three parts: "biological process", "cellular component" and "molecular function" (Supplementary Fig. S2). Most of the unigenes were in the terms "cellular process", "cell", and "binding". KOG is a tool to identify orthologous and paralogous proteins in eukaryotes and assign them to functional categories. In P. clarkii, the most enriched terms were in "General function prediction only" (Supplementary Fig. S3). In the KEGG database, the most of unigenes were assigned to the "Global and overview O n l i n e

F i r s t A r t i c l e
maps" pathways ( Supplementary Fig. S4).

Identification of candidate ionotropic receptors (IRs, IGluRs and VIGluRs)
Twenty-one candidate ionotropic receptors were identified including 9 PcIRs, 2 PcIGluRs, and 10 PcVIGluRs in P. clarkii antennules transcriptome. All PcIRs contained at least one characteristic structural domain (LCD-domain, LBD-domain) and 1-2 TMDs except PcIR93a with 4 TMDs ( Table Ⅲ). To understand the phylogenetic relationships of IRs in Decapoda, NJ tree and ML tree were constructed for 9 PcIRs and IRs of additional Decapoda species (Fig. 1). Although two different methods were adopted, the results tended to be consistent in this study. The phylogenetic trees (NJ tree and ML tree) indicated that PcIR93a was clustered with highly-conserved IR93a subfamilies. Then, IR coreceptors (IR25a and IR8a) of Decapoda formed other branch. Moreover, the majority of PcIRs and PargIR7 were clustered into one branch implying the closely orthologous relationships with P. argus.   Two PcIGluRs (PcIGluR1 and PcIGluR2) exhibited four structural domains (LCD-domain, LBD-domain, ATD-domain: PF01094, SBD-domain: PF00497) were identified ( Table Ⅲ). According to the structural analyses, complete ORFs and 2-3 TMDs were observed in PcIGluRs. Furthermore, we identified 10 PcVIGluRs denoted as PcVIGluR1-10. All PcVIGluRs contained 1-4 TMDs and LCD-domain. To explore the relationships between IRs, IGluRs, and VIGluRs, the phylogenetic trees for the analysis of ionotropic receptors in P. clarkii, L. vannamei, E. j. sinensis, and other Decapoda were constructed with ML and NJ methods (Fig. 2). Due to diverse algorithms, litter differences in locations in the NJ tree and ML tree were found, but the taxa are generally consistent. Results showed that co-receptors (IR8a and IR25a) were clustered into one branch, and IR93a clustered together alone. In addition, all co-receptors formed clusters together with IGluRs suggesting that these co-receptors had distinct orthologous relationships with IGluRs (Fig. 2). VIGluRs were divided into two branches which clustered with IRs in one large branch demonstrating a closer relationship between the VIGluRs and IRs (Fig. 2).

Identification of CXEs and CYPs
One putative PcCXE and four PcCYPs were selected from the P. clarkii antennules transcriptomes ( Table Ⅲ).

O n l i n e F i r s t A r t i c l e
Z. Wang et al.

Analysis of differential gene expression
To identify the key genes that are responsible for olfactory perception in P. clarkii, the unigenes of different libraries were subjected to a comparative analysis. A total of 4,224 genes were significantly differentially expressed (2,393 up-regulated and 1,831 down-regulated) in the HG vs. CG; 7,520 DEGs (3,683 up-regulated and 3,837 downregulated) in the HG vs. WG; and 6,114 DEGs (2,331 upregulated and 3,783 down-regulated) in the CG vs. WG (Fig. 3). Remarkably, The KEGG functional enrichment analysis performed in the DEGs revealed variations in the "olfactory transduction" pathway (ko04740) mechanisms for the HG, CG, and WG (Fig. 4). In the CG vs. WG, Adenylyl cyclase Ⅲ (AC3) was down-regulated in the olfactory transduction pathway, which might have adverse effects on the cAMP signaling pathway. In the HG vs. WG, three DEGs, Cyclic nucleotide-gated cation channel beta-1 (CNGB1), Calcium/calmodulin-dependent protein kinase II (CAMK2), and Na + /Ca 2+ exchanger (NCX), were upregulated in the olfactory pathway. These data indicated the key genes that were responsible for the response of P. clarkii in olfactory transduction. Fig. 3. Volcano plot showing the gene expression differences between the WG group, HG group, and CG group. Each dot in the graph represents a specific gene or transcript. Red dots represent genes that are significantly up-regulated, blue dots represent genes that are significantly down-regulated, and gray dots represent genes that are not significantly different points on the left and top represent more significant expression differences.

DISCUSSION
Olfactory perception is crucial for shrimps. However, the molecular basis of olfactory perception in shrimps is still relatively few. Therefore, we focused on the antennules of P. clarkii to analyze its internal olfactory perception mechanism. A suite of chemosensory genes from the P. clarkii antennules were identified by transcriptome analysis, including nine PcIRs, two PcIGluRs, ten PcVIGluRs, four PcCYPs, and one PcCXE.
IRs play a vital role in the olfactory transduction of odor molecules in crustaceans (Charles et al., 2016). Nine PcIRs were found in the antennules of P. clarkii, including IR93a. The co-receptor IR93a is an essential component of the olfactory detection process, assisting other IRs to detect odor (Rytz et al., 2013). In phylogenetic trees (NJ tree and ML tree), PcIR93a and other Decapoda IR93a were all clustered in the same branch, and it may play a key role in the P. clarkii olfactory perception. Then other IR co-receptors (IR25a and IR8a) of Decapoda species formed one branch. The results showed that co-receptor families were highly conserved in Decapoda. Additionally, PcIRs had relatively closely orthologous relationships with PargIR7, which suggested the olfactory mechanisms correlation between P. clarkii and P. argus. Unfortunately, the number of IRs in P. clarkii was limited. We speculated that the reasons can be attributed to two aspects: (i) the time of odor stimulation was too short to respond; (ii) P. clarkii was not hungry enough to enhance their desire for ingestion. Fig. 4. Annotation diagram of up and down regulated gene KEGG pathway in the HG and CG group. The genes with red/blue borders in the figure belong to the differential genes detected by this sequencing, where red represents up-regulated genes and blue represents down-regulated genes.
In previous studies, the IR evolutionary process was indistinct. It was reported that IRs are a new family of IGluRs related to excitatory neurotransmission . According to the analysis of the phylogenetic tree, we found most IR co-receptors had distinct orthologous relationships with IGluRs. This phenomenon proved that IR co-receptors were probably most closely resemble the ancestral IRs derived from IGluRs (Benton et al., 2009;Rytz et al., 2013). However, IRs are different from IGluRs O n l i n e

F i r s t A r t i c l e
Candidate Chemosensory Genes Identified in the Procambarus clarkii 7 in the olfactory perception of crustaceans. The majority of IRs were clustered in one branch and formed a large branch with VIGluRs. It indicated that IRs have distinct orthologous relationships with VIGluRs, and they might play a role together in the Decapoda olfactory mechanism. Rapid catabolism of redundant odorant molecules by ODEs including CXEs and CYPs might regulate odorant/pheromone concentration, participating in signal termination to maintain the sensitivity of olfactory organs (Vosshall, 2008). SlCXE7 works in pheromone signal termination and reduction of surplus odorant in Spodoptera littoralis (Durand et al., 2011). It was reported that treatment of the sensilla by a P450 inhibitor caused anosmia, suggesting that the P450-mediated metabolic clearance of odor is required for olfactory detection (Pottier et al., 2012). In our study, PcCXE and PcCYPs may protect P. clarkii olfactory receptor neurons (ORNs) against disturbance from unnecessary odors.
The transduction of odorant signal was reported to occur primarily through a cAMP-mediated pathway in most olfactory neurons (Juilfs et al., 1997), and Cyclic nucleotide-gated ion channels (CNGs) can be directly activated by cAMP (Gutièrrez-Mecinas et al., 2008). In this study, CNGB1 was present in the CNGs of the olfactory sensory neurons and was up-regulated in the HG, suggesting CNGB1 was involved in the cAMP-mediated olfactory transduction pathway. Furthermore, AC3 is an important component of the cAMP-signaling pathway (Bakalyar and Reed, 1990;Monneron and D'Alayer, 1980) and can transform ATP into cAMP causing cellular signal response (Taussig and Gilman, 1995). AC3 was significantly differentially expressed in the CG, implying that the importance of the cAMP-mediated olfactory transduction pathway in the olfactory perception of P. clarkii. Then Calcium ion balance was beneficial to olfactory signal transduction (Chen et al., 2006). In our study, NCX and CAMK2 were up-regulated in the HG. NCX plays an important quantitative role in the Ca 2+ transport of neurons as a relatively new family of transporters (Hassan and Lytton, 2019), and the up-regulation of NCX revealed that NCX promoted the dynamic adaptation of Ca 2+ signaling in sensory conduction neurons of the olfactory system to maintain calcium ion balance in P. clarkii. CAMK2 is a multifunctional serine/threonineprotein kinase that delivers calcium signals in a variety of cellular processes (Wong et al., 2019), and it was reported that calmodulin can inhibit G protein-coupled receptor kinases to indirectly activate the cAMP signaling pathway (Jones and Reed, 1989). These results suggested that the cAMP-mediated pathway might play an important role in the olfactory transduction of P. clarkii.

CONCLUSION
In summary, our study provided the first report on antennal transcriptome analysis in P. clarkii. A suite of chemosensory genes was identified including nine PcIRs, two PcIGluRs, ten PcVIGluRs, four PcCYPs, and one PcCXE. IRs originated from IGluRs are the main olfactory receptors of P. clarkii, and closely related to VIGluRs. IRs and VIGluRs might play a key role in the olfactory transduction of P. clarkii. Moreover, we found one type of olfactory transduction pathway in P. clarkii: cAMPsignaling pathway. In the olfactory transduction pathway, we report four genes (CNGB1, AC3, NCX, and CAMK2) related to olfactory detection in P. clarkii. This study will facilitate functional research on olfactory mechanisms in P. clarkii and other Decapoda. In the figure, the abscissa represents the secondary classification term of GO, the left ordinate represents the percentage of isogenic contained in the secondary classification, the right ordinate represents the number of isogenic compared with the secondary classification, and the three colors represent the three categories, in which green represents the biological process, blue represents a cellular component, and red represents the molecular function.
Supplementary Fig. S3. Histogram presentation of the results from the classification using EuKaryotic Orthologous Groups (KOG). In the figure, each color column represents A functional classification of KOG (represented by capital letters A~Z, the specific meaning is marked on the right), and the height of the column, namely the vertical coordinate, represents the transcript number with this kind of function.

O n l i n e F i r s t A r t i c l e
Supplementary Fig. S4. Distribution of Procambarus clarkii differentially expressed genes among Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways. The ordinate is the name of the KEGG metabolic pathway, and the abscissa is the number of genes annotated into the pathway and the proportion of the number of genes in the total number of genes annotated. According to the KEGG metabolic pathways involved, genes are divided into five branches: Metabolism (A), Genetic Information Processing (B), Environmental Information Processing (C), Cellular Processes (D), and organic Systems (E).
Candidate Chemosensory Genes Identified in the Procambarus clarkii