ZHANG Xuekai ,FAN Chao ,ZHANG Xingzhi ,LI Qiongzhen ,LI Yangchun ,MA Peizhen ,YUE Shaobo,HUO Zhongming,and WANG Zhaoping,

1) Fisheries College, Ocean University of China,Qingdao 266003, China

2) Guangxi Academy of Fisheries Sciences, Guangxi Key Laboratory of Aquatic Genetic Breeding and Healthy Aquaculture,Nanning 530021,China

Abstract Crossbreeding is an effective approach to manage the genetic decline in aquaculture.One-way hybrids of Crassostrea sikamea (♀) and Crassostrea gigas (♂) have advantages in growth traits and adaptation to high temperature.Here,we used highthroughput sequencing to analyze the molecular processes in the hybrids under and after thermal stress.The hybrids were cultured in the seawater with an increasing temperature from 25℃ to 40℃ during 10 hours,which is regarded as the thermal stress stage.Then the temperature decreased from 40℃ to 25℃ within 2 h,which is regarded as the recovery stage.In this study,1293 significant differentially expressed genes (DEGs) were obtained under thermal stress,of which 576 were upregulated and 717 were downregulated,and 740 significant differentially expressed genes (DEGs) were obtained in the recovery stage,of which the number of upregulated and downregulated genes was 409 and 331,respectively.The antigen processing and presentation,NOD-like,and NF-kappa B pathways were significantly enriched during the thermal stress stage.The MAPK and PPAR signaling pathways were significantly enriched during the recovery stage.The HSP70,HSP90,and CANX genes were strongly and rapidly upregulated in the control/thermal stress groups but were slightly less upregulated in the thermal stress/recovery group.These results indicate that the innate immune system or nonspecific immunity was deployed to protect interior tissues from thermal stress.In addition,85% of the mutual DEGs were involved in bidirectional regulation (up/down or down/up) when the oysters were removed from the thermal stress to recover.This study provides preliminary insight into the molecular response of C.sikamea (♀) and C.gigas (♂) hybrids to thermal stress and provides a basis for future studies on temperature-adaptation and the possible expansion of hybrid breeding.

Key words oyster hybrids;gill;thermal stress;recovery;transcriptome

1 Introduction

Oysters are representative bivalve mollusks that are widely distributed in the oceans around the world (Guoet al.,2015) and are also important global commercial aquaculture species.Questions have arisen about summer mortality,inbreeding depression,parasites,and disease in the oyster industry (Proestouet al.,2016).Hybridization can be an effective breeding method to resolve these problems(Gjedrem and Baranski,2010).Studies on oyster interspecific crossbreeding have been well documented.In recent years,Crassostrea gigas×C.nippona(Xuet al.,2019),C.sikamea×C.angulata(Yanet al.,2017,2018),C.gigas×C.hongkongensis(Zhanget al.,2017),andC.hongkongensis×C.rivularis(Huoet al.,2014) hybrids have shown different trophic values,and different adaptabilities to salinity,temperature,or growth traits.However,these studies have mostly focused on hereditary capacity without documenting expression at the transcriptional level under a specific stress.

The natural distribution of oysters is heavily affected by temperature.Oysters change body temperature by up to 10– 20℃ within a few hours during the diurnal/tidal cycles and even more dramatic changes (from 0 to 35– 40℃) occur over longer (seasonal) time-spans (Ivaninaet al.,2009).These changes increase with the tidal height from 37℃ to 43℃ (Hamdounet al.,2003).The temperature adaptability ofC.gigashas been well documented.Their optimum temperature range is 11– 34℃,but they can thrive in a wide range of water temperatures from 5℃ to 35℃ (Mannet al.,1994).The highest reported surviving temperature is 43℃(Shamseldinet al.,1997).However,wildC.gigashave never been found on the southern coast of China,as they are usually considered a coldwater species.Due to the influence of the warm Pacific current,C.sikameaare widely distributed in Kumamoto,Japan,but are seldom found in the Yellow Sea of China at the same latitude (Elstonet al.,2012).The earliest artificial crossing ofC.gigasandC.sikameawas reported in 1994 (Bankset al.,1994).The heterosis of survival and maternal comparable growth traits have been reported inC.gigasandC.sikameahybrids (Honget al.,2012).According to our previous study,the hybrids have lower mortality than both parents,particularly the sire-sibs,during culture in southern China,where the annual water temperature is high (Yueet al.,2021).Therefore,we have determined that it is possible to selectively breed culturedC.sikameaandC.gigashybrids in the southern China coastal area.

The temperature of seawater changes naturally over time.The body temperature of oysters is greatly affected by the environment.Previous studies have focused on the effect of a particular temperature on gene transcription and expression in oysters,but the survival of oysters under high temperature is limited by both the thermal degree and duration (Limet al.,2016),while the effect of gradual thermal stimulation and recovery processes after thermal stimulation have been rarely discussed (Zhanget al.,2012).Sustained molecular regulation is needed to repair metabolic damage after thermal stress in oysters.Although closing the shell may help the oyster withstand transient adverse environmental changes,such as salinity stress,it contributes little to protect the oysters from temperature change(Zhuet al.,2016).The gill is an early sensor of physical environmental changes,and thus,it has been assumed to be sensitive to thermal stress.Oysters trigger the expressions of stress response-related genes to adapt to harsh environmental conditions (Guoet al.,2015).Of the various organs,the gill plays a key role in oyster innate immunity.

The goal of this study was to explore the molecular mechanism of temperature tolerance in oyster hybrids.Highthroughput RNA-seq was used to analyze the hybrids under and after thermal stress.The results will provide a solid foundation for expanding the culture of a new hybrid species.

2 Materials and Methods

2.1 Hybrid Collection,Thermal Challenge,and Tissue Sampling

Hybrid collectionThe one-way hybrid ofC.sikamea(♀) ×C.gigas(♂) were accomplished in April 2018 (Hedgecocket al.,1999).After one year of culture in Beihai,China (109.4?N,21.5?W),we obtained thousands of hybrids with a shell length of 32.69 mm ± 3.84 mm and weight of 11.50 g ± 3.51 g (mean ± SD).From these hybrids,100 healthy non-injured individuals were screened and used in the experiment.The hybrids were reared in filtered,aerated seawater with temperature of 25℃± 0.6℃,dissolved oxygen of (6.12 ± 0.11) mg L-1,salinity of 30.6 ± 0.5,and pH 8.21 ±0.09 for two weeks of acclimatization and fed a mixture ofIsochrysis galbanaandChaetoceros muelleri.

Thermal stress and recoveryThe temperature of the seawater was increased from 25℃ to 40℃ within 10 h (heating efficiency 1.5℃ h-1) (Hamdounet al.,2003),then decreased to 25℃ over 2 h.The heating process was based on the daily surface temperature change that occurs during oyster culture.Thus,the oysters were deployed from 5 a.m.to 5 p.m.in an experimental day.

Tissue samplingAs oyster gills are a perfect thermal stress material (Meistertzheimet al.,2007;Langet al.,2009),gills of nine oysters in each group were randomly selected as biological replicates (three biological replicates× three oysters per biological replicate) for transcriptome analysis.Three independent thermal challenge experiments were conducted,and gills were collected at the initial 25℃(control group,CG),40℃ (thermal stress group,TSG),and the 25℃ of recovery (recovery group,RG),respectively.In addition,other 30 hybrid siblings were deployed and monitored for another 24 h to ensure that the thermal stress was nonlethal.

2.2 RNA Extraction,Library Construction,and Sequencing

Total RNA was extracted from the gills of the control,thermal stress,and recovery groups using TRIzol (Invitrogen,Carlsbad,CA,USA) according to the manufacturer’s instructions.RNA samples were qualified and quantified using a NanoDrop spectrophotometer (Nanodrop Technologies,Wilmington,DE,USA) and the Agilent Technologies 2100 Bioanalyzer (Thermo Fisher Scientific,Waltham,MA,USA),respectively.

The mRNA was purified by removing any potential genetic DNA pollution with DNase I,and removing any possible rRNA pollution with RNase H,and finally recovered with Oligo(dT)-attached magnetic beads (Illumina,San Diego,CA,USA).Purified mRNA was fragmented with a fragment buffer at an appropriate temperature.Then,first/second-strand cDNA was generated by polymerase chain reaction (PCR).The A-Tailing Mix and RNA index adapters were added and incubated to carry out end-repair.The cDNA fragments with adapters were amplified by PCR,and the products were purified with Ampure XP beads.The library was validated on the Agilent Technologies 2100 Bioanalyzer for quality control.The final library was amplified with phi29 to create DNA nanoballs (DNB),containing more than 300 copies of the original template.The DNBs were loaded into the patterned nanoarray and pair-end 100 base reads were generated on a BGISEQ500 platform (BGI,Shenzhen,China).

2.3 Transcriptome Assembly,Prediction of the Coding Sequences (CDSs),and Identification of Simple Sequence Repeats (SSRs)

The filtering software SOAPnuke v1.4.0 developed by BGI was used for statistic analysis.Trimmomatic v0.36 software was used for filtering.Raw reads with joint contamination and high numbers of unknown N bases were removed before the data analysis to ensure the reliability of the results.Then,the Q20/Q30 (the percentage of bases with a Phred value over 20/30) and GC-content of the clean data were measured.Clean reads were assembled with Trinity v2.0.6 software (Grabherret al.,2011),and the transcripts were clustered with Tgicl v2.1 to generate unigenes.The quality of the assembled transcripts was assessed using Benchmarking Universal Single-Copy Orthologs (BUSCO),a single-copy lineal homology database,and the integrity of the transcriptome assembly was demonstrated by comparison with conserved genes.TransDecoder v3.0.1 software was used to identify candidate coding regions in the unigenes by extracting the longest open reading frame,and then using Blast to compare the Swiss-Prot database and Hmmscan to search Pfam protein homologous sequences to predict the coding regions (Kimet al.,2015).According to the assembly results,we also tested the SSRs of the unigenes and designed primers using MISA v1.0 and Primer3 v2.2.2 software (Untergasseret al.,2012).

2.4 Functional Annotation of the Assembled Unigenes

The fragments per kb per million reads (FPKM) method was used to calculate the expression levels.Bowtie2 v2.2.5 software was used to compare clean reads with the reference gene sequence,and then RSEM was used to calculate the expression levels of the genes and transcripts.The prediction and annotation data of all unigenes were compared with protein databases,including NCBI nonredundant protein sequences and nucleotide sequences (Nr&Nt) (ftp://ftp.ncbi.nlm.nih.gov/blast/db),Swiss-Prot (http://www.ebi.ac.uk/uniprot),Gene Ontology (GO) (http://www.geneontology.org/),Clusters of Orthologous Groups of proteins/eukaryotic Orthologous Groups (COG/KOG) (https://www.ncbi.nlm.nih.gov/COG/),Kyoto Encyclopedia of Genes and Genomes (KEGG) (http://www.genome.jp/kegg/),and Protein family (Pfam) (http://pfam.xfam.org).Hmmscan v3.0,Blast v2.2.23,and Blast2GO software were used to obtain the annotation information for the unigenes.Priority order of alignments from the Nr,Nt,KEGG,Swiss-Prot,GO,COG,and Pfam databases were followed when the results of different databases conflicted.

2.5 Analysis of Differentially Expressed Genes(DEGs) and Functional Enrichment

The FPKM method was used to calculate the expression of unigenes in this study.Calculated gene expression can be directly used to compare gene expression levels between samples.The false discovery rate (FDR) was applied to identify the threshold of the P-value in multiple tests (Mortazaviet al.,2008).We set our threshold for significantly expressed genes at FPKM ≥ 1,|log2(Fold change)| ≥ 2,FDR<0.05.GO and KEGG functional pathway enrichment analyses were performed to identify the main biochemical and signal transduction pathways of the DEGs using FDR <0.05 as a cut-off value for significantly enriched DEGs.

2.6 Quantitative Real-Time Polymerase Chain Reaction (qRT-PCR)

Total RNA was reverse-transcribed into cDNA with the PrimeScript RT Reagent Kit (Takara,Beijing,China).Eight randomly selected DEGs,oligoadenylate synthase (OAS),aquaporin 8,interferon regulatory factor 1,meprin A subunit beta,heat shock protein 90a,tumor necrosis factor(TNF) superfamily member 10 (TNFSF10),TNFsuperfamily member 14 (LIGHT/ TNFSF14),DnaJ heat shock protein family memberB5,and a reference gene (tubulin α)were selected for verification by qRT-PCR.Gene-specific primers were designed with Primer Premier 5.0 software.All real-time RT-PCR experiments were carried out in 96-well PCR plates.Each 20 μL amplification reaction was conducted using the StepOne Plus system (ABI,Foster City,CA,USA) in triplicate.The PCR process included 95℃ for 3 min;45 cycles of 95℃ for 5 s,and 60℃ for 30 s;followed by a melting curve.The differences in gene expression were analyzed with the 2-ΔΔCTmethod (Livak and Schmittgen,2001).The expression fold change (2-ΔΔCT)of the qRT-PCR data and log2(Fold change) of the RNAseq data were used in the correlation analysis.

3 Results

3.1 RNA Sequencing and de novo Assembly

The average number of total raw reads generated from each library was 67.40 million.After filtering the low quality reads,the average clean reads accounted for 64.54 million,and the average Q30 value was 89.33%.The total mapped rate of the library was 83.67%– 86.92%.After removing the redundant data,98209 unigenes were assembled into nine libraries with an average unigene number of 53197,length of 947 bp,GC percentage of 39.76%,and N50 value of 1640 bp (Table 1).An evaluation of the transcript assembly showed that more than 80% of the unigenes were matched or partially matched to the BUSCO database and unmatched genes only accounted for 3% (Fig.1).A total of 52778 CDSs were detected in our study.In total,13776 SSRs were detected in 10663 unigenes,and 7826 unigenes encoding transcription factors were predicted (Table 2).

Fig.1 Evaluation of transcript assembly by Benchmarking Universal Single-Copy Orthologs (BUSCO).C,unigenes matched the sequence in the BUSCO database;F,only partial sequences were compared to the BUSCO database;D,multiple genes were aligned to the same BUSCO;M,unigenes were filtered out of the sequence.

Table 1 Sequencing data statistics for the control,treated and restored samples

Table 2 Quantity and quality indicators of the coding sequences (CDSs),and simple sequence repeats (SSRs)

These results show that the data were of high quality and that the unigenes could be used for the annotation analysis.All raw reads were deposited in the NCBI Short Read Archive database with Accession number PRJNA 559659.

3.2 Functional Annotation and Classification

We identified unigenes in seven functional databases.Of these unigenes,67108 in the NR protein sequence database (68.33%),79570 in the NT database (81.02%),39880 in the Swiss-Prot (40.61%),47470 in the KEGG (48.34%),36573 in the KOG (37.24%),46663 in the Pfam (47.51%),and 35500 in the GO (36.15%) were annotated.All unigenes were functionally annotated,with 18825 (19.17%)matched to all databases and 84874 (86.42%) matched to at least one database (Table 3).

Table 3 Annotation rates of the unigenes in the seven databases

There were five species in the NR annotations,includingC.gigas,C.virginica,Mizuhopecten yessoensis,Lottia gigantea,andExaiptasia pallida,with annotation rates of 89.77% (60238),6.36% (4267),0.69% (462),0.13% (88),and 0.12% (83),respectively.

The classification of 25 KOG gene homologs were identified,and the top three annotated terms were signal transduction (6698),general function prediction (6148),and posttranslational modifications,protein turnover,and chaperones (4087) (Fig.2).

Fig.2 Unigene sequence annotations submitted to the Eukaryotic Orthologous Groups (KOG) database to obtain classification information on the gene homologs.

3.3 Identification and Enrichment Analysis of the DEGs

The DEG screening results among the three stages of thermal stress (CG,TSG,and RG) are shown below.The GO classification was used to assess DEG function.The GO analysis assigns DEGs to three major functional categories:molecular functions (MF),cellular components (CC),and biological processes (BP).In the two comparisons of CG/TSG and TSG/RG,314 and 451 genes were grouped into MF,850 and 516 genes were grouped into the immune system,and 601 and 307 were related to BP,respectively.The top two annotated terms in each functional category were binding (313/158) and catalytic activity (182/89),membrane (205/149) and membrane part (199/147),cellular process (184/79),and metabolic process (113/46)(Fig.3).

Fig.3 Gene Ontology analysis of the differentially expressed genes in the control vs.thermal and thermal vs.recovery group.

In an analysis of gene expression among different groups,1293 DEGs and 740 DEGs were identified in the CG/TSG and TSG/RG groups,respectively.Of these genes,717 DEGs were downregulated and 576 were upregulated in the CG/TSG group,while 319 and 409 DEGs were downregulated and upregulated in the TSG/RG group.In addition,120 mutual DEGs were in the two comparisons.The mutual genes exhibited four modes of regulation:upregulation,up/downregulation,down/upregulation,and downregulation (Fig.4).

Fig.4 Differentially expressed genes (DEGs) between paired comparisons of the control vs.thermal groups (CG/TSG) and thermal vs.recovery groups (TSG/RG).A,venn diagram of the DEGs for comparing the CG/TSG and TSG/RG groups;B,four-quadrant scatterplot of the mutual DEGs in both comparisons;C,volcano plot of the DEGs for the CG/TSG comparison;D,volcano plot of the DEGs for the TSG/RG comparison.All points in the four-quadrant scatterplot represent expression log2 (Fold change) between the two comparisons;all points in the volcano plots represent the difference in expression log2 (Fold change) plotted against the level of statistical significance.

The DEGs in the CG/TSG and TSG/RG groups were analyzed with the KEGG to investigate the predicted functional information of DEGs in the hybrids under and after thermal stress.Our KEGG analysis mainly focused on environmental information processing,organismal systems,cellular processes,genetic information processing,and metabolism in KEGG classification level 1.The analysis of the KEGG pathways provides an overview of the DEGs that are involved in the thermal stress and recovery phases.

The DEGs in the CG/TSG group were mapped to the KEGG database and assigned to 13 KEGG pathways (FDR<0.05),among which 120,101,78,43,41,37,and 30 genes were grouped into signal transduction;folding,sorting and degradation;immune system;endocrine system;cell growth and death;aging;and transcription,respectively,in KEGG classification level 2.The main enriched pathways included the ko04064_NF-kappa B signaling pathway (e.g.,DDX58,BCL2L1,andVCAM1),ko04390/ko04391/ko043 92_hippo signaling pathway (e.g.,MOB1,ED,andFAT4),ko04141_protein processing in the endoplasmic reticulum(e.g.,DNAJC3,EIF2AK4,andHSP90B),ko03060_protein export (e.g.,HSPA5 andRGLG),ko04621_NOD-like receptor signaling pathway (e.g.,TRP,OAS,HSP90,andTBK1),ko04612_antigen processing and presentation (e.g.,HSP70/90 andCTSL),ko04918_thyroid hormone synthesis (e.g.,HSPA5,HSP90B,andCANX),ko04217_necroptosis (e.g.,CYLD,PKR,andPARP),and ko04213_ longevity regulating pathway (e.g.,CRYABandHSPA1s).Fourteen pathways in the TSG/RG group were significantly enriched (FDR <0.05),and 63,53,48,38,27,26,and 24 genes were grouped into the immune system;endocrine system;folding,sorting,and degradation;cellular community-eukaryotes;signal transduction;cell growth and death;and aging,respectively.The immune system,folding,sorting,and degradation system,and the endocrine system were the top three enriched level 2 KEGG,containing four,four,and two pathways,respectively.Five KEGG pathways were highly enriched,including ko04666_Fc gamma R-mediated phagocytosis (e.g.,RHOQ,vasodilator-stimulated phosphoprotein(VASP),andSLC30A10),ko04670_Leukocyte transendothelial migration (e.g.,ITGA8,MYL12,andVASP),ko04010_MAPK signaling pathway (e.g.,FCN,GRB2,andFGFR3),ko03320_PPAR signaling pathway(e.g.,PLIN2,ACDC,UBC,andPCK),ko04915_Estrogen signaling pathway (e.g.,GRB2,CREB3,and mitogen-activated protein kinase,MAPK1/3).In addition,six pathways,including ko04213,ko04915,ko04918,ko04141,ko03060,and ko04612 were identified.Among these pathways,ko-04141_protein processing in the endoplasmic reticulum was the most enriched,as 89 and 39 genes were enriched,respectively (Fig.5).

Fig.5 Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment of the DEGs between the control (CG),thermal(TSG),and recovery (RG) groups.The shape of the marker represents a different group;the color of the marker represents the richness factor (the ratio of the number of genes annotated to an item in the selected gene set to the number of genes annotated to the item by the species);the square area represents the number of DEGs enriched in each pathway.

The GO enrichment of the DEGs in the CG/TSG and TSG/RG groups was analyzed to clarify the changes in the hybrids from the molecular to the organismal level under thermal stress and during the recovery phase.The GO terms were grouped into three categories of BP,MF,and CC.Overall,20 and 13 GO terms were significantly enriched in both comparisons (FDR <0.05),and the distributions were distinct in the thermal stress and recovery phases.The order of the GO term enriching number in the CG/TSG group was MF (ten terms) >BP (eight terms) >CC (two terms),while the order changed to BP (eight terms) >MF (four terms)>CC (one term) in the TSGvs.RG group.GO:0005524_ATP binding,GO:0051082_unfolded protein binding,GO:0006457_protein folding,and GO:0003700_DNA-binding transcription factor activity were the top four enriched GO terms in the CG/TSG group,accounting for 77,29,28,and 21 DEGs.Similarly,GO:0051082 (7 DEGs) and GO:0006457 (8 DEGs) terms were also highly enriched in the TSG/RG group,and the other two highly enriched terms in the top four were GO:0005576_extracellular region (20 DEGs) and GO:0008191_metalloendopeptidase inhibitor activity (8 DEGs) (Fig.6).The lists of these representative DEGs and their functional results are shown in Supplementary Information.

Fig.6 The Gene Ontology enrichment of the differentially expressed genes (DEGs) between the control groups (CG),thermal groups (TSG),and recovery groups (RG).The shape of the marker represents a different group;the color of the marker represents the richness factor (the ratio of the number of genes annotated to an item in the selected gene set to the number of genes annotated to the item by the species);the square area represents the number of DEGs enriched in each pathway.

3.4 Validation of Transcriptome Data by qRT-PCR

The qRT-PCR gene expression results are presented in Fig.7.The qRT-PCR expression pattern was consistent with the RNA-seq results,which not only validated the expression of the genes but also verified the reliability and accuracy of our transcriptome analysis.

Fig.7 Comparison between qRT-PCR and the RNA-seq data set;significant difference between a treated sample and a blank sample is indicated by an asterisk.log2 (Fold change) of RNA-seq,FPKM of TSG/FPKM of CG(t/c) and FPKM of RG/FPKM of TSG(r/t);log2 (Fold change) of qRT-PCR,ratio of relative gene expression of TSG/CG and RG/TSG.

4 Discussion

The threat of summer heat toC.gigashas been of great concern to the aquaculture industry.In contrast,C.sikameahas a higher temperature tolerance because it is commonly distributed along the southern China coast.In our previous study,hybrids ofC.sikamea(♀)×C.gigas(♂) suffered lower mortality during the second farming summer in a south China oyster farm (Yueet al.,2021).Thus,we performed a thermal challenge transcriptome experiment to assess their responses to thermal stress and recovery.In our study,553 DEGs were involved in the regulation of thermal stress than the recovery stage.According to GO analysis,the main annotated terms were similar between the CG/TSG and TSG/RG groups,indicating that similar regulatory processes may be occurring in the oyster during and after thermal stress.A transcriptome study ofC.gigasshowed that the CC terms of GO increase over time and become the most frequent GO characteristic after 24 h (Yanget al.,2017).Similarly,many DEGs were grouped into CC terms in the thermal stress group in this study.In this study,the proportion of up regulated and downregulated DEGs changed from 576:717 to 409:331 when oysters were removed from the thermal stress and placed in recovery process.This dramatic change indicates that although similar pathways or processes were deployed (according to the GO analysis),specific genes may have been regulated in the opposite direction.This conclusion was confirmed by the analysis of mutual genes;77 and 25 DEGs had bidirectional regulation (up/down and down/up),accounting for 85% of the total mutual DEGs.

Oysters rely heavily on the innate immune system for protection from the ambient environment (Guoet al.,2015).An inflammatory immune response requires the recruitment of leukocytes to the inflammation site upon foreign insult.Chemokines are small chemoattractant peptides that provide directional cues for cell trafficking and thus are vital for a protective host response (Wong and Fish,2003).In the present study,the immune system and signal transduction signaling pathways were both highly enriched with DEGs in the two comparisons,indicating that the regulation of innate immunity relies heavily on these two pathways.This pattern of innate immunity regulation under thermal stress has been reported in mollusks by other studies,which supports our findings (Chapmanet al.,2011).The expression of different anti-stress proteins appears to be important to ease cellular stress,inflammation,and stimulation of immune function components during the response to thermal stress (Meistertzheimet al.,2007).The transporter associated with antigen processing is essential for the adaptive immune system (Abele and Tampe,2004).In our study,the antigen processing and presentation pathways were significantly enriched in the two comparisons.Calnexin (CANX) is an endoplasmic reticulum chaperone system gene that ensures proper folding and quality control of newly synthesized glycoproteins (Hamdounet al.,2003),and plays a crucial role in recognizing antigens.In this study,it was strongly upregulated in the CG/TSG but slightly less in the TSG/RG group.Previous studies have identified the function of theNOD-like receptor signaling pathway,which is involved in detecting various pathogens and generating the innate immune responses (Vandenabeele and Bertrand,2012).An important part of the invertebrate immune system participates in the immune regulation ofC.sikameaandC.agulatahybrids under low temperature and salinity stress (Yanet al.,2018),which further supports our findings.Genes,such as B-cell linker protein,VASP,and actin beta/gamma 1 (ACTB_G1),were strongly up/downregulated,and are related to the inflammatory response (FcepsilonRIsignaling pathway),phagocytosis (Fcgamma Rmediated phagocytosis),and cellular migration (leukocyte transendothelial migration).The differences in molecular expression between the thermal stress and recovery stages indicate that various immune responses were deployed after thermal stress.

Heat shock proteins play an essential role in maintaining homeostasis by interacting with stress-denatured proteins during exposure to proteotoxic stressors (Lindquist,1986).In the present study,theHSP70Dk proteins 1/2/5/6/8 were upregulated in the CG/TSG and TSG/TR groups,indicating thatHSP70s play a role in maintaining cell homeostasis during thermal stress and recovery.The cytoskeleton is the foundation for maintaining normal physiological activities in cells.Maintaining the stability of the skeletal structure is a prerequisite for cell homeostasis during changes in temperature (Kotas and Medzhitov,2015).According to a previous study,the expression level of the constitutive form increases after thermal stress,whereas the inducible form is expressed only after exposure to 32℃(Clegget al.,1998) and overexpression ofHSP70 may correspond to the constitutive form (Meistertzheimet al.,2007).A similar result was found in our study,as a greater number of downregulatedHSP70 genes were identified in the TSG/TR group.HSF1 is a regulatory stress-induced protein that binds to the heat shock element.In our observations,the downregulation of heat shock factor 1 (HSF1) indicated that the oyster immunoreaction was more complex under thermal stress.Previous studies have reported that the expressions ofHSP70 andHSP90 inhibit the expression ofHSF1 under thermal stress (Jaegeret al.,2014),which is consistent with our results.The molecular chaperoneHtpGpresented the same expression pattern as the other chaperones in response to heat stress in this study.HtpGis a class IV heat shock gene induced by heat that is involved in unfolded protein binding,ATPbinding,and protein folding (Craig and Schlesinger,1985).In the absence of theHtpGprotein,recovery of cells from 53℃ was restarted,and this delay was eliminated by overproduction ofHtpG(Versteeget al.,1999).A study of cold acclimation in cyanobacteria demonstrated that theHtpGprotein contributes significantly to the ability to acclimate to low temperatures (Hossain and Nakamoto,2002).Similar results were found in our study,as theHtpGgene was upregulated compared to the TSG/TR.All of these results suggest physiological flexibility or adaptation during heat stress regardless of the controversy concerning repression or expression of protein metabolic genes (Limet al.,2016).The upregulation of all chaperones confirmed the severity of the thermal stress in our experiment.The impact on the hybrid oyster immune system remains unclear and such topics will be considered in future studies.

In vertebrates,hormones exert rapid,nongenomic effects,which are mediated by plasma membrane-associated receptors.The estrogen signaling (ko04915) and thyroid hormone synthesis (ko04918) pathways are essential for metabolic and cellular homeostasis (Meyeret al.,2009).Among these two pathways,many genes involved in cytoskeleton structural adjustment were identified in our study,such asACTB_G1,calmodulin (CALM),CANX,and adiponectin (ACDC).ACTB_G1 encodes for a protein involved in the organization of the actin filament and is expressed in response to environmental changes (Krauseet al.,1989).Similarly,ACTB_G1 was rapidly upregulated in the CG/TSG group,suggesting that cytoskeletal reorganization occurs in response to thermal stress,as reported forC.gigasand fish (Buckleyet al.,2006).Among other signaling pathways,three Hippo pathways,ko04390,ko04391,and ko-04392,were strongly enriched during the recovery phase,and previous studies have shown that Hippo pathways play an important role in controlling animal organ size by regulating cell proliferation and apoptosis rates in the sea-anemoneNematostella(Hilman and Gat,2011).These findings indicate that thermal stress can be a stimulus for hybrid oysters to deploy apoptosis or proliferative processes.

Apoptosis is an evolution-conserved process used by multicellular organisms,which has long been reported as an important process in the molluscan immune system.Several studies have shown that apoptosis-related genes are upregulated during recovery from thermal stress (35℃) inC.gigas(Kimet al.,2015).This observation corresponds with our data that the apoptosis pathway was significantly enriched by DEGs in the hybrids during the recovery stage,such as baculoviralIAPrepeat-containing protein (BIRC),TNFSF,and cathepsin L (CTSL).BIRCs are inhibitors of cell death,which act by binding to active caspases (Silke and Vaux,2001).Our results showed thatBIRCs were upregulated in the hybrids during the recovery phase,suggesting that suppressing apoptosis led to better cell survival after thermal stress and could be a key reason for the adaptability of the hybrid.Our findings also agreed with previous studies on osmotic stress strategies inL.vannameiandC.sikamea(♀) ×C.gigas(♂) hybrids (Penget al.,2015).Apoptosis is the sole form of programmed cell death during the life cycle,whereas necrosis is regarded as an unregulated process.Recent studies on cell death have been revitalized by understanding that necrosis can occur in a highly regulated and genetically controlled manner (Bergheet al.,2014).Necroptosis is initiated by different stimuli,such asTNFandTNF-related apoptosis-inducing ligand(TRAIL),requiring the kinase activity of receptor-interacting protein 1 (RIPK1) andRIPK3.Necroptosis contributes to the innate immune response by killing pathogen-infected cells and by alerting the immune system through the release of danger signals (Dondelingeret al.,2016).HSP90 is a key regulator that prevents necroptosis under physiological and stress conditions in eukaryotic cells (Kimet al.,2012).Our data showed that necroptosis was apparent in the hybrids during thermal stress,while theTRAILgene was downregulated.Upregulated apoptosis and necroptosis were also found inSinopotamon yangtsekienseexposed to cadmium (Liuet al.,2011).We speculate thatHSP90 genes were deployed when necroptosis was activated by theTRAILgene under thermal stress in the hybrids to protect against unpredicted apoptosis.Although several studies have described the suppression of necroptosis under thermal stress(Yang and He,2016),this hypothesis needs further verification.

Ubiquitin is an important aspect of the innate immune response to stress by inducing the degradation of irreparably damaged proteins (Pickart and Eddins,2004).TheE3 ubiquitin-protein ligase is upregulated in hybrid oysters when ambient osmotic pressure changes.Another study suggested that ATP-dependent proteolysis for the survival ofC.gigasis suppressed by high temperatures (Limet al.,2016).In our study,the regulation ofE3 ubiquitin-related genes was expressed during thermal stress,as the casitas b-lineage lymphoma and zinc and ring finger 3 (ZNRF3-like) genes were upregulated,but the ring-box 1 (RBX1),tripartite motif-containing (TRIM71),and ring finger protein 115 (RNF115) genes were downregulated.Additionally,an alternative isoform (X1/X2) of theHECTD1 gene displayed opposite regulation during thermal stress,which further supports its multifunctional role in innate immunity.The transcriptome data also showed thatE3 ubiquitin plays an important role during the thermal stress recovery phase,asTRIM13-like,TRIM36-like,and x-linked inhibitor of apoptosis proteins were upregulated.Therefore,it is presumed that the regulation of protein degradation in response to thermal stress was related to the survival and temperature-adaptation of the organisms.Ubiquitination is not only involved in protein degradation but is also closely related to energy metabolism.The cellular stress response has an energetic cost and controls the balance between ATP supply and demand in the gill during thermal stress (Meistertzheimet al.,2007).The rapid induction of lactate dehydrogenase (LDH) in the gill indicates that anaerobic metabolism is required during the early phase of warming(Long,1976).This is in line with our results,as a number of genes involved in metabolic regulation were also regulated in response to temperature changes.

5 Conclusions

Our study represents the first transcriptome report ofC.sikamea♀ ×C.gigas♂ hybrids under thermal stress and recovery using next-generation sequencing technology.A total of 1293 and 740 significant DEGs were detected during and after thermal stress.In our study,several up-and down-regulated genes were identified in diverse pathways,including immune stress,apoptosis,and energy metabolism.Our findings highlight a complex network of immunological and metabolic pathways in the gill of hybrid oyster.In conclusion,these results facilitate future research on immune mechanisms and metabolism in hybrid oyster during thermal stress and recovery.

Acknowledgement

This study was supported by the National Natural Science Foundation of China (No.31172403).