This study was approved by the Animal Care and Use Committee of Jiangsu Ocean University (protocol no. 2020-37; approval date: September 1, 2019). All procedures involving animals were performed in accordance with guidelines for the Care and Use of Laboratory Animals in China.

Healthy M. japonicus individuals were obtained from the Jiangsu Key Laboratory of Marine Biotechnology of Jiangsu Ocean University (Lianyungang, China). All individuals were acclimated for two weeks in environmentally controlled breeding tanks. The average body length and body weight of M. japonicus were 6.30 ± 1.12 cm and 1.48 ± 0.75 g, respectively. The experimental shrimp were euthanized with the anesthetic alcohol: eugenol = 10:1. The carapaces of 20 individuals were sampled to detect the contents of different carotenoids, including astaxanthin, lutein, zeaxanthin, and beta-carotene. Three individuals from each of the two varieties were used to observe the composition of pigment cells in the carapace using an optical microscope. A total of 60 individuals, including 30 variety I individuals and 30 variety II individuals, were sampled for comparative transcriptome analysis. The carapace of these prawns was sampled and immediately preserved in liquid nitrogen for total RNA extraction.

The obtained carapace was dried at 60°C for 48 hours, crushed with a high-speed grinder, and stored at 20°C for later use. The astaxanthin extraction was slightly modified from the method of Hu et al. (Hu et al., 2018). Approximately 1 g of the sample was accurately weighed and placed in a 50 mL round-bottomed centrifuge tube with 4 mL of 0.1% BHT ethanol solution. The homogenizer (IKA-T 18 basic, Germany) was fully homogenized for 2 min. The mixed homogenate was ultrasonically extracted for 10 min below 15°C, and centrifuged for 5 min at 8000 r/min. The extract was dehydrated through anhydrous sodium sulfate filtration and filtered into a 100 mL brown rotary evaporation bottle. Depressurized concentration at 40°C ± 2°C was performed on the rotary evaporator, nearly dry. After drying with nitrogen, 1.0 mL of 0.1% BHT ethanol solution was accurately added and the extract was fully dissolved. It was filtered through a 0.22 μm organic membrane and detected by the high-performance liquid chromatography (HPLC) (LC 2030C 3D, Shimadzu, Japan).

The extraction methods of the other three carotenoids have been slightly modified. Approximately 4 g of the sample was accurately weighed and placed in a 50 mL centrifugal tube, and the 20 mL of 0.1% BHT ethanol solution was added. The mixed homogenate was extracted for 10 min below 15°C, added to 10 mL of 30% KOH solution, placed in a thermostatic oscillating water bath preheated to 50°C ± 2°C, and saponified for 20 min. The solution was removed and allowed to cool to room temperature. Then, 50 mL extraction solvent was added and extracted for 3 min in vortex oscillation against light. Ten milliliters extraction was dried with nitrogen. For example, 1.0 mL of 0.1% BHT ethanol solution was added accurately to fully dissolve the extract. The solution was filtered through a 0.22 μm organic membrane and examined by the HPLC.

Total RNA was extracted from the carapace of 30 individuals of types I and II. TRIzol reagent (Takara, Dalian, China) was used to extract RNA from the carapace tissue samples. The degradation and contamination of total RNA were assessed by 1% agarose gel electrophoresis, RNA purity and concentration were determined by spectrophotometry, and total RNA integrity was determined by a bioanalyzer. The total RNA from each group of 30 shrimp was randomly divided into 6 portions. Total RNA from each sample (containing 5 individuals) was mixed equally. Finally, each group had 6 RNA samples, of which 3 RNA samples were used for transcriptome analysis, and the remaining samples were used for fluorescence quantitative PCR analysis.

Approximately 3 micrograms of RNA per sample was used to synthesize first- and second-strand cDNA. Six cDNA libraries were generated by an Illumina kit and sequenced on the sequencing platform. High-quality data are obtained by deleting joint sequences and low-quality data. Reassembled high-quality data were obtained by the Trinity software (Grabherr et al., 2013) with the recommended parameters. The resulting reference sequence was used for subsequent analysis. HISAT2 software (https://daehwankimlab.github.io/hisa t2/) was used to align the high-quality reads with the reference genome of Marsupenaeus japonicus (unpublished), and the location information of reads on the reference genome was obtained. The assembled transcripts were annotated against the Nr, GO and Swiss-Prot databases.

The feature v1.5 tool (Yang et al., 2014) of the subread software was used to determine the number of reads corresponding to each gene. RSEM software (Simon and Huber, 2010) was used to normalize the expression levels of transcripts from different groups to transcripts per kilobase per million fragments. Gene function annotation was performed using the NCBI and KEGG databases. The DESeq2 method (Love et al., 2014) was used to screen differentially expressed genes between the two groups (variety I and variety II). When the p value < 0.05 and log₂(fold change) ≥ 1, the gene was considered to be significantly differentially expressed. GO enrichment and KEGG analysis were performed using ClusterProfile software (Wu et al., 2021) with a q-value < 0.05.

The fluorescence quantitative PCR was used to detect the expression levels of ten differentially expressed genes in the two groups. Quantitative primers were designed using Primer software based on transcript sequences ( Table S1 ). The elongation factor 1-α (EF1-α) was served as the house-keeping gene. The PCR efficiency was determined by calculating the slopes and regression curves of the standard curves. The SYBR kit (Takara, Dalian, China) was used for quantitative PCR. Three biological replicates and three technical replicates were performed for all samples.

Total RNA was extracted using the TRIzol reagent (Takara, Dalian, China). The purity and integrity of RNA were also determined. First strand cDNA was synthesized using the TransScript kit (Takara, Dalian, China) as instructed. Primers were designed based on our transcripts and the sequences of related species in NCBI. The PCR amplification product was connected to the vector and double-ended sequencing was performed. The gene’s open reading frame was predicted using the ORF finder tool. Expasy software (https://web.expasy.org/compute_pi/) was used to predict the molecular weight and isoelectric point of the protein sequence. The amino acid sequence of the gene was deduced by EMBOSS software. SOPMA (Geourjon and Deleage, 1995) and PROSITE software (Sigrist et al., 2010) were used to predict the secondary structure of protein sequences. ClustalW (Larking et al., 2007) was used for the multiple sequence alignment and MEGA software (Kumar et al., 2016) was used for phylogenetic analysis.

We observed the composition of pigment cells in the carapace of the two varieties using an optical microscope. The setae structure on the carapace was clearly observed. It is possible that the manipulation of sampling makes the setae disorganized. The results showed that there were red pigment cells in the stripe pattern and yellow pigment cells in the other parts ( Figure 2 ). Both red pigment cells and yellow pigment cells showed dendritic morphology. The red pigment cells are larger than the yellow pigment cells.

We determined the contents of different carotenoids, including astaxanthin, lutein, zeaxanthin, and beta-carotene in the carapace using the HPLC ( Figure 3 ). The results showed that the content of astaxanthin was the highest, which was significantly (P < 0.05) higher than that of other carotenoids. The contents of zeaxanthin and beta-carotene were comparable, and both were significantly (P < 0.05) higher than the content of lutein. The lutein content was significantly (P < 0.05) lower than that of the other carotenoids.

In the variety I group (VI-A, VI-B, VI-C), Illumina sequencing obtained 134.5 million raw data and 132.4 million clean data, with a total of 19.86 Gb of data ( Table S2 ). The variety II group (VII-A, VII-B, VII-C) obtained 134.7 million raw data points, 131.6 million clean data points, and a total of 19.74 Gb of data. The mapping rate of 6 cDNA libraries to the reference genome of M. japonicus ranged from 72.85% to 89.2%. A total of 23629 unigenes were obtained, with an average length of 1630.17 bp, and 2457 new genes were predicted. The Pearson correlation coefficient R² between groups calculated based on the FRKM value was greater than 0.86.

A total of 1984 differentially expressed genes were identified in this study, of which 1125 genes were upregulated and 859 genes were downregulated ( Figure 4 ). Among these differentially expressed genes, we observed some genes that may be associated with body color formation, such as CRCN, apolipoprotein D (ApoD), tubulin alpha-1 chain, cuticle protein, and low-density lipoprotein receptor ( Table S3 ).

To understand which biological processes and pathways are involved in the body color formation or material transport of M. japonicus, we performed GO and KEGG enrichment analyses of differentially expressed genes. In the biological process category, most DEGs were annotated with involvement in proteolysis, transmembrane transport, and oxidation-reduction processes ( Figure 5 ). Among the DEGs categorized as cellular components, intracellular membrane-bounded organelle, membrane-bounded organelle, and extracellular regions were the best represented. Additionally, most of the molecular function-related DEGs were associated with transmembrane transporter activity, peptidase activity, and oxidoreductase activity. In general, “anion transport”, “inorganic anion transport”, “anion transmembrane transporter activity”, and “inorganic anion transmembrane transporter activity” were the terms with the most differentially expressed unigenes.

The obtained differentially expressed genes were annotated in the KEGG database. The DEGs were mapped to 108 KEGG pathways. Figure 6 showed the 20 significantly enriched pathways of the upregulated and downregulated genes. More upregulated genes were enriched in the lysosome (dpx04142), metabolism of xenobiotics by cytochrome P450 (dpx00983), glycolysis/gluconeogenesis (dpx00010), and pentose and glucuronate interconversions (dpx00040). More downregulated genes were enriched in the endocytosis (dpx04144), purine metabolism (dpx00230), glycosaminoglycan biosynthesis (dpx00532), and the Wnt signaling pathway (dpx04310).

Ten differentially expressed genes of interest were screened from the carapace tissue of M. japonicus for real-time quantitative PCR analysis to verify the expression patterns of unigenes in transcriptome data. The fluorescence quantitative PCR results of candidate genes showed good correlation with the expression level from RNA-seq data, and the expression trend of the same gene was consistent between the two groups ( Figure 7 ).

According to the transcriptome analysis, a total of 15 unigenes encoding cuticular proteins or related proteins exhibited higher expression levels in variety II than in variety I ( Table S3 ). These unigenes were classified into two groups based on domains. The unigene Evm.Hic29.46 encoded a protein containing a carbohydrate-binding domain (CBM) 14 domain, and other unigenes encoded proteins containing a chitin_bind_4 domain.

GO analysis also showed that the binding process was the term with the most differentially expressed unigenes. In this study, we identified several genes with binding functions, such as tubulin- and fibronectin-related genes, which were significantly higher in the variety I than in the variety II ( Figure 8 ). As carotenoid-binding proteins, CRCN-A2 and CRCN-C1 subunit also had higher expression levels in the variety I individuals. The results of sequence alignment showed that crustacyanin A2 in the variety I and variety II shared 98.42% sequence identity, with three amino acid differences (18:S & G, 112:A & V, 113:L & F). The results of sequence alignment showed that crustacyanin C1 in the variety I and variety II shared 99.49% sequence identity, with one amino acid difference (28:A & S). The lipocalin consensus region marked in orange showed no amino acid residue differences ( Figure 7 ). There were two obvious differences in the secondary structures of the CRCN-A2 and CRCN-C1 ( Figure S1 ). Phylogenetic analysis showed that the CRCN-A2 and CRCN-C1 of M. japonicus were highly similar to those of its genetically close species ( Figure 9 ).

Unlike the crustacyanin gene, apolipoprotein D, a lipid transfer protein, had a higher expression level in the variety II individuals ( Figure 7 ). The results of sequence alignment showed that apolipoprotein D in the variety I and variety II shared 98.47% sequence identity, with three amino acid differences (17: F & S, 31: S & P, 92: V & A). The consensus region marked in orange showed no amino acid residue differences ( Figure 8 ). Similarly, there are two obvious differences in the secondary structure of these two amino acid sequences, which may be related to their functions ( Figure S1 ).

The type and number of pigment cells determine the color and vividness of the animal (Border et al., 2019; Vissio et al., 2021). The species of pigment cells in crustaceans include melanocytes, yellow pigment cells, red pigment cells and mixed pigment cells (Qin et al., 2021). N. denticulata sinensis has five basic colors, including red, yellow, blue, black and white, and the pigment forms include dot, branch, snowflake, line, ribbon, etc (Lu et al., 2022; Zhang et al., 2022). The red pigment cells were mainly punctate and dendritic, and the yellow pigment cells were mainly linear or banded (Lindsay-Mosher and Pearson, 2019). Most of the time, there are five common pigment cells on N. denticulata sinensis, some of which overlap. In this study, the setae structure on the carapace was clearly observed, and there were red pigment cells in the stripe pattern and yellow pigment cells in the other parts. The distribution of pigmented cells directly affects the pattern characteristics of the carapace, which has been studied in similar ways in other animals (Bullara and De Decker, 2015; Djurdjevič et al., 2019). Red pigment cells are the main pigments in the stripe pattern of M. japonicus. The pigment cells in crustaceans are regulated by the actions of pigment-dispersing and pigment-concentrating peptide hormones on epidermal chromatophores (Gaus et al., 1990; Wei et al., 2021).

The crustacean exoskeleton perfectly reflects the pattern and color of the underlying epithelial tissue (Ghidalia et al., 1985). Carotenoids, particularly astaxanthin, are the primary pigment in crustacean shell color (Wade et al., 2005). Compared with other carotenoids, the content of astaxanthin was the highest in the carapace of M. japonicus. In Panulirus cygnus, Wade et al. found a 2.4-fold difference in the amount of total carotenoids present in the shell extracts of reds compared to whites, and a correlated twofold difference in the amount of esterified astaxanthin present in the epithelium of red versus white individuals (Wade et al., 2005). The effect of esterified astaxanthin on body color was more significant, which has been studied in Oncorhynchus mykiss, Amphiprion ocellaris, and P. monodon (White et al., 2002; Tume et al., 2009; Ho et al., 2013).

In crustaceans, carotenoid-binding protein is important for stabilizing carotenoids, regulating color, and participating in development and antioxidant activities (Wade et al., 2017; Maoka, 2020; Tan et al., 2020). The shell color of crustaceans is primarily determined by incorporated of the carotenoid astaxanthin into a macromolecular protein complex known as crustacyanin (Chayen et al., 2003). Crustacyanin covalently binds to ingested astaxanthin to form a complex that results in a conformational change in astaxanthin that determines the species-specific color and pattern of the crustacean exoskeleton (Wade et al., 2012; Jin et al., 2021; Zhao et al., 2021; Chen et al., 2022). In this study, we found that the expression levels of CRCN A2 and CRCN C1 in variety I individuals were higher than those in variety I individuals. A lower expression level of crustacyanin might limit the binding and subsequent utilization of free astaxanthin (Jin et al., 2021). Whether low expression levels are directly or indirectly related to the formation of markings still needs further investigation.

Once loaded with carotenoids, apolipoproteins can help absorbed carotenoids target various tissues (Bhosale and Bernstein, 2007). As a lipid transfer protein, apolipoprotein D (ApoD) had higher expression levels in the variety II individuals. Jin et al. found that ApoD had a higher expression level in the wild-type E. carinicauda (Jin et al., 2021). In addition, the ABC transporter member 14 and solute carrier family 45 member 3 showed higher expression in variety II individuals. The pigmentation transport function of ABC transporters has been verified in silkworm, Drosophila and Harmonia axyridis (Tatematsu et al., 2011; Tsuji et al., 2018; Zuber et al., 2018; Shirk et al., 2023). Solute carriers (SLCs) are a group of membrane transport proteins, comprising over 300 members whose primary function is to facilitate the transport and absorption of various substrates across biofilms (Saikrithi et al., 2020; Girardin, 2022). Multiple SLC family genes were differentially expressed in different developmental stages of yellow mutant rainbow trout (Wu et al., 2022).