Genome-Wide Identification of heat shock protein 10/60 Genes in Japanese Flounder (Paralichthys olivaceus) and Their Regulated Expression After Bacterial Infection

YAN Weijie, QIAO Yingjie, LIU Yuxiang, LIU Xiumei, ZHANG Quanqi,and WANG Xubo,*

1) Key Laboratory of Marine Genetics and Breeding, Ministry of Education, Ocean University of China, Qingdao 266003, China

2) College of Life Sciences, Yantai University, Yantai 264005, China

Abstract Heat shock proteins 10/60 (hsp10/60) are a family of conserved ubiquitously expressed heat shock proteins which are produced by cells in response to exposure to stressful conditions. Besides the chaperone and housekeeping functions, they are also known to be involved in immune response during bacterial infection. In this study, we identified and annotated 10 hsp10/60 genes through bioinformatic analysis in Japanese flounder (Paralichthys olivaceus). Among them one member of hsp10 (hspe) family and nine members of hsp60 (hspd) family were identified. Phylogenetic and selection pressure analysis showed that the hsp10/60 genes were evolutionarily constrained and their function was conserved. Besides, hsp10/60 genes were involved in different embryonic and larval stages and acted as the sentinel role in an unchallenged organism. In addition, we also observed the expression patterns of hsp10/60 genes after Edwardsiella tarda infection, for the first time in Japanese flounder. Eight out of 10 genes were differentially expressed after bacterial challenges, the significantly regulated expressions of flounder hsp10/60 genes after bacterial infections suggested their involvement in immune response in flounder. Our results provide valuable information for clarifying the evolutionary relationship, and early insights of the immune functions of hsp10/60 genes in Japanese flounder.

Key words heat shock protein; Hsp10/60; Japanese flounder; Edwardsiella tarda; immune response

1 Introduction

Heat shock proteins (Hsps), first reported inDrosophilaunder heat stress, are a suite of highly conserved proteins found in most organisms from bacteria to humans(Ritossa, 1962; Whitleyet al., 1999). Although first discovered in response to heat shock, hsps were also found to participated in a great number of stress conditions such as radiation (UV), heavy metals, pesticides, hypoxia, oxygen radicals, anti-inflammatory drugs, malignant transformation, bacterial and viral infection (Fulleret al., 1994; Sørensenet al., 2003; Gehrmannet al., 2004; Multhoff, 2006;Akiraet al., 2008). Based on the molecular weight, sequence homology as well as domain structures, hsps were classified into several subfamilies: hsp90, hsp70/hsp110,hsp10/hsp60, hsp40, and hsp20 (shsp) families (Gething,1997).

The hsp60 family is a group of proteins with distinct ringshaped, or toroid (double doughnut) quaternary structures(Quintana and Cohen, 2005) which possess a conserved Cpn60_TCP1 domain (PF00118). While hsp10 is a subfamily of 10 kDa, highly conserved, mitochondrion-resident protein, which co-chaperones with hsp60 by acting as a dome-like cover on the ATP active form of hsp60 for protein folding as well as the assembly and disassembly of protein complexes. Besides,hsp10/60play important roles in transporting the mitochondrial proteins into the mitochondrial matrix from the cytoplasm (Kollet al., 1992;Ranfordet al., 2000; Jiaet al., 2011). In addition,hsp10/60were the major hsps involved in autoimmune diseases,innate immunity and cancer (Itohet al., 2002; Ranford and Henderson, 2002; Lebretet al., 2003; Tsan and Gao,2004; Urushibaraet al., 2007). Mammalianhsp10/60have been demonstrated to act as endogenous stress signal molecules, with implicated involvement in autoimmune disease (Liang and MacRae, 1997; Vabulaset al.,2001; Tsan and Gao, 2004), which hinted thathsp10/60play central roles in mammalian defenses against pathogenic infection and responses to damage or stress in addition to normal cell functions (Vabulaset al., 2001). Moreover,hsp10/60was proposed to interact with immune cells as a ligand for a variety of cell-surface receptors such as Toll-like receptors to stimulate adaptive immune responses(Ohashiet al., 2000; Vabulaset al., 2001).

Japanese flounder is an economically important flatfish in Asian countries, especially Japan, Korea and China, due to a great number of advantages such as rapid growth rate and delicious taste (Fujiet al., 2006). Nevertheless, in recent years, increased industrial farming has resulted in the susceptibility of Japanese flounder to various pathogens including bacteria, parasites, and viruses, which in turns has led to multiple infectious diseases and serious economic losses (Isshikiet al., 2001; Moustafaet al., 2010). Edwardsiellosis is a kind of gram-negative bacterial disease induced byEdwardsiella tardathat causes considerable economic losses worldwide in a number of important economic fishes (Hoshina, 1962; Meyer and Bullock, 1973; Banget al., 1992). Japanese flounder infected withE. tardadisplay spiral movement, loss of pigmentation, cutaneous lesions, and a swollen abdomen, and this disease has become the major bottleneck in flounder fishery industry (Banget al., 1992; Moonet al., 2014). Moreover, other vertebrates including amphibians, reptiles and mammals can be infected byE. tarda. Remarkably, bacteraemia caused byE. tardacan be a fatal disease in humans (Golubet al., 2010;Hiraiet al., 2015).

Compared to plants, the systemic identification and analysis of howhspsinvolved in economic fishes’ pathogen infection is far from complete. Considering the huge economic losses incurred byE. tardain cultured fishes including Japanese flounder, a deeper understanding of the participation of stress proteins in the response toE. tardain flatfish especially Japanese flounder is necessary. Previous studies had demonstrated that the expression ofhsp10,hsp40,hsp60andhsp70genes of Japanese flounder could be influenced by the response toStreptococcus parauberisand several viruses (Donget al., 2006; Chenet al.,2010; Sung and MacRae, 2011; Chaet al., 2013; Weiet al.,2013), latest report also observed the differential expression ofhsp70genes in response toE. tardain Japanese flounder (Kimet al., 2020). However, to the best of our knowledge, comprehensive information on the roles ofhspsin the immune response in Japanese flounder is limited.Herein, we report the genome-wide identification of a full set of 10hsp10/60genes in Japanese flounder, including their sequence features, phylogeny, and selective pressures.We also determined their expression patterns in embryonic development, challenged and healthy tissues, to provide insights into the roles ofhsp10/60genes in the immune response and disease defense of flatfish.

2 Material and Methods

2.1 Identification of hsp10/60 Genes

All available non-redundanthsp10/60sequences of 7 teleosts including zebrafish (Danio rerio), medaka (Oryzias latipes), tilapia (Oreochromis niloticus),channel catfish (Ictalurus punetaus), yellow catfish (Tachysurus fulvidraco),spotted gar (Lepisosteus oculatus) andfugu (Takifugu rubripes) were downloaded from ‘Heat Shock Protein Database Information Resource’ (Sinhaet al., 2012). The latest Japanese floundergenome was obtained from NCBI(http://www.ncbi.nlm.nih.gov) under the accession number MPLB00000000.1. Further, thehsp10/60genes of zebrafish were used as queries to search against all available floundergenomic resources by TBLASTN and BLASTP,with an e-value threshod of 1e-10 to acquire the candidate genes. Then, the output putativehsp10/60sequences were submitted to Pfam (El-Gebaliet al., 2018) and SMART(Letunic and Bork, 2017) databases to confirm the conservedhsp10domain (PF00166) andhsp60domain (PF00118),and sequences without ahsp10/60domain were deleted.Allhsp10/60genes retrieved from the floundergenome were named according to ZFIN (Zebrafish Nomenclature Guidelines) and the Guidelines for the nomenclature of the human heat shock proteins (Kampingaet al., 2009). Detailed information is given in Table 1.

Table 1 Comparison of copy numbers of hsp10/60 genes among selected teleosts genomes

2.2 Phylogenetic Analysis

The full-length amino acid sequences of hsp10/60 proteins derived from zebrafish,medaka,tilapia,channel catfish,yellow catfish,spotted gar,fugu and newly identified hsp10/60 proteins of Japanese flounder were used for phylogenetic analysis. The maximum-likelihood method were used to construct a phylogenetic tree using MEGA7 (Kumaret al., 2016) software with the WAG model. Then, Evolview (Zhanget al., 2012) was used to visualise the phylogenetic tree of the divergent species above.

2.3 Gene Structure and Motif Comparison

The MEME programme (version 4.11.2, http://alternate.meme-suite.org/tools/meme) was used to identify the conserved motifs in thehsp10/60sequences, with the following parameters: any number of repetitions, maximum of 8 motifs and an optimum motifwidth of 6-200 amino acid residues. The exon-intron structures of thehsp10/60genes were determined by DNA sequences on the Gene Structure Display Server (GSDS, http://gsds.cbi.pku.edu.cn/).Then the TBtools were used to visualise the results (Chen,et al., 2018). For eachhsp10/60gene sequence, sequencespecific molecular weight, protein length, isoelectronic value (pI) and estimated variability of their respective proteins were calculated using the ProtParam Tool (Walker,2005) database.

2.4 Calculation of Selective Pressure

To estimate the selective pressure ofhsp10/60genes, the relative rates of non-synonymous substitutions (dN) and synonymous substitutions (dS) were used to represent the natural selective pressure of eight different teleosts, including Japanese flounder, zebrafish, medaka, tilapia, channel catfish, yellow catfish, spotted gar and fugu (Kryazhimskiy and Plotkin, 2008). Codon-based alignment was employed utilising ClustalW and the terminator was manually removed. To make inferences concerning selective pressure on individual codons (sites) within the coding sequence of thehsp10/60genes, the Single Likelihood Ancestor Counting (SLAC) method (Kosakovsky Pond and Frost, 2005) in Datamonkey (Delportet al., 2010) (http://www.datamonkey.org/) was used based on previous reports. SLAC can detect non-neutral evolution in relatively larger sequence alignments (＞ 50 sequences) and can process an alignment with 100 sequences and 400 codons in approximately one minute, using likelihood-based branch lengths, nucleotide and codon substitution parameters and ancestral sequence reconstructions (Pond and Frost, 2005).

2.5 Expression Analysis Using the Available RNA-Seq Datasets

2.5.1 Expression profiles of hsp10/60 genes in different tissues and developmental stages

Expression patterns ofhsp10/60genes in healthy Japanese flounder were analysed. A total of eleven tissues including heart, liver, spleen, kidney, brain, gill, muscle, intestines, stomach, testis and ovary, and six embryonic and larval stages, which were defined as stages 1-6, of healthy Japanese flounder were used to analysehsp10/60genes expression profiles.

2.5.2 Expression analysis of hsp10/60 genes in response to E. tarda infection

The RNA-seq datasets from our previous research were obtained to calculate the expression levels ofhsp10/60mRNA in Japanese flounder blood, gill, and kidney duringE. tardainfection. All fish specimens were acclimatised in aerated seawater at 19℃ for a week before bacterial infection. TheE. tardastrain EIB202 used for bacterial challenge was a highly virulent strain isolated from an outbreak in farmed turbot in Shandong province of China and was acquired from the Key Laboratory of Microbial Oceanography, Ocean University of China. In addition, we ensured that the strain was chloramphenicol, tetracycline,rifampicin, and streptomycin (Xiaoet al., 2008; Wanget al.,2009). The Japanese flounder was tested by using primers based on the specific esaV gene of pathogenicE. tardabefore the experiment and the pathogenicE. tarda, which was not detected (Tanet al., 2005). TheE. tardastrain was incubated in Luria-Bertani medium to mid-logarithmic stage,and then harvested by centrifugation, and re-suspended to a final concentration of 2×107colony forming units (CFU)mL-1in Ringer’s solution for marine teleosts.

After we initially confirmed thatE. tardainjections caused disease and even death healthy fish, we conducted a formal experiment. Healthy experimental fish were randomly divided into three groups: 60 individuals in the bacteria-challenge experiment group (BCEG), 60 individuals as Ringer’s solution control group (RSCG), and 10 individuals as a blank control group (BCG). Japanese flounder in BCEG were injected intraperitoneally with 1 mLE. tardasuspension, while fish in RSCG were injected with 1 mL Ringer’s solution. Ten fishes randomly chosen as BCG were not given any treatment before the injection experiment (0 h), and their blood, gill and kidney samples were collected from the caudal veins and immediately stored in liquid nitrogen until RNA extraction. Two individuals were randomly selected for RNA extraction from blood, gill,and kidney at each time point: BCG for 0 h (Bl-0h), BCEG for 8 h (Bl-8hE), RSCG for 8 h (Bl-8hC), BCEG for 48 h(Bl-48hE), and RSCG for 48 h (Bl-48hC). At each time point, we selected four individuals randomly for RNA extraction and equal molar ratios of RNA from two individuals were pooled as one replicate to provide templates for RNA-Seq library construction. Finally, we calculated the average value of transcriptome data in unchallenged and challenged groups respectively to better understand the dynamic process ofhsp10/60gene expression in response to bacterial infection using the normalized FPKM (Fragments Per Kilobase per Million) value (Liuet al., 2017; Liet al., 2018a, 2018b). The R package pheatmap was then applied to illustrate the expression patterns at different developmental stages (Kolde, 2018).

3 Results

3.1 Identification and Nomenclature of the hsp10/60 Genes

A total of 2 non-redundanthsp10sequences and 13hsp60sequences from zebrafish were collected from Heat Shock Protein Database Information Resource, which were used as a query to screen their orthologous genes in the Japanese flounder genome. Then, Pfam and SMART databases were used to confirm the candidate genes whose proteins contain the conservedhsp10domain (Pfam accession number PF00166) andhsp60domain (Pfam accession number PF00118). Sequences that were repetitive or lacked a typicalhsp10/60domain were removed. Finally, onehsp10sequence (hspe1-mob4) and ninehsp60genes (hspd1,tcp1,cct4,cct6a,cct7,cct8,bbs12,mkksandpikfyve) identified from the Japanese flounder genome were named according to the Zebrafish Nomenclature Guidelines. The copy numbers of eight teleosts are listed in Table 1.

3.2 Phylogenetic Analysis of the hsp10/60 Genes from Eight Teleosts

To analyse the evolutionary relationship of the eight teleosts’hsp10/60genes, an unrooted phylogenetic tree was constructed using 103 complete protein sequences from zebrafish, medaka, tilapia, channel catfish, yellow catfish,spotted gar, fugu and Japanese flounder. The evolutionary history was inferred using the Maximum Likelihood method based on the WAG model (Fig.1). The 103hsp10/60genes were divided into 14 distinct subfamilies: 10hspe1,7hspd1, 8tcp1, 6bbs10, 8bbs12, 6cct2, 6cct3, 8cct4, 6cct5, 7cct6a, 8cct7, 8cct8, 8mkksand 7pikfyve. Each of thehsp10/60genes was classified into the same cluster as the homologous genes. As is shown in Fig.1, different genes with similar protein structures were more closely related in the phylogeny tree.

3.3 Gene Structure and Motif Analysis of the hsp10/60 Genes

The exon-intron organisations and motif patterns ofhsp10/60genes were examined to investigate their similarities and divergences and gain greater insights into their evolution(Fig.2). Among thehsp10/60genes, except for one gene(bbs12)belonged to ‘one intron’ group, which possessed only one intron in the gene structure, other genes were classified into ‘multiple introns’ that possessed more than two introns andpikfyveowned the maximum number of introns. Eight conserved motifs ofhsp10/60genes were identified using the MEME website. Genescct4,cct6aandtcp1had six motifs, andmkkscontained the minimum number of motifs. Furthermore, Gene IDs, amino acid numbers, pIs and molecular weights are listed in Table 2. The lengths of hsp10/60 proteins ranged from 174 (hspe1-mob4) to 2094 (pikfyve) amino acids; molecular masses were between 19.36047 kDa (hspe1-mob4) and 235.60819 kDa (pikfyve), and the predicted pI values ranged from 5.44 (cct4)to 10.04 (hspe1-mob4).

Fig.1 The phylogenetic tree was conducted using ML (maximum-likelihood) method with WAG model. The fourteen subfamilies were distinguished in different colors, and the hsp10/60 genes of Japanese flounder were labeled using asterisk. Ol, Oryzias latipes; On, Oreochromis niloticus; Ip, Ictalurus punetaus; Tf, Tachysurus fulvidraco; Lo, Lepisosteus oculatus; Tr, Takifugu rubripes; Dr, Danio rerio.

Fig.2 Phylogenetic relationship, gene structure and conserved motif analysis of hsp10/60 genes in Japanese flounder. (a),Phylogenetic tree of 10 hsp10/60 proteins. The unrooted maximum-likelihood phylogenetic tree was constructed with MEGA7 using full-length amino acid sequences of 10 hsp10/60 proteins with WAG model. (b), Distributions of conserved motifs in hsp10/60 genes of Japanese flounder. Eight putative motifs are indicated in different colored boxes. (c),Exon/intron organization of hsp10/60 genes. Green boxes represent UTRs and yellow lines with same length represent CDSs. The upstream/downstream region of Japanese flounder hsp10/60 genes are indicated in blue boxes. The length of exons can be inferred by the scale at the bottom.

Table 2 Features of Pohsp10/60 proteins identified

3.4 Pronounced Negative Selection of the hsp10/60 Genes by dN/dS Analysis

For a better understanding of the high degree of interspecies evolutionary constraints ofhsp10/60genes, we used the synonymous rate as a benchmark to infer whether fixation of non-synonymous mutations is aided or hindered by natural selection. We calculated the dN/dS ratio using coding sequences of the 10hsp10/60genes in eight teleosts respectively (Table 3). A dN/dS ratio of less than 1 indicates negative selective pressure, also known as purifyingselection; greater than 1 indicates positive selection, selection for diversity; equal to 1 indicates neutral selection (Nei and Gojobori, 1986). In the result, the dN/dS ratio of all the 10hsp10/60genes was less than 1, indicating pronounced negative selection. The positive selection sites of all thehsp10/60genes were zero. The dN/dS ratio ofmkksandbbs12was greater than 0.2, which is higher than that of other genes, indicating that they may have a relatively high evolutionary dynamic.

Table 3 Selection pressure of hsp10/60 genes (dN/dS)in teleosts

3.5 Expression Patterns of the hsp10/60 Genes

3.5.1 Expression profiles of hsp10/60 genes in unchallenged tissues and developmental stages

The expression patterns ofhsp10/60genes at different developmental stages and under normal conditions were analysed using available RNA-seq datasets. Heat maps with the phylogenetic tree were constructed to display the expression patterns ofhsp10/60genes more comprehensively(Fig.3). The results indicated thathsp10/60genes participated in all stages of embryonic and larval development,but different genes had diverse expression patterns, there were several genes which had high expression levels such ascct4andbbs12, but several genes such asmkksandpikfyvehad relatively low expression levels during the development of embryo and larva. Besides, the expression of a specific gene changed in different stages, for instance,the expression ofmkksdecreased from stage 1 to stage 6 in general and the expression level ofpikfyveincreased from stage 1 to stage 2 but decreased from stage 3. Under unchallenged conditions, allhsp10/60genes were expressed to maintain homeostasis and the expression ofhsp10/60genes in different tissues had diverse patterns. Interestingly, although we did not find a mainhsp10/60expression site, three tissues (intestines, heart, and stomach) had the highest expression levels ofhsp10/60genes. Moreover, several genes showed tissue-specific expression profiles. For instance, there were severalhsp10/60genes highly expressed in two tissues, but poorly expressed or not expressed in other tissues, namelyhspe1-mob4in the brain and the intestines,cct4in the heart and the stomach. Regarding sexually dimorphic gene expression patterns, mosthsp10/60genes had higher expression levels in the testis than in the ovary excepthspd1andpikfyve. Furthermore,we also observed the weak expression levels ofmkksandpikfyvein most tissues above.

Fig.3 Expression profiles of hsp10/60 genes in different tissues and embryonic developmental stages. Each cell in the heat map corresponds to an expression level. The number in cell are FPKM values. Stage 1 (from two cells to morula); stage 2(from early gastrula to late somites); stage 3 (from hatching stage to 2 d after hatching); stage 4 (before metamorphosis);stage 5 (metamorphosis stages 1 to 2); stage 6 (metamorphosis stages 3 to 5).

3.5.2 Expression analysis of hsp10/60 genes in response to E. tarda infection

To explain howhsp10/60family members were regulated in response toE. tarda, we obtained the expression profiles of the 10hsp10/60genes from blood, gill, and kidney samplesfrom Japanese flounder at 8 h and 48 h afterE. tardainfection respectively (Fig.4). In summary, 9 out of a total of 10hsp10/60genes were involved in disease defence responses of Japanese flounder. AfterE. tardaadministration in the blood, expression of six genes (tcp1,cct4,cct6a,cct7,cct8andbbs12) was significantly downregulated from 8 h to 48 h, whereas that ofpikfyveexpressed higher at 8 h but dropped significantly to a low level at 48 h. Gene expression ofhspd1showed a decrease at 8 h and then rose to a high level at 48 h. Notably, of the 10 genes analysed, two (mkksandhspe1-mob4) did not exhibit differential regulation throughout, while the expression of nine genes (tcp1,cct4,cct6a,cct7,cct8,bbs12,hspe1-mob4,hspd1andpikfyve)were influenced by Ringer’s solution significantly. AfterE. tardaadministration in the gill, expression of eight genes (hspe1-mob4,tcp1,cct4,cct6a,cct7,cct8,bbs12andhspd1) were prominently upregulated from 8 h to 48 h, while the expression ofpikfyvedownregulated at 8 h and showed a significant increase at 48 h. The expression ofmkksdid not show differential expression from 8 h to 48 h. In addition, the expression of nine genes (hspe1-mob4,tcp1,cct4,cct6a,cct7,cct8,bbs12,pikfyveandhspd1)were influenced prominently by Ringer’s solution. AfterE. tardaadministration in the kidney, expression of eight genes (hspe1-mob4,tcp1,cct4,cct6a,cct7,cct8,hspd1and

bbs12) were prominently upregulated from 8 h to 48 h, and gene expression ofpikfyvedownregulated from 8 h to 48 h.Besides, the expression of eight genes (hspe1-mob4,tcp1,cct4,cct6a,cct7,cct8,bbs12andhspd1) were influenced significantly by Ringer’s solution of marine teleosts. In summary, mosthsp10/60genes in blood, gill and kidney samples were found to participate in response toE. tardainfection, except formkks.

Fig.4 Expression profiles of hsp10/60 genes under immune stress. Each column represents a time point, and each row represents a gene. The relative expression level is indicated by the color bar on the top right. G represents gill, K represents kidney and Bl represents blood. 0 h represents the blank control group at the beginning of the experiment, C 8 h and C 48 h indicates Ringer’s solution control group, whereas E 8 h and 48 h indicates bacteria-challenged experimental group,BC represents blank control.

4 Discussion

E. tardais one of the serious fish pathogens, infecting many fish species including channel catfish (Meyer, Bullock, 1973), tilapia (Van Damme, Vandepitte, 1980), common carp (Sae-Ouiet al., 1984), European sea bass (Blanchet al., 1990) and turbot. In addition, isolation ofE. tardahas been reported in invertebrates (Wyattet al., 1979), amphibians (Wyattet al., 1979), reptiles, mammals (Leottaet al., 2009) and birds (Whiteet al., 1973). Previous studies have suggested the immunity responsibilities of heat shock proteins in a number of hosts including teleosts(Srivastava, 2002; Robertset al., 2010; Xieet al., 2015;Songet al., 2016). However, despite these illustrations,deeper understanding of the molecular roles of Japanese flounder immune responses is required. In this study, we identified 14hsp10/60sequences from tilapia, 12 fromchannel catfish,15 fromzebrafish, 9 from yellow catfish,10 from Japanese flounder, 14 from medaka,13 from spotted gar and 14 from fugu, indicating a similarhspgene number among the teleosts. As for the phylogenetic analysis,each of thehsp10/60genes were divided into 14 clades(hspe1,hspd1,cct2,cct3,cct4,cct5,cct6a,cct7,cct8,tcp1,bbs10,bbs12,mkksandpikfyve) with no obvious extension between species, and Japanese flounder harbored most genes, indicating that there was a degree of evolutionary conservation. Although exhaustive searches were performed with all Japanese flounder genomic resources,hspd1,cct2,cct3andcct5were not found in the Japanese flounder genome. Combined with selection pressure results,dN/dS analysis of 10hsp10/60genes displayed a pronounced negative selection, indicating that there were no non-synonymous nucleotide changes at that codon. Previous researchers have found that new genes often evolve with rapid changes both in their structure and sequence,and mutation is the initial condition in evolution. In addition, positive Darwinian selection may be another important force driving the development of new genes. Furthermore, selection pressure may only affect specific sites rather than the whole gene, which is called site-specific selection (Koesteret al., 2012). All thehsp10/60genes were evolutionarily constrained and their functions were conserved; thus, purifying selection may account for the gene loss in the Japanese flounder genome.

Hsp10/60genes have been shown to be involved in the embryonic development of teleosts. Wanget al. (2017) observed the relatively high expression level ofhsp60genes duringblastopore closure stages to 1 dpf inSiniperca chuatsiembryos, Xuet al. (2011) reported the expression ofhsp60genes in embryonic development of grass carp. In this study, we observed that mosthsp10/60genes were expressed in different embryonic and larval stages and two of them (cct4andbbs12) had high expression levels in one or several specific developmental stages, which suggests thathsp10/60genes may be involved in the embryonic and larval development of Japanese flounder. In a healthy organism, thehsp10/60genes may act in the sentinel role.The RNA-seq datasets ofhsp10/hsp60genes in 11 tissues of Japanese flounder showed that allhsp10/60genes were expressed in different tissues during normal conditions.This anticipatory expression of thehsp10/60genes may function to prevent viral infection, deepening the notion that,hspgenes are constitutively expressed in most tissues of healthy fish (Xieet al., 2015; Songet al., 2016).In addition, the highest expression levels ofhsp10/60genes were in the intestines, heart and stomach, which have been shown to have the responsibilities of immune response (Press and Evensen, 1999; Salinas and Parra, 2015;Soonthornchaiet al., 2016). Besides, a sexually dimorphic gene expression pattern of mosthsp10/60genes showing testis significantly higher than ovary could be observed in this study. Previous studies also demonstrated thathspswere crucial for rodent spermatogenesis (Sarge and Cullen, 1997). In aquatic animals, endogenoushspsin gonad changed significantly during gametogenesis, which hinted that the high level ofhspswere related to the reproduction of Pacific oyster (Meistertzheimet al., 2009). We hypothesize that thehsp10/60genes participate in the spermatogenesis of Japanese flounder, but further studies are required.

Previous studies have shown thathspgenes participated in a great number of immune reactions of aquatic animals includingMarcobrachium rosenbergii,Sebastes schlegelii,Channa striatus,Misgurnus anguillicaudatus,Larimichthys crocea,Labeo rohita,Pinctada martensii,Nibea albifloraandProcambarus clarkiiin response to bacterial and heavy metal’s challenges (Chaurasia, 2016; Giriet al.,2016; Heet al., 2016; Kim and Kang, 2016; Liet al.,2016; Sathyamoorthyet al., 2017; Chenget al., 2019; Daiet al., 2020; Xianget al., 2020; Yanet al., 2020). There were limited researches revealed the characters ofhspgenes in the immune response of Japanese flounder, detaillyhsp10/60genes againstStreptococcus parauberisinfection (Chaet al., 2013) andhsp70genes againstE. tardainjection (Kimet al., 2020). Besides, compared to other subfamilies ofhspsuperfamily, the reports abouthsp10/60genes’ immune responsibilities in other aquatic animals were also insufficient (Zhouet al., 2010; Xuet al., 2011;Xieet al., 2015; Wanget al., 2017). Interestingly, we also found that the expression of manyhsp10/60genes, 9 in total, was up- or downregulated after bacterial infection.Although the mechanisms of regulation are not known at present,hsp10/60genes were significantly regulated after bacterial challenges. This suggests that they were involved in disease responses, and that they may participated in disease defences against infectious bacteria. We found a dynamic response pattern in mosthsp10/60genes, namely different expression profiles at different time points, similar response patterns were reported inhsp70/hsp90genes in the reactions towardsAeromonashydrophila,Pseudomonas plecoglossicidaand poly I:C (Minget al., 2010;Xianget al., 2020). Specifically, different expression profiles were found in three tissues afterE. tardaadministration: nine genes were found up- or downregulated in gill;and 8 genes were found regulated in blood and 9 in kidney. Surprisingly, among these 9 regulatedhsp10/60genes,we found a tissue-specific response profiles; in other words,regulatedhsp10/60genes had diverse response patterns in different tissues. In addition, there were still eighthsp10/60genes (hspe1-mob4,tcp1,cct4,cct6a,cct7,cct8,hspd1andbbs12) that had the same response patterns in two tissues, namely in gill and kidney. Furthermore, a distinct spatial-specific expression profile of the tissue expression ofhsp10/60genes was observed:pikfyvehad relatively low expression level in all the tissues but responded strongly toE. tardainfection in kidney, gill, and blood. We assumed that there were two kinds ofhsp10/60 genes in Japanese flounder, similar to thehsp40genes (Ojimaet al.,2005; Donget al., 2006):pikfyve, which belong to the inducible form showing strongly stress-induced expression,whereas others belonged to the constitutive form that seemed weakly or not inducible. The constitutive expression ofhsp10/60genes was possibly gene-specific for different functions and requires further study. In addition, we cannot ignore the influence of Ringer’s solution. In brief, 9 out of 10hsp10/60genes exceptmkkswere influenced significantly by Ringer’s solution by contrast between BC(blank control) and C (control), and the influence ofE.tardachallenge were evaluated by contrast between C (control) and E (experiment), so we confirm that although Ringer’s solution influenced the expression ofhsp10/60genes prominently,hsp10/60genes responded toE. tardainfection in Japanese flounder actively.

Overall, 8hsp10/60genes were likely to be activated in response toE. tardainfections. However, the regulatory mechanism underlying such a wide spectrum of transcriptional activation has not yet been fully elucidated. Further studies are warranted to explore the mechanisms of regulation and to understand the roles ofhsp10/60genes in host defences against infectious diseases of flatfish.

5 Conclusions

Hsp10/60sis a major low-molecular-weight subfamily ofhspsuperfamily which involves in many physiological processes. In this study, 10hsp10/60genes from Japanese flounder were systematically identified and characterised for the first time. The hsp10/60proteinswere divided into 14 distinct subfamilies after phylogenetic analysis of eight teleosts, and dN/dS analysis showed that all thehsp10/60genes had pronounced negative selection. In addition, we observed the expression patterns of thehsp10/60genes in challenged and unchallenged organisms. In unchallenged organisms,hsp10/60genes are likely to expressed as the sentinel role. In challenged organisms, 8 out of 10hsp10/60genes responded to infection withE. tarda, indicating that they may be involved as a part of the disease response genes while several of these genes may be a part of disease defence against the bacterial pathogens. Our results provide useful information for further investigation of the role ofhsp10/60genes in the immune response in Japanese flounder.

Acknowledgements

This work was supported by the National Key Research and Development Program of China (No. 2018YFD0900 601) and the Natural Science Foundation of Shandong Province (No. ZR2017MC072).