Transcriptomic Responses of Two Ecologically Divergent Populations of Japanese Mantis Shrimp (Oratosquilla oratoria) under Thermal Stress

Simple Summary Rising ocean temperature would change the seawater chemistry and affect the external and internal physiology of crustaceans due to their lack of certain efficient temperature regulators. In addition, the infraspecific populations of crustaceans might also have different response strategies to the rising of temperature. Therefore, we identified the transcriptomic variations to the same thermal stress between ecologically divergent populations of Oratosquilla oratoria. The aim of this study was to investigate the population-specific function genes and relevant pathways in response to thermal stress in O. oratoria. The results showed that gene-expressed variation was in a population-specific pattern, which indicated that the local environment could lead to the evolvement of changes in gene regulation, ultimately leading to adaptive divergences. Additionally, we found several genes with large pleiotropic effects in the Zhoushan population, which might indicate that the regulation mechanisms of the Zhoushan population were more efficient than those of the Qingdao population under same thermal stress. The results provided some novel insights into the local adaptive differences of the infraspecific populations of O. oratoria and other crustaceans. Abstract Crustaceans are generally considered more sensitive to ocean warming due to their lack of certain efficient regulators. However, the alterations in the physiology and behavior of crustaceans in response to thermal stress differ vastly even among the infraspecific populations of heterogeneous landscapes. Consequently, understanding the impact of temperature fluctuation on crustacean infraspecific populations might be essential for maintaining a sustainable persistence of populations at existing locations. In the present study, we chose the Japanese mantis shrimp (Oratosquilla oratoria) as the representative crustacean population, and conducted transcriptome analyses in two divergent O. oratoria populations (the Zhoushan and Qingdao populations) under same thermal stress (20–28 °C) to identify the population-specific expression response to thermal stress. The results showed significant differences in gene expressions, GO terms and metabolic pathways between the two populations. We hypothesized that intraspecific mutations in the same or different genes might lead to thermal adaptive divergences. Temperature increases from 20–28 °C produced significant enrichment in GO terms and altered the metabolic pathways in the Zhoushan population despite the lack of differentially expressed unigenes. Therefore, several functional genes with large pleiotropic effects may underlie the response to thermal stress in the Zhoushan population. Furthermore, the most significantly enriched biological processes of the Qingdao population were associated with the state or activity of cells and its significant enriched pathways with genetic information processing as well as immune and environmental information processing. In contrast, the differentially regulated unigenes of the Zhoushan population were primarily involved in the regulatory cellular and transcription processes and the most significant pathways found were metabolic and digestive. Consequently, the regulatory mechanisms of the Zhoushan population are probably more efficient than those of the Qingdao population under the same thermal stress.

significant differences in gene expressions, GO terms and metabolic pathways between the two populations. We hypothesized that intraspecific mutations in the same or different genes might lead to thermal adaptive divergences. Temperature increases from 20-28 • C produced significant enrichment in GO terms and altered the metabolic pathways in the Zhoushan population despite the lack of differentially expressed unigenes. Therefore, several functional genes with large pleiotropic effects may underlie the response to thermal stress in the Zhoushan population. Furthermore, the most significantly enriched biological processes of the Qingdao population were associated with the state or activity of cells and its significant enriched pathways with genetic information processing as well as immune and environmental information processing. In contrast, the differentially regulated unigenes of the Zhoushan population were primarily involved in the regulatory cellular and transcription processes and the most significant pathways found were metabolic and digestive. Consequently, the regulatory mechanisms of the Zhoushan population are probably more efficient than those of the Qingdao population under the same thermal stress.

Introduction
The current global warming trend is thought to cause widespread effects on the structure and functioning of global fauna and ecosystem processes. Although the ocean has a higher heat capacity and its temperature is relatively stable in comparison with the land, the world's oceans are also continuously warming and the average ocean temperature is predicted to increase approximately 3 • C by the end of the century [1]. It is likely that the rising ocean temperature will change the seawater chemistry and then affect the external and internal physiology of marine organisms [2]. Sunday et al. also considered that most marine organisms belong to poikilotherms and therefore are particularly vulnerable to a rise in temperature of their surroundings [3]. Thus, a rapid regulation capacity and a higher physiological endurance to rising temperatures are very necessary because they can help organisms reduce maladaptive responses [4]. However, responses to rising temperatures might vary in different species. For instance, evidence has suggested that microhabitat species might encounter more pressure than macrohabitat species when habitat temperature continues to rise [5]. In addition, tropical species are expected to be particularly sensitive to increases in temperature because relatively stable thermal environments play an important role in their ecological evolutionary processes [6]. However, the infraspecific populations of heterogeneous landscapes might also have different response strategies to rising temperatures. This may be due to their physiological characteristics and biological interactions being directly affected by their local habitat environment and ultimately leading to genetic adaptation differences [7,8]. Consequently, as oceans become warmer, we can expect to see different pressure on the infraspecific populations that live in different environments. Elucidation of these variations at the transcriptome level would facilitate our understanding of whether these infraspecific organisms have the same ability to cope with ocean warming.
Although numerous studies have revealed that differential habitat temperatures may lead to evolutionarily significant units (ESUs) of many crustacean species [9,10], these studies were constrained by neutral markers that became an impediment for the precise evaluation of the infraspecific adaptability divergence of crustaceans. Recently-developed next-generation sequencing technology has provided a convenient and highly effective solution for biological study. As a result, several studies have more accurately estimated the genetic divergence based on reduced-representation sequencing [11,12], but they could not explain how the temperature fluctuations affected the regulatory mechanisms of marine organisms. However, RNA sequencing (RNA-seq) has clear advantages over the aforementioned approaches and can reveal the complex dynamics of the regulatory processes with both accuracy and sensitivity. To date, the RNA-seq approaches have been successfully used to analyze the regulatory mechanisms of interspecific organisms [13][14][15], and some functional genes and physiological pathways were also identified [16][17][18]. Therefore, with respect to infraspecific organisms in their heterogeneous environments, RNA-seq can also serve as a useful technique for identifying the genetic divergences associated with local environments, which can then be subjected to a further examination of the thermal adaptability of infraspecific organisms [19,20].
It is still worth noting that most marine calcifying organisms (e.g., crustaceans) are generally considered more sensitive to ocean warming because they lack certain efficient regulators [21]. Thus, it is necessary to discuss the long term ability of the crustacean infraspecific organisms to respond to ocean warming. As a crustacean representative, Oratosquilla oratoria (De Haan, 1844) spans a wide geographical area throughout the tropical, subtropical and temperate coastal waters, and its suitable survival temperature ranges from 20-27 • C [22]. Across this wide distribution, the O. oratoria might encounter highly heterogeneous environmental elements. In fact, Du et al. divided two O. oratoria geographical populations (the Yellow Sea and East China Sea populations) [23] and suspected that their genetic differentiation might be associated with temperature variability (Figure 1; [24]). Furthermore, since their habitat differences might ultimately lead to genetic divergence, it is not surprising that two O. oratoria populations vary in thermotolerance. Therefore, O. oratoria provided an ideal research source for investigating the thermal response variation of crustacean infraspecific populations. Furthermore, since their habitat differences might ultimately lead to genetic divergence, it is not surprising that two O. oratoria populations vary in thermotolerance. Therefore, O. oratoria provided an ideal research source for investigating the thermal response variation of crustacean infraspecific populations. The muscle is the largest energy and amino acid pool in the maintenance of homeostasis in ocean animals and because it exhibits a high surface area directly suspended in the water. Therefore, the muscle might be associated with thermal response. However, whether the O. oratoria infraspecific organisms have the same ability to cope with thermal stress is not well understood. In the present study, we determine the transcriptome changes that occur in response to thermal stress in the muscle of two divergent O. oratoria populations. Specifically, in order to avoid population mixture, the mitochondrial DNA Cytochrome oxidase subunit I (COI) and the control region sequences were used to distinguish between the O. oratoria populations. Then, 20 °C and 28 °C were selected as the control and heat stress treatment temperatures, respectively. The aim of this research was to identify the population-specific function genes and relevant pathways in response to thermal stress in O. oratoria. Furthermore, these results can further enhance our understanding of how environmental heterogeneity and climate change have contributed to adaptive diversification within lineages.  The muscle is the largest energy and amino acid pool in the maintenance of homeostasis in ocean animals and because it exhibits a high surface area directly suspended in the water. Therefore, the muscle might be associated with thermal response. However, whether the O. oratoria infraspecific organisms have the same ability to cope with thermal stress is not well understood. In the present study, we determine the transcriptome changes that occur in response to thermal stress in the muscle of two divergent O. oratoria populations. Specifically, in order to avoid population mixture, the mitochondrial DNA Cytochrome oxidase subunit I (COI) and the control region sequences were used to distinguish between the O. oratoria populations. Then, 20 • C and 28 • C were selected as the control and heat stress treatment temperatures, respectively. The aim of this research was to identify the population-specific function genes and relevant pathways in response to thermal stress in O. oratoria. Furthermore, these results can further enhance our understanding of how environmental heterogeneity and climate change have contributed to adaptive diversification within lineages.  30 September 2006). It should be noted that frost anesthesia was made to minimise the suffering of all animals.

O. oratoria Maintenance and Heat Exposure
The O. oratoria were collected from the coastal waters of Qingdao and Zhoushan regions in China, which belong to the Yellow Sea and the East China Sea, respectively ( Figure 1). All O. oratoria were dispatched into separated aquariums with recirculating, aerated seawater (28 salinity) and acclimated over a 24 h period at 20 • C. The heating schemes of the two populations were similar, involving water temperature graded increases from 20 • C to 28 • C, slowly, at a rate of 1 • C per day and acclimated over a 24 h period at 28 • C. Each population was divided into a control group (20 • C) and a treatment group (28 • C). Muscle tissues for transcriptome analyses were collected in four testing groups: Qingdao control group (QD control ), Qingdao treatment group (QD treatment ), Zhoushan control group (ZS control ), and Zhoushan treatment group (ZS treatment ). Three adult females per testing group were euthanised and the muscles were then extracted using sterilized scissors and forceps [25]. Subsequently, 12 individual muscles (2 populations × 2 temperature stages × 3 biological replicates) were separately snap-frozen in liquid nitrogen and stored at −80 • C prior to the following experiments.

Total RNA Extraction and Illumina Sequencing
Total RNA of individual muscles was extracted using a standard Trizol Reagent Kit (Huayueyang Biotech Co. Ltd., Beijing, China), in accordance with the manufacturer's protocol, and quantified using Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). We then purified the mRNA from the total RNA (4 µg) using the RNA Purification Beads (Illumina, San Diego, CA, USA). The remaining RNA was cleaned three times using the Beads Binding Buffer (Illumina, San Diego, CA, USA) and the eluted RNA were incubated at 94 • C for 8 min.

De Novo Assembly, Gene Expression Variation, and Gene Annotation
All raw reads in the FASTQ format were filtered by removing the reads with sequencing adaptors, unknown nucleotides (N ratio > 10%) and low quality (quality scores ≤ 5). We then used the Trinity package (version 2.0.6; [28]) on the de novo assembled all remaining high-quality reads from 12 samples and the redundancy sequences were removed using the Tgicl (Linux x86) software package, additionally splicing the longest unigenes (universal genes) for further analyses. In order to analyze the gene expression variation, we mapped clean reads of individual sample to a multi-fasta file of all unigenes using BWA-mem [29] and the expression level of the unigenes overall was normalised to determine the FPKM (Fragments per kilobase of exon model per million mapped fragments) using RSEM and Bowtie2 at default settings [30]. The number of differentially expressed unigenes between the two populations was quantified by generating Venn diagrams. Furthermore, we performed homology searches to analyze the functional classifications by comparing all unigenes against the NR, NT, Swiss-Prot, KEGG, COG, and GO databases using the Blastx alignment (E-value < 0.00001).

Testing for Population-Specific Genetic Differentiation
We hypothesized that intraspecific mutations in the same or different genes might lead to local adaptive divergences. To test this hypothesis, we first included a population as an independent variable and assessed whether there was genetic differentiation between the two populations. To quantify the differential gene expressions between populations, we treated the FPKM of Qingdao population as a control value and conducted pairwise comparisons with the FPKM of Zhoushan population inhabiting same temperatures (QD control -vs-ZS control and QD treatment -vs-ZS treatment ) using edgeR package with the following parameters: FDR ≤ 0.01 and |log 2 FC| ≥ 2. Furthermore, we assessed which physiological functions might cause genetic differentiation, ultimately leading to adaptability divergence in different O. oratoria populations. To do so, we first calculated and compared the enriched functional categories of each differentially expressed unigene between two experimental pairs. Then, pathway enrichment analyses for each population were performed based on the KEGG database [31]. Finally, we constructed pictures using R software to visualize the variations between the two populations.

Identifying Expression Responses to Thermal Stress
To test for local adaptation, expression responses to thermal stress were identified separately for each population. We included temperature factors as independent variables and then differentially expressed unigenes of two experiment pairs (QO control -vs-QD treatment and ZS control -vs-ZS treatment ) were identified using edgeR package (http://bioconductor.org), and FDR ≤ 0.01 and |log 2 FC| ≥ 2 were used as the filtering thresholds. Venn diagram was applied to quantitatively analyze the number of differentially regulated unigenes of two experiment pairs. Furthermore, we assessed the biological functions of each population to explore the potential functional consequences associated with local temperature adaptation. GO term and pathway enrichment analyses were conducted to evaluate the differential evolutionary adaptation between the two populations. We also constructed pictures using R software (https://CRAN.R-project.org/) to visualize the differential expression responses between them.

Quantitative Reverse Transcription PCR (qRT-PCR) Validation
qRT-PCR was applied to validate the transcriptomic data.
Within each of the categories for up-and down-regulated unigenes, five randomly selected unigenes or contigs were applied to the qRT-PCR analysis and the gene-specific primers were designed using the Primer Premier 5.0 (Table 1).
In addition, the β-actin (Forward, 5 -ATCGTTCGTGACATTAAGGA-3 ; Reverse, 5 -CAAGGAATGAAGGCTGGAA-3 ) and 18S rRNA (Forward, 5 -GAAGGATTGACAGATTGAGAG-3 ; Reverse, 5 -GTAGCGACGGACACATAT-3 ) were chosen as reference genes for internal standardization. Furthermore, standard curves were constructed to identify the ideal dilution times of the cDNA samples and were used as calibrators. A total of 12 cDNA samples were diluted 20-fold using nuclease-free water and were used as templates for PCR. Furthermore, the qRT-PCR analysis was designed following the manufacturer's instructions for the SYBR®Premix Ex TaqTM (Tli RNaseH Plus) RR420A (TaKaRa Biotech Co., Ltd., Dalian, China). A reaction system of 25 µL was amplified using the ABI PRISM 7300 Real-Time PCR System (Applied Biosystems, Thermofisher Scientific, MA, USA). Three parallel experiments for every cDNA template were performed to increase the veracity of the result. After the PCR program, the data were analyzed with ABI7300 SDS software (Applied Biosystems, Thermofisher Scientific, MA, USA). The relative expression levels of all target unigenes, or contigs, were calculated by the log22-∆∆CT analysis method (∆CT = CT target unigene − CT reference gene , ∆∆CT = ∆CT treatment − ∆CT control ).

Phylogenetic Relationship, Illumina Sequencing, and Annotation of the O. oratoria Muscle Transcriptome
Mitochondrial DNA COI and control region sequences were used to identify the phylogenetic relationship of O. oratoria in this experiment. The sequencing results of the COI and control region were shown in Supplementary Files S1 and S2, respectively. The Neighbor-Joining (NJ) trees were constructed using the complete data set of 12 individuals. Results identified that the Qingdao and Zhoushan samples, in the present study, belonged to two distinct lineages ( Figure 2).

Phylogenetic Relationship, Illumina Sequencing, and Annotation of the O. oratoria Muscle Transcriptome
Mitochondrial DNA COI and control region sequences were used to identify the phylogenetic relationship of O. oratoria in this experiment. The sequencing results of the COI and control region were shown in Supplementary Files S1 and S2, respectively. The Neighbor-Joining (NJ) trees were constructed using the complete data set of 12 individuals. Results identified that the Qingdao and Zhoushan samples, in the present study, belonged to two distinct lineages ( Figure 2). We obtained the sequencing information of 12 individuals in the present study and these details are listed in Table 2. The statistics for the de novo assembly showed that the number of merged transcriptomes was 130,102,343 bp, with an N50 of 2760 bp, and all unigenes were composed of 30,539 distinct clusters and 56,816 distinct singletons. The transcriptomic raw reads in this publication are archived on the NCBI Short Read Archive (SRA, SRR7346280; SRR7346279; SRR7346288; SRR7346287; SRR7346282; SRR7346281; SRR7346284; SRR7346283; SRR7346286; SRR7346285; SRR7346278; and SRR7346277) under BioProject PRJNA475657. The assembled and annotated transcript has been deposited at the DDBJ/EMBL/GenBank under the accession number GGQQ00000000. The version described in this paper represents the first version, GGQQ01000000. Furthermore, a total of 60,013 unigenes were annotated based on protein databases using the Blastx alignment. Of all annotated unigenes, 50,668, 39,785, 39,340, 37,790, 22,840, and 13,041 unigenes had significant matched with the sequences in the non-redundant protein sequences (NR), nucleotide sequences (NT), swiss prot protein sequence (Swiss-Prot), kyoto encyclopedia of genes and genomes (KEGG), clusters of orthologous groups of proteins (COG) and gene ontology (GO) databases, respectively.  We obtained the sequencing information of 12 individuals in the present study and these details are listed in Table 2. The statistics for the de novo assembly showed that the number of merged transcriptomes was 130,102,343 bp, with an N50 of 2760 bp, and all unigenes were composed of 30,539 distinct clusters and 56,816 distinct singletons. The transcriptomic raw reads in this publication are archived on the NCBI Short Read Archive (SRA, SRR7346280; SRR7346279; SRR7346288; SRR7346287; SRR7346282; SRR7346281; SRR7346284; SRR7346283; SRR7346286; SRR7346285; SRR7346278; and SRR7346277) under BioProject PRJNA475657. The assembled and annotated transcript has been deposited at the DDBJ/EMBL/GenBank under the accession number GGQQ00000000. The version described in this paper represents the first version, GGQQ01000000. Furthermore, a total of 60,013 unigenes were annotated based on protein databases using the Blastx alignment. Of all annotated unigenes, 50,668, 39,785, 39,340, 37,790, 22,840, and 13,041 unigenes had significant matched with the sequences in the non-redundant protein sequences (NR), nucleotide sequences (NT), swiss prot protein sequence (Swiss-Prot), kyoto encyclopedia of genes and genomes (KEGG), clusters of orthologous groups of proteins (COG) and gene ontology (GO) databases, respectively.

Quantifying Gene Expression Variation
In the present study, the differentially expressed unigenes among four experiment pairs (QD control -vs-ZS control , QD treatment -vs-ZS treatment , QD control -vs-QD treatment , and ZS control -vs-ZS treatment ) were identified by setting the criterions of |Log2FC| ≥ 2 and FDR < 0.05. A Venn diagram was applied to quantitatively reveal the number of differentially expressed unigenes ( Figure 3). The results showed that a substantial overlap of differentially expressed unigenes was observed in different pairs, although just three unigenes (CL15652.Contig1_All, Unigene38717_All, Unigene5742_All) were significantly differentially expressed. According to the annotation information, CL15652.Contig1_All, Unigene38717_All, Unigene5742_All may be the homologous sequence of Craniofacial development protein 2, C-type lectin 2 and beta-lactamase hcpA.

Quantifying Gene Expression Variation
In the present study, the differentially expressed unigenes among four experiment pairs (QDcontrol-vs-ZScontrol, QDtreatment-vs-ZStreatment, QDcontrol-vs-QDtreatment, and ZScontrol-vs-ZStreatment) were identified by setting the criterions of |Log2FC| ≥ 2 and FDR < 0.05. A Venn diagram was applied to quantitatively reveal the number of differentially expressed unigenes ( Figure 3). The results showed that a substantial overlap of differentially expressed unigenes was observed in different pairs, although just three unigenes (CL15652.Contig1_All, Unigene38717_All, Unigene5742_All) were significantly differentially expressed. According to the annotation information, CL15652.Contig1_All, Unigene38717_All, Unigene5742_All may be the homologous sequence of Craniofacial development protein 2, C-type lectin 2 and beta-lactamase hcpA.

Population-Specific Gene Expression
In the present study, we compared the gene expression between populations. The results indicated that the number of differentially expressed unigenes between the two populations studied were highly significant and that this number decreased with the degree of thermal stress increase (11,391 at 20 °C and 9395 at 28 °C ). Furthermore, in comparison with the Qingdao population, a total of 6477 (56.86%) and 6188 (65.86%) upregulated unigenes were obtained from the Zhoushan

Population-Specific Gene Expression
In the present study, we compared the gene expression between populations. The results indicated that the number of differentially expressed unigenes between the two populations studied were highly significant and that this number decreased with the degree of thermal stress increase (11,391 at 20 • C and 9395 at 28 • C). Furthermore, in comparison with the Qingdao population, a total of 6477 (56.86%) and 6188 (65.86%) upregulated unigenes were obtained from the Zhoushan population at 20 • C and 28 • C, respectively. As predicted, we also found that significant differences existed in the biological processes between the two populations (Figure 4). At 20 • C, only oxidoreductase activity (GO:0016491) and acyl-CoA dehydrogenase activity (GO:0003995) were significantly enriched in the functional responses of the two populations (Table 3). In contrast, the amounts of GO terms were much higher at 28 • C, with about 194 significantly enriched terms that appeared in the comparison of the two populations. For brevity, we only discuss the top 10 enriched terms (Table 3).  We recorded the networks of molecular interactions in the cells and the variants specific to particular organisms by comparing the differentially expressed unigenes to the KEGG pathway. At 20 °C , a total of 40 pathways were significantly enriched in the KEGG database (Q ≤ 0.05) and the top 20 statistically significant KEGG (https://www.kegg.jp/) classifications are shown in Figure 5. The results showed some of the digestive-and metabolism-related pathways which were predicted in the KEGG database (ko04972, ko01100, and associated pathways). However, the population-specific expression responses and molecular interactions were lower at 28 °C , with 27 differentially regulated pathways found to be significantly enriched in the two populations (see Figure 5 for the top 20 pathway enrichment statistics). These included the pathways associated with genetic information processing and immunity (ko03010, ko05410, and other associated pathways). For brevity, we only  We recorded the networks of molecular interactions in the cells and the variants specific to particular organisms by comparing the differentially expressed unigenes to the KEGG pathway. At 20 • C, a total of 40 pathways were significantly enriched in the KEGG database (Q ≤ 0.05) and the top 20 statistically significant KEGG (https://www.kegg.jp/) classifications are shown in Figure 5. The results showed some of the digestive-and metabolism-related pathways which were predicted in the KEGG database (ko04972, ko01100, and associated pathways). However, the population-specific expression responses and molecular interactions were lower at 28 • C, with 27 differentially regulated pathways found to be significantly enriched in the two populations (see Figure 5 for the top 20 pathway enrichment statistics). These included the pathways associated with genetic information processing and immunity (ko03010, ko05410, and other associated pathways). For brevity, we only discuss the top 10 statistically significant KEGG classifications of the two populations.

Differentially Expressed Responses of the Two Populations Exposed to the Same Thermal Scheme
To elucidate the pattern of gene expression under thermal stress, we first compared the number of differentially expressed unigenes of each population. The results indicated that the gene expressed differences were closely related to the temperature variation. When the temperature increased from 20 °C to 28 °C , a total of 361 and 1217 differentially expressed unigenes were obtained from the Zhoushan and Qingdao populations, respectively. However, 69.23% (243 up-and 118 downexpressed) and 20.14% (246 up-and 971 down-expressed) up-expressed unigenes were obtained from the Zhoushan and Qingdao populations, respectively. Therefore, we must take into account both the number and the degree of differentially expressed unigenes when conducting transcriptome research.
To identify the functional changes potentially associated with temperature adaptation of the two

Differentially Expressed Responses of the Two Populations Exposed to the Same Thermal Scheme
To elucidate the pattern of gene expression under thermal stress, we first compared the number of differentially expressed unigenes of each population. The results indicated that the gene expressed differences were closely related to the temperature variation. When the temperature increased from 20 • C to 28 • C, a total of 361 and 1217 differentially expressed unigenes were obtained from the Zhoushan and Qingdao populations, respectively. However, 69.23% (243 up-and 118 down-expressed) and 20.14% (246 up-and 971 down-expressed) up-expressed unigenes were obtained from the Zhoushan and Qingdao populations, respectively. Therefore, we must take into account both the number and the degree of differentially expressed unigenes when conducting transcriptome research.
To identify the functional changes potentially associated with temperature adaptation of the two populations, we conducted a GO term enrichment analysis for differentially expressed unigenes of each population when exposed to the same thermal stress. The results showed that 221 GO terms were enriched in the Zhoushan population and about 26.24% of the terms were statistically significant ( Figure 6). For brevity, we only discuss the top 10 enriched terms (Table 4). While the instances of the GO terms were much higher (504 terms), there were only 19 GO terms that were significantly enriched in the Qingdao population ( Figure 6; see Table 4 for the top 10 enriched terms). Differentially expressed unigenes in the Qingdao population were predominately associated with the state or activity of cells. These included the processes of cellular movement, cellular secretion, enzyme production, gene expression and cellular apoptosis (GO: 0031435, GO: 0030968, GO: 0034620, GO: 0035967 and other associated terms). The differentially expressed unigenes in each pair (ZScontrol-VS-ZStreatment, QDcontrol-VS-OQtreatment) were then used to identify the enriched KEGG pathways. There were 154 and 174 pathways found to be differentially regulated between the Zhoushan and Qingdao populations exposed to the same thermal stress. However, the results also indicated that 18 and 4 pathways were significantly enriched in the Zhoushan and Qingdao populations, respectively (Q ≤ 0.05; see Figure 7 for the top 20 pathway  The differentially expressed unigenes in each pair (ZS control -VS-ZS treatment , QD control -VS-OQ treatment ) were then used to identify the enriched KEGG pathways. There were 154 and 174 pathways found to be differentially regulated between the Zhoushan and Qingdao populations exposed to the same thermal stress. However, the results also indicated that 18 and 4 pathways were significantly enriched in the Zhoushan and Qingdao populations, respectively (Q ≤ 0.05; see Figure 7 for the top 20 pathway enrichment statistics). Differentially expressed unigenes in the Zhoushan population were predominately associated with the pathways of the metabolism (ko00520 and ko00531) and the digestive system (ko04972 and ko04974). In contrast, the differentially expressed unigenes in the Qingdao population were enriched for pathways associated with genetic information processing (ko04141), immune (ko05202 and ko05110) and environmental information processing (ko04010).

Validation of the Transcriptome Data by qRT-PCR
We selected five target unigenes, or contigs, to evaluate the transcriptome data of each experimental pair based on the qRT-PCR analysis. For these candidate unigenes, or contigs, the variation trend in expression was concordant between the qRT-PCR data and the transcriptome data, although the values derived from both analytical methods did not perfectly match ( Figure 8). Consequently, the results indicated that the transcriptome data were credible.

Validation of the Transcriptome Data by qRT-PCR
We selected five target unigenes, or contigs, to evaluate the transcriptome data of each experimental pair based on the qRT-PCR analysis. For these candidate unigenes, or contigs, the variation trend in expression was concordant between the qRT-PCR data and the transcriptome data, although the values derived from both analytical methods did not perfectly match ( Figure 8). Consequently, the results indicated that the transcriptome data were credible.

Discussion
Crustaceans are under increasing pressure due to elevated water temperatures [21]. Their ability to cope with ocean warming is critical for our understanding of the probability of crustacean sustainability [32]. Nevertheless, there are variations in the rate of and extent to which different populations of crustaceans respond to ocean warming as a result of their long-term habitat

Discussion
Crustaceans are under increasing pressure due to elevated water temperatures [21]. Their ability to cope with ocean warming is critical for our understanding of the probability of crustacean sustainability [32]. Nevertheless, there are variations in the rate of and extent to which different populations of crustaceans respond to ocean warming as a result of their long-term habitat differences [7]. As a crustacean representative, O. oratoria (De Haan, 1844) spans a wide geographical area and thus different habitats in which temperatures might affect the competition and resource predation regimes and the distribution of the O. oratoria populations [22]. Transcriptome sequencing has been widely employed as an effective and accessible approach for understanding many fundamental evolutionary questions when no genome sequencing data are available [33][34][35]. Consequently, this study examined the thermal response differences of two O. oratoria populations to better investigate the infraspecific temperature adaption differences in O. oratoria. Furthermore, the present study also provided a foundation from which to further our understanding of how different crustacean populations adapt to ocean warming.
It is undeniable that although there is a conservative transcriptome response in all tissues of the organism under thermal stress, actually each tissue will show highly different responses that are associated with their physiological functions. Muscle tissue is the largest energy and amino acid pool in the maintenance of homeostasis, and it best demonstrates the effect of temperature stress on marine organisms. Additionally, the muscle may exhibit remodeling in response to temperature stress, thereby ensuring that the organism adapts to locomotory activity and load. Therefore, muscle tissues were collected for transcriptomic analyses [36]. As predicted, numerous differentially expressed genes were associated with population structures and thermal stress. We also confirmed that several important pathways, associated with metabolic processes, immunity response, genetic information processing processes and digestive processes, were involved in the responses expressed between the two populations. Earlier work has already identified these pathways as critical physiological components and as involved in the responses of the aquatic animals to elevated water temperatures [37,38]. Therefore, such evidence demonstrated that gene expression variation has an interactive effect between the environmental factors and the habitat types. At the same time, the putative functional consequences of expression variation were closely related to different organisational analyses levels (transcript, gene and physiological pathway) [39].

Gene Expression Variation between Geographic Populations
Numerous differentially expressed unigenes were detected between the two studied O. oratoria populations, with the more upregulated unigenes obtained from the Zhoushan population (56.86% at 20 • C and 65.86% at 28 • C). Considering that two populations came from different habitat environments, it was not surprising that the majority of the gene expressed variations were observed between them. Such a phenomenon might demonstrate that adaptive differences to temperature had evolved in the ecologically divergent populations of O. oratoria.
Functional annotations of differentially expressed genes can provide insights about the differences of infraspecific adaptive evolution. At 20 • C, there were only two enriched GO terms that were significantly different between the two O. oratoria populations, which included oxidoreductase activity and acyl-CoA dehydrogenase activity. In addition, there were significant differences in some important digestive-and metabolism-related pathways between the two populations. Previous studies demonstrated that oxidative stress increases when environmental factors move away from the optimum [40,41]. Oxidative stress might reduce lipid beta-oxidation with consequent lipid accumulation, leading to hepatic steatosis [42,43]. Specifically, oxidoreductase and acyl-CoA dehydrogenase activity played important roles in lipid beta-oxidation, potentially shielding the organisms from stress effects. This process might also affect the digestive system and the metabolism differently, ultimately leading to adaptability divergences between organisms. However, the genes involved in the oxidative stress response might over-express towards a threshold when the temperature is far from optimum. Thus, many differential GO terms were observed when the two populations were exposed to 28 • C and the terms were mainly associated with cellular components. Furthermore, there were significant differences in some genetic information and immunity processes between the two populations. Most cellular activities, such as metabolic pathways, cell division and others, always occurred in the cytoplasm. The movement of Ca 2+ in the cytoplasm represented a metabolic process signalling activity, and the ribosomes affected many cellular functions, such as repairing damage or directing chemical processes through a synthesis of proteins [44][45][46]. The significant difference in cellular component terms implied that thermal stress may cause cellular damages in the structure. The cellular stress response is a complex mechanism that aims at preventing the damage of functional proteins. Remarkably, stress responses might compromise the protein function directly through damaged cellular structures and thus present a stress signal for organisms [47,48]. In the present study, the differentially expressed genes might be associated with protein processing, antigen processing, presentation and other cellular activities and aimed at ultimately preventing cellular death. Consequently, we hypothesized that the intraspecific mutations in the same gene or different genes could lead to local adaptive divergences [49] and that the adaptability variation of different populations might correlate with habitat temperatures.
Furthermore, we also discovered that the number of differentially expressed unigenes between the populations decreased when the temperature increased (11,391 at 20 • C and 9395 at 28 • C). This trend confirmed that some functional genes of two populations might over-express towards a threshold due to the increase of the thermal stress degree [50][51][52]. On the other hand, organisms might moderate the gene expressed response over the stress degree and the gene responses of immediate survival are a priority [53].

Differentially Expressed Responses of the Two Populations Exposed to Thermal Stress
When the temperature increased from 20 • C to 28 • C, 1216 functional genes were differentially expressed and a lower proportion of up-regulated genes (20.14%) were found in the Qingdao population. In addition, several GO terms and pathways were significantly enriched in the Qingdao population when the temperature increased. The most significantly enriched biological processes were associated with the state or activity of cells, and the significantly enriched pathways were associated with genetic information, immune and environmental information processing. The most interesting finding in the present study was that the mitogen-activated protein kinase (MAPK) and the unfolded protein response (UPR) played an important role in the response to heat. The UPR presented a cellular stress response related to the endoplasmic reticulum (ER) stress and aimed to restore the normal function of the cell, degrade misfolded proteins and activate the signalling pathways associated with protein folding [47]. As a protective response, UPR might restore cell homeostasis by preventing ER overload with misfolded or damaged proteins. However, Han et al. considered that the functional protein synthesis and adenosine triphosphate (ATP) depletion will increase with the enhancement of stress [54]. In that case, UPR is unable to relieve the ER overload and ultimately leads to cell death [48]. On the other hand, MAPK enzymes phosphorylated and then activated the mitogen-activated protein kinase kinase 1 (MKK1) and the mitogen-activated protein kinase kinase 2 (MKK2). The latter's phosphorylation was important for cellular proliferation, cellular cycle progression, cellular division and differentiation [55]. In addition, some MAPK enzymes might activate both the p38 and c-Jun N-terminal kinase (JNK) pathways [56]. Both the JNK and p38 signalling pathways responded to stress stimuli and were involved in cellular apoptosis or cellular differentiation. Therefore, such results indicated that thermal stress might have a more significant influence on the Qingdao population and damaging the protein function through protein denaturation. This process also represented a stress signal for UPR and MAPK. UPR and MAPK formed a complex regulatory mechanism aimed at prevented cellular damaged or apoptosis. However, we also suspected that prolonged over-expression of UPR and MAPK might ultimately lead to cellular apoptosis with the increase of thermal stress.
Our analyses also indicated that only 351 differentially expressed unigenes were obtained in the Zhoushan population due to rising temperatures, however, the stress response genes were highly upregulated (69.23%) when exposed to heat stress. At 28 • C, the differentially regulated unigenes of the Zhoushan population were primarily involved in the cellular and transcription regulation process, and the most significant pathways were those of the metabolism and the digestive system. It should be noted that the protein kinase-C (PKC) was a critical enzyme involved in the thermal response processes of the Zhoushan population. PKC is a family of protein kinase enzymes that were involved in controlling the function of other proteins [57]. This included the functions associated with receptor desensitisation, modulating membrane structure, regulating transcription, mediating immune responses, and regulating cellular activity. However, the effects of PKC were cell-type-specific [58]. Higher temperatures might lead to metabolic and digestive depressions in O. oratoria. In order to ensure the survival of organisms, PCKs, as functional proteins, were activated to regulate the metabolic and digestive pathways. Ibarz et al. found that prolonged fasting also caused a depletion of glycogen and such loss of energy might also activate the regulatory mechanism of the metabolic and digestive pathways [59].
We tested how thermal stress actually affected the expressed responses of two populations of O. oratoria. As mentioned earlier, the adaptability variation of different populations might be population-specific. Consequently, we first identified the number of differentially expressed unigenes of each population exposed to same thermal stress, respectively. Although most differentially expressed unigenes were obtained from the Qingdao population with the temperature increase, the stress response genes were highly upregulated in the Zhoushan population. Narum and Campbell also investigated the adaptive response of the functional genes in redband trout and emphasized that a lower gene variation was observed in warmth adapted populations under conditions of thermal stress [60]. Therefore, we suspected that the Qingdao population might sustain enormous thermal stress and that it should mobilize more specific genes to resist the rising of temperatures. The functional annotations of the differentially expressed genes provided insights about potential regulatory and physiological mechanisms for mediating the adaptation to thermal stress. Despite the fact that the number of differentially expressed genes was lower than in the Qingdao population, it is important to note that there were more significantly enriched GO terms and pathways in the Zhoushan population. Such results might confirm that several functional genes with large pleiotropic effects are applied in response to thermal stress in the Zhoushan population. In addition, the results of the GO terms and pathways also indicated that the degree of organism damage in the Qingdao population may be more serious. Thus, we considered that different populations can embark on different solutions to cope with heat stressors [61,62] and that a thermal adaptation of the functional genes had evolved in the Zhoushan population. In short, the higher annual temperature of the East China Sea might have generated more thermal adaptive evolutionary responses in the Zhoushan population.

Conclusions
Overall, our analyses have confirmed that the infraspecific gene expression variations of the O. oratoria were closely related with specific populations and their environmental factors. The results showed a gene expressed variation in a population-specific pattern, which indicates that the local environment could lead to evolutionary changes in gene regulation, ultimately leading to adaptive divergences. Therefore, we considered that different populations might use different solutions to cope with environmental stressors. Higher temperatures apparently generated a series of stress reactions in O. oratoria, such as oxidative stress, metabolic disorder, and cell apoptosis. We identified and compared the potential regulatory variations between two divergent populations under the same thermal stress conditions, a few genes with large pleiotropic effects in Zhoushan population. The discrepancy of regulated degree also indicated that thermal stress had a much stronger influence on the survival of the Qingdao population. In addition, this study revealed adaptive patterns in different organization levels (from transcript to gene to physiological pathway) of the two populations, and we suggest that the regulation mechanisms of the Zhoushan population was more efficient than that of the Qingdao population under thermal stress conditions. Finally, the present study provided some novel insights into the discrepancy between temperature adaptive differences of two O. oratoria populations. The results can also serve as a foundation for further studies aiming to explain the genetic basis of local adaptation in O. oratoria and other crustaceans. Accordingly, future studies will need to identify the potential functional consequences of these transcriptome variations and to test how functional genes affect the distribution of crustaceans. In addition, multiple environmental factors might cause synergistic and generate different local adaption for crustaceans. From a transcriptome perspective, future research needs to characterise transcriptome variation among crustacean populations under multiple environmental stressors.