Transcriptome analysis in hepatopancreases reveals the response of domesticated common carp to a high-temperature environment in the agricultural heritage rice–fish system

Qingtian paddy field carp (PF-carp) is a local carp cultivated in the paddy field of Qingtian, Zhejiang. This rice–fish co-culture system has been recognized as one of the Globally Important Agriculture Heritage Systems (GIAHS). PF-carp has been acclimatized to the high-temperature environment of shallow paddy fields after several centuries of domestication. To reveal the physiological and molecular regulatory mechanisms of PF-carp, we chose to use 28°C as the control group and 34°C as the treatment group. We measured biochemical parameters in their serum and hepatopancreases and also performed transcriptome sequencing analysis. Compared with the control group, the serum levels of malondialdehyde (MDA), glucose (GLU), glutathione peroxidase (GSH-Px), catalase (CAT), alanine aminotransferase (ALT), and aspartate aminotransferase (AST) show no significant change. In addition, superoxide dismutase (SOD), GSH-Px, and CAT also show no significant change in hepatopancreases. We identified 1,253 differentially expressed genes (DEGs), and their pathway analysis revealed that heat stress affected AMPK signaling pathway, protein export, and other biological processes. It is worth noting that protein processing in the endoplasmic reticulum (ER) was the most significantly enriched pathway identified by the Kyoto Encyclopedia of Genes and Genomes (KEGG) and gene set enrichment analysis (GSEA). Significantly higher levels of HSP40, HSP70, HSP90, and other ubiquitin ligase-related genes were upregulated. In summary, heat stress did not lead to tissue damage, inflammation, oxidative stress, and ER stress in the hepatopancreases of PF-carp. This study provides valuable insights into the adaptation mechanism of this species to the high-temperature environment of paddy fields.


Introduction
Modern human civilization has developed based on the successful domestication of various plants and animals (Diamond, 2002).The earliest truly domesticated fish is the common carp (Cyprinus carpio) (Balon, 2004).Influenced by geographic, cultural, and other factors, the common carp has been domesticated into a variety of local varieties that are widely distributed in different farming systems (ponds or paddy fields).Due to the disparity in altitude and the few plains in Qingtian, Zhejiang, humans have domesticated the most successful paddy-farmed fish in this region and have been named Qingtian paddy field carp (PFcarp).In Qingtian, a system of rice-fish symbiosis has been created by the integration of PF-carp and rice farming.As one of the first Globally Important Agricultural Heritage Systems (GIAHS), this system was acknowledged by the Food and Agriculture Organization (FAO) in 2005 (Lu and Li, 2006).
From rivers and lakes to paddy fields, the environment in which common carp live has changed dramatically.Common carp is a natural demersal fish inhabiting in the lower layer of rivers and lakes, where the temperature is more stable and less susceptible to external environmental influences.However, paddy fields have a small water body and are shallow (5-25 cm) compared to rivers and lakes (Chen, X. et al., 2021).The shallow paddy water environment is often hot in the summer afternoons due to the intense sunlight.According to our monitoring, the summer temperature of the paddy field where PFcarp lives is approximately 34 °C.Yet, PF-carp has been domesticated The values presented are the sum of the means and standard deviations (mean ± SD) of three replicates.Values in the same column with different lowercase letters indicate significant differences (p < 0.05).G, G0, G2, G6, G12, and G24 represent the control group and heat stress at 0, 2, 6, 12, and 24 h, respectively.The values presented are the sum of means and standard deviations (mean ± SD) of three replicates.Values in the same column with different lowercase letters indicate significant differences (p < 0.05).G, G0, G2, G6, G12, and G24 represent the control group and heat stress at 0, 2, 6, 12, and 24 h, respectively.The values presented are the sum of means and standard deviations (mean ± SD) of three replicates.Values in the same column with different lowercase letters indicate significant differences (p < 0.05).G, G0, G2, G6, G12, and G24 represent the control group and heat stress at 0, 2, 6, 12, and 24 h, respectively.
in this environment for more than 12 centuries (Xie et al., 2011).Therefore, we hypothesize that they have been domesticated with special physiological and molecular regulatory mechanisms to adapt to the high-temperature environment of shallow paddy fields in Qingtian.
In order to investigate the adaptation mechanism of this species to the high-temperature environment of shallow paddy fields, we measured its biochemical parameters and performed RNA-seq analyses.The main objectives of our study were to characterize 1) the physiological changes in PF-carp in response to the hightemperature environment of shallow paddy fields and 2) the major signaling pathways and genes involved in the adaptation of PF-carp to the high-temperature environment of shallow paddy fields.

Ethical statement
The experiments were conducted in accordance with the Guidelines for the Care and Use of Laboratory Animals in China.The animals used in this study were cultured and euthanized following the terms approved by the Institutional Animal Care and Use Committee at Shanghai Ocean University (Shanghai, China) with approval number: SHOU-DW-2018-026.

Animals
Fifty-four healthy PF-carp juveniles (weighing 104.69 ± 3.08 g and measuring 14.65 ± 0.46 cm in length) used in this experiment were procured from Yugong Ecological Agricultural Technology Co., Ltd (Qingtian, Lishui, Zhejiang, China).They were then transported to the PF-Carp Research Center (Qingtian, Lishui, Zhejiang, China) for a 7-day acclimation period.They were randomly divided into three circular tanks (18 fish per tank).Furthermore, they were acclimatized in laboratory settings with the aerating water maintained at 28 °C ± 0.5 °C and dissolved oxygen levels of approximately 7 mg/L.The water was changed daily, and the fish were given artificial food twice daily.

Experimental design and sample collection
We randomly selected nine individuals in three tanks after 7 °days of acclimatization to serve as the control group.Then, the other experimental fish were elevated from 28 °C to 34 °C at a rate of 1 °C per hour and maintained at 34 °C for 24 h.
After maintaining the experimental temperature for 0, 2, 6, 12, and 24 h, samples were collected at each time point.Three PF-carp were randomly selected from each of the three tanks and anesthetized with MS-222 (300 mg/L) prior to sampling.A syringe was used to draw blood from the caudal vessel.Then, the fish were immediately dissected, and their hepatopancreases were collected for examination.The obtained samples were immediately frozen in liquid nitrogen and stored at −80 C until subsequent use.The same sampling procedure was applied to the control group.

Biochemical parameter determination
The changes in superoxide dismutase (SOD), malondialdehyde (MDA), glutathione peroxidase (GSH-Px), catalase (CAT), glucose (GLU), triglyceride (TG), alanine aminotransferase (ALT), and aspartate aminotransferase (AST) levels in the serum were determined.The supernatants of hepatopancreatic tissue homogenate were used for oxidative stress analysis, including SOD, MDA, GSH-Px, and CAT.According to the standard protocols, all the biochemical parameters were determined using reagent kits (Jiancheng Institute of Biotechnology, Nanjing, China).

Transcriptome sequencing
In this study, we selected hepatopancreatic tissues subjected to 6 h of heat stress for RNA-seq as the GM group.The TRIzol reagent (Invitrogen, Carlsbad, CA, United States) was used to isolate RNA of hepatopancreases (n = 3 per group); genomic RNA was removed using RNase I (Takara, Shanghai, China).Bioanalyzer 2100 (Agilent Technologies, United States) was used to determine RNA quality, and then ND-2000 (NanoDrop Technologies, United States) was used to quantify RNA.OD 260/280 ≥ 1.8 and OD 260/230 ≥ 1 were used for sequencing libraries.
These RNAs were reversed into cDNAs after a quality control process, and sequencing was carried out by Major Bioinformatics Technology Co., Ltd.(Shanghai, China) on Illumina NovaSeq 6000 (Illumina, United States).All transcriptome datasets can be found in the Sequence Read Archive (SRA) of the National Center for Biotechnology Information (NCBI).The accession number of the GM group was PRJNA1002641.RNA-seq results from a concurrent experiment were used as a control group (GC group), and its accession number was PRJNA971384.

Transcriptome analysis
SeqPrep and Sickle software programs were used to eliminate lowquality raw reads, and then the clean reads were mapped to the genome of common carp (accession number: ERS541549) (Xu et al., 2014) using HISAT 2 software and then aligned using StringTie.The levels of gene expression were estimated using transcripts per million reads (FPKM) in the present study.Genes with |log 2 FoldChange| ≥ 2 and p-adjust <0.05 were regarded as differentially expressed genes (DEGs) using DESeq2.Furthermore, GOATOOLS software and R package were used to perform GO and KEGG enrichment analyses, respectively.The gene sets of the KEGG pathway were used to perform further GSEA using GSEA 3.0 software.

Real-time quantitative PCR validation
We randomly chose nine DEGs for real-time quantitative PCR (RT-qPCR) verification to assess the accuracy of the transcriptome sequencing results.Primer Premier v.6.0 was used to design primers for DEGs (Supplementary Table S1).RT-qPCR was performed using the ChamQ SYBR Color qPCR Master Mix (2X) (Novozymes Bio, Nanjing, China) on a real-time fluorescence quantitative PCR instrument (ABI 7300, United States).The PCR conditions were as follows: 95 °C for 5 min, followed by 40 cycles of 95 °C for 5 s, 55 °C for 30 s, and 72 °C for 40 s.Gene expression levels were standardized relative to β-actin and calculated using Ct values (2 −ΔΔCT ).

Statistical analysis
All statistical analyses were carried out using SPSS 19.0.All results were analyzed using one-way ANOVA.Values are presented as mean ± standard deviation (mean ± SD). p < 0.05 was considered significant.In addition, GraphPad Prism 9 was used to display the results.

Changes in biochemical indicators after heat stress
When PF-carp was maintained at 34 °C, significant differences in the activity of SOD in serum (p < 0.05)were observed between the control and treatment groups at 6 h and 12 h.For the TG content in serum, significant changes were also observed to be higher in the treatment group than the control group at 12 h and 24 h (p < 0.05).However, there were no significant changes in the contents of MDA and GLU and the activities of GSH-Px, CAT, ALT, and AST in serum (p > 0.05) (Tables 1, 2).
The MDA content in hepatopancreases of the treatment group was significantly higher than that of the control group at 24 h (p < 0.05).In contrast, the activities of SOD, GSH-Px, and CAT were not significantly changed (p > 0.05) (Table 3).

Overview of RNA-seq
We used six samples for transcriptome sequencing with three replicates for the hepatopancreatic tissue in the two distinct groups.After filtering low-quality reads, 272,933,850 clean reads were obtained.Parameter statistics of clean reads among two groups were Q20: 97.84%-98.21%;Q30: 93.88%-94.62%;GC content: 48.4%-50.43%;and error rate: 0.0246%-0.0254%.The mapping of clean reads to the reference genome was also carried out, and the total mapped ratio was in the range of 78.07%-79.87%(Table 4).

Identification of DEGs
In our study, 1,253 DEGs were identified, including 446 upregulated genes and 807 downregulated genes (Figure 1).Clustering analysis indicated that the GC and GM groups were clustered separately, with quite different expression patterns in each group (Figure 2).

Verification of RNA-seq using RT-qPCR
In order to determine the accuracy and reliability of RNA-seq, nine genes were randomly chosen for RT-qPCR validation.In both RNA-seq and RT-qPCR, the chosen genes had consistent expression patterns (Figure 6).

Discussion
Domestication is a widely known example of artificial selection and has helped understand some of the most extreme within-species phenotypic variations over the years (Hoglund et al., 2020).PF-carp has been domesticated by humans for more than 1,200 years in paddy fields, developing phenotypic traits and genetic structure to adapt to the environment of the paddy fields (Ren et al., 2018;Qi et al., 2020).The liver of fish plays a crucial role in metabolic and immunological processes (Nakamura and Nishina, 2009).Many studies have reported the effect of heat stress on the structure and function of the liver (Li B et al., 2019;Dettleff et al., 2022;Yang et al., 2022).However, underlying physiological and molecular regulatory mechanisms of hepatopancreases in PF-carp responses to the hightemperature environment of shallow paddy fields remain elusive.Therefore, to better understand the adaptation mechanism of PFcarp to the paddy field environment after extensive domestication, we assessed the pertinent physiological parameters and performed RNA-seq on PF-carp that were under heat stress.
In our current study, only SOD activity and TG content showed significant changes in serum, and MDA content in the hepatopancreas also showed significant changes after undergoing 24 h of hightemperature stress at 34 °C, whereas the rest of the various enzyme activity indexes measured in serum and hepatopancreas did not show any significant changes.Fish have been reported to produce large amounts of ROS at high temperatures (Messina et al., 2023).In fish, SOD is the most important system that is used as the first line of defense against oxidative stress to remove the ROS that is released as a response to heat stress (Wang et al., 2019).Free radicals attack unsaturated fatty acids in the cell membrane to produce MDA (Papadimitriou and Loumbourdis, 2002).However, only SOD activity in the serum antioxidant defense system showed significantly changes, and MDA content in the hepatopancreas was significantly upregulated after 24 h of heat stress, suggesting that the body of PF-carp does not undergo stress when experiencing heat stress in rice paddy fields.In contrast, the antioxidant systems of Acipenser baerii (A.baerii) (Yang et al., 2021), rainbow trout (Oncorhynchus mykiss) (Li et al., 2022), and pikeperch (Sander lucioperca) (Li, C. et al., 2019) are significantly activated after thermal stress to maintain homeostasis in their internal environment.
TG is the main form of intracellular fat in fish and is critical for the storage and supply of energy (Komprda et al., 2014).GLU is the primary source of energy in fish that is produced by the digestion and absorption of glucose from meals and the breakdown of glycogen in the liver (Mommsen et al., 1999).TG levels increased significantly during heat stress in PF-carp, while the GLU content did not significantly change, indicating that PF-carp does not need to consume energy during heat stress.Certain particular enzymes may be released into the blood by damaged tissues or organs (Islam et al., 2021).It is reported that cytolysis and enzyme leakage into the bloodstream could lead to an increase in ALT and AST levels, which suggests liver damage (Bacchetta et al., 2014).In our findings, serum levels of ALT and AST in PF-carp during heat stress showed no significant changes, suggesting that heat stress did not affect the dysfunction of the hepatopancreas.In addition, Japanese flounder (Paralichthys olivaceus) (Han et al., 2023), pikeperch (Sander lucioperca) (Chen, Y. et al., 2021), and largemouth bass (Micropterus salmoides) (Zhao et al., 2022) often show damage to their tissues after experiencing heat stress.
RNA-seq has been well established and used to study the impact of high temperatures on fish, such as channel catfish (Ictalurus punctatus) (Tan et al., 2019), turbot (Scophthalmus maximus) (Huang, Z. et al., 2022), and blunt snout bream (Megalobrama amblycephala) (Li, B. et al., 2019).In this novel study, RNA-seq was used to determine the molecular changes in hepatopancreases of PF-carp.Through GO enrichment analysis, we reported that protein processing was significantly enriched.Furthermore, protein processing in the endoplasmic reticulum was the most significantly enriched pathway identified via KEGG enrichment analysis.According to the GSEA results, it was demonstrated that protein processing in the endoplasmic reticulum was the most significantly enriched and upregulated gene set.The endoplasmic reticulum is essential for intracellular calcium homeostasis, modification, and transport, as well as protein synthesis and folding (Kang and Jeon, 2021).The change in temperature could cause ER stress and protein misfolding (Tang et al., 2023).Heat shock proteins (HSPs) aid the folding and function of many proteins and could prevent protein misfolding (Lanneau et al., 2010).In addition, misfolded proteins through endoplasmic reticulum-associated protein degradation (ERAD) can be promoted by the coordination between the HSPs and ubiquitin ligase (Bozaykut et al., 2014;Kang and Jeon, 2021).In our study, HSP40, HSP70, HSP90, and various genes related to ubiquitin ligase were significantly upregulated.Therefore, we hypothesized that PF-carp removes abnormal proteins from the body mainly through protein processing in the endoplasmic reticulum during acclimatization to the high-temperature environment of shallow paddy fields.This is consistent with the way in which Atlantic salmon (Salmo salar) GC1-GC3 are the three parallel experimental group samples of the control group.GM1-GM3 are the three parallel experimental group samples of the heat stress group.
Frontiers in Physiology frontiersin.org(Shi et al., 2019), rainbow trout (Oncorhynchus mykiss) (Li et al., 2017), and largemouth bass (Micropterus salmoides) (Zhao et al., 2022) maintain cellular homeostasis during heat stress.Furthermore, GSEA revealed that the insulin, AMPK, FoxO , Tolllike receptor, and TNF signaling pathways were downregulated and insulin, AMPK, and FoxO signaling pathways could regulate body energy metabolism (Glauser and Schlegel, 2007;Wang et al., 2015;Schell et al., 2021).The insulin signaling pathways could be activated and take part in regulating glucose production in the liver.Additionally, it plays a role in the absorption of glucose into fat and muscle cells, helping maintain the body's glucose balance.(Suren Garg et al., 2023).A variety of physiological stimuli could activate the AMPK signaling pathway, such as glucose deprivation and oxidative stress, which may lead to a reduction in the cellular energy level and an increase in the AMP/ATP ratio (Schultze et al., 2012).Many genes in the FoxO signaling pathway are involved in the production of fat and glucose, and their upregulated expression can stimulate the production of these substances (Gross et al., 2009).It is reported that when fish are subjected to stress, the organism spends a large amount of energy to protect itself from external stresses (Petitjean et al., 2019).For example, gilthead sea bream (Sparus aurata L.) will improve mitochondrial metabolism when it faces stress (Bermejo-Nogales et al., 2014).However, these signaling pathways associated with metabolism were downregulated in our results, suggesting that PF-carp does not need to consume a lot of energy to adapt to the high-temperature environment of shallow paddy fields, which is also consistent with the results of our physiological parameters.Therefore, we believe that it is well adapted to paddy fields.
When the organism is exposed to heat stress, it activates a variety of immunomodulatory pathways, such as the Toll-like receptor and TNF signaling pathways.The activation of these immune pathways mediates the inflammatory response and reduces the damage caused by heat stress (Basu et al., 2015;Huang, T. et al., 2022).However, these immune pathways were downregulated when PF-carp was subjected to heat stress, suggesting that the organism did not initiate immune regulation.Pikeperch (Sander lucioperca) (Liu et al., 2022), rainbow trout (Oncorhynchus mykiss) (Guo et al., 2023), and grass carp (Ctenopharyngodon idella) (Zhang et al., 2022) activate these immune pathways to fight against damage when exposed to heat stress.
In a word, to maintain cellular homeostasis in PF-carp, protein processing in the endoplasmic reticulum plays a key role when it is subjected to high temperature stress in shallow paddy fields.In addition, the organism does not produce a stress response, nor does it consume a large amount of energy and trigger an inflammatory response to withstand any harm caused by heat stress.Instead, it adapts well to the high-temperature environment of the paddy field.

FIGURE 1
FIGURE 1Volcano plot of DEGs in GM vs. GC groups.Each dot represents one gene, red color indicates significantly upregulated genes, blue color indicates significantly downregulated genes, and gray color indicates non-significantly differentially expressed genes.

FIGURE 3
FIGURE 3GO enrichment analysis of DEGs in GM vs. GC groups.Colors ranging from blue to red indicate GO term significance from low to high.The size of the dot represents the number of genes enriched from more to less.

FIGURE 4
FIGURE 4 KEGG enrichment analysis of DEGs in GM vs. GC groups.Colors ranging from blue to red indicate the KEGG pathway significance from low to high.The size of the dot represents the number of genes enriched from more to less.

FIGURE 6
FIGURE 6Comparison of gene expression between RT-qPCR and RNAseq, with gene names on the X-axis and log 2 value of relative gene expression on the Y-axis.

TABLE 1
Effects of the oxidative index of the enzymatic type in the serum of PF-carp during heat stress.

TABLE 2
Effects of metabolism and hepatopancreatic injury in the serum of PF-carp during heat stress.

TABLE 3
Effects of the oxidative index of the enzymatic type in the hepatopancreas of PF-carp during heat stress.

TABLE 4
RNA-seq library sequencing data statistics.