Transcriptome Analysis Reveals the Effect of Stocking Density on Energy Metabolism in the Gills of Cherax quadricarinatus under Rice-Crayfish Co-Culture

Stocking density is a crucial factor affecting productivity in aquaculture, and high stocking density is a stressor for aquatic animals. In this study, we aimed to investigate the effects of stocking densities on oxidative stress and energy metabolism in the gills of Cherax quadricarinatus under rice-crayfish farming. The C. quadricarinatus were reared at low density (LD), medium density (MD), and high density (HD) for 90 days. The results showed that the superoxide dismutase (SOD), catalase (CAT), glutathione (GSH), and malondialdehyde (MDA) levels were higher in the HD group than those in the LD group. Transcriptomic analysis revealed 1944 upregulated and 1157 downregulated genes in the gills of the HD group compared to the LD group. Gene ontology (GO) enrichment analysis indicated that these differentially expressed genes (DEGs) were significantly associated with ATP metabolism. KEGG (Kyoto Encyclopedia of Genes and Genomes) analysis also showed that high stocking density resulted in the dysregulation of oxidative phosphorylation. Furthermore, high stocking density upregulated six lipid metabolism-related pathways. Overall, our findings, despite the limited number of samples, suggested that high stocking density led to oxidative stress and dysregulation of energy metabolism in the gills of C. quadricarinatus under rice–crayfish co-culture. Alteration in energy metabolism may be an adaptive response to adverse farming conditions.


Introduction
Integrated rice-fish farming is a promising ecological production model that combines rice planting and aquaculture in the same space, with the potential to enhance efficiency and reduce the environmental impact of agriculture [1]. This approach can be used to effectively utilize the unique spatial environment and water resources of rice fields, significantly reducing the use of pesticides and fertilizers [2]. Furthermore, the circular flow of energy and material resources in rice field ecosystems promotes mutual benefits and safe production of both rice and aquatic animals [3]. Notably, in China, integrated rice-fish farming practices have remarkably expanded, covering a vast area of 39.6612 million mu and producing nearly 20 million tons of rice and 3.5569 million tons of aquatic products in 2021 [4]. Several models of integrated rice-fish farming have been developed in China, including rice-carp, rice-turtle, rice-crayfish, rice-crab, rice-loach, and rice-snail [5]. Among these, the rice-crayfish co-culture is the most extensively applied and highest-yielding integrated rice-fish farming model in China [4,6]. Currently, research on rice-crayfish co-culture primarily focuses on its ecological benefits. For instance, rice-crayfish co-culture In addition, the CAT activity in the medium stocking density (MD) group was significantly higher than that in the LD group (p < 0.05). However, there was no significant difference in glutathione peroxidase (Gpx) activity in the gills of crayfish among the different density groups ( Figure 1D). The total antioxidant capacity (T-AOC) level was increased in the MD group, while it was decreased in the HD and significantly lower than that in the MD group (p < 0.05; Figure 1E). Compared to the LD group, malondialdehyde (MDA) content was markedly increased in the MD and HD groups (p < 0.05; Figure 1F).

Changes in Oxidative Stress Parameters
As shown in Figure 1A-C, the activities of superoxide dismutase (SOD) and catalase (CAT) and reduced glutathione (GSH) content in the gills of the high stocking density (HD) group were significantly higher than those in the low stocking density (LD) group (p < 0.05). In addition, the CAT activity in the medium stocking density (MD) group was significantly higher than that in the LD group (p < 0.05). However, there was no significant difference in glutathione peroxidase (Gpx) activity in the gills of crayfish among the different density groups ( Figure 1D). The total antioxidant capacity (T-AOC) level was increased in the MD group, while it was decreased in the HD and significantly lower than that in the MD group (p < 0.05; Figure 1E). Compared to the LD group, malondialdehyde (MDA) content was markedly increased in the MD and HD groups (p < 0.05; Figure 1F).

Transcriptome Sequencing and Analysis of Differently Expressed Genes (DEGs) in Gills
RNA-seq data showed that, after filtering, 39,311,670 (99.76%) to 46,816,068 (99.43%) clean reads were obtained, with Q20 and Q30 values ranging from 97.45% to 97.89% and 92.75% to 93.87%, respectively, and GC content ranging from 37.27% to 41.33% (Table 1). The total mapping ratio ranged from 90.45% to 95.16% (Table 1). These results indicated that the sequencing data of the gills from the LD and HD groups of C. quadricarinatus are reliable.
The principal component analysis (PCA) result showed that the LD and HD groups were obviously divided into two categories (Figure 2A), indicating that different stocking densities had a significant effect on gene expression patterns in the gill tissue. Analysis of DEGs between the LD and HD groups revealed a total of 3101 genes, with 1944 upregulated genes and 1157 downregulated genes in the gill tissue ( Figure 2B (Table 1). The total mapping ratio ranged from 90.45% to 95.16% (Table 1). These results indicated that the sequencing data of the gills from the LD and HD groups of C. quadricarinatus are reliable. The principal component analysis (PCA) result showed that the LD and HD groups were obviously divided into two categories (Figure 2A), indicating that different stocking densities had a significant effect on gene expression patterns in the gill tissue. Analysis of DEGs between the LD and HD groups revealed a total of 3101 genes, with 1944 upregulated genes and 1157 downregulated genes in the gill tissue ( Figure 2B,C). Radar plots depicting the data further demonstrated that, among the top 20 DEGs, thirteen genes were upregulated, and seven genes were downregulated in response to changes in stocking density ( Figure 2D). upregulated, and seven genes were downregulated in response to changes in stocking density ( Figure 2D). Note: Q20 and Q30, the base quality scores, were no less than 20 and 30, respectively, for clean reads. GC, GC content of clean reads.

Kyoto Encyclopedia of Genes and Genomes (KEGG) Enrichment Analysis of DEGs
The KEGG enrichment analysis showed that 874 DEGs were enriched in 136 pathways of the five KEGG A classes: metabolism, organismal systems, cellular processes, genetic information processing, and environmental information processing ( Figure 4A). The top 20 enriched pathways are listed in Figure 4B, with many pathways related to metabolic functions, such as oxidative phosphorylation, arachidonic acid metabolism, and biosynthesis of unsaturated fatty acids. Protein synthesis and processing, including pathways related to ribosomes and protein processing in the endoplasmic reticulum, were also affected to varying degrees by different stocking densities.
The KEGG enrichment analysis showed that 874 DEGs were enriched in 136 pathways of the five KEGG A classes: metabolism, organismal systems, cellular processes, genetic information processing, and environmental information processing ( Figure 4A). The top 20 enriched pathways are listed in Figure 4B, with many pathways related to metabolic functions, such as oxidative phosphorylation, arachidonic acid metabolism, and biosynthesis of unsaturated fatty acids. Protein synthesis and processing, including pathways related to ribosomes and protein processing in the endoplasmic reticulum, were also affected to varying degrees by different stocking densities.
The ATP metabolic process was significantly affected by stocking density (p.adj = 0.00005; Figure 5B). HD treatment resulted in 112 DEGs, including 58 upregulated and 54 downregulated genes.

Changes in Lipid Metabolism-Related Pathways
Gene set enrichment analysis (GSEA) results showed that stocking density notably impacted lipid metabolism in the gills of C. quadricarinatus. In the HD group, six lipid metabolism-related pathways were upregulated, including arachidonic acid metabolism, biosynthesis of unsaturated fatty acids, ether lipid metabolism, fatty acid elongation, glycerolipid metabolism, and glycerophospholipid metabolism ( Figure 6). The ATP metabolic process was significantly affected by stocking density (p.adj = 0.00005; Figure 5B). HD treatment resulted in 112 DEGs, including 58 upregulated and 54 downregulated genes.

Changes in Lipid Metabolism-Related Pathways
Gene set enrichment analysis (GSEA) results showed that stocking density notably impacted lipid metabolism in the gills of C. quadricarinatus. In the HD group, six lipid metabolism-related pathways were upregulated, including arachidonic acid metabolism, biosynthesis of unsaturated fatty acids, ether lipid metabolism, fatty acid elongation, glycerolipid metabolism, and glycerophospholipid metabolism ( Figure 6).

Changes in Lipid Metabolism-Related Pathways
Gene set enrichment analysis (GSEA) results showed that stocking density notably impacted lipid metabolism in the gills of C. quadricarinatus. In the HD group, six lipid metabolism-related pathways were upregulated, including arachidonic acid metabolism, biosynthesis of unsaturated fatty acids, ether lipid metabolism, fatty acid elongation, glycerolipid metabolism, and glycerophospholipid metabolism ( Figure 6).

Validation via RT-qPCR
Five genes were randomly selected to validate the accuracy of the transcriptome results via RT-qPCR analysis. The results showed that the RNA-seq data were significantly consistent with the RT-qPCR data (r = 0.908, p = 0.033; Figure 7A), further supporting the reliability of the transcriptome sequencing results.

Validation via RT-qPCR
Five genes were randomly selected to validate the accuracy of the transcriptome results via RT-qPCR analysis. The results showed that the RNA-seq data were significantly consistent with the RT-qPCR data (r = 0.908, p = 0.033; Figure 7A), further supporting the reliability of the transcriptome sequencing results.

Discussion
Stocking density is a key factor in aquaculture that determines production, economic benefits, and resource utilization. Consequently, high stocking densities are common in aquaculture practice. However, high stocking density may cause water quality deterioration, inhibition of aquatic animal growth, and decreased quality. Furthermore, high stocking density as a stressor may activate the stress response system and impair the physiological function and health of farmed aquatic animal [28]. Previous findings have indicated that aquatic animals under high stocking conditions experience oxidative stress, which often manifests as alterations in enzymatic antioxidants (e.g., SOD and CAT) and non-enzymatic antioxidants (e.g., GSH) [29,30]. For example, studies on O. mykiss and S.

Discussion
Stocking density is a key factor in aquaculture that determines production, economic benefits, and resource utilization. Consequently, high stocking densities are common in aquaculture practice. However, high stocking density may cause water quality deterioration, inhibition of aquatic animal growth, and decreased quality. Furthermore, high stocking density as a stressor may activate the stress response system and impair the physiological function and health of farmed aquatic animal [28]. Previous findings have indicated that aquatic animals under high stocking conditions experience oxidative stress, which often manifests as alterations in enzymatic antioxidants (e.g., SOD and CAT) and non-enzymatic antioxidants (e.g., GSH) [29,30]. For example, studies on O. mykiss and S. maxima have shown that high stocking densities reduced the levels of antioxidant parameters, such as SOD, GSH, and CAT, resulting in impaired antioxidant function [30,31]. Nevertheless, other studies on L. vannameiin, P. olivaceus, and Clarias gariepinus have demonstrated that antioxidants (e.g., SOD, CAT, and GSH) tend to increase with increased stocking density, indicating a positive response of the antioxidant system to stress [32][33][34]. Similarly, we found a significant increase in SOD and CAT activities and GSH content in the gills of C. quadricarinatus after being farmed for 90 days, suggesting that C. quadricarinatus responded to stress by activating its antioxidant defense system. Exposure to oxidative stress can lead to excessive amounts of reactive oxygen species (ROS) that attack lipids on the cell membrane and result in lipid peroxidation [35]. MDA is a well known lipid peroxidation product and serves as an indicator of oxidative stress levels [36]. Previous studies have shown that stress induced by high density results in an increased MDA level in various tissues [37]. In this study, the MDA contents in the gills of the MD and HD groups of C. quadricarinatus were significantly higher than those in the LD group. These findings are in agreement with previous studies conducted on O. mykiss [38] and Fenneropenaeus chinensis [39], which also suggested that high density stocking induced lipid peroxidation in the gills of C. quadricarinatus.
In general, aquatic animals can trigger physiological responses to adapt and to cope with changes in culture conditions, which require additional energy [40]. ATP plays a critical role as a direct energy source in many cellular process, including cell proliferation, metabolism, and survival [41]. Environmental fluctuations can modify the ATP metabolism in aquatic animals, ultimately affecting their cellular functions and physiological processes [42]. The gills of fish and crustaceans are pivotal organs responsible for gas exchange and osmoregulation, which are limited by the energy supply and are particularly sensitive to environmental changes [43]. Many previous studies have indicated that aquatic animals experience changes in their energy metabolism, owing to environmental stressors. For example, an early study suggested that ATP production in O. mossambicus gills was boosted in response to salinity challenges, allowing them to acclimate to adverse environments [44]. In metabolomic analysis, ATP synthase was found to increase in the gills of Cynoglossus semilaevis with increasing salinity [45]. Chronic exposure to tributyltin (TBT) decreased ATP levels by altering ATP enzymatic systems in gills of the common carp (Cyprinus carpio) [46]. The gills of Matrinxa (Brycon amazonicus) exhibit decreased ATP availability, owing to lipid peroxidation and ROS overproduction when exposed to air, leading to a change in their energy supply [47]. In the gills of Penaeus vannamei, ammonia stress reduces ATP content and cellular energy levels, which may alter physiological processes in gills [48]. In line with previous findings, in this study, transcriptome analysis revealed that the high stocking density significantly affected ATP metabolic process and ATP biosynthetic process, resulting in the upregulation of 58 genes and the downregulation of 54 genes in the gills of C. quadricarinatus. We speculate that the alteration in the ATP metabolic process may be attributed to oxidative stress induced by high stocking density, which could potentially compromise energy supply, gas exchange, and osmoregulation in the gills.
It is a well established fact that the production of cellular ATP is reliant on the process of oxidative phosphorylation that occurs within the mitochondrial electron transport chain [49]. This process is also accompanied by the generation of ROS, which contributes to both homeostatic signaling, as well as oxidative stress during pathological conditions [50]. In crustacean gills, several studies have shown that environmental stress can significantly impact oxidative phosphorylation. For instance, salinity stress alters the oxidative phosphorylation pathway in the gills of Eriocheir sinensis, leading to changes in energy-metabolism homeostasis [51]. Similarly, L. vannamei responded to nitrite stress by upregulating genes related to oxidative phosphorylation, resulting in enhanced energy metabolism in gills [52]. Chronic hypoxic stress affects the oxidative phosphorylation pathway in the gills of Macrobrachium nipponense, potentially coping with increased energy demand [53]. In our study, KEGG analysis demonstrated that oxidative phosphorylation was the most enriched pathway in the gills, with 34 upregulated and 41 downregulated genes. Meanwhile, cytochrome c oxidase, F-type ATPase, and V-type ATPase in the oxidative phosphorylation pathways were significantly upregulated in the HD group. These changes may represent an adaptation to supply more energy to cope with a crowded environment. Moreover, this may explain the alterations in ATP metabolic processes and oxidative stress under high stocking density.
Lipid metabolism plays a crucial role in the energy supply of organisms because fatty acids serve as a significant source of ATP [54]. In aquatic animals, this process is highly susceptible to environmental stress, particularly at high stocking density. An early study found that Sparus aurata mobilized its triglyceride content to fulfill the increased energy demand under high rearing density [55]. Similarly, O. niloticus adapts to high density conditions by altering fatty acid composition and stimulating lipid metabolic enzymes [56]. Transcriptomic analysis of Ctenopharyngodon idella also revealed that lipid metabolism is accelerated at high stocking density [57]. However, conflicting findings have been reported in the literature. For example, high stocking density in Brachymystax lenok led to the activation of energy metabolism while repressing lipid metabolism [58]. Similarly, in Micropterus salmoides, high stocking density tends to suppress fatty acid metabolism, biosynthesis of unsaturated fatty acids, and steroid biosynthesis [59]. It is worth noting that our data were consistent with former studies, as we found that six lipid metabolism-related pathways were upregulated in the gills of C. quadricarinatus under high stocking density conditions. These findings suggest that enhanced lipid metabolism may provide a coping mechanism to counteract stress caused by high stocking density.
To minimize experimental errors caused by individual variability, a pooled sample was employed for transcriptome sequencing in this study, which was a common approach in studies of aquatic animals, such as Chlamys farreri [60,61], Rachycentron canadum [62], Danio rerio [63], Plectropomus leopardus [64], E. sinensis [65], and Oryzias melastigma [66]. Despite the limited sample size of three (each sample contains fifteen individuals) in the study, it satisfied the requirements for transcriptome analysis, and numerous published reports have demonstrated that three samples are acceptable for transcriptome analysis [67][68][69][70][71][72][73][74][75]. Furthermore, the primary objective of this study was to investigate the impact of stocking density on C. quadricarinatus populations, and, therefore, individual effects were disregarded.

Crayfish, Experimental Design and Sampling
The experiment was conducted at the Jingjiang Farming Base of the Freshwater Fisheries Research Center (Jingjiang, China). The red-clawed crayfish with a 14.29 ± 1.05 g average initial weight was provided by the Zhejiang Freshwater Fisheries Research Institute (Huzhou, China). Nine rice-crayfish co-culture systems were used (Figure 8) in this experiment. The area of each rice-crayfish co-culture system was 400 m 2 , including a 360 m 2 rice field and a 40 m 2 canal refuge (ditch, 0.8 m in depth). In each co-culture system, the rice variety (Nangeng 5055) was planted in each rice field. density conditions. These findings suggest that enhanced lipid metabolism may provide a coping mechanism to counteract stress caused by high stocking density.
To minimize experimental errors caused by individual variability, a pooled sample was employed for transcriptome sequencing in this study, which was a common approach in studies of aquatic animals, such as Chlamys farreri [60,61], Rachycentron canadum [62], Danio rerio [63], Plectropomus leopardus [64], E. sinensis [65], and Oryzias melastigma [66]. Despite the limited sample size of three (each sample contains fifteen individuals) in the study, it satisfied the requirements for transcriptome analysis, and numerous published reports have demonstrated that three samples are acceptable for transcriptome analysis [67][68][69][70][71][72][73][74][75]. Furthermore, the primary objective of this study was to investigate the impact of stocking density on C. quadricarinatus populations, and, therefore, individual effects were disregarded.

Crayfish, Experimental Design and Sampling
The experiment was conducted at the Jingjiang Farming Base of the Freshwater Fisheries Research Center (Jingjiang, China). The red-clawed crayfish with a 14.29 ± 1.05 g average initial weight was provided by the Zhejiang Freshwater Fisheries Research Institute (Huzhou, China). Nine rice-crayfish co-culture systems were used (Figure 8) in this experiment. The area of each rice-crayfish co-culture system was 400 m 2 , including a 360m 2 rice field and a 40 m 2 canal refuge (ditch, 0.8 m in depth). In each co-culture system, the rice variety (Nangeng 5055) was planted in each rice field. In this experiment, the crayfish were randomly assigned to three groups: low stocking density (LD), medium stocking density (MD), and high stocking density (HD). Each group consisted of 3000, 6000, and 9000 crayfish, respectively, distributed among the three replicates. Throughout the 90-day experimental period, crayfish were fed once a day with a commercial feed containing 30% crude protein, 3% crude fat, 8% crude fiber, 18% crude ash, 1% total phosphorus, and 1 to 3.5% calcium. The paddy field was managed following conventional local agricultural practices. In this experiment, the crayfish were randomly assigned to three groups: low stocking density (LD), medium stocking density (MD), and high stocking density (HD). Each group consisted of 3000, 6000, and 9000 crayfish, respectively, distributed among the three replicates. Throughout the 90-day experimental period, crayfish were fed once a day with a commercial feed containing 30% crude protein, 3% crude fat, 8% crude fiber, 18% crude ash, 1% total phosphorus, and 1 to 3.5% calcium. The paddy field was managed following conventional local agricultural practices.
At the end of the experiment, 45 individuals from each group were randomly selected, and the gills were sampled under anesthesia on ice. For the analysis of oxidative stress parameters, the gills of five individuals were pooled into a single sample, whereas the gills of 15 individuals were combined into a sample for transcriptome sequencing. All samples were temporarily stored in liquid nitrogen and transferred to long-term storage at −80 • C. It is important to note that the use of crayfish in this study was approved by the Freshwater Fisheries Research Centre (Wuxi, China), and all experimental procedures were performed with proper regard for animal welfare requirements.

Measurement of Oxidative Stress Parameters
The frozen gill samples were thawed on an ice plate, and approximately 0.1g of each sample was homogenized in normal saline (0.86%). The homogenate was centrifuged at 3600 r/min and 4 • C for 10 min, and the supernatant was collected to determine the oxidative stress parameters. Oxidative stress parameters, including T-AOC, SOD, GSH, Gpx, CAT, and MDA, were measured using commercial kits. The T-AOC level was measured using the Ferric Reducing Ability of Plasma (FRAP) method [76], while SOD activity was tested via the water-soluble tetrazolium-1 (WST-1) method [77]. GSH and MDA contents were determined using the 5,5 -dithiobis-(2-nitrobenzoic acid) (DTNB) method [78,79] and the thiobarbituric acid (TBA) method, respectively. CAT activity was measured by analyzing the decomposition of H 2 O 2 [80], and total protein in the gills was measured using the bicinchoninic acid (BCA) method [81].

Transcriptome Sequencing and Analysis
Total RNA was isolated from the gill tissue using a commercial kit (Trizol reagent, Invitrogen, Carlsbad, CA, USA), following the manufacturer's protocol. After assessing the quality of the RNA, we enriched the mRNA by removing the rRNA. The enriched mRNA was fragmented into short fragments and reverse-transcribed into first-strand cDNA. Second-strand cDNA was synthesized using DNA polymerase I, RNase H, dNTP, and buffer. The cDNA was subjected to end-repair, A-base addition, adapter ligation, and PCR amplification to obtain the final library. The library was subsequently sequenced using an Illumina Novaseq 6000 by Gene Denovo Biotechnology Co. (Guangzhou, China), and the raw data were deposited in the open database (Sequence Read Archive database, No. PRJNA970695).
Raw reads from the sequencing machines were filtered using Fastp (version 0.18.0) [82] to obtain high-quality clean reads. The clean reads were utilized for subsequent assembly using the Trinity software (v2.1.1) and evaluated for integrity using BUSCO [83]. Annotations were assigned to the resulting unigenes using established databases, including Nr, KEGG, KOG, and SwissProt. The expression of unigenes was calculated and normalized to RPKM (Reads Per kb per Million reads).
Based on the gene expression in each sample, principal component analysis (PCA) was conducted to reveal correlations among the samples. The differentially expressed genes (DEGs) between LD and HD were analyzed using DESeq2 software (version 3.0) with a significance threshold p < 0.05 and |log 2 Fold Change (FC)| > 1 [84]. To gain insight into the functions of the DEGs, we performed functional enrichment analysis, referring to GO and KEGG databases. GO and KEGG enrichment analysis were performed using Goseq software (version 2.12) and KOBAS 2.0 software, respectively, and significantly enriched GO terms and KEGG pathways were defined by a hypergeometric test with a significance threshold q value < 0.05 (corrected by the Benjamin-Hochberg method). Furthermore, we employed gene set enrichment analysis (GSEA) to identify specific KEGG pathways between the LD and HD groups, and |Normalized enrichment score (ES)| > 1,Nominal p-value < 0.05, and FDR < 0.25 in each gene set were set as threshold values. All analyses and transcriptome charts generated using the OmicShare tools (www.omicshare.com, accessed on 1-23 April 2023)

Quantitative Real-Time PCR (RT-qPCR) Analysis
To validate the accuracy of the transcriptome sequencing, we performed RT-qPCR analysis. Total RNA was extracted from gill tissue using RNAiso Plus reagent (Takara, Beijing, China), along with chloroform and isopropyl alcohol (Sinopharm Chemical Reagent Co., Ltd., Shanghai, China). The RNA concentration was measured using a Nanodrop 2000 spectrophotometer (Thermo Scientific, Waltham, MA, USA), and quality was evaluated using the Agilent 4200 system (Agilent Technologies, Santa Clara, CA, USA) with RNA integrity number (RIN) > 7. After quality assessment, RNA was reverse-transcribed into cDNA using the primeScript™ RT reagent kit (Takara). In brief, 1 µg of RNA was mixed into gDNA Eraser (Takara) to remove the genomic DNA at 42 • C for 2 min. After that, the reaction mixture was supplemented with 1 µL of PrimeScript RT Enzyme Mix I, 1 µL of RT Primer Mix, 4 µL of 5× PrimeScript Buffer, and 4 µL of RNase Free dH 2 O. The mixture was then incubated at 37 • C for 15 min and 85 • C for 5 s.
The cDNA was then amplified using a commercial kit (Takara, RR820A) to determine the Cq value of the target gene. Briefly, the reaction mixture for RT-qPCR (total volume 25 µL) included 12.5 µL of TB Green Premix Ex Taq II (Takara), 2 µL of cDNA, 1 µL of forward and reverse specific primers each, and 8.5 µL of Rnase-free water. The RT-qPCR reaction was carried out under the following conditions: initial denaturation at 95 • C for 30 s, followed by 40 cycles of denaturation at 95 • C for 5 s, and extension at 60 • C for 1 min. After 40 cycles, the specificity of the PCR product was confirmed by analyzing the melting curve. Relative expression of the gene was calculated using the 2 −∆∆Cq method, with β-actin serving as the reference gene [85]. The specific primers for the target genes analysis in this study can be found in Table 2.

Statistical Analysis
Data in the study were analyzed using SPSS (version 25.0) and EXCL software (version 2016). All results are expressed as mean ± SE (standard error). The normal distribution of the data was analyzed using the Shapiro-Wilk test, and the variance homogeneity of the data was analyzed using the Levene test. Differences in oxidative stress parameters among different groups were analyzed using one-way analysis of variance (ANOVA), followed by the least significant difference (LSD) test. The correlation between RT-qPCR and RNA-seq was assessed using the Pearson test. Statistical significance was established when the p-value was less than 0.05.

Conclusions
This study identified, for the first time, changes in energy metabolism in the gills of crustaceans under high stocking density. High stocking density induced oxidative stress and lipid peroxidation in the gills of C. quadricarinatus. Further, we conducted transcriptome analysis to examine the molecular mechanisms underlying the stress caused by high stocking density in the gills. After 90 days of farming, we identified 3101 DEGs in the gill tissue between the LD and HD groups, with 1944 upregulated genes and 1157 downregulated genes. GO and KEGG enrichment revealed that the DEGs significantly enriched ATP metabolic process and oxidative phosphorylation pathway. Meanwhile, high stocking density caused the upregulation trend in six lipid metabolism-related pathways. The changes in ATP metabolic process, oxidative phosphorylation, and lipid metabolism may be an adaptive response to stress induced by high stocking density. Our study contributed to understanding the mechanism that C. quadricarinatus responded to the adverse conditions in the rice-crayfish farming system.