Noninvasive Epidermal Mucus RNA-Seq Analysis for Developing Methods to Evaluate Environmental Stress in the Japanese Eel (Anguilla japonica)

TL;DR

This study explores noninvasive methods to assess environmental stress in Japanese eels by analyzing RNA from their epidermal mucus.

Contribution

The study introduces a noninvasive RNA-Seq approach using epidermal mucus to evaluate environmental stress in wild Japanese eels.

Findings

SS stress elevated genes related to cellular immunity and inflammation.

Low temperature increased lipid metabolism genes, while high temperature increased heat-shock response genes.

Epidermal mucus RNA-Seq can identify stressor-specific gene responses in Japanese eels.

Abstract

Environmental changes in rivers caused by fluctuations in water temperature and increasingly intense rainfall linked to climate change are likely to influence the physiological condition (environmental stress) of wild Japanese eels (Anguilla japonica). This study examined the epidermal mucus of A. japonica, which can be collected noninvasively in the field, and evaluated whether RNA-Seq analysis could characterize the physiological state of this species in response to such environmental stressors. We conducted controlled tank trials simulating (1) suspended solids (SS) stress resulting from sediment runoff during heavy rainfall (using kaolin clay addition) and (2) water temperature variations (10 °C, 20 °C, and 30 °C) as physical stresses experienced by wild A. japonica under climate-driven changes. Subsequently, we performed RNA-Seq analysis to determine genes in the epidermal mucus responsive to each experimental stressor. As a result, under (1) SS stress, genes associated primarily with cellular immunity and inflammatory pathways showed elevated expression. In the (2) water temperature trials, relative to the control (20 °C), the low-temperature group (10 °C) displayed increased expression of lipid-metabolism-related genes such as scd, prom2, and tmem64, whereas the high-temperature group (30 °C) exhibited heightened expression of genes linked to the heat-shock response. Thus, because this study enabled the identification of genes specifically expressed in reaction to each stressor from the epidermal mucus of A. japonica, these findings are expected to support the future development of noninvasive approaches for evaluating the physiological state (health condition) of wild A. japonica.

The online version contains supplementary material available at 10.1007/s10126-026-10582-w.

Figures

Principal component analysis (PCA) of RNA-seq data from epidermal mucus and four other tissues (fin, gill, muscle, and skin) under normal conditions

Estimation of tissue-specific gene expression in A. japonica epidermal mucus and other tissues under normal conditions. (a) Number of genes specifically expressed in each tissue based on the tissue specificity index τ (Yanai et al. 2005). (b) Heatmap of the top 20 enriched terms across tissue-specific expressed genes in epidermal mucus and three other tissues (skin, fin, and gills) in direct contact with the external environment. Heatmap colors represent p-values, with darker colors indicating lower p-values; gray indicates that the term was not enriched in that tissue

Number of differentially expressed genes (DEGs) in response to (a) suspended sediment (SS) stress, (b) low-temperature stress, and (c) high-temperature stress (n = 5; FDR p-value <0.05 with a twofold cutoff)

Venn diagram illustrating the numbers of upregulated genes in the (a) suspended sediment (SS) stress group and (b) low-temperature (LT) and (c) high-temperature (HT) stress groups (n = 5; false discovery rate (FDR) p < 0.05 with a twofold expression cutoff)

Heatmap of the top 20 terms identified by enrichment analysis of upregulated genes in the suspended sediment (SS) stress group. Darker colors indicate lower p-values, whereas gray indicates that the term was not enriched in the corresponding tissue

RT-qPCR of the highly expressed genes in the suspended sediment (SS) stress group (n = 3). Relative gene expression was calculated by 2−ΔΔCT method using the β-actin gene (actb) as an internal control. The vertical axis represents the log 10 scale. Asterisks (*) indicate a statistically significant difference compared with the control group (Student’s t-test, p < 0.05)

Heatmap of the top 20 terms identified by enrichment analysis of upregulated genes in the (a) low-temperature (LT) and (b) high-temperature (HT) stress groups. Darker colors indicate lower p-values, whereas gray indicates that the term was not enriched in the corresponding tissue

RT-qPCR of the highly expressed genes in the low temperature (LT) stress and high temperature (HT) stress groups (n = 3). Relative gene expression was calculated by 2−ΔΔCT method using the β-actin gene (actb) as an internal control. The vertical axis represents the log 10 scale. Asterisks (*) indicate a statistically significant difference compared with the control group (Student’s t-test, p < 0.05)

Tables

Table 1. Top 10 genes upregulated in epidermal mucus by suspended sediment stress (ranked by FDR p-value)

Gene ID	Putative gene product	Gene name	FDR p-value	Fold change
Aja4_Angja6g001740	Tubulin alpha-1 chain	tuba.1	1.2E-10	31.0
Aja4_Angja1g003940	Galactose-binding lectin l-1	gbl	1.5E-07	82.6
Aja4_Angja80548s000010	Mast cell protease 1 A	mcpt1a	1.5E-07	26.6
Aja4_Angja8g008180	Interleukin-1 family member A	il1fa	9.8E-07	4.6
Aja4_Angja7g017710	Tubulin beta chain	tubb	2.1E-06	17.2
Aja4_Angja7g016300	Neutrophil cytosol factor 2	ncf2	2.1E-06	7.3
Aja4_Angja7976s000010	GTPase IMAP family member 9	gimap9	2.3E-06	10.9
Aja4_Angja80754s000010	Mast cell protease 4	mcpt4	2.3E-06	9.1
Aja4_Angja6g001760	Tubulin alpha-1 A chain	tuba.2	5.8E-06	121.1
Aja4_Angja12g010690	Tubulin alpha chain	tuba.3	1.8E-05	12.4

Table 2. Top 10 genes upregulated in epidermal mucus by low-temperature stress (ranked by FDR p-value)

Gene ID	Putative gene product	Gene name	FDR p-value	Fold change
Aja4_Angja1g016530	Stearoyl-CoA desaturase	scd	5.8E-20	11.6
Aja4_Angja16g008220	Prominin-2	prom2	1.6E-18	53.0
Aja4_Angja4g006630	Cartilage matrix protein	cmp	4.0E-16	28.5
Aja4_Angja19g000940	Microtubule-associated protein, RP/EB family, member 3b	mapre3b	3.8E-14	8.5
Aja4_Angja4g010010	protein Jumonji	jarid2b	6.5E-14	6.7
Aja4_Angja4g009780	Transmembrane protein 64	tmem64	1.8E-13	5.8
Aja4_Angja18g010160	Mucin-2	muc2	6.5E-13	49.2
Aja4_Angja15g001060	High mobility group protein B1	hmgb1	2.2E-12	6.1
Aja4_Angja2g000310	RNA-binding protein with multiple splicing	rbpms	9.0E-12	23.8
Aja4_Angja1g004880	Immunoglobulin mu heavy chain	ighm	1.5E-11	32.1

Table 3. Top 10 genes upregulated in epidermal mucus by high-temperature stress (ranked by FDR p-value)

Gene ID	Putative gene product	Gene name	FDR p-value	Fold change
Aja4_Angja15g004970	Keratin, type I cytoskeletal 13	krt13	4.4E-17	36.6
Aja4_Angja7g002830	Thioredoxin-interacting protein	txnip	4.9E-14	6.8
Aja4_Angja19g008080	RING finger protein 208	rnf208	1.1E-09	28.0
Aja4_Angja6g008190	Keratin, type II cytoskeletal 8	krt8	2.0E-09	31.6
Aja4_Angja10g002520	C-terminal-binding protein 2	ctbp2	4.5E-08	5.1
Aja4_Angja13g001750	Zona pellucida-like domain-containing protein 1	zpld1	5.5E-08	12391.4
Aja4_Angja6g011050	Period circadian protein homolog 3	per3	6.7E-07	4.6
Aja4_Angja12g007580	Cadherin-like protein 26	cdh26	1.3E-06	16.4
Aja4_Angja6g011060	Period circadian protein homolog 3	per3	1.7E-06	4.0
Aja4_Angja8g012040	SLIT and NTRK-like protein 6	slitrk6	1.7E-06	4.7

Full text

Introduction

The Japanese eel (Anguilla japonica) is a catadromous fish that spawns west of the Mariana Ridge, where the hatched leaf-like larvae (leptocephali) are transported by the North Equatorial and Kuroshio currents (Tsukamoto 1992, 2006; Tsukamoto et al. 2011). After traveling thousands of kilometers to the East Asian coast, they metamorphose into transparent juvenile A. japonica (glass eels) and migrate into estuaries and rivers for further growth (Fukuda and Shinoda 2023). Analyses of strontium (Sr) and calcium (Ca) concentrations in otoliths have revealed that some A. japonica populations do not enter freshwater and instead spend their entire lives in marine habitats (Tsukamoto et al. 1998), while others repeatedly shift between oceanic and freshwater environments throughout their growth phases (Tsukamoto and Arai 2001). A. japonica is an economically important species for inland water fisheries in Japan and a key resource for aquaculture industries across East Asian countries (Sakurai and Shibusawa 2021). However, natural populations of A. japonica in Japanese rivers have declined sharply in recent decades, and the species is now listed on the Red List of the Ministry of the Environment of Japan (Japanese Ministry of the Environment 2020) as well as the International Union for Conservation of Nature Red List of Threatened Species (Pike et al. 2020).

It has been noted that factors driving the decline in the population of A. japonica include overfishing, degradation of freshwater habitats, and alterations in the marine environment (Hakoyama et al. 2023; Tsukamoto et al. 2009). To counter this downward trend, management actions such as stocking and habitat restoration have been carried out in rivers across Japan (Kaifu et al. 2014; Tsukamoto et al. 2009). However, there is concern that warming and increasingly extreme rainfall associated with future climate change will alter river environments and impose physiological stress on freshwater fish (Ficke et al. 2007; Nimma et al. 2025). When exposed to such environmental shifts, fish can adjust physiologically to some extent, but once a threshold is exceeded, mortality increases, or individuals migrate to more suitable habitats. Therefore, monitoring the physiological condition of this species represents an important approach for advancing effective conservation management.

In fish conservation, it has been proposed that transcriptional response profiles provide an effective means of assessing the intensity of environmental stress and the extent of physiological impact (i.e., health status) (Connon et al. 2018). Typically, tissues or organs are used for transcriptome analysis, but because this requires sacrificing specimens, noninvasive sampling approaches are preferable for conservation purposes. Fish epidermal mucus functions as both a physical and chemical barrier between fish and their environment, contributes to homeostasis, and is therefore gaining attention as a medium for noninvasive monitoring of stress and health (Brinchmann 2016; Dash et al. 2018; Ekman et al. 2023). Gene expression profiles derived from epidermal mucus have already been reported as an effective noninvasive tool for examining immune responses to bacterial infection (Ren et al. 2015) and for evaluating the effects of diluted bitumen exposure (Andrzejczyk et al. 2022).

In this study, we first compared the transcriptome (RNA-Seq) of A. japonica epidermal mucus under normal conditions with that of tissues directly exposed to the external environment—the skin, muscle, caudal fin, and gills—to examine the characteristics of genes expressed in epidermal mucus. We then conducted laboratory trials simulating two physical stressors associated with climate-driven environmental change in wild A. japonica: (1) increased “turbidity” from sediment runoff following heavy rainfall and (2) changes in water temperature. For the turbidity trial, kaolin, a natural clay mineral commonly used to investigate the effects of suspended solids (SS) from sediment loading on fish, was applied (Awata et al. 2011; Redding et al. 1987). Using these samples, we analyzed RNA-Seq data from epidermal mucus, skin, gills, and liver to evaluate whether epidermal mucus is a suitable medium for noninvasive assessment of the physiological condition of A. japonica under each environmental stressor.

Materials and Methods

Ethical Statement

All experimental fish were handled in accordance with the Guidelines for Animal Experimentation of the Japan Fisheries Research and Education Agency (FRA).

Sample Collection Under Normal Conditions

Six A. japonica individuals weighing approximately 190 g were obtained from a commercial eel supplier in Shizuoka Prefecture, Japan. The specimens obtained from this facility that did not undergo stress trials or other tank trials was considered “normal condition.” Epidermal mucus samples were collected and preserved using DNA/RNA Shield Collection Tubes with Swabs (ZYMO RESEARCH, USA). Eels secrete large amounts of epidermal mucus, making it possible to collect samples without significantly damaging the epidermis. Thus, eels were captured with a net, and the epidermal mucus from the central part of their body was immediately collected using a swab, taking care to avoid damaging the epidermal surface.

Additionally, skin, muscle, caudal fin, and gill tissues were dissected and preserved in DNA/RNA Shield (ZYMO RESEARCH, USA). These samples were incubated overnight at 4 °C and subsequently stored at − 20 °C until use.

Sample Collection of Stress Conditions

Cultured A. japonica of approximately 300 g body weight were purchased from a commercial eel supplier in Kagoshima Prefecture, Japan. About 50 individuals were maintained in a fiberglass-reinforced plastic tank containing six polyvinyl chloride (PVC) pipe shelters (Φ100 mm × 1 m) for 8 months under freshwater, closed-circulation conditions at 16℃–17℃ and natural sunlight at the Nikko Field Station, FRA.

Suspended Sediment (SS) Stress

The experimental fish were divided into two groups: an SS-stress group exposed to kaolin clay and a control group with no additions. Two acrylic tanks (W 89 cm × D 29 cm × H 30 cm) were prepared, each equipped with a 300-W heater and five PVC shelters (Φ50 mm × 80 cm), and five A. japonica were held in aerated freshwater (water depth 21 cm). The water temperature was gradually increased to 20 °C over a period of 1 day, and the fish were acclimated for 4 days. In European eel (Anguilla anguilla), it has been reported that although plasma cortisol levels return to baseline within 6 h following transport stress, energy metabolism remains elevated for at least 72 h (Boerrigter et al. 2015). Therefore, the acclimation phase of this trial was set to more than 72 h. For the SS-stress group, kaolin clay (Fujifilm Wako Pure Chemicals, Japan) was added every six hours to maintain turbidity between 50 and 100 nephelometric turbidity units until sampling 1 day later. Epidermal mucus, skin, gills, and liver were collected and stored as described above (Sample Collection under normal conditions).

Temperature Stress

The experimental fish were divided into three groups: 10 °C (low temperature, LT), 20 °C (control), and 30 °C (high temperature, HT). Three acrylic tanks (W 88 cm × D 44 cm × H 44.5 cm) were prepared, each filled with aerated freshwater and five PVC shelters (Φ50 mm × 70 cm), and covered with a lid containing eight holes (Φ30 mm) to allow natural light to enter. The 10 °C group was supplied with spring water, whereas the 20 °C group was heated using a 300-W heater and the 30 °C group with a 1,000-W heater. The experiments were conducted under natural daylight conditions. Five A. japonica were placed in each tank and gradually acclimated to their respective temperature settings over 2 days, after which the temperatures were maintained until sampling 5 days later. Epidermal mucus, skin, gills, and liver were collected and stored as described above (Sample Collection under normal conditions).

RNA-Seq Analysis

Total RNA was extracted using the Quick-RNA Microprep kit (ZYMO RESEARCH, USA) for epidermal mucus and the Direct-zol RNA MicroPrep kit with TRI-Reagent (ZYMO RESEARCH, USA) for the other tissues. Library preparation for RNA sequencing (RNA-Seq) was carried out using the Stranded mRNA Prep, Ligation Kit (Illumina, USA) following the manufacturer’s protocol. Library quality and concentration were assessed using the TapeStation 4150 (Agilent Technologies, USA) and a Qubit Fluorometer (Thermo Fisher Scientific, USA), respectively. Sequencing was performed on a NextSeq 500 platform (Illumina, USA) using a paired-end (2 × 74 bp) strategy. RNA-Seq data were analyzed with CLC Genomics Workbench 22.0.2 (CLC Bio, Denmark) as follows: low-quality sequences were removed using default settings (quality limit = 0.05, ambiguous limit = 2), and filtered reads were mapped to the A. japonica reference genome (Sudo et al. 2024) using default parameters (mismatch cost = 2, insertion cost = 3, deletion cost = 3, length fraction = 0.5, similarity fraction = 0.8). Transcripts per million (TPM) values were calculated for each gene (Wagner et al. 2012). To estimate tissue-specific gene expression in A. japonica epidermal mucus and the four other tissues under normal conditions, the tissue-specificity index τ (Yanai et al. 2005) was applied as described previously (Nakamura et al. 2024). Principal component analysis (PCA) was conducted using iDEP 2.01 (Ge et al. 2018) to classify gene expression patterns across samples. Differentially expressed genes (DEGs) were identified using a ≥ 2-fold expression change and a false discovery rate (FDR) p-value < 0.05. Venn diagrams were generated with Venny 2.1 (Oliveros 2007–2015) to identify DEGs unique to or shared among the experiments. To incorporate zebrafish (Danio rerio) annotation data, a well-annotated model species, homologous genes between zebrafish (genome assembly GRCz11: GCA_000002035.4) and A. japonica were identified via BLASTP (Altschul et al. 1990) using an E-value threshold of 1.0E − 4. Functional enrichment analysis was then performed in Metascape (Zhou et al. 2019) using the zebrafish gene list.

RT-qPCR Assay

To validate the RNA-Seq results of epidermal mucus, we performed RT-qPCR (n = 3) for several selected genes, including eight genes highly expressed in the SS-stress group (tuba.2, tubb, gbl, mcpt1a, il1fa, ncf2, gimap9, and mcpt4), four genes highly expressed in the LT-stress group (scd, prom2, rbpms, and muc2), and four genes highly expressed in the HT-stress group (krt8, krt13, txnipa, and hps90). cDNA was synthesized from 5 ng of total RNA using the PrimeScript™ II 1st Strand cDNA Synthesis Kit (Takara, Japan). PCR primers are listed in Supplementary Table S1. PCR amplification was performed in a 20-µL reaction mixture containing 1 µL of cDNA, 1× PowerTrack™ SYBR Green Master Mix (Thermo Fisher Scientific, USA), and 0.4 µM of each primer. qPCR was conducted using a QuantStudio 3 Real-Time PCR System (Thermo Fisher Scientific, USA) under the following cycling conditions: initial denaturation at 95 °C for 1 min, followed by 40 cycles of 95 °C for 5 s and 60 °C for 20 s. All samples, along with a no-template control, were analyzed in duplicate. The relative expression of target genes was calculated using the 2^−ΔΔCT^ method (Livak and Schmittgen 2001) and normalized to the expression of the β-actin gene (actb). qPCR data were analyzed using Student’s t-test, with statistical significance set at p < 0.05.

Results

RNA-Seq Analysis of Epidermal Mucus Under Normal Condition

RNA-Seq of samples collected under normal conditions (epidermal mucus, skin, muscle, caudal fin, and gills) generated an average of 16.4 million clean paired-end reads per sample, with approximately 90.2% mapping to the reference sequence (Supplementary Table S2).

PCA of tissues under normal conditions showed that gene expression patterns separated into two major clusters: muscle and the tissues in direct contact with the external environment (epidermal mucus, skin, fins, and gills) (Fig. 1). Within the latter group, epidermal mucus displayed a clearly distinct expression pattern from the other three tissues (skin, fins, and gills) (Fig. 1). In contrast, fin and gill expression patterns were highly similar and overlapped substantially (Fig. 1).

Principal component analysis (PCA) of RNA-seq data from epidermal mucus and four other tissues (fin, gill, muscle, and skin) under normal conditions

In this study, the normal condition experiment used approximately 190 g of A. japonica from a farm in Shizuoka Prefecture, whereas the stress trials used approximately 300 g of A. japonica from a farm in Kagoshima Prefecture. Consequently, we investigated whether differences in body weight or rearing location influenced gene expression patterns in the epidermal mucus, skin, and gills common to both trials under the normal condition and the A. japonica specimens used as controls in stress trials. PCA revealed clear separation of gene expression patterns by tissue type, with epidermal mucus displaying distinct differences from other tissues (Supplementary Fig. S1). Therefore, it was suggested that gene expression does not vary significantly across tissues based on body weight or rearing location.

Analysis of tissue-specific expressed genes identified 102 genes specifically expressed in epidermal mucus. Muscle showed the highest number of tissue-specific genes (871), followed by gills (489), fins (156), and skin (118) (Fig. 2a). Enrichment analysis of these tissue-specific gene sets revealed that the extremely large number of muscle-specific genes generated a predominance of muscle-related terms, making it difficult to characterize the other tissues (Supplementary Fig. S2). Consistent with the PCA results demonstrating a marked separation between muscle and the other tissues (Fig. 1), we performed enrichment analysis focusing on epidermal mucus and the three other tissues in direct contact with the environment (skin, fins, and gills). This analysis highlighted “ribosome” and “adipocytokine signaling pathway” as enriched terms associated with epidermal mucus–specific expressed genes (Fig. 2b). Note that the epidermal mucus-specific genes included in the “adipocytokine signaling pathway” were pck1, nfkbie, socs3b, irf6, and tradd (Supplementary Table. S3). The top 10 epidermal mucus-specific expressed genes based on TPM values are presented in Supplementary Table S4. Of these 10 genes, eight were reported to play roles in biological defense, including five genes involved in inflammatory responses (s100p, mfap4, trim39, epx and mep1b), two antimicrobial peptides (lgals1 and lgals4), and one gene from the perforin superfamily of pore-forming immune effector molecules (sntxa).

Estimation of tissue-specific gene expression in A. japonica epidermal mucus and other tissues under normal conditions. (a) Number of genes specifically expressed in each tissue based on the tissue specificity index τ (Yanai et al. 2005). (b) Heatmap of the top 20 enriched terms across tissue-specific expressed genes in epidermal mucus and three other tissues (skin, fin, and gills) in direct contact with the external environment. Heatmap colors represent p-values, with darker colors indicating lower p-values; gray indicates that the term was not enriched in that tissue

Overview of RNA-Seq Analysis of Stress Samples

RNA-Seq of SS and temperature stress samples (epidermal mucus, skin, gills, and liver) generated an average of 21.5 million clean paired-end reads per sample, with approximately 88.8% mapping to the reference sequence (Supplementary Table S5).

Figure 3 summarizes the number of DEGs detected across the three stress experiments (SS, LT, and HT). In the SS experiment, epidermal mucus showed the highest number of DEGs (674), followed by liver (314), gills (125), and skin (116) (Fig. 3a). In the LT stress experiment, skin exhibited the largest number of DEGs (3,909), followed by liver (3,708), epidermal mucus (2,890), and gills (2,006) (Fig. 3b). In the HT stress experiment, gills had the greatest number of DEGs (1,554), with liver (870), skin (800), and epidermal mucus (515) showing fewer (Fig. 3c). These findings suggest that epidermal mucus exhibited the strongest physiological response to SS stress. Furthermore, the LT experiment generated the highest number of DEGs across all tissues, indicating that low temperature induced the greatest overall physiological changes in A. japonica.

Number of differentially expressed genes (DEGs) in response to (a) suspended sediment (SS) stress, (b) low-temperature stress, and (c) high-temperature stress (n = 5; FDR p-value <0.05 with a twofold cutoff)

We then examined whether the genes upregulated by each stress condition overlapped among the four tissues (Fig. 4). Under SS stress, only 4 genes (0.8%) were commonly upregulated across all tissues, and epidermal mucus exhibited the largest number of tissue-specific upregulated genes, with 242 genes (50.6%) (Fig. 4a). Under LT stress, 92 genes (2.5%) were shared among the four tissues, while the liver showed the highest number of tissue-specific induced genes with 925 genes (24.6%), followed by epidermal mucus with 876 genes (23.3%) (Fig. 4b). In the HT stress experiment, 16 genes (0.9%) overlapped among the four tissues, with the gills displaying the greatest number of tissue-specific upregulated genes (689 genes, 39.6%) and epidermal mucus the fewest (131 genes, 7.5%) (Fig. 4c). Overall, these findings indicate limited overlap in upregulated genes among tissues under each stress condition, suggesting that distinct physiological responses occur in epidermal mucus and other tissues.

Venn diagram illustrating the numbers of upregulated genes in the (a) suspended sediment (SS) stress group and (b) low-temperature (LT) and (c) high-temperature (HT) stress groups (n = 5; false discovery rate (FDR) p < 0.05 with a twofold expression cutoff)

RNA-Seq Analysis of Epidermal Mucus Under SS Stress

The top 20 terms identified by the Metascape enrichment analysis of upregulated genes in the SS stress group are shown in Fig. 5. Numerous terms were significantly enriched in epidermal mucus compared with the other three tissues, and most were associated with biological defense processes (GO:0006952: defense response; GO:0009617: response to bacterium; R-DRE-877312: regulation of IFNG signaling). In particular, many terms were related to cellular immunity, including phagocytosis (dre04145: phagosome), leukocyte migration and regulation (GO:0071674: mononuclear cell migration; GO:0002685: regulation of leukocyte migration; R-DRE-114608: platelet degranulation), and inflammatory responses (R-DRE-6798695: neutrophil degranulation; dre03320: PPAR signaling pathway). Among these, “R-DRE-6798695: neutrophil degranulation” was also enriched in gills and skin, while “dre03320: PPAR signaling pathway” was enriched in liver. The top 10 genes upregulated in epidermal mucus under SS stress are listed in Table 1. These genes could be grouped into two functional categories: (1) tubulin genes (three tuba and one tubb) and (2) immune-related genes (gbl, mcpt1a, mcpt4, il1fa, ncf2, and gimap9). Note that the three tuba genes are located at distinct loci on the reference genome. Among these 10 genes, RT-qPCR was performed on eight genes (excluding two of the three tuba genes), confirming that their expression was elevated in the SS-stress group (Fig. 6).

Heatmap of the top 20 terms identified by enrichment analysis of upregulated genes in the suspended sediment (SS) stress group. Darker colors indicate lower p-values, whereas gray indicates that the term was not enriched in the corresponding tissue

RT-qPCR of the highly expressed genes in the suspended sediment (SS) stress group (n = 3). Relative gene expression was calculated by 2−ΔΔCT method using the β-actin gene (actb) as an internal control. The vertical axis represents the log 10 scale. Asterisks () indicate a statistically significant difference compared with the control group (Student’s t-test, p < 0.05)*

RNA-Seq Analysis of Epidermal Mucus Under Temperature Stress

The LT stress enrichment analysis showed that immune-related terms (R-DRE-1222556: ROS and RNS production in phagocytes; R-DRE-168256: immune system), terms associated with protein synthesis (GO:0042274: ribosomal small subunit biogenesis; GO:0022613: ribonucleoprotein complex biogenesis; GO:0042273: ribosomal large subunit biogenesis), and protein metabolic processes (GO:0006457: protein folding) were enriched across all tissues, including epidermal mucus (Fig. 7a). In addition, the immune-related terms “GO:0007229: integrin-mediated signaling pathway” and “GO:0050853: B cell receptor signaling pathway” were enriched exclusively in the epidermal mucus dataset (Fig. 7a). In contrast, the stress-related term “GO:0033554: cellular response to stress” was enriched only in the liver dataset (Fig. 7a). Among the top 10 genes upregulated in epidermal mucus under LT stress, the identified genes were related to lipid metabolism (scd, prom2, tmem64), immune response (muc2, hmgb1, ighm), collagen secretion and fibrillogenesis (cmp), histone demethylase activity (jarid2b), regulation of pre-mRNA alternative splicing (rbpms), and regulation of dynamics of the microtubule cytoskeleton (mapre3b; Table 2). Among these genes, scd, prom2, rbpms, and muc2 were also confirmed to be upregulated in the LT-stress group by RT-qPCR (Fig. 8).

Heatmap of the top 20 terms identified by enrichment analysis of upregulated genes in the (a) low-temperature (LT) and (b) high-temperature (HT) stress groups. Darker colors indicate lower p-values, whereas gray indicates that the term was not enriched in the corresponding tissue

RT-qPCR of the highly expressed genes in the low temperature (LT) stress and high temperature (HT) stress groups (n = 3). Relative gene expression was calculated by 2−ΔΔCT method using the β-actin gene (actb) as an internal control. The vertical axis represents the log 10 scale. Asterisks () indicate a statistically significant difference compared with the control group (Student’s t-test, p < 0.05)*

In the HT stress experiment, enrichment analysis across all four tissues showed that key heat shock–related terms—“GO:0009266: response to temperature stimulus,” “GO:0042026: protein refolding,” and “dre04141: protein processing in endoplasmic reticulum”—were enriched (Fig. 7b). Additional stress–associated terms, including “GO:0070482: response to oxygen levels,” “GO:0009628: response to abiotic stimulus,” and “R-DRE-2262752: cellular responses to stress,” were enriched in three tissues but not in epidermal mucus (Fig. 7b). In epidermal mucus, muscle, and liver, the enriched terms included “R-DRE-168249: innate immune system” and “WP1331: adipogenesis” (Fig. 7b). The only term enriched exclusively in epidermal mucus was “GO:0031076: embryonic camera-type eye development.” The genes included in “embryonic camera-type eye development” are listed in Supplementary Table S6. There were 19 genes involved, and eight genes were found only in the epidermal mucus group under HT stress. These 19 genes included two heat shock protein genes (hsp8 and hsp70), suggesting that genes in this GO category might also be related to heat shock. In contrast, the gills showed a particularly large number of enriched terms, especially those associated with the cell cycle, such as “GO:0051726: regulation of cell cycle,” “GO:0007052: mitotic spindle organization,” “dre04110: cell cycle,” “GO:1903046: meiotic cell cycle process,” and “GO:0000278: mitotic cell cycle” (Fig. 7b).

Among the top 10 upregulated genes in epidermal mucus under HT stress, we identified genes associated with oxidative stress (txnip, ctbp2), regulation of inflammation (cdh26), two keratin genes (krt8, krt13), two circadian clock–related genes (per3 × 2), two sensory-function genes (zpld1, slitrk6), and one gene involved in protein ubiquitination (rnf208) (Table 3). Among these genes, krt8, krt13, and txnipa were also confirmed to be upregulated in the HT-stress group by RT-qPCR (Fig. 8). Additionally, hsp90, commonly used as an indicator of HT stress, was also upregulated in the epidermal mucus of the HT-stress group according to both RNA-Seq (fold change = 9.2, FDR p-value = 2.46E − 5) and RT-qPCR (Fig. 8).

Discussion

In this study, we focused on epidermal mucus as a noninvasive sampling medium and examined the potential of RNA-Seq to assess the physiological state of A. japonica in response to turbidity (SS stress) and temperature stresses. We first performed RNA-Seq analyses of epidermal mucus together with reference tissues—skin, muscle, fins, and gills—under normal conditions. The results showed that the gene expression profile of epidermal mucus was distinct from that of the other tissues, including the skin (Fig. 1). Enrichment analysis of the 102 epidermal mucus–specific genes highlighted “ribosome” and “adipocytokine signaling pathway” as enriched terms (Fig. 2b). Although ribosomes are organelles responsible for protein synthesis, many ribosomal proteins are present in fish epidermal mucus and have documented antimicrobial activity (Bergsson et al. 2005; Dawar et al. 2024a). This suggests that A. japonica may also express ribosomal proteins in their epidermal mucus as part of a defense mechanism against pathogenic microorganisms. The “adipocytokine signaling pathway”; however, has been proposed to contribute to inflammatory responses in grass carp (Ctenopharyngodon idella) (Sun et al. 2019). Furthermore, genes involved in inflammatory responses, antimicrobial peptides, and the perforin superfamily of pore-forming immune effectors were highly expressed in epidermal mucus (Supplementary Table S4). Fish epidermal mucus is widely recognized as playing a critical role in innate immunity against pathogens (Guardiola et al. 2014; Ren et al. 2015; Subramanian et al. 2007). Taken together, these findings imply that the epidermal mucus of A. japonica functions as an important first-line barrier against pathogenic microorganisms.

In fact, RNA-Seq analysis in this study showed that epidermal mucus exhibited the highest number of DEGs among all analyzed tissues under SS stress, and enrichment analysis revealed that many terms associated with biological defense were included (Fig. 5). Moreover, the top 10 upregulated genes in epidermal mucus consisted of immune-related genes and tubulin genes (Table 1). Tubulin is known to play an important role in the phagocytosis process during bacterial infection (Duan et al. 2020; Zhang et al. 2019), and proteome analysis of Atlantic cod (Gadus morhua) epidermal mucus has similarly reported increased beta 2-tubulin expression following Vibrio anguillarum infection (Rajan et al. 2013). These findings suggest that kaolin was recognized as a foreign substance in A. japonica epidermal mucus, triggering biological defense mechanisms, and further support that A. japonica epidermal mucus functions as a primary barrier against pathogenic microorganisms. Conversely, mechanical abrasion of the epidermis by kaolin might also directly induce an immune response. In grass carp, mechanical abrasion was found to increase the activity of immune-related enzymes in the epidermal mucus (Wang et al. 2024). Further histological studies are required to confirm whether kaolin causes mechanical damage to the epidermal surface of A. japonica. In summary, it is inferred that turbidity (SS) stress caused by sediment discharge in rivers activates biological defense responses in A. japonica, likely resulting in additional energy expenditure.

The LT experiment generated the highest number of DEGs across all tissues (Fig. 3). Histopathological analysis of A. japonica found no significant changes in individuals held at warm temperatures (25 °C), whereas exposure to “cold water stress” at < 15 °C caused damage to A. japonica (Kobayashi et al. 2000). Therefore, A. japonica is more sensitive to cold temperatures than to warm temperatures, suggesting that stronger physiological responses are induced and more DEGs are detected. RNA-Seq analysis of LT stress revealed activation of immune-related terms and those associated with ribosomal biogenesis and protein folding in all analyzed tissues, including epidermal mucus (Fig. 7a). Interactome analysis of gilthead sea bream (Sparus aurata) epidermal mucus has likewise shown LT stress–induced upregulation of proteins related to ribosomal activity and protein folding (Sanahuja et al. 2019). As proposed by Sanahuja et al. (2019), these responses may function to block pathogen invasion and regulate epidermal cell turnover under low-temperature conditions. Among the top 10 genes upregulated under LT stress in epidermal mucus, three were involved in lipid metabolism (scd, prom2, tmem64) (Table 2). Of these, scd is well documented to be induced by cold stress in fish (Tiku et al. 1996; Xu et al. 2015; Zerai et al. 2010). Scd encodes stearoyl-CoA desaturase, the rate-limiting enzyme in monounsaturated fatty acid synthesis, and its activation under cold stress increases membrane fluidity, thereby enhancing cold adaptation capacity (Hsieh and Kuo 2005; Xu et al. 2015). The top 10 also included an RNA-binding protein gene (rbpms) and a mucin gene (muc2), which encodes a major glycoprotein component of epidermal mucus (Table 2). In human epithelial cells, cold stimulation has been shown to elevate RNA-binding protein expression, which in turn promotes increased mucin expression (Ran et al. 2016). Because epidermal mucus secretion is known to increase in fish under cold stress (Refaey et al. 2022), the increased expression of rbpms and muc2 in A. japonica epidermal mucus suggests that they may contribute to enhanced mucus secretion during LT stress. Collectively, these results indicate that RNA-Seq analysis of A. japonica epidermal mucus can effectively characterize physiological responses to LT stress.

In HT stress, the number of DEGs was highest in the gills (Fig. 3). Fish gills represent the largest internal organ directly interacting with the external environment. HT exposure accelerates respiration rates to increase oxygen uptake and support basal metabolic needs, and this process has been reported to induce changes in gill phenotype (Mohamad et al. 2021). Thus, HT exposure might have caused sufficient changes to affect the gill phenotype, resulting in the observation of dynamic gene expression changes in the gills. However, enrichment analysis of epidermal mucus and all analyzed tissues showed that heat shock–related terms were enriched (Fig. 7b). Furthermore, hsp90, a commonly used indicator of HT stress (Akbarzadeh et al. 2018), was upregulated in epidermal mucilage under HT stress as determined by both RNA-Seq and RT-qPCR (Fig. 8). These findings indicate that the epidermal mucus RNA-Seq data successfully captured the physiological state of A. japonica under HT stress. Among the top 10 genes upregulated under HT stress, txnip and ctbp2, both linked to oxidative stress (Choi and Park 2023; Kainoh et al. 2021), were identified (Table 3). Because heat stress has been shown to trigger oxidative stress in fish (Banh et al. 2016; Hao et al. 2024; Singh et al. 2025), elevated txnip and ctbp2 expression indicates that HT stress–associated oxidative stress arose in A. japonica epidermal mucus under study. Increased expression of per3 has been reported under heat stress in chickens (Jastrebski et al. 2017) and rats (Kotlarz et al. 2022), and heat stress is known to influence circadian rhythms across diverse animals and plants (Laosuntisuk and Doherty 2022; Long et al. 2012; Okamoto-Mizuno and Mizuno 2012; Shehab-El-Deen et al. 2010; Wu et al. 2021). Therefore, altered per3 expression suggests that HT stress may also affect the circadian rhythm of A. japonica. Two keratin genes (ktr8 and krt13) also showed increased expression under HT stress (Table 3). These keratin genes exhibited increased expression under HT stress, whereas their expression tended to be decreased under LT stress in RT-qPCR assays, suggesting their potential as indicators of thermal stress (Fig. 8). Keratin proteins are abundant components of fish epidermal mucus (Rajan et al. 2011; Sanahuja and Ibarz 2015), and their expression increases during bacterial (Dawar et al. 2024b) or parasitic (Easy and Ross 2009) infection. Keratins possess antimicrobial properties due to their pore-forming capabilities (Molle et al. 2008). In grass carp, keratin gene expression increases at high water temperatures when fish are exposed to formalin-killed bacteria (Yang et al. 2016), and elevated keratin 8 (krt8) protein levels have been reported in the epidermal mucus of gilthead sea bream under chronic stress (Pérez-Sánchez et al. 2017). In mammals, keratin 8 is thought to act as a phosphate “sponge” for stress-activated kinases, thereby mitigating harmful effects and protecting cells from injury (Ku and Omary 2006). Thus, the increased expression of keratin genes in A. japonica epidermal mucus under HT stress may represent a protective response aimed at reducing HT-induced epidermal cell damage. However, given the limited information regarding keratin involvement in HT stress in fish, further research will be necessary to clarify its functional role.

In conclusion, this study demonstrated that RNA-Seq analysis of A. japonica epidermal mucus enables noninvasive assessment of physiological responses to environmental stresses, including SS and temperature stress. Although the results of this study were based on a controlled single-replication tank experiment, RNA-Seq of A. japonica epidermal mucus enables the acquisition of gene expression profiles for stress responses. In the future, it will be necessary to collect epidermal mucus samples from wild A. japonica and verify the applicability of this method. Thus, this study is expected to contribute to the future development of methods for non-invasively monitoring the physiological condition (health status) of farmed and wild A. japonica.

Supplementary Information

Below is the link to the electronic supplementary material.

Supplementary File 1 (XLSX 10.4 KB)

Supplementary File 2 (XLSX 10.9 KB)

Supplementary File 3 (XLSX 10.1 KB)

Supplementary File 4 (XLSX 13.0 KB)

Supplementary File 5 (XLSX 15.2 KB)

Supplementary File 6 (XLSX 12.3 KB)

Supplementary File 7 (PDF 212 KB)

Supplementary File 8 (PDF 284 KB)