The Phenotypic and the Genetic Response to the Extreme High Temperature Provides New Insight Into Thermal Tolerance for the Pacific Oyster Crassostrea gigas

Investigating responses of organisms to stressful or new environments with selection pressure is one of the crucial problems in evolutionary biology, and it is of importance to understand the phenotypic and molecular mechanism underlying thermal tolerance under the context of the climate change. The Pacific oyster, Crassostrea gigas, inhabiting the environment with high variation in temperature, is a worldwide aquaculture species. However, summer mortality relevant to the high temperature is one of the problems challenging the oyster industry. An artificial selective breeding program was initiated to select for the thermal tolerance of oysters in an attempt to increase the summer survival rates since 2017 in our study. Parents of thermotolerance oyster selection is are based on acute thermal tolerance under controlled heat stress to strengthen the selection intensity. Furthermore, the phenotypic and the genotypic response to extreme high temperature were investigated based on the comparison of the F1 progeny of the selected and natural populations in growth, physiology, mortality rate post heat stress, genetic structure, and gene expression. The parameter of growth showed no significant change between the selected and natural populations for the 6-month-old oysters. The selected population exhibited a higher survival rate after exposure to heat stress in the laboratory, which is in line with result of the filed experience that summer mortality of selected population was significantly lower than that of natural population. Further, the respiration rate of the selected population increased at 38°C, while it increased at 35°C in the natural population. Simultaneously, metabolism-related enzymes (PK, SOD) showed higher activity levels in the selected population. Furthermore, phylogenetic analysis, population structure, and principal component analysis (PCA) suggested that the selected and natural populations exhibited genetic divergence, with eight genes (IF4A2, IF6, EIF3A, MANBA, DDX43, RECS, CAT2, and BAG4) in the selected regions showing differential expression patterns in response to heat stress in the two populations. This study suggests that artificial selection has a significant effect on phenotype and genome structure for the oyster, our study providing an alternative way to reveal the mechanism underlying thermotolerance that plays an important role to predict the potential adaptation to the climate change.

Investigating responses of organisms to stressful or new environments with selection pressure is one of the crucial problems in evolutionary biology, and it is of importance to understand the phenotypic and molecular mechanism underlying thermal tolerance under the context of the climate change. The Pacific oyster, Crassostrea gigas, inhabiting the environment with high variation in temperature, is a worldwide aquaculture species. However, summer mortality relevant to the high temperature is one of the problems challenging the oyster industry. An artificial selective breeding program was initiated to select for the thermal tolerance of oysters in an attempt to increase the summer survival rates since 2017 in our study. Parents of thermotolerance oyster selection is are based on acute thermal tolerance under controlled heat stress to strengthen the selection intensity. Furthermore, the phenotypic and the genotypic response to extreme high temperature were investigated based on the comparison of the F 1 progeny of the selected and natural populations in growth, physiology, mortality rate post heat stress, genetic structure, and gene expression. The parameter of growth showed no significant change between the selected and natural populations for the 6-month-old oysters. The selected population exhibited a higher survival rate after exposure to heat stress in the laboratory, which is in line with result of the filed experience that summer mortality of selected population was significantly lower than that of natural population. Further, the respiration rate of the selected population increased at 38 • C, while it increased at 35 • C in the natural population. Simultaneously, metabolismrelated enzymes (PK, SOD) showed higher activity levels in the selected population. Furthermore, phylogenetic analysis, population structure, and principal component analysis (PCA) suggested that the selected and natural populations exhibited genetic divergence, with eight genes (IF4A2, IF6, EIF3A, MANBA, DDX43, RECS, CAT2, and

INTRODUCTION
Investigating organisms' responses to stressful or new environments with selection pressure is one of the important problems in evolutionary biology. Moreover, it is of importance to understand the mechanism of the thermal adaptation not only because the temperature has profound effects on organisms' physiological and biochemical reaction, but also the ability of organisms to cope with thermal stress is a significant predictor of adaptive potential to climate change. Local adaptation and phenotypic plasticity of organisms may play vital roles in forming intraspecific thermal adaptation divergence under the selection for the tolerance to suboptimal temperature (Sanford and Kelly, 2011;Yampolsky et al., 2014). Although more and more evidences appeared to support the fact that the natural selection is pervasive leading to the divergence of thermotolerance even in the marine invertebrate with long larval dispersal, it is still challenging to infer the adaptive divergence among the natural populations with high connection due to the slight genetic divergence . The thermal response to acute heat stress was proved to be highly relevant to the capability to adapt to chronic stress (Yampolsky et al., 2014). The artificial selection for the thermal tolerance to the acute heat stress may provide an effective strategy to intensify the adaptive response to the high temperature under the controlled condition.
As marine molluscs distributing in the intertidal and shallow subtidal zone, oysters experience daily and seasonal variation in temperature. As such it provides an ideal model to study the thermal adaption (Przeslawski et al., 2005). The Pacific oyster Crassostrea gigas originates in East Asia and has become the primary oyster species supporting the worldwide shellfish industry, due to its rapid growth rate and high productivity (Guo, 2009). However, mass mortality is a major problem challenging the oyster industry. Significant recurrent summer losses of the C. gigas have occurred globally for over five decades and the concern over its impact on oyster aquaculture has increased in recent years (Glude, 1975;Cheney and Macdonald, 2000;Soletchnik et al., 2005;Paul-Pont et al., 2013;Alfaro et al., 2018). Elevating temperature is considered as one of the most important abiotic causing the summer mortality of oysters. It has been reported that summer mortality occurs regularly in culture farms as well as on natural oyster beds during the period of summer when temperature increased (Soletchnik et al., 2005). Elevated temperature have been proposed to increase the spread of Vibrio sp. from infected to non-infected oysters (Petton et al., 2013). Therefore, exposure to high temperatures was one of most important causes linked to summer mortality of Pacific oyster.
The artificial selection for the survival rate during summer was proved to be effective to improve the survival rate of the oyster in summer. A study on the West Coast of the United States increased survival rate among the resulting progeny. Average mortality rates of selected third generation families ranged from 13 to 23% compared with 63% of the natural group (Beattie et al., 1980;Hershberger et al., 1984). The studies in France observed a high genetic basis for survival of young C. gigas (under 1 year of age) in field trials during the summer (Samain et al., 2007). The high heritability (0.47-1.8) of this trait has also been demonstrated by divergent selection over successive years and generations (Dégremont et al., 2010). All these studies focused on the morphological traits and survival rate while the genetic and molecular mechanism underpinning the selection was missing. An artificial selective breeding program aimed at increasing the survival rate during the summer since 2017 was initiated in our study. Parents of thermotolerance oyster selection are based on acute thermal tolerance under controlled heat stress to strengthen the selection intensity. The comparison of the selected population with the natural population would provide new insight into the thermal adaptation by investigating the phenotypic and the genetic response to the artificial selection for the thermotolerance.
Most of the studies of the response to heat stress in marine mollusks focused on the physiological aspects at the organismal, cellar and molecular level (Gagnaire et al., 2006;Farcy et al., 2008;Zhang et al., 2016). Thermal stress has a pronounced impact on the physiology of many intertidal species (Yao and Somero, 2014), including reef fishes (Mora and Ospína, 2001), marine snail Chlorostoma funebralis (Gleason and Burton, 2013), geoduck clams Panopea globose (Juarez et al., 2018). Variation in temperature influences the activity of enzymes that are involved in the metabolic processes that support their growth and survival (Tang and Liu, 2005;Alter et al., 2017;Wang et al., 2018). Animals will increase metabolism to ensure adequate energy supply for survival when temperature is elevated (Mora and Maya, 2006;Sokolova et al., 2012). Beyond the thermal optimum, aerobic respiration may transition to anaerobic metabolism because of insufficient oxygen supply and decreased metabolism (Portner, 2001;Anestis et al., 2010). Transcriptomic studies revealed an extensive set of heat stress responsive genes related to energy metabolism, anti-apoptotic pathways and mRNA processing (Yang et al., 2017). There have been reports on dramatic difference in expression level of the HSP70 gene family and genes associated with oxygen binding between susceptible and thermotolerance oysters (Kim et al., 2017). Divergent transcriptional profiles of genes associated with intracellular stress mechanisms of antioxidant defense, energy metabolism and the cytoskeleton were observed in two population of Sydney rock oysters Saccostrea glomerata (one that had been selectively bred over seven generations for fast growth and disease resistance and one wild type) on the hottest pavers (McAfee et al., 2018).
The studies related genetic responses to the natural and artificial are relatively rare, especially for identification of genetic variation linked to the adaptation. A small proportion of adaptive loci were inferred based on the natural populations with slight genetic divergence . Genetic analysis based on neutral marker such as microsatellite has been used to evaluate genetic differentiation of the populations of cultured oysters subject to artificial selection (Li et al., 2006) but has not produced sufficient information to understand the mechanism of selection. Next-generation sequencing technologies have promoted the generation of massive SNP markers, which may serve as a useful tool for analyzing population genetics and phylogenomics and identifying varieties (Gupta et al., 2015). Specific length amplified fragment sequencing (SLAF-seq) was developed based on highthroughput sequencing technology, and it has been widely used to elevate understanding of genetic populations in plants and animals (Li Z. et al., 2017;Ge et al., 2019). In summary, the integration of the physiological and genetic response may provide an efficient way to reveal the underpinning of the thermotolerance.
Therefore, the motivation behind this study was to further understand the mechanism in the selection for thermotolerance to acute high temperature. Specifically, we aimed to collect phenomics and genomics evidence to identify discrepant responses to elevated temperature between thermotolerance oysters after one-generation of selection and unselected oysters. The effects of thermal stress on survival rate, respiration rate, the activities of key metabolic enzymes, the expression levels of candidate genes were assessed. The genetic difference was identified by SLAF-seq approach. This study will be useful to understand the potential molecular mechanisms for heat tolerance in oysters.

Sample Collection
In October 2017, wild populations of Pacific oyster (Crassostrea gigas) were collected from the Shentanggou (36 • 21 N), Qingdao, Shangdong Province, China. Oysters were cleaned to eliminate attaching organisms and acclimated in an aquarium for 14 days with sand-filtered and aerated seawater before experimentation. During the acclimation period, individuals were fed spirulina powder, and seawater was changed daily. Oysters were randomly divided into two populations (the number of individuals in each population is greater than 300). One population served as the parent of the natural population (hereinafter referred to as control population), and continued to be cultured without any treatment. The other population was subjected to heat shock at 42 • C (LT 50 ) for 1 h . Specifically, oysters were submerged in preheated, aerated filtered seawater baths of 42 • C for 1 h, then returned to sea water of ambient temperature for 14 days, before removing and counting dead individuals. Finally, the surviving individuals were selected as parents to produce the next generation of the thermal tolerance population (intensity of selection = 37%).
In April 2018, the progeny of selected parents and control population, were used as the experimental F 1 generation population and the control population, respectively. We use 30 males and 30 females to reproduce progeny. Specifically, the eggs of the 30 females were mixed and equally divided into 30 beakers, each of which was fertilized by sperm from one of the 30 males. The larvae were reared in the hatchery and juvenile oysters were cultured in the same sea area using standard protocol (Guo et al., 2012). Six months old F 1 generations were collected from the sea and acclimatized to the laboratory aquarium for 2 weeks for follow-up experiments. The treatments of experiment material and processes are shown in Figure 1.

Growth Measurements
The individuals of the selected F 1 generation were used as the selected population and the individuals of control F 1 generation were used as the control population. A total of 150 oysters were collected from each population. The wet weight and shell height were measured after cleaning the epibionts.

Analysis of Thermal Tolerance
Thermotolerance of the populations was assessed by detecting the survival rate post-acute thermal treatment. Six-month-old oysters were used to examined thermal tolerance (n = 300). The oysters from each population were subjected to heat shock treatment at 42 • C (LT 50 ) for 1 h and then removed to ambient temperature. During the recovery process, oysters were fed spirulina powder, and the seawater was changed daily. Visual monitoring of oyster shell activity (open or closed) was performed. The individuals whose shell remained open even after being challenged were considered dead and immediately removed. The survival rate after heat shock for 20 days was recorded.

Respiration Rate Analysis
Six-month-old oysters of similar size (n = 6, shell height: 42.76 ± 3.01 mm [selected population], 41.33 ± 3.56 mm [control population], p = 0.215) were randomly selected to estimate respirate rate at different temperature. One oyster was measured at a time. Individuals were cleaned using seawater to avoid the influence of periphyton. The metabolic rate was measured at 20 • C to simulate the natural environment and at 35, 37, 38, and 39 • C to simulate acute heat stress, respectively. Individuals were placed in a 330 mL glass bottle filled with sand-filtered, air-saturated seawater at the designed temperature, which was regulated using a water bath. The seawater was mixed using a magnetic stirrer bar beneath the chamber. Needle-type fiber-optic microsensor oxygen sensor and temperature probe were connected to the oxygen transmitter that simultaneously recorded the temperature and oxygen concentration of water every 3 s for 1 h . The slope of the decrease in oxygen concentration was calculated as the respiration rate (mg/mL/h). FIGURE 1 | Flow chart of the experiment. (A) Wild C. gigas was randomly divided into two population (n > 300). One population was subjected to heat shock at 42 • C for 1 h, and surviving individuals served as the parent of the next generation of selection for thermotolerance population, while the other population served as the control population without any treatment. The individuals from selected F 1 population and control population were exposed at 37 • C seawater. Gills from 15 individuals were sampled at each time points of 0, 6, 24 h, which used to determine physiological parameters and gene expression. (B) After heat shock at 42 • C for 1 h, selected F 1 population and control population were kept at ambient temperature restoration for 2-weeks recovery, and mortality was recorded. The gills of dead individuals from control population during the recovery period and that of the survival individuals from selected F 1 population after 2 weeks were sampled for SLAF-sequencing (n = 50).

Assessment Impact of Heat Stress
The gills of oysters after heat shock were used to examine the physiological parameters and expression level of heat responsive genes. Oysters from each population were randomly divided into three replicates and subjected to heat stress of seawater at 37 • C for 24 h (shell height: 41.7 ± 5.63 mm [selected population], 40.38 ± 6.64 mm [control population], p = 0.14). The gills of five oysters were individually sampled from each replicate per population at 0, 6, 24 h (n = 15 per sampling time per population). The gills were immediately frozen in the liquid nitrogen and stored at -80 • C for follow-up analysis of physiology and gene expression.

Physiological Assessment
According to the ratio of weight: volume = 1: 9, 30 ∼ 40 mg of frozen gill tissue from each sample of every sample time of each population was added to an excess volume (9X) of precooled saline and ground to homogenate in an ice bath. The supernatant was diluted to a 1% homogenate by adding ninefold saline after centrifuging for 10 min at 2500 rpm and 4 • C. The 1% homogenate was then stored at -80 • C until required for further analysis.
The fundamental physiological indexes, including the content of total protein and malondialdehyde (MDA), and the enzymatic activity of pyruvate kinase (PK), and superoxide dismutase (SOD), were chosen to evaluate the dynamic physiological responses to high temperature in oysters . The concentration of protein was used to normalize other enzyme activity indexes. The first three indexes were measured using the corresponding kit from Nanjing Jiancheng Bioengineering Institute (Nanjing, China). SOD levels were determined using the Total Superoxide Dismutase Assay kit with NBT (Beyotime, Shanghai, China). Absorbance values were measured using the Varioskan Flash (Thermo Fisher Scientific, Waltham, MA, United States).

Field Experiment
Oysters were cultured in Shentanggou, Qingdao. The monthly average sea surface temperature (SST) of Qingdao in 2019 was showed in Figure 2 which is not significantly different from the past 18 years (Ghaffari et al., 2019). The highest average sea surface temperature was observed in August is 25.57 ± 0.71 • C. The salinity of seawater was 30.0 ± 0.5 . When the F 1 oysters were 19-month-old, the survival rate and shell height of the selected and control population were recorded to investigate the thermal resistance of the selected population in the field.

Sample Collection and DNA Extraction
After heat shock at 42 • C for 1 h, the F 1 selected and control populations were placed to restore at room temperature for 2 weeks. The gills of individuals from F 1 control population that died during the recovery period, and those of the live individuals from the F 1 selected population were sampled after 2 weeks, before freezing in liquid nitrogen and storing at -80 • C (n = 50 per population). Thus, the oysters with the potential largest differences in thermal tolerance were used for genomic analysis to identify SNPs and candidate genes associated with thermal tolerance. Total DNA was extracted from the gill tissue using the TIANamp Marine Animals DNA Kit (Tiangen, Beijing, China) according to the manufacturer's protocol. DNA quality and concentration were measured by 1.0% agarose gel electrophoresis and UV spectrometry on a NanoDrop 2000 device, respectively.

SLAF Library Construction and Sequencing
Simulated restriction enzyme digestion was performed on the oyster genome in order to avoid repetitive SLAF tags and to obtain a relatively uniform distribution of restriction fragments in the genome for maximum SLAF-seq efficiency. As a result, genomic DNA was digested with RsaI and HaeIII restriction enzyme combination. To evaluate the efficiency and accuracy of the experiment procedure, Oryza sativa ssp. Japonica DNA (Ohyanagi, 2006) was served as control to evaluate the availability of the enzyme digestion. Specific-locus amplified fragment library construction was carried out as previously described (Sun et al., 2013). DNA fragments of 264-314 bp were defined as SLAF tags and diluted for pair-end sequencing on an Illumina High-seq 2500 system (Illumina, Inc., San Diego, CA, United States) at Beijing Biomarker Technologies Corporation.

SLAF-seq Data Processing and SNP Marker Development
The raw paired-end reads were mapped to the oyster reference genome (oyster_v9, GenBank accession number GCA_000297895.1) using BWA software (Li and Durbin, 2009). The position of clean reads was located by comparison with the reference genome, the sequencing depth of each sample was counted, and mutation were detected. The development of SNP markers was performed by both GATK (McKenna et al., 2010) and SAMtools  analysis, the intersection of SNP markers obtained by the two methods as the final reliable SNPs.

Population Genetic Structure
Filtering high-quality SNPs for downstream analysis according to the following criteria: minor allele frequency ≥ 5% and the integrity of each SNP ≥ 80%. A total of 379,002 SNPs from 100 individuals were development for analysis of phylogenetic relationships and population structure. MEGA X (Kumar et al., 2018) was employed to construct phylogenetic trees of each sample with Kimura 2-parameter based on the neighbor joining (NJ) method. The bootstrap repeated 1,000 times was applied to assure the reliable of the phylogenetic trees. The population structure was analyzed using the ADMIXTURE program (Alexander et al., 2009). We conducted a K-means clustering analysis supposing the K-value ranged from 1 to 10. The clustering results were cross-validated, and the optimal group number was determined according to the valley value of the cross-validation error rate. A good K-value will exhibit a low cross-validation error compared to other K-values. Principal component analysis (PCA) was performed with EIGENSOFT program (Price et al., 2006). We can access the relationship among samples through PCA analysis which can assist evolutionary analysis.

Identification of Differentially Selected Regions and Functional Enrichment Analysis
The divergence index, F-statistics (F st ) and the index of nucleotide polymorphism (θπ) were calculated via the PopGenome package (Pfeifer et al., 2014) based on 100kb sliding windows in 10 -kb steps. We calculated and used the θπ ratio (θπ, control population/θπ, selected population) and F st values between the two populations to detect genomic regions with signatures of selection. The region with significantly low and high θπ ratios (low 5% regions and top 5% regions) and significantly high F st values (top 5% regions) was identified as selected regions with strong selection signals. Gene Ontology (GO) enrichment analysis of target genes were performed in the topGO R package. GO terms with a p-value < 0.05 were considered significantly enriched and genes on these terms were considered as candidate genes.

Expression Level of Candidate Genes
The materials we used here were described in above-mentioned assessment impact of heat stress. GO terms with a p-value < 0.05 were considered significantly enriched and the genes in these terms were identified as candidate genes. Although not enriched in the GO analysis, BAG4 was selected due to its potential role in temperature adaptation, and was also found in the region of top 5% F st in our study. Therefore, BAG4 was considered a candidate gene.
The gills were used for total RNA extraction with the RNA prep Pure Tissue kit (Tiangen, China) following manufacturer's protocol. All RNAs were treated with DNase I to avoid genomic contamination in the RNA. The quality and concentration of RNA were assessed using 1.2% agarose gel electrophoresis and a NanoDrop 2000 spectrophotomer. Equal amounts (1 µg) of RNA from five oysters at each sample time of each population were pooled together to reduce the effect of error (three biological replicates × three technical replicates). Reverse transcription was performed using PrimeScript RT reagent Kit (TaKaRa Bio, Shiga, Japan) on 1 µg of total RNA to synthesize cDNA following the manufacture's protocol. Synthetic cDNA was diluted 20-fold and its expression was determined.
The 18 candidate genes enriched in GO terms with p < 0.05 and BAG4 were selected for real-time PCR analysis. The gene of elongation factor 1 alpha (EF1-α) was used as a reference gene due to its low expressional variability during heat stress . Primers designed for each gene with Primer Premier 5.0. qPCR was performed using the ABI7500 Fast Real-Time Detection System (Applied Biosystems, Foster City, CA, United States). The total volume of the reaction system was 20 µL and consisted of 10 µL SYBR Green 2X Supermix (TaKaRa), 6.8 µL DEPC H2O, 0.4 µL of each primer pair and ROX DyeII, 2 µL diluted cDNA. The program started with a 30 s activation of DNA polymerase at 95 • C, followed by 40 cycles of 5 s at 95 • C and 30 s at 60 • C. The melt curve stage was conducted as: 15 s at 95 • C, 1 min at 60 • C, 30 s at 95 • C, and 15 s at 60 • C. The expression level of gene of induced transcripts (6 and 24 h) was calculated using the Livak 2 − CT method (Livak and Schmittgen, 2001). The expression of the control population was used as a control to calculated basal expression level.

Statistical Analysis
All data were expressed as mean ± standard error of mean (SEM). Statistics analysis was performed with SPSS 22.0 statistical package (IBM, NY, United States). The differences in shell height and wet weight between the two populations were analyzed using one-way analysis of variance (ANOVA), followed by Duncan's multiple range tests for post hoc comparison of means. The same method was used to calculate the variation in respiration rate, gene expression, and physiological parameters between the two populations. Survival time in the laboratory survival module was determined by the Kaplan-Meier analysis with long-rank test. The survival rate in field was obtained by counting the number of deaths and survival on each of attachment substance. The level of statistically significant differences was set at p ≤ 0.05.

Growth
The shell height and wet height of selected population were 41.70 ± 5.63 mm, 8.88 ± 2.87 g respectively, and those of control population were 40.38 ± 6.64 mm, 8.71 ± 3.58 g ( Figure 3A). There was no significant difference in shell height and wet weight between the selected and control populations (p > 0.05).

Survival
The survival rate after heat shock at 42 • C for 1 h was used to evaluate the thermal resistance of oysters. The selected oysters exhibited a significant higher survival rate (91.67%) during the recovery period compared control oysters (83.33%) following post-acute heat stress (p < 0.05) ( Figure 3B). The first mortality of oysters under post treatment was observed in both selected and control populations on day 2. The highest mortality was observed on day 4 in the selected population while on day 7 in the control population. The continuous low and average mortality was observed after highest mortality.

Respiration Rate
The effects of temperature on oxygen consumption rate in two populations were not significant at 20 • C (p > 0.05). Specifically, the average oxygen consumption of the selected population was 0.085 mg/mL/h while that of control population was 0.077 mg/mL/h. However, the differences were significant when temperature increased to 35 • C (p < 0.01), and the average oxygen consumption of the control population reached 0.142 mg/mL/h. In contrast oxygen consumption in the selected population was 0.087 mg/mL/h. The respiration rate of the selected oysters was significantly lower than that of unselected oysters at 35 • C, but the opposite result was shown when the temperature increased from 35 to 38 • C (p < 0.05), when the average oxygen consumption of the selected population reached 0.166 mg/mL/h compared with 0.120 mg/mL/h in the control population ( Figure 3C).

Physiological Assessment
Significant effects of time on MDA concentration were observed in the selected population which displayed a tendency to accumulate during exposure to high temperature. Moreover, we observed significant differences at 6 h between the selected and control populations (p < 0.05). Total protein content of the FIGURE 3 | Comparison of growth and heat tolerance between selected and control populations, respectively, and comparison of the effect of elevating temperature on respiration between two populations. The error bar denotes the standard error of the mean. * indicated p < 0.05, * * indicates p < 0.01. (A) The shell height and wet height of selected and control populations, respectively (n = 150 per population). (B) The survival rates of selected oysters and control oysters after heat shock at 42 • C for 1 h (n = 300 per population). (C) The effect of increasing temperature on the respiration between selected oysters and control oysters (n = 6, six oysters from each population were used for each temperature). selected population was significantly higher than that of the control population during heat stress (p < 0.01). Furthermore, there was significant variation in activity of enzymes between the two populations. The activity of SOD and PK in both populations showed a trend of increasing in the beginning and then decreasing during heat shock. Moreover, the enzyme activities of SOD and PK were significantly higher in the selected population than in the control population during heat stress (p < 0.05) (Figure 4).

Field Experiment
The mortality of the 19-month-old selected population was lower than the control population during summer under field conditions, which is consistent with our laboratory result. The survival rate of the selected population was 74.19% which is significantly higher than the value of 62.33% (p < 0.05) for control population. The average shell height of the selected and control populations was 69.53 and 73.38 mm, respectively. The shell height of the selected population was not significantly different from the control population (p > 0.05).

SLAF Sequencing and SNP Discovery
A total of 286.82 million reads were obtained from 100 individuals, which showed the average Q30 was 94.40%, and the average GC content was 34.55%, indicating the high quality of the generated data ( Table 1). A total of 593,103 SLAF tags were obtained, including 454,594 polymorphic SLAF tags, with an average sequencing depth of 14.83x. In addition, Oryza sativa ssp. Japonica was used as a control in this project. Our results revealed that the percentage of paired-end mapped reads and normally digestion of control was 96.31 and 91.08%, respectively, indicating that the SLAF-seq process was normal and reliable. After genomic mapping and SNP calling with GATK and SAMtools, a total of 7,354,783 SNPs were revealed using all individuals, of which the per individual averages of 198,765 were intergenic, 300,145 were intron and 368 were exons. A total of 1,076,631 indels were identified in this program, of which 9,095 were in coding regions.

Population Structure and Phylogenetic Relationship
A strict quality control filtering of SNPs was performed to identify 379,002 SNPs used for further analysis. Phylogenetic analysis was conducted to determine the relationships between the two populations (100 individuals). Neighbor-joining cluster analysis divided 100 individuals into three groups ( Figure 5A). The result clearly showed that the clustering of most selected individuals separated from the control population, and a few control individuals mixed with the selected population, which was consistent with the assignments made by PCA. The PCA  analysis separated the 100 individuals into three clusters. The first three PCs accounted for 15.7% of the population-wide SNP variation, with PC1 accounting for 7.78% ( Figure 5B). The estimated membership fractions for different values of K ranged from 1 to 10, and the maximum likelihood revealed by the population structure showed an optimum value of 3 (K = 3) (Figures 5C,D), which indicated that the entire population could be categorized into three groups: group 1, group 2, and group 3. Group 1 contained 27 individuals from the control population. Group 2 contained 23 individuals, two of which were from the control population, while the remainder were from the selected population. Group 3 contained 50 individuals, 29 of which were from the selected population, and 21 of which were from the control population ( Table 2).

Identification of Selected Regions and Candidate Genes
To obtain the genomic difference between the selected and control populations, a SNP-based F st and θπ estimate was performed to identify selected regions. The value of the top 5% paired-F st ranged from 0.0851 to 0.3253. As a result, analysis identified 472 genes from 115 genomic regions showing selection signals between the two populations. GO enrichment analysis of all 472 genes under selection showed that seven GO terms covering 18 genes were significantly over-represented (p < 0.05) ( Table 3). The previously identified adaptive gene, BAG4 , which is in the selected regions of top 5% F st value was also selected as a candidate gene for subsequent gene expression investigation.

The Expression Level of Candidate Genes
The expression of the 19 genes in the selected region were further investigated by real-time PCR. The relative induced expression of eight genes (IF4A2, IF6, EIF3A, MANBA, DDX43, RECS, CAT2, and BAG4) were significantly different between the selected and control populations (p < 0.05), and the primers of each gene are shown in Table 4. The relative basal expression under ambient environment of six genes (IF6, EIF3A, DDX43, RECS, CAT2, and BAG4) in the control population exhibited significantly higher expression than in the selected population (1.42-2.64fold change, p < 0.05) (Figure 6). In the selected population, the expression level of IF4A2, DDX43, and BAG4 showed upregulation under heat stress, whereas induced expression of IF6, EIF3A, MANBA, RECS, and CAT2 showed down-regulation. The changed range of genes expression in selected population from 0.41 to 35.41-fold. However, most genes in the control population showed more moderate expression patterns under heat stress, the induced genes expression changed range from 0.55 to 2.27fold (Figure 7).

DISCUSSION
In this study, the genetic and phenotypic effect of one-generation selection oysters targeting to improve the thermal tolerance was evaluated. The selection population showed higher survival rate, increased oxygen consumption at higher temperatures, and higher enzyme activity than control population (p < 0.05). In addition, the two populations exhibited difference in both genome structure and gene expression. No significant differentiation of growth was observed between the selected and control populations in our study. This result is consistent with studies which have shown that selection on survival during the summer period after three generations does not have any impact on growth (Dégremont et al., 2010). Moreover, selection for growth was conducted but no effect on survival was observed in Sydney rock oyster Saccostrea glomerata (Hand et al., 2004). Our results show no significant difference of growth between the selected and control populations, which suggested that the selection for thermotolerance would not affect the growth at least at juvenile stage (<1-year-old) Pacific oysters.
We used the survival rate after acute heat stress to assess the thermotolerance of oysters. High heritability for survival during summer was observed, which suggests that selective breeding could effectively improve the thermotolerance of juvenile Pacific oysters . In this study, the survival rate of the progeny post exposure to acute heat stress at 42 • C in the selected population was significantly higher than that in the control population even only after one-generation selection, suggesting the thermotolerance of selected oysters was improved. This result was in line with our expectations, which demonstrated that our methodology for thermotolerance selection of oysters was effective. Moreover, it showed that the mortality of selected population was lower than that of the control population under field conditions (p < 0.05). The selected population demonstrated a better survival rate than the control population, indicating a positive response to selection for survival. This result is consistent with the results of a previous study that significant higher survival was observed in the "high survival" selected population during summer at three different sites for F2 and F3 generations in 6-month-old Pacific oysters (Dégremont et al., 2010).
Temperature is one of the most important abiotic factors affecting the physiological response in marine organisms (Windisch et al., 2011). In this study, the oxygen consumption rate of Pacific oyster was greatly affected by temperature and showed a fluctuation with increase in temperature. The respiration rate obtained in this study increased with temperature, up to a maximum or optimum limit beyond which they rapidly decreased. This was consistent with the result of previous studies for other bivalve species, the temperature has a significant effect on oxygen consumption and respiratory energy loss of Calafia mother-of-pearl oyster, Pinctada mazatlanica, which suggested the existence of compensatory mechanisms that facilitate the organism to adapt to the temperature changes (Bayne and Newell, 1983;Saucedo et al., 2004). The individuals from the control population exposed to acute thermal stress (35 • C) showed an increment in their metabolic rates compared to the condition (20 • C), which is a trend also observed in different species of marine bivalves such as the flat oyster Ostrea edulis and mother-of-pearl oyster Pinctada mazatlanica (Haure, 1998;Saucedo et al., 2004). The oxygen consumption rate of selective oysters increased at higher temperature (38 • C) compared to the control population, indicating the respiration of selective oysters reacted less sensitively than unselected oysters at these temperatures by increasing rapidly and the selective oysters had a higher tolerance to high temperature than the control population. Furthermore, 35 and 38 • C represent different stress strength, the selected population showed higher oxygen consumption rate under higher stress strength than control population, indicating that the selected population has greater capacity of oxygen supply that decides the wider scope for metabolic adjustment, which is limited to a particular thermal tolerance window. Thus, the wider thermotolerance window of selected population decides higher survival at extreme high temperature compared with control population. This phenomenon is common in aquatic animals (Portner, 2001;Portner et al., 2006;Portner and Knust, 2007). The greater capacity of oxygen supply in selected population may be the reason that selected population is more thermotolerance. Low metabolic sensitivity of selected population to fluctuations of temperature may translate into lower long-term metabolic costs and higher heat tolerance. The metabolic rate increases as temperature increases, generating higher production of reactive oxygen species or ROS (Lushchak and Bagnyukova, 2006). Marine invertebrates are equipped with antioxidant enzymes 4 | The genes with differential expression level between selected and control oysters and the primers used for q-PCR quantification. such as SOD, CAT, and glutathione peroxidase to prevent cellular damage from ROS accumulation, such as lipid peroxidation, protein degradation and DNA breaking (Lesser, 2006). The activity of enzymes PK and SOD first exhibited an increasing pattern and then a decreasing pattern during the period of heat stress, which indicated the transition from aerobic to anaerobic metabolism (Anestis et al., 2010;Sussarellu et al., 2012). Enzyme activity and total protein of the selective population were significantly higher than those of the control population at most time points (p < 0.05). However, the concentration of MDA did not exhibit divergence between the two populations. This result was consistent with our previous study , which showed southern oysters inhabiting warmer inshore locations had higher enzyme activity than northern oysters. The higher metabolic, enzymes' activities could be explained as protective mechanism to alleviate oxidative damage from elevated amounts of ROS for organism's cope with extremely high temperature (Abele et al., 2002).
FIGURE 7 | Induced expression of candidate genes at 6 and 24 h in response to heat stress in selected oysters and control oysters (n = 15) and significance analysis was conducted. * p < 0.05; * * p < 0.01.
In the present study, we used SLAF-seq to analyze the genetic structure and differentially selected regions between the selected and control populations using a reducing representation library sequencing method (Sun et al., 2013). It is clearly indicated that the surviving oysters from selected population and dead oyster from the control population exhibited divergence in genome structure after the acute stress selection for thermotolerance. This suggested the selection of heat stress is directional rather than random. Moreover, 115 differentially selected regions including 472 genes were identified, which provide a comprehensive view of the genome-wide basis for genetic difference between the two populations.
Population structure analysis identified three subgroups instead of two. The result is consistent with the phylogenetic analysis and PCA. This may be due to that fact that artificial selection is at the early stage and the restricted genetic differentiation between the oysters with one-generation selection and the unselected oysters. As the quantitative trait, thermotolerance of the oyster shows the continuous distribution. Furthermore, the third subgroup may represent oysters with moderate thermotolerance. In another oyster selection program for growth improvement, the population structure analysis showed that the selected lines and wild populations were not divided into two groups even after continuous selective breeding of six-generations (Zhong et al., 2016).
In this study, the top 5% F st and θπ values between selected and control populations were used to determine the selected regions, in which 18 genes were enriched in p < 0.05 GO terms and were considered to be candidate for gene expression detection. In addition, BAG4 in the region with top 5% F st also was selected because of its potential function in temperature adaptation . Surprisingly, no canonical heat shock protein family were enriched in our study. This is consistent with the result that no significant difference of expression of HSP genes between selected and control population (unpublished).
Of the 19 genes showing genomic differentiation between the two populations, eight genes (IF4A2, IF6, EIF3A, MANBA, DDX43, RECS, CAT2, and BAG4) also diverged in transcriptional expression response to heat stress, which suggested divergence of genomic structure may regulate the expression of those genes and the expression differentiation may mediate the superior thermotolerance in the selected oysters. IF4A2, IF6, and EIF3A are related with protein synthesis. The studies of the three genes are related with diseases in human, such as diabetes, gallbladder cancer and nasopharyngeal carcinoma. IF6 and EIF3A play regulatory roles at translation initiation process as a rate-limiting step (Cheyssac et al., 2006;Liu et al., 2011;Golob-Schwarzl et al., 2019). The proteins of DDX43 and RECS are helicase. The overexpression of DDX43 was found in human cancer (Ambrosini et al., 2014). The proteins of MANBA represent one of the most thermostable versions and probably hydrolyze extracellular galactomannans to monosaccharides as energy sources in hyperthermophilic eubacterium (Duffaud et al., 1997). CAT2 is involved in dehydrated process of the fermentation of 4-aminobutyrate to ammonia, acetate, and butyrate in anaerobe (Peter and Wolfgang, 1990). The studies of these genes are rare in invertebrates. The potential effects of these genes in responding to high temperature require further investigations for oysters. It is worth noting that expression pattern of BAG4 in this study is consistent with previous study. BAG4, in intertidal populations and southern populations whose habitat exposed to higher temperature and more fluctuating in temperature exhibited higher plasticity compared with subtidal populations and northern populations respectively under heat stress . This suggested BAG4 may play an important role in thermal adaptation.
The selected population exhibited higher expression flexibility in candidate genes (p < 0.05). Specifically, IF4A2, DDX43, and BAG4 showed a high extent of up-regulated in response to heat shock in the selected population. The other five genes in the selected population exhibited high extent down-regulation in response to elevated temperature, which indicated the selected population has a higher gene expression plasticity under thermal stress. This was consistent with our previous study that plasticity in gene expression positively correlated with evolved divergence, high plasticity in gene expression may be adapt and favored by organisms under thermal regime . Some studies have showed that in order to adapt to fluctuating environments organisms have evolved higher phenotypic plasticity than those habited in the stable environment, which suggested that plasticity may be adaptive (Kenkel and Matz, 2016).

CONCLUSION
We explored divergence of phenotype and genome between the selected and control populations after one-generation selection for thermotolerance in oysters. From analyses of different levels of biological organization (e.g., growth, thermotolerance, oxygen consumption rate, enzyme activities, gene expression, and genetic structure), we can propose that selected population has evolved higher phenotypic plasticity ensuring the oysters have potential adaptation to face extreme high temperature. Furthermore, 115 differentially selected regions including 472 genes were identified indicated the divergence of phenotype between selected and natural population maybe have a genetic background. GO analysis identified and q-PCR result indicated eight genes (IF4A2, IF6, EIF3A, MANBA, DDX43, RECS, CAT2, and BAG4) may be critical in response to thermal stress for oysters. Further, the result of BAG4 was consistent with our previous studies, which suggested its vital function in thermal adaptation.

DATA AVAILABILITY STATEMENT
The raw sequence reads have been submitted to the Sequence Read Archive (SRA) in the NCBI, and accession number PRJNA604406.

ETHICS STATEMENT
This study was conducted on oysters which are not included in the Animals (Scientific Procedures) Act. All experimental works was reviewed and approved by the Experimental Animal Ethics Committee, Institute of Oceanology, Chinese Academy of Sciences, China.

AUTHOR CONTRIBUTIONS
GZ, HQ, and LL conceived the experiment. WW and RC collected oysters and performed the spawning. AL and XW contributed to the experiment in parental oysters. FD performed the experiment, collected and analyzed the data. FD, LL, AL, and GZ wrote and edited the manuscript.