Transcriptomic Networks Reveal the Tissue-Specific Cold Shock Responses in Japanese Flounder (Paralichthys olivaceus)

Simple Summary Low temperature is an often overlooked stress that many fish face due to both natural and anthropogenic causes. Japanese flounder (Paralichthys olivaceus) is an economically important aquaculture fish species broadly cultivated in east Asia, mainly along the coast of the Bohai and Yellow Seas in China, as well as the coast of Korea and east Japan. Natural and cultivated P. olivaceus may suffer from cold stress during the winter months. In this study, modulated transcriptomic responses to 10 °C acute cold stress were investigated in the gills, hearts, livers, and spleens of P. olivaceus. Based on transcriptome and weighted gene coexpression network analysis, tissue-specific cold responsive modules (CRMs) were identified, which revealed a cascade of specific cellular responses to cold stress in different tissues. Our results illustrate the diverse and modulated regulation of the cellular process and stress response to low temperature, which provide essential insights for the conservation and cultivation of P. olivaceus in cold water. Abstract Low temperature is among the important factors affecting the distribution, survival, growth, and physiology of aquatic animals. In this study, coordinated transcriptomic responses to 10 °C acute cold stress were investigated in the gills, hearts, livers, and spleens of Japanese flounder (Paralichthys olivaceus), an important aquaculture species in east Asia. Histological examination suggested different levels of injury among P. olivaceus tissues after cold shock, mainly in the gills and livers. Based on transcriptome and weighted gene coexpression network analysis, 10 tissue-specific cold responsive modules (CRMs) were identified, revealing a cascade of cellular responses to cold stress. Specifically, five upregulated CRMs were enriched with induced differentially expressed genes (DEGs), mainly corresponding to the functions of “extracellular matrix”, “cytoskeleton”, and “oxidoreductase activity”, indicating the induced cellular response to cold shock. The “cell cycle/division” and “DNA complex” functions were enriched in the downregulated CRMs for all four tissues, which comprised inhibited DEGs, suggesting that even with tissue-specific responses, cold shock may induce severely disrupted cellular functions in all tissues, reducing aquaculture productivity. Therefore, our results revealed the tissue-specific regulation of the cellular response to low-temperature stress, which warrants further investigation and provides more comprehensive insights for the conservation and cultivation of P. olivaceus in cold water.


Introduction
Temperature, as the major environmental stimulating factor, has essential effects on almost all biological processes in organisms, affecting their physiological and metabolic activities, immune and antioxidative responses, growth and behavior, as well as migration, especially for poikilothermic teleost species [1,2]. In the context of global climate change, a large number of studies on thermal stress in fish have mainly focused on high-temperature challenge in conservation and fishery management [1]; compared with heat stress, lowtemperature exposure can also frequently occur in aquatic environments as the result of nature and/or anthropogenic causes [2].
A sharp reduction in ambient temperature is defined as cold shock, which has the potential to cause a rapid decline in body temperature, resulting in a cascade of physiological and behavioral responses [3]. Teleosts have evolved biochemical and physiological adaptations to low-temperature challenges. A conservative responsive mechanism exists in different fish species under chronic long-term low-temperature and acute short-term cold stress. For example, cold stress can result in the increase in endogenous reactive oxygen species (ROS) [4]. ROS generation may further contribute to DNA damage and apoptosis [5], which can help to maintain cellular homeostasis and play vital roles in the immune response [6,7]. Moreover, cold stress can also affect the survival of teleosts, cause mass mortality in aquaculture, and has important impacts on reproductive, metabolic, and developmental functions in teleosts [3].
In addition, many studies have pointed out that tissues involved in the physiological response and acclimation to low-temperature stress can exhibit tissue-specific characteristics [8][9][10]. For example, fish gills are directly exposed to aquatic environment, are the target for the temperature stress response, and are involved in many physiological processes, including respiration, ion-and osmo-regulation, immune response, and the acid-base balance [11]. The heart supplies power for the circulatory system to maintain the survival of the organism [12]. The liver plays an essential role in the metabolism, storage, and distribution of carbohydrates, proteins, and lipids in fish [13]. The spleen is mainly composed of lymph nodes and plays important roles in the hematopoietic, immune, and blood storage of fish [13]. Moreover, there are many functional pathways, including mitochondrial function, antioxidation, apoptosis, carbohydrate/lipid metabolism, DNA/RNA processing, and protein catabolism, that show temperature-dependent regulation [14][15][16][17]. Therefore, investigating the molecular mechanisms among the various tissues of fish that are affected by low temperature could provide a better understanding of how organisms adapt or respond to aquatic environmental challenges.
Recently, with the increasing abundance of sequencing data, weighted gene coexpression network analysis (WGCNA) is a powerful method that has been used to provide a comprehensive understanding of gene interaction networks. WGCNA focuses on the modulated regulation of the full-complement gene sets that comprise complex networks, transitioning data analysis from being gene-centric to systematically network-centric [18]. Therefore, the combination of transcriptome and network analysis will provide insightful implications about the tissue-specific responses of fish to dynamic environmental stresses through the functioning of cellular signaling pathways [19].
Japanese flounder (Paralichthys olivaceus) is an economically important aquaculture teleost species broadly cultivated in east Asia, mainly along the coast of the Bohai and Yellow Seas in China, as well as along the coast of Korea and east Japan [14]. Many studies have illustrated the effect of heat stress on the mechanisms of sexual differentiation, energy metabolism, neurosecretion, and stress-related gene regulation in P. olivaceus [14,[20][21][22][23]. However, as a temperate water fish, nowadays, temperature decline has also influenced the survival, growth, distribution, and reproduction of this species [24]. For example, in the natural habitat of P. olivaceus, there are several months of cold winter, with the lowest temperature being 0 • C [25], and a seasonal water mass, the Yellow Sea Cold Water Mass, appears during the summer half of the year in the central bottom water with low temperature (<10 • C) [26], which makes the P. olivaceus growth rate significantly decrease and reduces aquaculture productivity [25,27]. To illustrate the tissue-specific regulation under cold stress, in this study, P. olivaceus specimens were subjected to 10 • C acute cold shock stress to observe the cellular damage to the gill, heart, liver, and spleen tissues, and transcriptome data were comprehensively analyzed with WGCNA. The results present a fundamental reference for cold shock damage in different P. olivaceus tissues and provide vital insights for the conservation and cultivation of P. olivaceus in cold water.

Acute Cold Shock Challenge for Adult P. olivaceus
A total of 18 adult P. olivaceus specimens (average body weight 866 ± 166 g, average total length 33 ± 4.5 cm) were obtained from Nanshan market, Qingdao, China, and acclimated in seawater in a laboratory aquarium (18 • C, 14:10 h light:dark, 30 ppt salinity) for one week. The seawater was aerated and refreshed once per day. Then, the P. olivaceus specimens were randomly divided into two groups: a control group under a temperature of 18 • C, and an acute cold shock group, which were abruptly transferred to 10 • C sea water for a short term of 6 h. There were 9 specimens in each group with 3 individuals in each of the 3 tanks as replicates, and no fish died during the experiment. Circulating water refrigerators (RESUN, Shenzhen, China, CL650, 650 W, 1/4HP) were employed to maintain the cold temperature. After the cold shock challenge, P. olivaceus specimens were euthanized with MS-222, and the gill, heart, liver and spleen tissues from each specimen were sampled and stored at −80 • C with lipid nitrogen for RNA isolation and in Bouin's fluid for histological observation.

Histological Examination of P. olivaceus Tissues
To observe the histological changes under cold stress, the P. olivaceus tissues from Bouin's fluid were dehydrated following a successive ethanol gradient of 50%, 70%, 90%, 95%, and 100%. Then, the tissue samples were transparentized with xylene and ethanol mixture, embedded into paraffin, and sliced into a thickness of 5 µm. The samples were further dewaxed and stained with hematoxylin and eosin (H&E) (Solarbio, Beijing, China) following the traditional method [14]. Finally, the samples, which were sealed with neutral gum, were observed and photographed with a Nikon Eclipse TiU microscope (Nikon, Tokyo, Japan) for histological examination.

RNA Isolation, Library Preparation, and Sequencing
The gill, heart, liver, and spleen tissues of P. olivaceus were sampled and frozen with liquid nitrogen. Total RNA was isolated using a method with SparkZol Reagent (SparkJade, Jinan, China). The genomic DNA contamination was removed by RNase-free DNase I (TaKaRa, Beijing, China) treatment. The total RNA was then qualified and quantified by 1.5% agarose gel electrophoresis and spectrophotometry, respectively. The RNA samples were further employed for RNA sequencing library preparation and high-throughput sequencing. The RNA sequencing was conducted on an Illumina Hiseq 4000 platform at the Beijing Novogene company (Novogene, Beijing, China).

Transcriptome Analysis of P. olivaceus Tissues in Response to Cold Shock Stress
To investigate the gene expression profiles of P. olivaceus under cold shock stress, 24 transcriptome datasets from four tissues (gill, heart, liver, and spleen) of both the control and cold shock groups were analyzed. After trimming adaptors and removing low-quality reads, the clean reads were mapped against the P. olivaceus genome using the Hisat and StringTie pipeline [28] with default parameters. The fragments per kilobase of exon per million mapped reads (FPKM) values were used to estimate gene expression levels with StringTie [28]. Differentially expressed genes (DEGs) between the control group and the cold shock group were analyzed using edgeR, with a threshold q-value of <0.05 and |log 2 FoldChange(FC)| > 1.5. The R package clusterProfiler 4.0 was employed for Gene Ontology (GO) functional enrichment analysis. TBtools [29] was utilized to draw heatmaps with the log 2 FC values.

Gene Coexpression Network Construction and Functional Characterization
WGCNA [18] was conducted to characterize the modulated gene interaction patterns across different P. olivaceus tissue samples with WGCNA R library [18]. A gene dendrogram was employed to identify modules using the dynamic tree cut method (minimum module size = 200, cutting height = 0.99, and deepSplit = F). The intramodular connectivity (Kwithin) was then used to characterize the hub genes in each module, which represents the strong gene connection to other genes in the module. Each node (gene) usually connected to many other nodes (genes) as the connection edges with different weight values. Moreover, the cold responsive modules (CRMs) were identified based on the over-representation of tissue-specific DEGs using a hypergeometric test (p < 0.05). Furthermore, GO terms were enriched for genes in each module as the functional annotation by EnrichPipeline [30]. Cytoscape was used to visualize the coexpression network, which filtered the edges and focused on the nodes (genes) with the strongest edges [31].

Histological Observation of P. olivaceus Tissues after Cold Shock Challenge
After cold shock, there was a slight swelling of the lamellae and epithelial cells in P. olivaceus gills compared with those of the control group, and the supporting structure of the gill filament shrunk and the branchial lamellae were curlier (Figure 1a,b). This result was similar to that of four-finger threadfin (Eleutheronema tetradactylum) gills under cold shock, in which the branchial lamellae of E. tetradactylum shrunk, swelled, and curved [32]. There was no significant difference in the P. olivaceus hearts between the cold shock group and the control group (Figure 1c,d), which has also rarely been described in other fish species. In the livers, the hepatic sinusoids contracted and the gap between hepatocytes expanded under cold shock, which was possibly caused by the contraction of the hepatocytes (Figure 1e,f). There was also obvious shrinkage of and reduction in the lipid drops in the hepatic stellate cells, which suggested the consumption of fat under cold shock stress (Figure 1e,f). Under cold stress, morphological changes such as cavitation of the hepatocytes, contraction of the liver and hepatic blood sinuses, deepening of hepatocyte staining, severe damage and deformation of hepatocytes, as well as atrophy or disappearance of some hepatocyte nuclei, were also reported in the liver of E. tetradactylum [32]. Moreover, there was no significant difference in P. olivaceus spleens between the cold shock group and the control group (Figure 1g,h). Therefore, compared with at a medium temperature (18 • C), acute cold shock (10 • C) induced varying degrees of P. olivaceus tissue structure changes, especially in the gills and livers.

Tissue-Specific Transcriptomic Response of P. olivaceus under Cold Shock Stress
To characterize the transcriptomic response of P. olivaceus tissues under acute cold shock stress, 24 RNA-seq libraries from the gill, heart, liver, and spleen tissues were analyzed. Firstly, compared with the control group, a total of 276 (105 up-and 171 downregulated), 270 (80 up-and 190 downregulated), 140 (57 up-and 83 downregulated), and 457 (234 up-and 223 downregulated) differentially expressed genes (DEGs, |log 2 FC| ≥ 1.5 and q < 0.05) were identified in the gill, heart, liver, and spleen tissues, respectively (Figures 2a and S1). There were more tissue-specific DEGs than the few overlapped DEGs among the challenged tissues ( Figure 2b), with the gills and spleens sharing 20 upregulated and 53 downregulated DEGs (Figure 2b). For the overlapped DEGs, the single upregulated DEG for all the four tissues was period circadian protein homolog 2 (PER2) and the four downregulated DEGs shared by the four tissues were interleukin-10 receptor subunit beta (Il10rb), C-X-C motif chemokine 10 (CXCL10), putative nuclease HARBI1 (HARBI1), and serine/threonineprotein phosphatase (PPP3CA). Moreover, the abundant tissue-specific DEGs indicated the particular gene sets participating in the diverse cellular processes in P. olivaceus tissues.
Specifically, the tissue-specific upregulated DEGs were enriched in GO terms including "cell growth", "demethylase activity", "extracellular matrix," and "immune response" in the gills; "circulatory system", "oxidation", "extracellular region", and "cytoskeleton" in the hearts; "heme oxidation", "lipase inhibitor activity", and "cellular metabolic process" in the livers; and "transmembrane transport", "cell cycle", and "cell growth" in the spleens ( Figure 2c and Table S1). The tissue-specific downregulated DEGs were mainly enriched in GO terms such as "cell cycle", "response to stimulus", "cell division", and "chromatin" in the gills; "transport", "cytoskeleton", and "oxidation" in the hearts; "oxidation" and "lipid metabolism" in the livers; "cell cycle", "response to stress", "chromosome segregation", and "microtubule-based process" in the spleens (Figure 2c and Table S1). These regulated DEGs indicated that under the cold shock challenge, significant tissue-specific cellular process alterations and stress effects occur in P. olivaceus.

Tissue-Specific Transcriptomic Response of P. olivaceus under Cold Shock Stress
To characterize the transcriptomic response of P. olivaceus tissues under acute cold shock stress, 24 RNA-seq libraries from the gill, heart, liver, and spleen tissues were analyzed. Firstly, compared with the control group, a total of 276 (105 up-and 171 downregulated), 270 (80 up-and 190 downregulated), 140 (57 up-and 83 downregulated), and 457 (234 up-and 223 downregulated) differentially expressed genes (DEGs, |log2FC| ≥ 1.5 and q < 0.05) were identified in the gill, heart, liver, and spleen tissues, respectively (Figures 2a and S1). There were more tissue-specific DEGs than the few overlapped DEGs among the challenged tissues ( Figure 2b), with the gills and spleens sharing 20 upregulated and 53 downregulated DEGs (Figure 2b). For the overlapped DEGs, the single upregulated DEG for all the four tissues was period circadian protein homolog 2 (PER2) and the four downregulated DEGs shared by the four tissues were interleukin-10 receptor subunit beta (Il10rb), C-X-C motif chemokine 10 (CXCL10), putative nuclease HARBI1 (HARBI1), and serine/threonine-protein phosphatase (PPP3CA). Moreover, the abundant tissue-specific DEGs indicated the particular gene sets participating in the diverse cellular processes in P. olivaceus tissues. Specifically, the tissue-specific upregulated DEGs were enriched in GO terms including "cell growth", "demethylase activity", "extracellular matrix," and "immune response" in the gills; "circulatory system", "oxidation", "extracellular region", and "cytoskeleton" in the hearts; "heme oxidation", "lipase inhibitor activity", and "cellular metabolic process" in the livers; and "transmembrane transport", "cell cycle", and "cell growth" in the spleens (Figure 2c and Table S1). The tissue-specific downregulated DEGs were mainly enriched in GO terms such as "cell cycle", "response to stimulus", "cell division", and "chromatin" in the gills; "transport", "cytoskeleton", and "oxidation" in the hearts; "oxidation" and "lipid metabolism" in the livers; "cell cycle", "response to stress", "chromosome segregation", and "microtubule-based process" in the spleens (Figure 2c

Tissue Specific Coexpression Network of P. olivaceus DEGs under Cold Shock Stress
To further understand the modulated gene regulation and interaction among tissues of P. olivaceus under cold shock stress, 20,877 expressed genes from the 24 tissue transcriptome samples were used to perform WGCNA, and they were assigned into 18 modules with a size ranging from 288 to 3613 genes (Figures 3 and S2 and Table S2). The cold responsive modules (CRMs) were identified based on the over-representation of DEGs. As a result, 10 CRMs were identified (FDR < 0.05), including 5 modules enriching activated and inhibited DEGs each ( Figure 3 and Table S3). Among the CRMs, two, one, two, and two modules showed enriched upregulated DEGs (FDR < 0.05) for gills (turquoise and pink), hearts (yellow), livers (turquoise and cyan), and spleens (pink and red), respectively; while two, three, two, and one modules were enriched in downregulated DEGs (FDR < 0.05) for gills (green-yellow and purple), hearts (green-yellow, magenta, and green), livers (green-yellow and blue), and spleens (green-yellow), respectively ( Figure 3 and Table S3). Interestingly, among the 10 SRMs, 3 modules represented overlapped responsive pattern between the tissues, with the turquoise (gill and liver) and pink (gill and spleen) modules being upregulated, whereas the green-yellow module being downregulated for all four tissues ( Figure 3). This could confirm the diversified regulatory patterns as revealed by the transcriptome analysis ( Figure 2) among the P. olivaceus tissues under cold shock stress.

Tissue Specific Coexpression Network of P. olivaceus DEGs under Cold Shock Stress
To further understand the modulated gene regulation and interaction among tissues of P. olivaceus under cold shock stress, 20,877 expressed genes from the 24 tissue transcriptome samples were used to perform WGCNA, and they were assigned into 18 modules with a size ranging from 288 to 3613 genes (Figures 3 and S2 and Table S2). The cold responsive modules (CRMs) were identified based on the over-representation of DEGs. As a result, 10 CRMs were identified (FDR < 0.05), including 5 modules enriching activated and inhibited DEGs each ( Figure 3 and Table S3). Among the CRMs, two, one, two, and two modules showed enriched upregulated DEGs (FDR < 0.05) for gills (turquoise and pink), hearts (yellow), livers (turquoise and cyan), and spleens (pink and red), respectively; while two, three, two, and one modules were enriched in downregulated DEGs (FDR < 0.05) for gills (green-yellow and purple), hearts (green-yellow, magenta, and green), livers (green-yellow and blue), and spleens (green-yellow), respectively ( Figure 3 and Table S3). Interestingly, among the 10 SRMs, 3 modules represented overlapped responsive pattern between the tissues, with the turquoise (gill and liver) and pink (gill and spleen) modules being upregulated, whereas the green-yellow module being downregulated for all four tissues (Figure 3). This could confirm the diversified regulatory patterns as revealed by the transcriptome analysis ( Figure 2) among the P. olivaceus tissues under cold shock stress.

Coordinated Regulation of the Tissue-Specific P. olivaceus DEGs
In the CRMs, many DEGs represented tissue-specific regulation under cold shock challenge (Figure 2b), which resulted in diverse cellular responses in P. olivaceus tissues. These DEGs from enriched GO functions (Table S4) were interlinked in the networks (Figure 4). Therefore, the CRMs and DEGs as key functional nodes were investigated and discussed to illustrate the core gene interactions in the tissue-specific responses of P. olivaceus. Firstly, several gill-specific DEGs were interlinked in the networks. For example,

Coordinated Regulation of the Tissue-Specific P. olivaceus DEGs
In the CRMs, many DEGs represented tissue-specific regulation under cold shock challenge (Figure 2b), which resulted in diverse cellular responses in P. olivaceus tissues. These DEGs from enriched GO functions (Table S4) were interlinked in the networks (Figure 4). Therefore, the CRMs and DEGs as key functional nodes were investigated and discussed to illustrate the core gene interactions in the tissue-specific responses of P. olivaceus. Firstly, several gill-specific DEGs were interlinked in the networks. For example, as an upregulated DEG, neuropilin-2 (NRP2) could act as a receptor for entry into the epithelial and endothelial cells of the gills; integrin alpha-11 (ITGA11) is the receptor for collagen; and cell antigen CD34 is the adhesion molecule mediating the attachment of stem cells to the extracellular matrix (ECM) [33][34][35]. All these DEGs are from the "extracellular region" and "cell adhesion" functions ( Figures 4 and 5), as key mediators of cell-matrix and cell-cell adhesion, suggesting their essential roles in the ECM structure maintenance and cell matrix adhesion in the cellular homeostasis of P. olivaceus gill filaments and lamella under the cold shock stress. Moreover, the downregulated DEGs, mainly from the "cell cycle" and "cell division" functions, were interlinked in the gill-specific networks ( Figure 4). For instance, DNA polymerase (POLA2 and POLE2) and histone (H2B and H3) play essential roles in transcription regulation, DNA repair, DNA replication, and chromosomal stability; cyclin-dependent kinase inhibitor 1 (CDKN1A) may be involved in the p53/TP53-mediated inhibition of cellular proliferation in response to DNA damage; and protein mis12 is required for normal chromosome alignment and segregation and for kinetochore formation during mitosis [36][37][38]. These interlinked gill-specific DEGs suggested cellular damage in the gill filaments and lamina under cold shock challenges, which could also be observed from the H&E staining (Figure 1a,b). Similar regulated functions were also identified in the gills of P. olivaceus under heat shock, in which the upregulated functions were mainly enriched in the cellular response to stimuli, protein refolding, and regulation of apoptotic process, while the downregulated functions were mostly enriched in DNA repair and replication, as well as the cell cycle [14]. Interestingly, cyclin-dependent protein kinase activity was identified in P. olivaceus gills under both heat [14] and cold stress, which may suggest that the downregulation of the cell-cycle function may be a common response in P. olivaceus gills under both heat and cold shock stresses.
In the heart-specific upregulated DEGs, cellular-structure-related functions such as "actin filament" and "cytoskeleton" were enriched. For example, xin actin-binding repeatcontaining protein 2 (XIRP2) can protect actin filaments from depolymerization; leiomodin-3 (LMOD3) increases the rate of actin polymerization; and microtubule-associated protein 6 (MAP6) specifically shows microtubule cold-stabilizing activity [39][40][41][42] (Figures 4 and 5). All these DEGs indicated the activation of intracellular movements and membrane trafficking in the P. olivaceus heart, which could interact with the circulatory functional gene endothelin (EDN1, EDN2) to cope with cold stress (Figure 4). Moreover, superoxide dismutase (SOD3), which can protect the extracellular space from the toxic effects of reactive oxygen intermediates [43], was also interlinked in the network, suggesting antioxidative protection in the P. olivaceus cardiovascular system (Figures 4 and 5). In addition, many transport and cytoskeleton-related downregulated DEGs were interlinked in the networks. For example, tubulins (TUBA and TUBB) are the major constituents of microtubules; myosin (Myo16) is an actin-based motor molecules with ATPase activity, which can serve in intracellular movements; homeodomaininteracting protein kinase 2 (HIPK2) is involved in transcription regulation, p53/TP53-mediated cellular apoptosis, and regulation of the cell cycle [37,[44][45][46] (Figures 4 and 5). These heartspecific DEGs indicated that even with fine protection, cellular injury could still be found in P. olivaceus heart under cold stress. A sharp temperature decline can lead to impaired cardiac contractile function in fish, as the thin filaments in the cardiac muscle are less sensitive to Ca 2+ at lower temperatures [47]. For example, when compared at the respective physiological temperature and pH, rainbow trout (Oncorhynchus mykiss) cardiac actin-myosin ATPase activity showed more Ca 2+ sensitivity than that of rats, which may help O. mykiss to offset the cardioplegic effects of cold [47]. Moreover, after cold acclimation, myocardial fibrillary collagen content and/or connective tissue content increases, which may protect the myocardium from the increased hemodynamic stress due to pumping cold viscous blood [47]. In the heart-specific upregulated DEGs, cellular-structure-related functions such as "actin filament" and "cytoskeleton" were enriched. For example, xin actin-binding repeatcontaining protein 2 (XIRP2) can protect actin filaments from depolymerization; leiomodin-3 (LMOD3) increases the rate of actin polymerization; and microtubule-associated protein 6 (MAP6) specifically shows microtubule cold-stabilizing activity [39][40][41][42] (Figures 4 and 5). All these DEGs indicated the activation of intracellular movements and membrane trafficking in the P. olivaceus heart, which could interact with the circulatory functional gene endothelin (EDN1, EDN2) to cope with cold stress (Figure 4). Moreover, superoxide dismutase (SOD3), which can protect the extracellular space from the toxic effects of reactive oxygen intermediates [43], was also interlinked in the network, suggesting antioxidative protection in the P. olivaceus cardiovascular system (Figures 4 and 5). In addition, many transport and cytoskeleton-related downregulated DEGs were interlinked in the networks. For example, tubulins (TUBA and TUBB) are the major constituents of microtubules; myosin (Myo16) is an actin-based motor molecules with ATPase activity, which can serve in intracellular movements; homeodomain-interacting protein kinase 2 (HIPK2) is involved in transcription regulation, p53/TP53-mediated cellular apoptosis, and regulation of the cell cycle [37,[44][45][46] (Figures 4 and 5). These heart-specific DEGs indicated that even with fine protection, cellular injury could still be found in P. olivaceus heart under cold stress. A sharp temperature decline can lead to impaired cardiac contractile function in fish, as the thin filaments in the cardiac muscle are less sensitive to Ca 2+ at lower temperatures [47]. For example, when compared at the respective physiological temperature and There were fewer liver-specific DEGs in P. olivaceus livers than in other tissues from the network (Figure 2b and 4). For example, of the liver-specific upregulated DEGs, heme oxygenase (HMOX) catalyzes the oxidative cleavage of heme, which can protect against cold-injury-induced cell loss and damage [48], while in the downregulated DEGs, DNA polymerase (POLE) and cell division control protein (CDC45) are both required for the initiation of chromosomal DNA replication. Mitochondrial ribosome-associated GTPase 2 (GTPBP5) plays a role in the regulation of the mitochondrial ribosome assembly and of translational activity, all suggesting the inhibited function of the "cell cycle" and "cell division" [49][50][51][52]. In the livers of fish, high-temperature stress may affect apoptosis, DNA replication, protein metabolism, energy metabolism, lipid metabolism, immune, cell cycle, protein processing, and transport, as well as antioxidative responses [53][54][55], while low-temperature stress may affect signaling transduction, membrane fluidity, lipid metabolism, antioxidative responses, ion binding, macromolecule catabolism, mitochondrial enzymes related to transport, carbohydrate metabolism, and cell cycle endocrine system [56,57]. Interestingly, inositol polyphosphate-5-phosphatase A (INPP5A) and phosphatidylinositol 4,5-bisphosphate 5phosphatase A (INPP5J) were interlinked as the up-and downregulated DEGs in the network, respectively (Figures 4 and 5), which may suggest the effects of membrane fluidity and lipid metabolism involved in the modulation of inositol function and the influence of cell migration, adhesion, and polarity in P. olivaceus livers under cold shock [58].
Biology 2023, 12, x FOR PEER REVIEW 9 of 15 pH, rainbow trout (Oncorhynchus mykiss) cardiac actin-myosin ATPase activity showed more Ca 2+ sensitivity than that of rats, which may help O. mykiss to offset the cardioplegic effects of cold [47]. Moreover, after cold acclimation, myocardial fibrillary collagen content and/or connective tissue content increases, which may protect the myocardium from the increased hemodynamic stress due to pumping cold viscous blood [47]. There were fewer liver-specific DEGs in P. olivaceus livers than in other tissues from the network (Figure 2b and 4). For example, of the liver-specific upregulated DEGs, heme oxygenase (HMOX) catalyzes the oxidative cleavage of heme, which can protect against cold-injury-induced cell loss and damage [48], while in the downregulated DEGs, DNA polymerase (POLE) and cell division control protein (CDC45) are both required for the initiation of chromosomal DNA replication. Mitochondrial ribosome-associated GTPase 2 (GTPBP5) plays a role in the regulation of the mitochondrial ribosome assembly and of translational activity, all suggesting the inhibited function of the "cell cycle" and "cell division" [49][50][51][52]. In the livers of fish, high-temperature stress may affect apoptosis, DNA replication, protein metabolism, energy metabolism, lipid metabolism, immune, cell cycle, protein processing, and transport, as well as antioxidative responses [53][54][55], while lowtemperature stress may affect signaling transduction, membrane fluidity, lipid Moreover, several spleen-specific upregulated DEGs interlinked in the networks were membrane-and (ion)-transporter activity-related genes (Figure 4), such as solute carrier transporters (SLC6a13, SLC7A9, SLC9A5, SLC20a1a, SLC22A31, SLC35F2, SLC41a1, and SLC43A3) (Figures 4 and 5), which indicated the diverse functions of SLCs in transmembrane transport [59]. In addition, cell-cycle-and DNA-replication-related functions were inhibited in the downregulated spleen-specific DEGs. For example, DNA replication licensing factors (MCM2 and MCM3), G2/mitotic-specific cyclin-B (CCNB3), kinesin-like proteins (KIF18A, KIF20B, and KIF23), and mitotic spindle assembly checkpoint protein (MAD2A) can prevent the onset of anaphase until all chromosomes are properly aligned at the metaphase plate [60][61][62]. Proliferating cell nuclear antigen (PCNA) is also involved in the control of eukaryotic DNA replication by increasing the polymerase's processability during elongation of the leading strand [50,63] (Figures 4 and 5). There are only limited studies focused on fish spleen tissues under temperature challenges. For example, in the spleen of Nile tilapia (Oreochromis niloticus), immune and antioxidative process were affected, and xanthine oxidase, peroxidase, catalase, and superoxide dismutase played an important role [64]. Therefore, the strongly activated transport and inhibited cell cycle functions in the spleens revealed the severe effects of cold shock on P. olivaceus spleens, which may warrant further investigation of their specific functions.

Modulation of P. olivaceus CRMs with the Cold Shock Response
In addition to the DEGs functioning in coordination in the networks, many more non-DEGs, with known and novel functions, were also clustered in the CRMs (Table S4). To further investigate the tissue-specific modulated gene interaction in the given modules, the genes in each CRMs were annotated with GO (Table S4), and the functions of CRMs from overlapped tissues were characterized. For example, the turquoise module was mainly represented with gill-and liver upregulated genes and enriched functions including "membrane", "extracellular matrix", "cell junction", "response to stimulus", and "glycosylation" (Figure 6 and Table S4), indicating the activated ECM and stress gene regulation in P. olivaceus gills and livers under cold shock stress. More specifically, several claudins (CLDN1, CLDN4, CLDN6, CLDN7a, and CLDN8) and occludin (OCLN) genes were clustered as hub genes in the network (Figure 6), which may play an essential role in the formation and regulation of the tight junction through calcium-independent celladhesion activity [65]. It was reported that cold exposure increases intestinal paracellular permeability to nutrients in mice, which was associated with a transient increase in the expression of CLDN2 [66]. In addition, myosins (MYO1E, MYO5B, MYO6) are actin-based motor molecules serve in intracellular movements; epithelial splicing regulatory protein 1 (ESPR1) regulates the formation of epithelial-cell-specific isoforms; and syntaxin-19 (STX19) plays a role in endosomal trafficking of the epidermal growth factor receptor [67][68][69]. These hub genes were also correlated in the network from the cytoskeleton and anatomical structure functions (Figure 6), which, together with CLDNs and OCLNs, could contribute to the maintenance of the cellular structure of P. olivaceus gills and livers under cold stress ( Figure 1).
More interestingly, the green-yellow module, representing the downregulated DEGs from all four tissues, was enriched in functions mainly including "cell cycle", "cell division", "DNA package complex", "chromosomal region", and "response to stress" (Figure 7 and Table S4), indicating the severe cellular function inhibition in P. olivaceus tissues under cold shock stress. For example, several DNA polymerases (POLA1, POLD2/3, POLE/E2, and POLE2) and DNA replication complex GINS proteins (GINS2, GINS3), both playing essential roles in the initiation of DNA synthesis [36,50], were the top hub genes from the "cell cycle" and "DNA complex" functions ( Figures 7 and S3). In addition, G2/mitotic-specific cyclins (CCNA2 and CCNE2), interacting with cyclin-dependent kinase (CDK), are important for the control of the cell cycle at the G2/M (mitosis) transition [60]. Moreover, cell division cycle associated family proteins (CDCA5 and CDCA8) play a vital role in efficient DNA double-stranded break repair; DNA replication licensing factors (MCM2, MCM3, MCM4, and MCM5) are the replicative helicases essential for DNA replication initiation and elongation in eukaryotic cells; kinesin-like proteins (KIF4, KIF15, and KIF18) play a role in chromosome segregation during mitosis; and protein mises (MIS12, MIS18a, and MIS18BP) are required for normal chromosome alignment and segregation [38,61,[70][71][72]. All these genes were the key hub genes mainly from the "cell cycle process" and "DNA structure/polymerase" functions ( Figures 7 and S3), suggesting severe cellular function inhibition among P. olivaceus tissues under cold shock. It was reported that many types of mammalian cells failed to undergo the G(2)/M transition after cooling from 37 • C to 16-20 • C, while the progress at G(1)/S was not affected by the same temperatures [73]. Moreover, in some mammalian cultures, the exposure to 4-10 • C can cause a high abundance of mitotic synchrony (up to 80%), which may involve the cell cycle checkpoint in response to cold shock stress, inhibiting the G(1)/S transition [73]. In teleost fish under temperature changes, cell cycle progress is also regulated in P. olivaceus gills and livers under heat shock [14,53], and cold-responsive genes in zebrafish (Danio rerio) and common carp (Cyprinus carpio) under cold stress are mostly involved in oxidative phosphorylation, protein folding and degradation, RNA processing, and translation, which comprise the set of evolutionarily conserved cold-responsive mechanisms in teleosts [10]. Therefore, in this study, the CDC45-MCM-GINS (CMG) helicase, the molecular machine that unwinds template DNA during replication, could be strongly inhibited among P. olivaceus gill, heart, liver and spleen tissues during the cold shock, which warrants comprehensive investigation of their functions in further studies. More interestingly, the green-yellow module, representing the downregulated DEGs from all four tissues, was enriched in functions mainly including "cell cycle", "cell division", "DNA package complex", "chromosomal region", and "response to stress" ( Figure  7 and Table S4), indicating the severe cellular function inhibition in P. olivaceus tissues under cold shock stress. For example, several DNA polymerases (POLA1, POLD2/3, POLE/E2, and POLE2) and DNA replication complex GINS proteins (GINS2, GINS3), both playing essential roles in the initiation of DNA synthesis [36,50], were the top hub genes from the "cell cycle" and "DNA complex" functions ( Figures 7 and S3). In addition, G2/mitotic-specific cyclins (CCNA2 and CCNE2), interacting with cyclin-dependent kinase (CDK), are important for the control of the cell cycle at the G2/M (mitosis) transition [60]. Moreover, cell division cycle associated family proteins (CDCA5 and CDCA8) play a vital role in efficient DNA double-stranded break repair; DNA replication licensing factors (MCM2, MCM3, MCM4, and MCM5) are the replicative helicases essential for DNA replication initiation and elongation in eukaryotic cells; kinesin-like proteins (KIF4, KIF15, and KIF18) play a role in chromosome segregation during mitosis; and protein mises (MIS12, MIS18a, and MIS18BP) are required for normal chromosome alignment and segregation [38,61,[70][71][72]. All these genes were the key hub genes mainly from the "cell cycle process" and "DNA structure/polymerase" functions ( Figures 7 and S3), suggesting severe cellular function inhibition among P. olivaceus tissues under cold shock. It was reported that many types of mammalian cells failed to undergo the G(2)/M transition after cooling from 37 °C to 16-20 °C, while the progress at G(1)/S was not affected by the same temperatures [73]. Moreover, in some mammalian cultures, the exposure to 4-10 °C can cause a high abundance of mitotic synchrony (up to 80%), which may involve the cell cycle checkpoint in response to cold shock stress, inhibiting the G(1)/S transition [73]. In teleost fish under temperature changes, cell cycle progress is also regulated in P. olivaceus gills and livers under heat shock [14,53], and cold-responsive genes in zebrafish (Danio rerio) and common carp (Cyprinus carpio) under cold stress are mostly involved in oxidative phosphorylation, protein folding and degradation, RNA processing, and translation, which comprise the set of evolutionarily conserved cold-responsive mechanisms in teleosts [10]. Therefore, in this study, the CDC45-MCM-GINS (CMG) helicase, the molecular machine that unwinds template DNA during replication, could be strongly inhibited among P. olivaceus gill, heart, liver and spleen tissues during the cold shock, which warrants comprehensive investigation of their functions in further studies.

Conclusions
The poikilothermal teleost species present diverse temperature responses reflected in the tissue-specific cellular, antioxidative and immune statuses. This study performed a systematic transcriptomic and gene coexpression network survey in Japanese flounder (P. olivaceus) tissues under cold shock stress. The tissue-specific modulated expression regulation of genes was observed in response to acute cold shock, which was mainly found in the activated "extracellular matrix" and "cytoskeleton" functions, as well as in the inhibited "cell cycle" and "DNA complex" functions in all tested tissues. These findings indi-

Conclusions
The poikilothermal teleost species present diverse temperature responses reflected in the tissue-specific cellular, antioxidative and immune statuses. This study performed a systematic transcriptomic and gene coexpression network survey in Japanese flounder (P. olivaceus) tissues under cold shock stress. The tissue-specific modulated expression regulation of genes was observed in response to acute cold shock, which was mainly found in the activated "extracellular matrix" and "cytoskeleton" functions, as well as in the inhibited "cell cycle" and "DNA complex" functions in all tested tissues. These findings indicated that even with different levels of tissue-specific responses, cold shock may induce severely disrupted cellular functions in many P. olivaceus tissues, which reduce aquaculture productivity. Therefore, further investigation into the functions of specific CRM networks will provide more comprehensive insights into the cold adaptation mechanisms and will enable better conservation and cultivation of P. olivaceus in cold water.
Supplementary Materials: The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/biology12060784/s1. Figure S1: Volcano plots for the DEGs of the four P. olivaceus tissues under cold shock; Figure S2: Modulation from the WGCNA; Figure S3: Expression profiles of the genes enriched in the cell-cycle-related functions in the four P. olivaceus tissues under cold shock; Table S1: GO enrichment of the tissue-specific DEGs from the four P. olivaceus tissues under cold shock; Table S2: General information for the modules of WGCNA; Table S3: Identification of CRMs; Table S4: GO enrichment of each CRM.  Data Availability Statement: The P. olivaceus transcriptome datasets used in this study can be found in the NCBI Sequence Read Archive (SRA) BioProject PRJNA717098, PRJNA717103, PRJNA717106, and PRJNA717107.