Comparative transcriptional profiling analysis of the effect of heat waves during embryo incubation on the hatchlings of the Chinese soft‐shelled turtle (Pelodiscus sinensis)

Abstract Temperature is one of most the important environmental factors that affect the ontogenesis of organisms. In this study, we incubated Chinese soft‐shelled turtle eggs at 28°C (control temperature, C treatment), a temperature with a 16°C cold shock and a 36°C heat shock twice per week (S treatment) or a ramp‐programmed temperature of 29 ± 9°C (with 12 hr (+) and 12 hr (−) every day) (F treatment). The incubation period, hatching success, hatchling weight, and locomotor performance were significantly different between the controls and the different heat treatment groups. The pathogen challenge results illustrated that hatchlings from the S treatment group were more resistant to bacterial infection, whereas hatchlings from the F treatment group were more vulnerable. We used RNA‐seq quantification analysis to identify differentially expressed genes (DEGs) of hatchlings in the S treatment group. Based on the functional annotation results for the DEGs, 9 genes were chosen to verify the RNA‐seq results. The background expression of DEGs was also analyzed for the three treatments, as was the regulation of the pathogen challenge. The results showed that 8 DEGs were related to the immune response after pathogen challenge and that temperature was an important factor in differential regulation of the immunity pathways.

Animals can choose the most suitable temperature for their development, and temperature can modify their performance and life history traits as well (Angilletta, Sears, & Pringle, 2009;Parker & Andrews, 2007). Unpredictable extreme temperatures mean that organisms must respond rapidly for protection. The process is energetically costly, with the impairment of enzymatic function, membrane structure, and oxygen supply (Leicht, Jokela, & Seppala, 2013). Research based on experiments that were designed with simulated heat waves demonstrated that the survival rate, morphology, growth, and reproduction could be affected, and the levels of defense traits were ultimately reduced (Bauchinger, McWilliams, & Pinshow, 2011;Leicht et al., 2013). In Drosophila, genetic perturbation occurred during the spring 2011 heat wave in Western Europe, during which the genetic constitution of the population transiently shifted to summer-like frequencies (Rodriguez-Trelles, Tarrio, & Santos, 2013).
In ectotherms, the development and efficiency of the immune system, which is susceptible to the external temperature and is one of the key systems for survival, were markedly affected by temperature changes (Rijkers, Frederix-Wolters, & Van Muiswinkel, 1980;Truscott & White, 1990). The expression of immune parameters, which are related to innate and acquired immune functions, is sensitive and responds quickly along with changes in phenotypes. Genes related to the immune system, such as heat-shock proteins (Hsps), superoxide dismutase (SOD), and tumor necrosis factor (TNF), have been selected for gene expression profiling (Dittmar, Janssen, Kuske, Kurtz, & Scharsack, 2013;Zhang et al., 2014). Despite the growing appreciation of the effects of heat waves on organisms, research about ectotherms is rare. Recently reported results about ectotherms have mostly addressed the effects of heat waves on the phenotypes of animals, whereas molecular mechanisms and gene adaptations were rarely considered.
In this study, we used the Chinese soft-shelled turtle, Pelodiscus sinensis, which is genetically sex determined (GSD) (Ji, Chen, Du, & Chen, 2003). Therefore, we designed this experiment without regard to sex. Most of the research about Chinese soft-shelled turtle embryo incubation has been conducted at a constant temperature in the laboratory. The constant incubation temperature range of the Chinese soft-shelled turtle is 24-34°C, and temperature extremes will cause stress, which would lead to deformities and to disruptions of embryonic development . In the wild, ectotherms experience temperature fluctuation, but there are few studies addressing the effects of extreme temperature changes. In this study, we used two experimental incubation temperatures with fluctuations to mimic the natural extreme temperatures with climate change, and the lowest and highest temperatures were outside the range of constant temperatures used for incubation. Morphological parameters of hatchlings were recorded and analyzed. We used a pathogen isolated from a diseased Chinese soft-shelled turtle to challenge the hatchlings and evaluated the immunity of hatchlings from different treatments. The effects of temperature on phenotypes of ectotherms have been widely reported, but the underlying molecular mechanisms remain unclear. In this study, we intended to identify the genes responsive to heat stress, especially genes related to immunity and other pathways. RNA-seq quantification analysis is one of the RNA deep-sequencing technologies for genome-wide transcriptome analysis (Wang, Gerstein, & Snyder, 2009). We used RNA-seq quantification analysis to determine the complexity of gene expression regulation affected by heat waves. The gene data mining is the first step to illustrate the thermal adaptation of ectotherms and is also important for identifying the problematic effects of global climate change.

| Study species and pathogen strain
The Chinese soft-shelled turtle is one of the main aquaculture species in China and is distributed across China and Southeastern Asia. We collected 633 freshly laid fertilized eggs, which were identified by a white patch on the shell surface, from a farm located in Hangzhou, Zhejiang, in July 2014. The mean clutch size of the female turtles is 18, and the mean egg mass was determined to be 4.19 ± 0.58 g using an electronic balance (Mettler Toledo AB135-S, Germany). The pathogen used with the Chinese soft-shelled turtles was Aeromonas hydrophila TL1, which was isolated from the liver of a diseased turtle. A. hydrophila TL1 was cultured in Luria-Bertani (LB) liquid medium at 28°C.

| Egg incubation and hatching rate
The eggs were placed in plastic containers filled with moist vermiculite at −220 kPa (1 g water/1 g vermiculite). Each container contained eight eggs. There were three types of temperature settings for incubators (Ningbo Life Science and Technology Ltd., China). The control temperature was 28°C (the C treatment), which is the average land surface temperature from May to July to in Hangzhou over 10 years and the most suitable temperature for incubation of Chinese soft-shell turtles. The first experimental temperature condition was a temperature with a 16°C cold shock and a 36°C heat shock twice per week (the S treatment, Figure 1a), and the second was a ramp-programmed temperature of 29 ± 9°C, with 12 hr above 29°C (+) and 12 hr below 29°C (−) every day (the F treatment, Figure 1b). The eggs were randomly distributed into three incubators to minimize maternal effects. The containers in incubators were examined twice per week to promptly remove the eggs in which embryogenesis had terminated and to keep the vermiculite at the appropriate moisture level. The containers were moved to different shelves according to a predetermined schedule to minimize any effects of thermal gradients inside the incubators.

| Incubation periods and hatching success
Toward the end of incubation, the containers were monitored once per day for newly emerging hatchlings. The incubation periods recorded for the day lasted from the beginning of incubation to the emergence of the hatchlings, and the hatching rates are expressed as the percentage of eggs that hatched live turtles at the end of the experiment. After emergence, the hatchlings were maintained in the containers until the yolk was entirely absorbed, after which they were kept in a temperaturecontrolled room at 28 ± 1°C with a 12 hr light/12 hr dark cycle.

| Hatchling traits and locomotor performance
After the yolk had been entirely absorbed, the weights of hatchlings were measured with an electronic balance (Mettler Toledo AB135-S, Germany) to the nearest 1 mg. The lengths, widths, and heights of the turtle shells were recorded with digital Vernier calipers (Mitutoyo).
Then, each hatchling was chased along a 1.2-m racetrack filled with water at a depth of 40 mm and at a constant temperature of 28°C.
Swimming performance was recorded with a Panasonic video camera.
Videotapes were later analyzed for sprint speed in the fastest 30-mm interval.

| Experimental infection and virulence analysis
We used the median lethal dose (LD50) for hatchlings after pathogen infection to determine the immunocompetence of the turtles from different incubation temperatures. Briefly, the hatchlings from each incubation temperature were divided into four groups, each containing 20 hatchlings. The four groups were separately challenged with 5 × 10 3 colony-forming units (CFU), 5 × 10 4 CFU, 5 × 10 5 CFU, or 5 × 10 6 CFU of TL1. The turtles were monitored for 10 days, and the mortalities were recorded every day.

| Transcriptome data processing and analysis of differential gene expression
The original image data collected from the Ion Proton platform were converted into sequence data. The adaptor sequences and low-quality reads were filtered to obtain high-quality clean reads. After data quality statistics, clean reads were mapped against the genome and genes of P. sinensis, which are available at http://www.ncbi.nlm.nih.gov, by the Torrent Mapping Alignment Program (TMAP). The gene expression levels were quantified with the Sailfish software package. The RPKM method, which calculates the reads per kilobase per million mapped reads (RPKM = total exon reads/mapped reads in million X exon length in kb), was used for calculating the expression level. The Noiseq package method (Tarazona, García-Alcalde, Dopazo, Ferrer, & Conesa, 2011) was used to screen differentially expressed genes F I G U R E 1 The two heat wave treatments of the Chinese softshelled turtle eggs. (a) The S treatment, an incubation temperature with a 16°C cold shock and a 36°C heat shock twice per week; (b) The F treatment, a ramp-programmed temperature of 29 ± 9°C with 12 hr (+) and 12 hr (−) every day (DEGs) between two groups. The data were log2 transformed to identify the genes with significant differential expression during heat stress [log2-fold change (FC) ≥1 or log2 FC ≤−1 and p value < .05]. Gene ontology analysis was used to identify the main biological functions performed by the DEGs. The calculated p-value underwent a Bonferroni Correction (Abdi, 2007), and a corrected p-value ≤ .05 was used as a threshold. Reads showing significant homology against Chinese softshelled turtle genes were mapped onto metabolic pathways using the KEGG database (http://www.genome.jp/kegg/pathway.html).

| Changes in DEG expression after heat treatment and bacterial infection
According to the analysis using the KEGG database, we designed primers for 9 DEGs; these primers are shown in Table S1. At the same time, we also examined the changes in expression of these genes at 12 hr and 48 hr after bacterial infection. Hatchlings (5/group) from the three treatments without bacterial challenge and from the 5 × 10 5 CFU group at 12 hr and 48 hr postchallenge were sacrificed. Livers were removed and frozen in liquid nitrogen. Total RNA was extracted using TRIzol (Invitrogen, CA, USA), and the first-strand cDNA was obtained using the PrimeScript ™ RT Reagent Kit with gDNA Eraser (Perfect Real Time) (Takara, Dalian, China). The designed primers were confirmed to be specific through PCR and sequencing. Then, quantitative realtime PCR was conducted according to a previously described protocol (Dang, Lu, et al., 2015;Dang, Zhang, & Du, 2015). The specific qRT-PCR products were verified by electrophoresis.

| Data analysis
The software packages Statistica 6.0 (StatSoft, Tulsa, USA) and SPSS 13.0 were used for data analysis. A Kolmogorov-Smirnov test was used for normality of distributions, and Bartlett's test was used for homogeneity of variances. Analysis of variance (ANOVA) was used to determine the influence of incubation temperature on the incubation period, and analysis of covariance (ANCOVA) was used to determine the influence of incubation temperature on hatchling mass, with the initial egg mass as the covariate. When the hatchlings were exposed to bacterial infections, a G test was used to analyze the among-treatment differences in mortality, and a probit model was used to determine the median lethal dose (LD50). Values are presented as the mean ± SE and range, and the significance level was set at p = .05 throughout this study.

| Incubation period and hatching rate
There were significant differences among the treatments for incubation period (F 2, 433 = 3100.2, p < .0001). Hatchlings from the S treatment possessed the longest incubation period and highest hatching success (76.78%). Both the incubation period and hatching success (71.09%) of the C treatment hatchlings were intermediate. Hatchlings from the F treatment possessed the shortest incubation period and lowest hatching rate (58.77%). The specific data are all shown in Table 1.

| Hatchling traits and locomotor performance
The weights of hatchlings from the F treatment, which were the lightest, were significantly different from those in the other two treatments (F 2, 433 = 11.77, p < .001). There were no significant differences among treatments in terms of shell length (F 2, 433 = 0.55, p = .575), shell width (F 2, 433 = 1.71, p = .182), and body height (F 2, 433 = 2.05, p = .130). The hatchlings from the F treatment showed the fastest swimming speed, which was significantly faster than in the other two treatments (F 2, 55 = 8.88, p < .001), whereas there were no significant differences between hatchlings from the C and S treatments. The specific data are all shown in Table 2.

| RNA-seq mapping and analysis of differentially expressed genes
The data analyzed in this study were submitted to the NCBI Sequence Read Archive under accession no. SRP104379. We sequenced three samples from each group, and for each sample, we collected more than 12 million clean reads. Clean reads were mapped to a reference genome and to reference genes at mapping rates that were higher than 97% and 58%, respectively. There were more than 18,000 expressed genes for each sample. The specific number of clean reads yielded from each sample, the genome and gene mapping rates, and the expressed genes are shown in Table S2. As there were three replicates per group, the RPKM values of each gene in each replicate showed slight differences. Based on the different DEG screening methods, we obtained various statistics about the DEGs, which are shown in Table 3. In this study, we obtained 84 genes that were differentially expressed between two groups. These significantly upregulated or downregulated genes are listed in Table S3.

| Functional annotation of differentially expressed genes
The 84 DEGs were mapped to GO and KEGG ontology (KO) terms to assign gene functions. Among the component ontology terms, most DEGs were associated with cells, cell parts, and organelles. Among the function ontology terms, most DEGs were related to binding and catalytic activity. Among the process ontology terms, most DEGs were involved in muscle-associated processes and metabolism (Table   S4). Seventy-three DEGs with pathway annotations were identified.
Because there were no available previous reports of pathways for the Chinese soft-shelled turtle or other reptiles, the annotations mostly used human pathways as references. Most of the enriched pathways were involved in disease or infection (39.73%), signaling pathways (13.70%), and metabolism (10.96%) ( Table S5). Most of the DEGs involved in signaling pathways were related to disease or infection.

| Effects of heat treatment and bacterial challenge
We chose 9 DEGs to verify the results of the RNA-seq and analyze the gene expression among the three temperature treatments by qRT-PCR. The results showed that the expression patterns of 6 of the 9 DEGs were identical to the patterns determined by RNA-seq ( Figure 3). For the same 9 DEGs, we also analyzed the expression patterns in the S treatment, and we found that the expression patterns were not simply opposite between the F treatment and the S treatment and that some DEGs showed similar expression patterns in the two groups ( Figure 3). As we could not simply determine whether these 9 DEGs were related to the immunity of hatchlings, we examined the changes in expression of these 9 DEGs after bacterial challenge in the different treatment groups. The expression patterns of each gene were different in different treatment groups. After bacterial challenge, BCL1 was observed to be upregulated in the C and F treatments but not in the S treatment. Also after bacterial challenge, NFATF5 and VEGF were found to be downregulated in the C and S treatments but not in the F treatment. In the F and S treatments, we found that three genes, RanGAP, DLG1, and hnRNPK, showed the same expression pattern, namely, upregulation after bacterial challenge, whereas these genes were not differentially expressed in the C treatment (Figure 4).

| Temperature and parameters
Temperature plays an important role in the life histories of organ-  (Bowden et al., 2014). The experimental temperatures for the red-eared slider turtle were within the optimal developmental range. In this study, we mimicked the natural temperature fluctuation and designed two different heat waves to illustrate the effects of temperature, and the lowest and highest temperatures of these two fluctuating treatments were outside the range of 24-34°C , which is the constant incubation temperature range for the Chinese soft-shelled turtle. The mean temperatures of the two heat wave treatments were 26°C (S treatment) and 29°C (F treatment). Although the temperatures in the S treatment were outside the optimal development range for Chinese soft-shelled turtles, the hatchlings in the S treatment showed a longer incubation period and higher hatching success. We cannot currently determine whether the mean temperature or the fluctuation frequency induced these results.
In reptiles, incubation temperature can have many effects on posthatching traits, including fitness. Immunity, providing a generalized defense against infection, is one important component of fitness in vertebrates. In Chinese soft-shelled turtle and map turtle hatchlings, stronger immune responses were shown after incubation under a lower constant temperature (24°C) (Dang, Lu, et al., 2015;Dang, Zhang, et al., 2015;Freedberg et al., 2008). In painted turtle, hatchlings from the fluctuating treatments mounted a greater swelling response to an intraperitoneal priming injection of phytohemagglutinin than did hatchlings from the constant treatments (Les et al., 2009). In this study, the bacterial challenge results showed that fluctuations in incubation temperature definitely affected the turtle mortality after bacterial challenge. As described above, the mean temperature and the magnitude and frequency of the fluctuations in the incubation temperature should be all important considerations. In the future, we suggest refining the temperature design to test the effects of these factors separately. Here, we focused on the underlying molecular mechanism associated with immunity and temperature, and the DEGs were examined first. Research about SERCA1 in humans has shown that this protein is one of the most important contributors to removing calcium from the cytoplasm in skeletal muscle (Guglielmi et al., 2016). Defective function of SERCA1 protein was associated with a type of myopathy and with loss of muscle mass (Mázala et al., 2015). The results as verified F I G U R E 3 The expression of the 9 selected DEGs in liver of Chinese softshelled turtle hatchlings from the three temperature treatments. The expression levels of 9 DEGs in liver were determined by quantitative real-time reversetranscription PCR and normalized to those of β-actin mRNA. Levels in the C treatment were set as 1. Vertical bars indicate the means ± SE (N = 5) by qRT-PCR were not consistent with the RNA-seq results. The expression of SERCA1 was slightly but nonsignificantly upregulated in the different temperature treatments and after pathogen challenge. In this study, we speculated that SERCA1 was not related to long-term temperature treatment and immunity in the Chinese soft-shelled turtle. The RNA-seq analysis also identified other DEGs, such as Myosins, that are involved in muscle development, but these genes were not verified here. A new question is whether the heat stress affected muscle development, which would merit further investigation. Significant differences between hatchlings without pathogen challenge and hatchlings at 12 hr and 48 hr postchallenge are indicated with asterisks. *p < .05 rate of GTP hydrolysis (Bernards, 2003). Rho GTPases are found in various kinds of cells; work as important intracellular signaling switches; are involved in actin polymerization, focal adhesion complex assembly and function, cell cycle progression, membrane trafficking, cell adhesion, and cell polarity; and are reported to be associated with diseases such as cancer and neurological and cardiovascular disorders (McKerracher, Ferraro, & Fournier, 2012;Peck, Douglas, Wu, & Burbelo, 2002). As negative regulators of Rho GTPases, knockdown or overexpression of Rho GAPs was reported to result in morphological changes in tumor cells (Saito, Ozawa, Hibino, & Ohta, 2012). Ran GTPases participate in various nuclear events, and mutation of Ran GAP caused defects in RNA processing and nucleocytoplasmic mRNA transport in yeast (Bischoff, Krebber, Kempf, Hermes, & Ponstingl, 1995). Upregulated expression of Rho GAP and Rho GAP after pathogen challenge was only found in the heat treatment group, not in the control group. These results may be related to the higher background expression of these two genes in the control group.

| DEGs associated with heat and bacterial challenge
hnRNP-K is a component of the ribonucleoprotein particles, serves as a nucleic acid-linking docking platform, and is involved in the regulation of chromatin structure, transcription, pre-mRNA processing, splicing, mature mRNA transport to the cytoplasm, translation, nuclear transport, signal transduction, and DNA repair (Gao et al., 2013). hn-RNP-K was found to bind to elongation factor 1α and to the promoter of VEGF (Bomsztyk, Van Seuningen, Suzuki, Denisenko, & Ostrowski, 1997;Uribe, Guo, Shin, & Sun, 2011). VEGF-A was originally discovered to function in vasculogenesis and angiogenesis and participates in preventing apoptosis in neurons and in promoting their proliferation and survival (Zhang, Bao, Hambly, & Gillies, 2009). More recently, VEGF-A has been found to be produced in response to stressors, such as hypoxia, glucose deprivation, and extracellular pH (Latham, Molina-Paris, Homer-Vanniasinkam, & Ponnambalam, 2010). VEGF-A was found to be associated with disease development and has even been proposed to function as a tumor biomarker (Latham et al., 2010). After the pathogen challenge, we found similar expression patterns for hn-

RNP-K and VEGF-A.
BCL6 is a transcriptional repressor participating in the germinal center (GC) reaction of mature B cells and CD4 + T-follicular helper cells (Basso & Dalla-Favera, 2012;Bunting & Melnick, 2013). In the presence of pathogens, the expression of BCL6 was significantly upregulated, except in the S treatment, in which the background expression level of BCL6 was significantly higher than in the other two treatments. NFAT proteins are nuclear transcription factors that were first found in activated T cells, and they regulate various cytokines involved in the immune system (De Boer, Mordvinov, Thomas, & Sanderson, 1999). Unlike the other four members of the NFAT family, which are regulated by calcium-calcineurin, NFAT5 is found in all cell types and is regulated by osmotic stress (Macian, 2005).
DLG1 regulates Ag-dependent cytoskeletal and signaling events by localizing to the junction between T cells and antigen-presenting cells, also known as the immunological synapse (Crocetti, Silva, Humphries, Tibbs, & Miceli, 2014;Zanin-Zhorov et al., 2012). Tumor necrosis factor receptor-associated factor 3 (TRAF3) is an adaptor molecule that links upstream receptor signals to downstream gene activation in both adaptive immune cells and the innate immune system (Oganesyan et al., 2006;Yi, Lin, Stunz, & Bishop, 2014).
TRAF3IP1 interacts with TRAF3 and is also important in the regulation of the immune response (Cheng et al., 2016). TRAF3IP1 was also found to bind both actin and tubulin to regulate cytoskeletal dynamics (Berbari et al., 2011). In this study, we found that NFAT5 was downregulated both in the heat treatment and after pathogen challenge, but the other genes related to immunity cells were upregulated after pathogen challenge in the groups subjected to heat treatment.
These 9 selected DEGs did not all show significant changes in expression in either temperature treatment. Temperature must be an important factor affecting gene expression pathways. The environmental temperature affects the entire physiology of ectotherms, and genes are the fundamental inherited factors that determine animal development. As the global environment deteriorates, population decreases or even the extinction of these animals seems likely. We attempted to identify the underlying effect of the environmental change in order to protect these animals.
Heat-shock proteins are also associated with changes in temperature. In loggerhead sea turtle, the expression of genes in the heat-shock protein family was upregulated after heat-shock treatment (Bentley, Haas, Tedeschi, & Berry, 2017), and in invertebrates, heat-shock proteins were found to be differentially expressed after heat treatment (DeSalvo et al., 2008;Gleason & Burton, 2015;Huang et al., 2017). However, we did not observe differential expression of heat-shock proteins in Chinese soft-shelled turtles. In our studies, the embryos completed development during the heat treatment, and we later analyzed the difference in just one organ.
In our previous studies about the Chinese soft-shelled turtle, the different constant incubation temperatures did not affect the expression of Hsp70 in liver of the hatchlings (Dang, Lu, et al., 2015;Dang, Zhang, et al., 2015). Further, the heat treatments in the other studies were short. We suggest that the embryos of Chinese softshelled turtles acclimated to the incubation temperatures, and the expression of heat-shock proteins returned to normal levels because constant high expression levels would be at the cost of other protein functions.
In conclusion, we incubated the Chinese soft-shelled turtle eggs in two heat wave conditions and used 28°C as the control temperature.
The incubation period, hatching success, hatchling weight, and locomotor performance were significantly different between the control group and the different heat wave treatment groups. The bacterial challenge results demonstrated that hatchlings from the group treated with two heat shocks during incubation (S treatment) were more resistant to bacterial infection. We used RNA-seq quantification analysis to identify differentially expressed genes in hatchlings in the S treatment group. The functional annotation of DEGs with the KEGG ontology database implied that they include genes related to immunity and metabolism. The genes reported in this study will provide a fundamental basis for discovering novel genes in the responses of Chinese softshelled turtles to different temperatures.