A transcriptional and functional analysis of heat hardening in two invasive fruit fly species, Bactrocera dorsalis and Bactrocera correcta

Abstract Many insects have the capacity to increase their resistance to high temperatures by undergoing heat hardening at nonlethal temperatures. Although this response is well established, its molecular underpinnings have only been investigated in a few species where it seems to relate at least partly to the expression of heat shock protein (Hsp) genes. Here, we studied the mechanism of hardening and associated transcription responses in larvae of two invasive fruit fly species in China, Bactrocera dorsalis and Bactrocera correcta. Both species showed hardening which increased resistance to 45°C, although the more widespread B. dorsalis hardened better at higher temperatures compared to B. correcta which hardened better at lower temperatures. Transcriptional analyses highlighted expression changes in a number of genes representing different biochemical pathways, but these changes and pathways were inconsistent between the two species. Overall B. dorsalis showed expression changes in more genes than B. correcta. Hsp genes tended to be upregulated at a hardening temperature of 38°C in both species, while at 35°C many Hsp genes tended to be upregulated in B. correcta but not B. dorsalis. One candidate gene (the small heat shock protein gene, Hsp23) with a particularly high level of upregulation was investigated functionally using RNA interference (RNAi). We found that RNAi may be more efficient in B. dorsalis, in which suppression of the expression of this gene removed the hardening response, whereas in B. correcta RNAi did not decrease the hardening response. The different patterns of gene expression in these two species at the two hardening temperatures highlight the diverse mechanisms underlying hardening even in closely related species. These results may provide target genes for future control efforts against such pest species.


| INTRODUC TI ON
During the invasive and adaptive process, species often encounter novel environmental conditions that require adaptation through genetically based evolutionary changes, phenotypic plasticity or a combination of these processes (Gibert et al., 2016;Vázquez, Gianoli, Morris, & Bozinovic, 2017). Phenotypic plasticity or evolutionary changes can help buffer organisms from environmental changes and thereby help their establishment and expansion, and even be the target of selection (Wellband & Heath, 2017). Many studies that consider the ability of plastic changes to buffer environmental effects consider temperature extremes, which play an important role in the success of the invasive process (David et al., 1997;Delpuech et al., 1995;Klepsatel et al., 2013).
Bactrocera dorsalis (Hendel) and B. correcta (Bezzi) (Diptera: Tephritidae) are invasive pests that damage fruits and vegetables (Permpoon, Aketarawong, & Thanaphum, 2011). These species have been investigated because of their increasingly wide distributions and repeated invasions. The Bactrocera genus has a competitive advantage for oviposition over other fruit flies such as Ceratitis capitata which is one of the most devastating and invasive worldwide pests (Liu, Zhang, Hou, Ou-Yang, & Ma, 2017;Malacrida et al., 2007). The two Bactrocera species are similar in many biological attributes such as mating duration, oviposition preference, and offspring performance as well as sharing a common origin. However, the geographical distribution of B. dorsalis is now much wider than that of B. correcta both in China and elsewhere in the invasive range, which B. dorsalis has invaded including Hawaii, Kenya, and Tahiti and gradually displaced pre-established Ceratitis species in recent years ( Figure 1; Hu, Chen, & Li, 2014;Liu et al., 2014;Liu, Jin, & Ye, 2013;Lux, Copeland, White, Manrakhan, & Billah, 2003;Reitz & Trumble, 2002). Thus far, there are no records of B. dorsalis being displaced by other tephritid fly species .
However, mechanisms that underpin heat hardening in Bactrocera species and how they might contribute to differences among the species are poorly characterized. As in other insects, there are limited data on how species might differ in gene transcription under hardening and how any differences relate to temperature adaptation and organism performance (Clarke, 2003). In this paper, we considered whether differences in temperature adaptability associated with hardening are linked to the regulation of certain genes or pathways, which in turn has contributed to the different distribution and F I G U R E 1 Global distribution map of Bactrocera correcta and B. dorsalis. The global distribution of the two species based on the information from GBIF (https://www.gbif.org/) and CABI (https:// www.cabi.org/). The blue dots and areas represent the distribution of B. correcta, and the red dots and areas represent the distribution of B. dorsalis invasive potential of the two species. We undertook a transcriptional analysis to identify mechanisms underlying differences of heat hardening responses in B. dorsalis and B. correcta. We first characterized the survival of these species under heat stress after a series of hardening treatments. Then, we undertook a transcriptional analysis to identify key pathways and genes involved in hardening. Lastly, quantitative real-time PCR (qRT-PCR) and RNA interference (RNAi) were used in a functional analysis of these two species around one gene that seemed particularly important in the hardening response and contributed to the different hardening response of these species.

| Samples
Bactrocera dorsalis and B. correcta were sourced from their first invaded range in China, from Guangdong province (N 23.40, E 113.22) for B. dorsalis with the annual mean temperature 21.7°C and from Yunnan province (N 23.60, E 102) for B. correcta with the annual mean temperature 25.8°C (Li, Wu, Chen, Wu, & Li, 2012;Liu & Ye, 2009). Cultures of these species were maintained in the laboratory for approximately 10 generations, across 3 years, at a temperature of 25°C, a humidity of 70% and light period of 10D:14L in an environment-controlled incubator (Guo, Zhao, Liu, Li, & Shen, 2018;Nyamukondiwa, Terblanche, Marshall, & Sinclair, 2011;Weldon, Nyamukondiwa, Karsten, Chown, & Terblanche, 2018). Cultures were maintained by turning over flies across three cages (each 45 cm × 45 cm × 50 cm), with 200 individuals per cage. To reduce the risk of commonly observed inbreeding in fruit flies, we also regularly initiated new cultures of the species from the same provinces as above and added 30 new field fruit flies per cage from these cultures every half year (Hoffmann & Ross, 2018). Adult flies were given 25% sucrose and 75% peptone as their diet and immature stages were cultured on an artificial diet described by Yuan et al. (2006).
The most temperature-sensitive development stage (7-day-old 3rd early-instar larvae) of these two species was chosen for hardening experiments (Hu et al., 2014;Jang, 1991). The 3-day-old 1st instar larvae for hardening response study and 4-day-old 2nd instar larvae for heat tolerance study of B. dorsalis and B. correcta were used in dsRNA-feeding experiments.

| Survival analysis under hardening
Five replicates with 30 larvae of each species, B. dorsalis and B. correcta, were exposed to heat hardening temperatures ranging from 34 to 40°C as described in Hu et al. (2014) for 4 hr in a circulating water bath (PolyScience Programmable Temperature Controller, Total Temperature Instrumentation, Inc., USA), while the control group was maintained at 25°C. After pretreatment, all tubes were returned to 25°C for 1 hr before exposure to heat stress at 45°C for 1 hr in the water bath, which allowed various heat shock proteins and other genes to develop (Hoffmann, Sørensen, & Loeschcke, 2003;Sørensen & Loeschcke, 2001). We chose 45°C for 1 hr as an extreme thermal stress which caused mortality in the range 40%-50% without heat hardening (Hu et al., 2014;Nyamukondiwa et al., 2011). After exposure to the heat stress, the larvae were returned to 25°C and their survival rate was scored after 4 hr (Supporting Information Figure S1A). During the thermal treatments, samples were enclosed in 2 ml tubes with a hole in the lid, topped with 4 g diet to ensure the food supply and humidity. Tubes were suspended in a circulating water bath set to the desired temperature. The control 25°C flies were also handled in the same way as the flies in the hardening experiments.
Each replicate had 30 larvae in the 2 ml tube and each hardening treatment involved five biological replicates.

| RNA extraction and cDNA synthesis
Individuals for mRNA sequencing, transcriptome verification and the detection of RNA interference efficiency in each species were collected randomly for RNA extraction. RNA was extracted from whole body of third instar larvae using an RNA simple Total RNA kit (Tiangen, China). cDNA was synthesized from 1,000 ng total RNA using PrimeScript™ RT reagent kit with gDNA Eraser (Perfect Real Time; Takara, Japan) following the manufacturer's instructions.

| Transcriptome and transcriptome analysis of heat hardening
For mRNA sequencing, total RNA (5 μg) was isolated from three groups of 30 3rd early-instar larvae after exposing them to specific temperatures (25, 35 and 38°C) for 4 hr followed by 25°C for 1 hr for both B. dorsalis and B. correcta (Supporting Information Figure   S1B) using an RNA simple Total RNA kit as described in Section 2.3.
The quantity and integrity of the RNA were assessed with a Qubit ® RNA Assay kit in Qubit ® 2.0 Fluorometer (Life Technologies, CA, USA) and an RNA Nano 6000 Assay kit of the Agilent Bioanalyzer To remove adapter contamination, low-quality bases and bases artificially introduced during library construction, we trimmed all raw reads using Trimmomatic 0.32 (http://www.usadellab.org/ cms/index.php?page=trimmomatic) with the following parameters: Phred33 LEADING: 3 TRAILING: 3 SLIDING WINDOW: 1:10 MINLEN: 75 before transcript assembly, while the unpaired reads were discarded (Bolger, Lohse, & Usadel, 2014). The clean reads were mapped to the sequences in the rRNA database of all published insects downloaded from NCBI to discard rRNAs using SOAP (Li, Li, Kristiansen, & Wang, 2008). Only the clean reads with the standard of Q30 > 85% and processed with Trimmomatic 0.32 and SOAP were used for further analysis. As we wanted to compare the two species and only the genome for B. dorsalis was available (ASM78921v2; https://i5k.nal.usda.gov/Bactrocera_dorsalis), Trinity for research without the need for a genome sequence with set parameters (min_ kmer_cov: 2 min_contig_length: 200 group_pairs_distance: 500) was used for de novo transcriptome assembly to obtain corresponding transcripts (Guo et al., 2018;Haas et al., 2013). The two species were firstly assembled to get their own UniGene database, and then the general UniGene library was obtained by clustering the two individual databases through CD-Hit to compare the expression and analyze differences between the two species. Trinity was used to break down clean reads into short fragments (Grabherr et al., 2011). Reads of certain lengths of overlap were combined in contigs to form longer fragments. These assembled unigenes were further processed for sequence splicing and redundancy removal using sequence clustering software to acquire maximum length nonredundant unigenes (Fu, Niu, Zhu, Wu, & Li, 2012). Briefly, Trinity was used for assembly to obtain corresponding transcripts, and we retained the transcript with the highest read coverage and removed the transcript with the lowest read coverage for each subcomponent. The transcripts selected in the clustering united as unigenes using the De Bruijn graph algorithm and CD-HIT to reduce sequence redundancy and improve the performance of other sequence analyses (Fu et al., 2012;Yang & Smith, 2013).
For functional annotation, sequence alignments of unigenes to the protein databases NR (nonredundant RefSeq proteins-NCBI),
Bowtie was used to align the reads to the Trinity transcripts and estimate the number of RNA-Seq fragments (Langmead, Trapnell, Pop, & Salzberg, 2009). With Bowtie alignment, FPKM (fragments per kilobase of exon per million reads mapped) were quantified using RSEM software (Li & Dewey, 2011). To evaluate the correlation of gene expression across different replicates, we calculated Pearson's correlation coefficient (r) between the three biological replicates (Schulze, Kanwar, Gölzenleuchter, Therneau, & Beutler, 2012). The analysis of differentially expressed genes (DEGs) was conducted in three steps. Firstly, read counts were adjusted by DESeq package through scaling normalized factor for each sequenced library before analysis (Anders & Huber, 2010;Wang, Feng, Wang, Wang, & Zhang, 2009). Secondly, differential expression analysis of samples was performed using the DESeq (Anders & Huber, 2010;Wang et al., 2009). In this analysis, p-value adjusted by the Benjamini-Hochberg approach was used to decrease false positives (Ferreira & Zwinderman, 2006). Lastly, the standard fold change>|2| was set as the threshold for differential expression. To visualize transcription differences among treatments, a principal component analysis

| Quantitative Real-time PCR for transcriptome verification and the detection of RNA interference efficiency
In order to validate the results from the transcriptome analysis, we quantified the expression profiles of candidate genes Hsp23, Hsp70 and Hsp90 for B. dorsalis and B. correcta. Thirty individuals for transcriptome verification (two species with three temperature treatments: 25, 35 and 38°C for 4 hr followed by 25°C for 1 hr as described in Supporting Information Figure S1B) were collected randomly. For the detection of RNA interference efficiency, 10 individuals fed on dsRNA after 96 hr in each species were randomly selected for Hsp23 expression detection. RNA extraction and cDNA synthesis followed the description in Section 2.3. The RNAi efficiency was quantified by quantitative real-time PCR using SYBR ® Premix Ex Taq™ II (TliRNaseH Plus; Takara, Japan) on an ABI 7500 instrument (Applied Biosystems Europe, Belgium).
All RNA samples were performed with same methods and analyzed in triplicate. The reaction included 1 μl cDNA, 12.5 μl SYBR Green mix, 1 μl each of forward and reverse primers (10 p.m.), 0.5 μl ROX Reference Dye II and 9 μl ddH 2 O. The thermocycler conditions were 95°C for 30 s, followed by 40 cycles at 95°C for 5 s and 60°C for 34 s. Melting curve analysis was performed at the end of each expression analysis, using the following conditions: 95°C for 15 s, followed by 60°C for 60 s with the decreasing rate of 1°C/s from 95°C. All the primers used are described in Supporting Information Table S4, and the amplification efficiency of all the primers in these two species is close to 100% (Supporting Information Table S5; Hu et al., 2014;Shen, Huang, Jiang, Dou, & Wang, 2013). We evaluated the performance of five reference genes, including 18s ribosomal RNA (18s), ribosomal protein L13 (RPL13), succinate dehydrogenase (SD), α-tubulin (α-TUB) and β-tubulin (β-TUB) in B. dorsalis and B. correcta using three softwarebased approaches (BestKeeper, NormFinder and geNorm) and 18s was chosen as the reference gene in the experiment (Supporting Information Table S6). The qRT-PCR data were analyzed using the 2 −ΔCT method (Chen & Wagner, 2012). Five biological replicates were carried out for statistical analysis. All results from experimental replicates were analyzed with Student's t tests or one-way analyses of variance (ANOVAs) with SPSS 20 (IBM Corporation, USA).

| Cloning the open reading frame of Hsp23
Hsp23 was selected as a candidate hardening gene based on the analysis of the transcriptomic data of B. dorsalis and B. correcta.
To verify the open reading frame (ORF) of Hsp23 in each species, due to the high similarity on sequence, one pair of primers named Hsp23-dsRNA-F/R (Supporting Information

| Phylogenetic analysis
The integrity of homologous amino acid sequences of other species was retrieved from the NCBI server (https://www.ncbi.nlm.nih.gov/; Dou et al., 2017). Sequences were first aligned by the conserved sequences, and then, phylogenetic analysis was performed using the neighbor-joining method in the Molecular Evolutionary Genetics Analysis software (MEGA version 5.1; Tamura et al., 2011). All the positions that contained gaps and missing data were eliminated before alignment and phylogenetic analysis.

| Validation of Hsp23 role through dsRNA
Hsp23 was chosen for functional analyses because it is the most expressed gene in B. dorsalis, the species with the stronger invasion potential, which likely involved the ability of this species to adapt to increased temperature. While Hsp70 was the most ex- Green fluorescent protein (dsGFP) was used as a negative control.
The dsRNA was synthesized by the primers named Hsp23-dsR-NA-F/R and Hsp23-F/R-T7 in Supporting Information Table S4 and using the T7 RiboMAX Express RNAi system (Promega). For larval feeding, three grams artificial diet with 30 μl of a dsRNA solution (1,000 ng/μl) was put into a 50 ml tube. Forty 3-day-old 1st instar larvae for hardening response study or 4-day-old 2nd instar larvae for heat tolerance study of B. dorsalis and B. correcta were collected and moved into a 50 ml tube with artificial diet, which had three holes in its lid for air. Five replicates per treatment were carried out and each replicate contained 40 larvae. For the temperature study, larvae were fed with dsHsp23 and dsGFP for 48 hr and transferred to new artificial diet for another 48 hr. After 96 hr feeding, the larvae developed to 3rd early-or late-instar larvae.

| Extreme heat stress and hardening response after feeding dsHsp23
After 96 hr dsRNA feeding, 10 larvae of each species of B. dorsalis and B. correcta were killed for detecting Hsp23 gene expression. The rest of the 3rd late-instar larvae were transferred to one 2 ml tube with diet for high temperature exposure. We chose 3rd late-instar larvae because of their high survival rate under extreme temperature to study the function of Hsp23 on heat tolerance (Jang, 1991). The larvae were directly exposed to heat stress at 45°C for 1 hr followed by 4 hr at room temperature before scoring. As there was no response in the Hsp23-knocking down in B. correcta at extreme high temperature, hardening response was only studied in B. dorsalis. Five groups of 40 3rd early-instar larvae with low survival rate under extreme temperature were used to study the function of Hsp23 on heat hardening, which showed the recovery of Hsp23 expression could increase survival rate. After the feeding, the larvae were exposed to heat hardening (35 and 38°C) and control (25°C) temperatures separately for 4 hr and returned to 25°C for 1 hr (Supporting Information Figure   S1B). Ten of the flies were frozen in liquid nitrogen for the detection of Hsp23 expression. The other 30 flies were exposed to heat stress at 45°C for 1 hr in a water bath (Supporting Information Figure S1A). After heat stress, the larvae were moved back to 25°C and their survival rate scored after 4 hr. Each treatment was replicated five times.

| Hardening response of B. dorsalis and B. correcta
We calculated the survival rate of 3rd instar larvae in the two species ( Figure 2

| Transcriptome and bioinformatic analysis of B. dorsalis and B. correcta
The percentage of bases with a Q30 score was over 92.64%, indicating that the sequencing was reliable. We obtained 94,821 unigenes from larval transcriptome sequencing, and a total of 42,931 unigenes were annotated in the larvae transcriptome from databases in these two species (Supporting Information Tables S1 and S2). The length distribution of the unigenes and statistics for larval unigenes assessment are shown in Supporting Information Figure S2 and Table S1, and the number of unigenes with a length >1,000 bp was 16,114.
The mapped ratio for all the samples was over 76% (Supporting   Information Table S3). In the transcriptome analysis, the PCA, bar charts and heatmaps of DGEs analyses revealed differences between the Bactrocera species under different hardening temperatures (Figures 3 and 4). In the PCA analysis of these two species,  Figures S3D and S4D).
However, genes such as larval serum protein, dehydrogenase, cecropin and mucin were significantly downregulated. In B. correcta, 181 genes were upregulated and 201 genes were downregulated.
Genes annotated with endocytosis, protein processing in endoplasmic reticulum and spliceosome in KEGG analysis were upregulated (Supporting Information Figures S3A and S4A). Downregulated genes were annotated with peroxisome, glycolysis, pyruvate metabolism and biosynthesis of amino acids (Supporting Information Figures S3A and S4A). These patterns indicated substantial changes in gene expression affecting a diversity of pathways.
In the comparison of 38 and 25°C, a total of 1,329 and 826 genes showed significantly different expression in B. dorsalis and B. correcta, respectively. In B. dorsalis, 815 genes were upregulated and 514 genes downregulated (Figure 3b). Genes related to lysosome, glycolysis, carbon metabolism and biosynthesis of amino acids were upregulated at 38°C (Supporting Information Figures S3E and S4E).
Genes related to carbon metabolism, metabolism of xenobiotics by cytochrome P450, protein processing in endoplasmic reticulum and pentose and glucuronate interconversions were significantly downregulated (Supporting Information Figures S3E and S4E). These changes reflected energy metabolism such as glycolysis which may have provided energy for the hardening response, and protein processing in the endoplasmic reticulum which may affect the accumulation of unfolded proteins (Ron & Walter, 2007). In B. correcta, 435 genes were upregulated and 391 genes downregulated (Figure 3b).
Genes related to endocytosis, protein processing in endoplasmic reticulum and spliceosome were upregulated at 38°C (Supporting Information Figures S3B and S4B), which was similar to the response at 35°C. Genes such as peptidoglycan-recognition protein, cuticle protein and lambda-crystallin homolog were significantly downregulated. For the KEGG annotation analysis, the regulated genes mainly involved the peroxisome, protein processing in endoplasmic reticulum and metabolism of xenobiotics by cytochrome P450 (Supporting Information Figures S3B and S4B), this suggests that heat stress may affect a range of oxidative processes.  Figures S3F and S4F). There were also many differentially expressed genes with unknown functions. In the KEGG analysis, most B. correcta genes related to the peroxisome, protein processing in endoplasmic reticulum and fatty acid activities (Supporting Information Figures S3C and S4C).
There were more genes and processes affected at 38°C than at 35°C (Figure 4a). Genes that showed certain expression patterns in the two species formed 26 clusters (Figure 4b). Genes in cluster 2 F I G U R E 2 Heat hardening response of Bactrocera correcta and B. dorsalis. The survival rate of 7-day-old 3rd instar larvae exposed to high temperature after heat hardening treatments. The temperatures used for hardening are indicated on the x-axis. The blue line represents the fitted curve of larval survival rate using an asymmetric model for B. correcta. The red line shows the fitted curve for B. dorsalis using the same model. Different letters above the bars indicate significant differences at p < 0.05, as determined by a Tukey HSD test. The capital letters in blue refer to B. correcta and the lowercase letters in red refer to B. dorsalis. The error bars indicate 1 SE and 9 annotated with immune response such as defense response to bacterium and innate immune response, stress response and some monooxygenase activities were upregulated only in B. correcta at both hardening temperatures, while the related genes were downregulated in B. dorsalis. However, genes in cluster 13,16,17,18,20 and 24 related to many developmental processes were only upreg-  Table S7). Also, many other sHsps were associated with the hardening process, such as Hsp18.4 and Hsp20 (Supporting Information Figures S5 and S6 and Table S7).  Figure S7). In general, Hsp23 is one of the highest expressed genes in all DEGs including Hsps in both species ( Figure 5).

| Cloning and characterization of Hsp23
Hsp23 from both species was cloned from cDNA. Based on the phylogenetic analysis and sequencing results, BcHsp23 and BdHsp23 were similar to each other, with a difference in nucleic acid and amino acid sequences of ORFs of only 17 out of 513 base pairs and 7 out of 170 amino acids.

| Phylogenetic analysis of sHsp from different insect species
The NJ tree (Supporting Information Figure S8) showed two annotated Hsp23 genes in B. correcta and B. dorsalis were clustered together, with high sequence similarity to the Hsp23 group from other species.  Figure 7b). In these two species, the other two Hsps had the similar expression pattern; Hsp70 expression was increased at 38°C but there was no change in Hsp90 at the two hardening temperatures (Figure 7).

| Gene function
Bactrocera dorsalis has been suggested to have a stronger invasion ability compared with other tephritid fruit flies, and we chose Hsp23 in this species to study whether this gene influenced temperature adaptation and caused the difference between these two species Malacrida et al., 2007). Despite the similar expression pattern and sequence, the gene showed distinct functions between these two species. Therefore, an exposure of 3rd late-instar larvae was used to study whether low expression of Hsp23 could decrease survival. Compared with the dsGFP-feeding group, dsHsp23 exposure led to a 0.48-fold reduction in expression in B. dorsalis (Figure 8a) while survival rate under 45°C was decreased by 42.9% (Figure 8b).
The result suggests that constitutive expression of Hsp23 increases survival of B. dorsalis under heat. We used two hardening temperatures of 35 and 38°C to study whether suppressed Hsp23 expression influenced survival rate following hardening. In the dsHsp23-feeding treatments, compared to the 25°C group, fold change expression of

| Species differences in hardening
We showed that while both Bactrocera species were hardened by nonlethal temperatures, the temperatures range of 37-40°C led to a stronger hardening response for B. dorsalis, while B. correcta performed better following exposure to 34-36°C (Figure 2). The more narrowly distributed B. correcta seemed to show a heat F I G U R E 5 Relationship between mRNA levels for genes and heat shock protein genes (Hsps) with fold change more than 2 under hardening/control treatment comparisons (35, 38 and 25°C controls). DEGs are measured as a logarithmic value of mean FPKM for three replicates (log 10 mean FPKM ) in Bactrocera correcta (Bc) and B. dorsalis (Bd). Transcripts encoding Hsps are shown as purple triangles, Hsp23 is highlighted in red (arrow), and all other transcripts are represented by gray circles. The black line represents k = 1, and the red line and zone (±95% confidence bands) represent the least-squares regression fit of Hsp transcripts, for which Hsp23 was a significant outlier in B. dorsalis under the two hardening temperatures when compared to 25°C hardening response more rapidly than the more widely distributed species. These results highlight differences among species in hardening responses which has also been noted in other groups of related species, particularly in Drosophila (Malmendal et al., 2006;Nyamukondiwa et al., 2011;Overgaard, Sørensen, Com, & Colinet, 2014;Willot et al., 2017). They also appear to be consistent with previous studies on temperature responses of other life stages of these two species (Hu et al., 2014;Liu & Ye, 2009;Pieterse et al., 2017). For B. dorsalis, heat shock tolerance of larvae was significantly enhanced by exposing to heat hardening 37 and 39°C for 1 or 2 hr (Hu et al., 2014;Pieterse et al., 2017). In B. correcta, the temperatures from 30 to 33°C appeared to be the most suitable for egg, larva and pupa development, and the preoviposition time was shortened following exposure to 33-36°C (Liu & Ye, 2009).
Whether these differences in hardening responses contribute to differences in the distribution of the two species is unclear.
Other comparisons of insects also provide mixed evidence for associations between hardening and species distributions. In a comparison of the widely distributed C. capitata with its more narrowly (1 hr) which altered survival at 41°C (Nyamukondiwa, Kleynhans, & Terblanche, 2010). Rapid cold hardening ability has been linked to species distributions in Collembola, and heat hardening has been related to different distributions in two Cataglyphis ants (Bahrndorff, Loeschcke, Pertoldi, Beier, & Holmstrup, 2009;Willot et al., 2017), but there is often no link between hardening and species distributions in Drosophila (Overgaard, Kristensen, Mitchell, & Hoffmann, 2011). Apart from hardening responses, other environmental factors can also influence the distribution and abundance of tephritid species (Vayssières, Carel, Coubes, & Duyck, 2008).

Some factors affect the fitness of fruit flies and have an indirect
influence on their distribution, while others have a direct influence.
It was found that high temperatures, low humidity and the absence of suitable host fruits for oviposition could cause ovarian immaturity thus influence further spread and distribution of Bactrocera species (De Meyer et al., 2010). Due to the preference for certain hosts, Dacus ciliatus enhances its biotic potentialities and maintains its population at low levels especially at low altitudes to avoid competition with the melon fly B. cucurbitae (Vayssières et al., 2008). In the PCA analysis and the comparison of up-and downregulated genes (Figure 3), it appears that the transcriptome of B. dorsalis is more responsive to high temperature than B. correcta, especially under 38°C, which contributed to a more variable survival rate.
Interestingly, B. correcta shows a more rapid hardening response yet its transcriptome is not as responsive and the genes in B. correcta regulated more stably, which was consistent with survival rate. The different responses between species were also evident in our GO and KEGG analysis (Supporting Information Figures S3 and S4), with the regulated genes in B. correcta annotated with protein synthesis that could protect organisms against heat and oxidative stresses by eliminating hydrogen peroxide (Lu, Bai, Zheng, & Lu, 2017), whereas in B. dorsalis, changes to genes involved in a range of metabolic processes such as carbon metabolism were involved under hardening temperatures ( Figure 4).
When considering both species, the response to stimulus was one of the most significantly enriched GO terms considering the DEG number and ratio of DEGs of the inquiry set in all related genes of the reference set (Supporting Information Figure S3).
Hsps were the largest family to respond to stimulus among all the DEGs, with the highest levels of expression ( Figure 5)

| Hsp23 expression and heat tolerance
We showed that suppressing Hsp23  18s was used as the control gene. The letter "*" and "**" above the bars represent significant differences at p < 0.05 and p < 0.01, respectively, and "ns" represents no significant differences. Different letters above indicate significant differences at p < 0.05, as determined by a t or Tukey HSD test suggested being important mechanisms for protecting other cellular proteins from denaturation under thermal stress and other stresses, and overexpression of sHsps can enhance the tolerance of cells to temperature changes (Li et al., 2009;Wang et al., 2017). Among sHsps, Hsp23 shows a key role in temperature adaptation in many species. In Drosophila melanogaster, gene knockdown experiments have suggested that Hsp22 and Hsp23 genes contribute to adaptive responses to fluctuating thermal conditions and particularly in chill coma recovery (Colinet, Lee, & Hoffmann, 2010). Also, overexpression of Hsp23 muscle-specific has been suggested to promote proteostasis and protect muscle from heat stress (Kawasaki et al., 2016).
It is also possible that this gene has no impact on the hardening response of B. correcta, despite the substantial upregulation of this gene under the hardening treatment. We did not consider the association between this gene and the hardening response in B. correcta for three reasons. Firstly, despite a high level of sequence conservation between the species, the thermal responses did not decrease at 45°C even with successful knockdown as we mentioned in the methods, suggesting that Hsp23 has limited or no effect on heat tolerance ( Figure 8) (Figures 2 and 7a).
Ambient temperature is a key environmental factor influencing the ecology and evolution of ectotherms, and invasive species can be used as models for investigating issues related to adaptive processes under changing ambient conditions and bringing ecological and evolutionary processes into a common framework (Bennett & Lenski, 2007;Diamond, Chick, Perez, Strickler, & Martin, 2017;Ghalambor, McKay, Carroll, & Reznick, 2007). In our research, despite being similar, these two invasive species differ in many respects such as heat hardening profiles, transcriptome responses and Hsp expression. These differences likely mediate thermal hardening and involve genes such as Hsps and pathways such as energy metabolism and biosynthesis. As more genes and pathways are regulated in B. dorsalis, this species may be more adaptable under high temperature stress.

| CON CLUS IONS
It is important to understand the role of genes such as Hsps and pathways in relation to the hardening response and stress resistance, in order to develop a mechanistic basis to evolutionary change and plastic responses . In our study, we used two successful invasive species that have invaded different thermal environments to understand the role of hardening and plasticity in climate adaptation. Based on our results, the Bactrocera species have a varied ability to adapt to temperatures. Different hardening responses may be related to the regulation of certain genes like Hsp23 which in turn could contribute to species distributions and invasive potential. The widely distributed species, B. dorsalis, seems more sensitive to temperature change at the transcriptome level.
This sensitivity may help adaptive thermal responses, whereas B. correcta might be at an advantage in milder environments.
However, this study only provides a starting point for understanding the genomic basis of climate adaptation in invasive fruit flies and the functional studies in particular need to be expanded to other gene families. Eventually, these types of studies may indicate targets for genetic manipulation to eventually control these pests. For instance, the importance of Hsp23 in heat hardening in B. dorsalis may point to this gene as being a useful candidate for gene drive technology to suppress this pest at warm times of the year.

ACK N OWLED G EM ENTS
The authors would like to thank Mr. Guoping Zhan and Mr. Chen Ma in Chinese Academy of Inspection and Quarantine for providing samples of B. dorsalis and B. correcta. This study was financially supported by the National Natural Science Foundation of China (no. 31772230).

CO N FLI C T O F I NTE R E S T
None declared.

DATA ACCE SS I B I LIT Y
RNA-Seq data are available at National Center for Biotechnology