Comparative Transcriptome and Proteome Analysis of Heat Acclimation in Predatory Mite Neoseiulus barkeri

In our previous study, we reported a high temperature adapted strain (HTAS) of the predatory mite Neoseiulus barkeri was artificially selected via a long-term heat acclimation (35°C) and frequent heat hardenings. To understand the molecular basis of heat acclimation, ‘omics’ analyses were performed to compare the differences between HTAS female adults to conventional strain (CS) at transcriptional and translational levels. We obtained a total of 5,374 differentially expressed genes and 500 differentially expressed proteins. Among them, 119 transcripts had concurrent transcription and translation profiles. It’s conserved that some processes, such as high expression of heat shock protein (HSP) genes, involved in heat tolerance of transcriptome analyses, while many protective enzymes including glutathione S-transferase, superoxide dismutase, peroxidase, and cytochrome P450 displayed down-regulated expression. KEGG analysis mapped 4,979 and 348 differentially expressed genes and proteins, to 299 and 253 pathways, respectively. The mitogen-activated protein kinases (MAPK) signaling pathway may provide new insights for the investigation of the molecular mechanisms of heat tolerance. Correlation enriched pathways indicated that there were four pathways associated with heat acclimation involving in energy metabolism and immunity. In addition, the expression patterns of ten randomly selected genes including HSP were consistent with the transcriptome results obtained through quantitative real-time PCR. Comparisons between transcriptome and proteome results indicated the upregulation of HSPs and genes participated in ATP production, immunity and energy metabolism process. A majority of antioxidant-related genes and detoxication-related genes were down-regulated suggesting a fitness cost of heat acclimation. Our results demonstrated that heat tolerance during a long-time acclimation of N. barkeri is a fairly complicated process of physiological regulations. These findings also contribute to a better understanding of the mechanisms of thermal responses of phytoseiid mites which could provide useful information for biological control through natural enemies.


INTRODUCTION
Phytoseiid mites (Acari: Phytoseiidae), as the effective natural enemies of spider mites and small insect pests, are widely used for biological control on fruits, vegetables, and other crops worldwide (Nomikou et al., 2001;Arthurs et al., 2009;McMurtry et al., 2013;Kamel et al., 2018). Both phytoseiid mites and pest mites are tiny organisms, and they are prone to environmental challenges in the field, such as thermal stress (Veerman, 1992;Gotoh et al., 2004), pesticides (Bin Ibrahim and Yee, 2000;Li et al., 2006;Sato et al., 2011), ultraviolet-B radiation (Tachi and Osakabe, 2012;Ghazy et al., 2016) and drought stress (Ferrero et al., 2010). Most studies have proved that phytoseiid mites were more vulnerable than the spider mites under environmental stressors including high temperature (Stavrinides and Mills, 2011) and ultraviolet radiation (Tachi and Osakabe, 2012), while pest mites could endure and survive from them, which caused herbivore natural control disruptions and failures in biological control strategies (Guzman et al., 2016). Climate change may become a crucial limiting factor for biological control in agriculture, since the chances of pests to escape from predator control may increase because vulnerability to rising temperature usually increases with trophic level. Vulnerability to rising temperatures usually increases with trophic level which may cause increasing chances of pest to escape predation in climate change situation (Voigt et al., 2003).
Neoseiulus barkeri Hughes (Acari: Phytoseiidae) is an effective natural enemy and important biological control agent for spider mites and several small insect pests (Bonde, 1989;Wu et al., 2014) which has been widely distributed and applied in China (Niu et al., 2014). However, due to their poor adaptability in the field condition, its application in agro-ecosystems is often unsatisfactory. To improve this worrying situation, a high temperature adapted strain (HTAS) of N. barkeri was selected from a conventional strain (CS) that maintained at 25 ± 1 • C via a long-term heat acclimation and frequent heat hardenings in 2012 (Zhang G. H. et al., 2018). It has been proved that the heat acclimation greatly improved their survival probabilities under heat stress, e.g., the survival probabilities of adult females up to 90% when exposed at 45 • C for 2-6 h, indicating a strong tolerance of HTAS N. barkeri to heat stress was obtained or developed after heat acclimation. In addition, long-time heat acclimation has also resulted in accelerated growth and developmental rate and reduced total fecundity and longevity which indicated the potential costs on fitness of heat acclimation (Zhang G. H. et al., 2018). Acclimatory responses (short or long-time) are often considered a multistep process, including detecting environmental cues, transducing signals into a cellular response, and activating certain genes, proteins and metabolites that regulate physiological function to enhance thermotolerance. The physiological adaptation to thermal acclimation has been extensively studied in Drosophila and other insects (Kregel, 2002). Acclimation of Drosophila melanogaster could induce adaptive changes in metabolic rates and membrane lipids composition, indicating acclimation-related physiological adjustments (Berrigan, 1997;Overgaard et al., 2008). Analysis of proteomic responses induced by heat stress in insects showed many thermal tolerance related proteins involved in iron ion and cell redox homeostasis, carbohydrate and energy metabolism (Colinet et al., 2013), structural elements of the cytoskeleton (Nguyen et al., 2009) and stress-induced signal transduction (Li et al., 2012). However, in spite of the findings and extensive knowledge of thermotolerance mechanism in genetics among many insects, the understanding of acclimation-related physiological adjustments of phytoseiid mites remains poorly studied.
For a more detailed understanding of the thermal acclimatory responses of phytoseiid mites and the physiological difference between two strains of N. barkeri, the comparative transcriptome and proteome were analyzed. We speculated that genes related to heat acclimation might be associated with heat shock protein, energy metabolism and immune system, which may have a different expression pattern from sudden heat stimulation. Response to heat tolerance is a complex process and our study could be conducive to identify the mechanisms and molecular pathways that are responsible for the capacity of acclimated predatory mite N. barkeri to maintain metabolic homeostasis and survive under thermal stress and heat acclimation.

Mites Rearing and Sampling
The CS and HTAS of N. barkeri were kept at a 25 ± 1 • C and 35 ± 1 • C, respectively with 70-80% RH and L: D = 14:10 h photoperiod in a climate-controlled room (Zhang G. H. et al., 2018). The CS mites were originally obtained from the Institute of Plant Protection, Chinese Academy of Agricultural Sciences in 2010. The HTAS mites had been heat acclimated for at least 300 generations at 35 ± 1 • C and heat hardened at 42 • C for 70 times (every 15-25 days) before our study. Both colonies were separately maintained on plastic non-transparent sheets resting on wet cotton and fed daily with the two-spotted spider mite (Tetranychus urticae) of different life stages. Five hundred newly emerged females (1-2-day-old) from each colony were collected into a 1.5-mL centrifuge tube respectively, and stored as a replicate sample at -80 • C for total RNA or protein preparation.

Preparation of cDNA and Transcriptome Sequencing
For each strain, three replicate samples were prepared for the total RNA isolation using Trizol Reagent (Invitrogen, CA, United States) according to the manufacturer's protocol. The integrity and quantity of RNA were determined using the Agilent 2100 Bioanalyzer (Agilent Technologies, CA, United States).
Library preparations were sequenced on the BGIseq-500 platform (BGI, Shenzhen, China) (Zhu et al., 2018). The library was amplified by phi29 to make DNA nanoball (DNB) which has >300 copies of a molecular. DNBs were load into the patterned nanoarray, and the single end 50 bases reads were produced in the way of sequenced via synthesis. The DNB patterned array technology not only provides sequencing accuracy but also increases chip utilization and sample density.

Annotation and Identification of Differentially Expressed Genes
To identify the differentially expressed genes from N. barkeri that are regulated by heat acclimation at a transcriptional level, a FDR (false discovery rate) ≤ 0.001 and an absolute value of log 2 (fold-change HTAS/CS) ≥ 1 were used as the thresholds to screen DEGs (differentially expressed genes) (Love et al., 2014). For functional annotation, distinct sequences were searched via BLAST against the NCBI NR database with a cut-off E-value of 10 −5 . Blast2GO 1 was used to assign Gene Ontology terms (GO 2 ), while the database resource of Kyoto Encyclopedia of Genes and Genomes (KEGG) 3 which integrates genomic, chemical, and systemic functional information was adopted to annotate molecular networks (pathways). The transcriptome raw reads were deposited in the Sequence Read Archive (SRA) bioproject number PRJNA492849.

Protein Preparation, iTRAQ Labeling and SCX Fractionation
Six samples were grounded into powder in liquid nitrogen and extracted with lysis buffer (8 M Urea, 40 mM Tris-HCl, 2 mM EDTA (Ethylene diamine tetraacetic acid) (E1170, Solarbio, China) and 10 mM DTT (Dithiothreitol) (D1070, Solarbio, China) pH 8.5, respectively. The mixtures were placed into a Tissue Lyser for 2 min at 50 HZ to release proteins. After centrifugation with 25,000 g at 4 • C for 20 min, the supernatant was transferred into another new tube. To reduce disulfide bonds, 10 mM DTT was added and incubated at 56 • C for 1 h. Then, the mixtures were alkylated by 55 mM IAA (Iodoacetamide) (I8010, Solarbio, China) in the dark at room temperature for 45 min. Following centrifugation (25,000 g, 4 • C, 20 min), the supernatant containing proteins was quantified by Bradford assay (P0006, Beyotime, China).
The protein solution (100 µg) with 8 M urea was diluted four times with 100 mM TEAB (Triethylammonium Bicarbonate) (T7951, Solarbio, China)and digested with Trypsin Gold (Promega, Madison, WI, United States) with the ratio of protein: trypsin = 40:1 at 37 • C overnight. After trypsin digestion, peptides were desalted with a Strata X column (Phenomenex, United States) and vacuum-dried according to the manufacturer's protocol for 8plex iTRAQ reagent (Applied Biosystems, Foster City, CA, United States). Samples were labeled with iTRAQ reagent. The peptides were dissolved in 30 µL 0.5 M TEAB with vortex. After the iTRAQ labeling reagents were recovered to ambient temperature, they were transferred and combined with proper samples. The labeled peptides with different reagents were combined and desalted with the Strata X column (Phenomenex, UCA) and vacuum-dried. The peptides were separated on a Shimadzu LC-20AB HPLC Pump system coupled with a high pH RP column. The peptides were reconstituted with buffer A1 (5% acetonitrile, 95% H 2 O, adjust pH to 9.8 with ammonia) to 2 mL and loaded onto a column containing 5 µm particles (Phenomenex, United States). The peptides were separated at a flow rate of 1 mL/min with a gradient of 5% buffer B1 (5% H 2 O, 95% acetonitrile, adjust pH to 9.8 with ammonia) for 10 min, 5-35% buffer B1 for 40 min, 35-95% buffer B1 for 1 min. The system was then maintained in 95% buffer B1 for 3 min and decreased to 5% within 1 min before equilibrating with 5% buffer B1 for 10 min. Then, elution was monitored by measuring absorbance at 214 nm, and the fractions were collected every 1 min. The eluted peptides were pooled as 20 fractions and vacuum-dried. Each fraction was resuspended in buffer A2 (2% acetonitrile and 0.1% formic acid in water) and centrifuged at 20,000 g for 10 min. The supernatant was loaded onto a C18 trap column 5 µL/min for 8 min using a LC-20AD nano-HPLC instrument (Shimadzu, Kyoto, Japan) by the auto-sampler. Then, the peptides were eluted from trap column and separated by an analytical C18 column (inner diameter 75 µm) packed in-house. The gradient was run at 300 nL/min starting from 8 to 35% of buffer B2 (2% H 2 O and 0.1% formic acid in acetonitrile) in 35 min, then going up to 60% in 5 min, then maintenance at 80% buffer B2 for 5 min, and finally returned to 5% in 0.1 min and equilibrated for 10 min.

Mass Spectrometry (MS) Analysis
The peptides separated from nanoHPLC were subjected into the tandem mass spectrometry Q EXACTIVE (Thermo Fisher Scientific, San Jose, CA, United States) for DDA (data-dependent acquisition) detection by nano-electrospray ionization to achieve a better coverage and reliable statistics. Mass spectrometry data were acquired in peptide recognition mode using the 20 most abundant precursor ions from a survey scan of 350-1600 at a resolution of 70,000 in orbitrap. The target value (full MS target and MS2 target: 3e6 and 1e5) was determined by predictive Automatic Gain Control with dynamic exclusion for selected precursor ions 15 s, and under fill ratio 0.1%. The under fill ratio is the minimum percentage of the target value likely to be reached at the maximum fill time.

iTRAQ Protein Identification and Quantification
The raw MS/MS data was converted into Mascot Generic Format (MGF) by Thermo scientific tool Proteome Discoverer (Thermo Fisher Scientific, United States), and the exported MGF files were searched using Mascot version 2.3.02 (Matrix Science, London, United Kingdom). The mass spectrometry proteomics data have been deposited with the ProteomeXchange Consortium via the PRIDE partner repository with the data set identifier PXD011379. For protein identification, the peptide mass tolerance was 20 ppm and the fragment mass tolerance was 0.05 Da. Only peptides with a score above the probability-based Mascot identity threshold at the 95% confidence interval were counted as identified.
To identify a protein, the global FDR must be less than 1% and each protein involved at least one unique peptide. For protein quantization by an automated software called IQuant for quantitatively analyzing the labeled peptides with isobaric tags, a protein must contain at least two unique peptides.

GO and Pathway Enrichment Analysis
Functional annotation of proteins identified in N. barkeri samples was carried out using Blast2GO which assigns gene ontology through BLAST searches against protein databases. The differentially expressed proteins were defined to be significantly regulated if the P-value is less than 0.05. The hyper geometric test was used to get the target GO terms in the GO enrichment analysis. The formula is as follows: Where N is the number of all identified proteins which matched to GO terms; n is the number of differentially expressed proteins; M is the number of proteins which matched a certain GO term and m is the number of differentially expressed proteins which matched a certain GO term. If the P-value is less than 0.05, the GO term is significantly enriched in differentially expressed proteins. All identified proteins were mapped to a pathway in the Kyoto Encyclopedia of Genes and Genomes (KEGG) database. Significantly enriched pathways containing differentially expressed proteins were identified using the same formula as in GO analysis.

Correlation Between Protein and mRNA Expression
To obtain the correlation between transcriptomic and proteomic database, the cutoff values were chosen to select the subsets of genes and proteins with distinctive expression signals. For each protein, we queried the RNA-seq data for expression patterns of matching transcripts (P-value < 0.05).

qRT-PCR Analysis
The genes selected according to the differentially expressed genes (DEGs) and differentially expressed proteins (DEPs) were investigated by qRT-PCR at the transcriptional level. The housekeeping gene β-actin was used as a reference gene  for normalization purposes and the primers used in this study were listed in Supplementary Table S1. One µg of total RNA was used to synthesize the first-strand cDNA using the PrimeScript TM RT reagent Kit (TaKaRa, Dalian, China). The qRT-PCR was determined with a Qtower3 System (Analytik Jena, Germany) using GoTaq qPCR Master Mix (Promega, Madison, WI, United States) according to the following thermal program: 1 cycle of 95 • C for 2 min, followed by 40 cycles of 95 • C for 15 s FIGURE 1 | Expression levels in HTAS vs. CS female adults. The variability pattern was displayed by volcano plot. Red points, blue points, and gray points represent up-regulated genes, down-regulated genes, and non-differentially expressed genes, respectively. The horizontal coordinate indicates expression level in CS female adults, while the vertical coordinate indicates expression level in HTAS female adults.
Frontiers in Physiology | www.frontiersin.org and 60 • C for 30 s. The relative expression levels were calculated using the 2 − Ct method (Livak and Schmittgen, 2001). All data was given in terms of relative mRNA expression as mean ± SE.

RNA-seq and iTRAQ Data Analysis
Approximately a total of 66.14 Gb transcriptome data and 39,623 unigenes were obtained from six samples of N. barkeri. The total length, average length, N 50 , and GC content of unigenes were 56,432,972 bp, 1,412 bp, 2,592 bp and 49.12%, respectively (Supplementary  Figure S1). The total number of sequences identified by mass spectrometry of N. barkeri proteomes was 336,503, which represented 92,671 peptide spectra and 24,811 distinct peptides (Supplementary Table S3). Among the 24,811 identified peptides, only 20% (5,082) were assigned to a putative protein by homology search against the transcriptome database, leaving approximately 80% (19,729) of the peptides unidentified.

Differentially Expressed Genes (DEGs) Between CS and HTAS
After heat acclimation, a total of 2,571 and 2,803 up-and down-regulated transcripts were differentially expressed in HTAS, respectively (FDR ≤ 0.001 and log 2 Ratio ≥ 1) (Figure 1). The top ten up-regulated genes that respond to heat acclimation were as follows: fructose-2,6-bisphosphatase, stromal interaction molecule, heat shock factor protein, phosphatidylinositolbinding clathrin assembly protein, AH receptor-interacting protein, calcium-binding protein, microtubule affinityregulating kinase, mucolipin, zinc finger protein, single-stranded DNA-binding protein and DNA replication licensing factor (Supplementary Table S4).
Functional annotation showed that 2,979 transcripts were differentially expressed between two strains in GO classification (≥2-fold change, FDR ≤ 0.001). These 2,979 differentially expressed genes were assigned to 53 GO terms (Figure 2) and classified into three categories including biological process, cellular component and molecular function, and mostly responses to the cellular process, metabolic process, singleorganism process, cell, cell part, membrane, organelle, binding and catalytic activity. In KEGG analysis, 4,979 DEGs were mapped to 299 pathways. Among the DEGs, many pathways FIGURE 2 | GO enrichment analysis revealed the biological processes most associated with detected DEGs. Based on the GO results, "cellular process," "metabolic process," "membrane" and "catalytic activity" were the most enriched GO terms under heat acclimation.
including ribosome, RNA transport, MAPK signaling pathway, citrate cycle, glyoxylate and dicarboxylate metabolism and fructose mannose metabolism were enriched (Figure 3). More specifically, we found mostly up-regulated genes were linked to MAPK signaling pathway, RNA transport and cell cycle (Table 1 and Figure 4). It is well known that genes encoding the heat shock proteins and the protective enzymes play a vital role in thermal tolerance in N. barkeri and other insects. In our present study, we found 9 out of 15 HSPs upregulated in HTAS. However, the protective enzymes including glutathione S-transferase, superoxide dismutase, peroxidase and cytochrome P450 displayed down-regulated expression in HTAS strain ( Table 2).
Heat Tolerance-Dependent Protein Expression in Two Strains of N. barkeri At proteome, a total of 500 differentially expressed proteins were identified in HTAS in comparison with CS, with 225 proteins up-regulated and 275 down-regulated after heat acclimation (Supplementary Table S3). Isocitrate dehydrogenase, an essential enzyme in the mitochondrial antioxidant system (Choi et al., 2018), was up-regulated by 3.16-fold in HTAS. Other up-regulated proteins were also identified in HTAS strain, including Chymotrypsin B (2.45-fold) which displays positive regulation of apoptotic process (Orlowski, 1999); methyltransferase (2.14-fold) involves in environmentally induced epigenetic modification mechanism of DNA methylation (Bird, 2002); serpin B10 (1.71-fold) an immune-related protein, plays a role in the regulation of protease activities during hematopoiesis and apoptosis (Schleef and Chuang, 2000), as well as cytochrome P450 (1.64-fold) which is ubiquitous heme-containing monooxygenases involve in a number of vital processes including detoxication and fatty acid metabolism (Nebert and Gonzalez, 1987). Proteins down-regulated in HTAS included maleylacetoacetate isomerase (0.5-fold) which catalyzes the GSH-dependent isomerization of maleylacetoacetate and play vital role in tyrosine metabolism (Kim et al., 2011); PHD finger (0.52-fold) which performs essential roles in the regulation of histone methylation and is an essential factor for epigenetic regulation and genome maintenance ; calcineurin-like phosphoesterase domain-containing protein 1 (CPPED1) (0.53-fold) which improves adipocyte  glucose metabolism when downregulation of its expression (Vaittinen et al., 2013); Ras-like GTP-binding protein family (0.53-fold) which are components of signaling pathways that link extracellular signals via trans-membrane receptors to cytoplasm and touch on virtually every aspect of cell biology including growth, differentiation, morphogenesis, cell division and motility and cytokinesis ; vitellogenin (0.63-fold) which belongs to a family of lipid transport proteins and is associated with reproduction (Sappington and Raikhel, 1998) and glutathione S-transferase (0.64-fold) which is involved in detoxification and antioxidant defense machinery protects against oxidative stress (Chasseaud, 1979;Gill and Tuteja, 2010) ( Supplementary Tables S5, S6).
To correlate protein and mRNA expression profiles, genes with both quantifiable mRNA and protein data were extracted and compared with the annotated RNA-seq libraries. Correlation between DEGs and DEPs showed 119 genes/proteins related to heat acclimation (Figure 5). Supplementary Tables S5, S6 represented the correlation between mRNA and protein, and the correlation coefficients between the protein and gene expression profiles of the same trend and opposite trend were 0.8410 and -0.6703, respectively (Figure 6 and Supplementary Tables S7, S8).

Gene Ontology and Pathway Enriched Analysis
Among the 5082 proteins, 2415 were sub-categorized into 57 hierarchically structured GO classes, including 24 biological process, 19 cellular component, and 14 molecular function, which mostly response to cellular process (52.3%), metabolic process (48.8%), cell part (47.1%), and catalytic activity (56.3%) (Figure 7). Three hundred and forty eight differentially expressed proteins were allocated to the reference pathways in KEGG ( Table 3). As a result, 23 pathways were enriched (P-value ≤ 0.05, Table 3). Correlation of the enriched pathways for DEGs and DEPs showed there were four enriched pathways related to heat acclimation including amoebiasis, glyoxylate and dicarboxylate metabolism, citrate cycle (TCA cycle) and protein digestion and absorption (Tables 1, 3).

Validation of Data Through Real-Time PCR (qRT-PCR)
Thousands of genes showed significantly different expression levels in this study. As a result, 10 genes and 10 proteins FIGURE 4 | Heat map of the expression levels of different pathways. The expression levels of some genes in "MAPK signaling pathway," "RNA transport" and "Citrate cycle (TCA cycle)" were shown in (A-C), respectively. The color scale is shown at the upper right, which denotes FPKM value from the lowest (green) to the highest (red).
The changing trend of the qRT-PCR results was consistent with the DEG expression profiling (Figure 8).

DISCUSSION
Thermal acclimation can significantly alter the thermotolerance of invertebrates, such as Tribolium castaneum   (Lü and Huo, 2018), Drosophila melanogaster (Colinet et al., 2013) and Nilaparvata lugens (Piyaphongkul et al., 2014). For example, high temperature acclimation (42 • C) could significantly improve the heat tolerance of adult females of T. castaneum (Lü and Huo, 2018). Our recent study proved that a long-term high temperature adapted strain (HTAS) of N. barkeri (acclimated at 35 • C) could withstand more severe heat stress, with a significantly lower mortality caused by heat stress than CS mites did (Zhang G. H. et al., 2018). This intraspecific variation on thermal susceptibility between HTAS and CS suggested N. barkeri might had been shaped by a differentially physiological strategy in response to temperature changes. As expected, the physiological strategy of HTAS N. barkeri to heat acclimation was related to the expression of some unique genes and enriched pathways. In total, 5,374 DEGs (2,571 upregulated and 2,803 down-regulated) has been identified in response to heat acclimation. The GO analysis results suggested that most of DEGs were associated with the metabolic reaction, enzyme catalysis and membrane process. Similar results were observed in T. castaneum in response to heat stress with the membrane, metabolic process, and catalytic activity were remarkably enriched (Lü and Huo, 2018). Many DEGs associated with "membrane" suggesting the cells related to heat stress mechanically were injured since the membrane is concerned with the requirement for repair promotion (Howard et al., 2011). Moreover, some important signaling pathways were involved in heat adaptation on N. barkeri adults, such as "Ribosome, " "Citrate cycle (TCA cycle), " "Cell cycle, " "Glyoxylate and dicarboxylate metabolism, " "RNA transport" and "MAPK signaling pathway, " were enriched in HTAS mites, indicating the importance of these pathways in forming heat tolerance. In insects, the innate immune system is a major effector response system (Gillespie et al., 1997). MAPK signaling pathway plays an important role in the immune processes of a variety of organisms, particularly with FIGURE 5 | Correlation between differentially expressed proteins and genes. The numerical value in each circle represents the quantity of genes or proteins, including identified genes and proteins, and genes or proteins related to heat acclimation respectively, and genes/proteins related to heat acclimation together.  respect to the regulation of the innate immune response (Zhan et al., 2018). It was proved that MAPK signaling pathway was activated by high temperature to break pupal diapause in the flesh fly (Fujiwara and Denlinger, 2007) and involved in the responses of Dendranthema morifolium to low temperature (Lu et al., 2018).
Surprisingly, the majority of DEGs within the MAPK signaling pathway were up-regulated in HTAS adults, suggesting these humoral immune-related pathways were activated in response to the heat stress acclimation . However, the metabolism of HTAS for heat acclimation was weak since the pathways such as glyoxylate and dicarboxylate metabolism and citrate cycle (TCA cycle) were enriched a majority of down-regulated DEGs (Figure 4). A similar reduction was found in Glyphodes pyloalis responding to heat tolerance . The suppression of metabolism reaction may reflect cellular homeostasis and an energy-saving mechanism to manage heat stress. Heat shock proteins (HSPs), acting as highly conserved molecular chaperones, displayed essential roles in various cellular processes, including preventing protein aggregation and denaturation and maintaining protein homeostasis during periods of thermal stress (Iwama et al., 1999). Fifteen HSPs of N. barkeri, including HSP40, HSP60, HSP70, HSP90 and several small HSPs, were observed and most of them exhibited upregulated expression in HTAS according to the transcriptome analysis. The preferential induction of HSPs related to a thermal challenge has been documented in many previous studies. In N. lugens, HSP70 is constitutively responsive to heat shock but down-regulated to cold tolerance, while remarkably upregulated after insecticide exposure, which reveals that HSP70 is evolutionarily and functionally diversified and involved in responses to environmental stresses such as the resistance or tolerance (Lu et al., 2017). Moreover, HSP70 could inhibit JNK and Caspase-3 activation to regulate innate immune responses (Shukla et al., 2014) and HSP60 could activate the stressactivated protein kinases p38 to regulate Toll signaling pathway in innate immune cells (Vabulas et al., 2001). In our previous study, HSP70 and HSP90 of CS N. barkeri female adults were significantly more intense in high temperature at 40 • C for just 1.5 h, which suggested HTAS was more tolerant under heat stress . Our transcriptome analysis showed the induced expression of four HSP70 may accelerate the refolding of damaged proteins and prevent protein aggregation during heat acclimation. Nine out of 15 HSPs of HTAS, including 1 HSP10, 1 HSP40, 1 HSP60, 4 HSP70, and 2 HSP90, displayed upregulated expression related to heat acclimation which served as an evolutionarily conserved mechanism.
Arthropods, like other aerobic organisms, face a constant oxidative challenge from the production of ROS during cellular metabolism when exposed to heat stress. Peroxidase genes have been confirmed on potential roles in protecting the organism from the oxidative damage challenge under three environmental stresses (thermal, UV and pesticide) in N. barkeri (Tian et al., 2017). It has been proved that the HTAS N. barkeri exhibited lower fitness under UV-B stress compared with the CS N. barkeri due to the lower levels of activity for some antioxidant enzymes in HTAS N. barkeri (Tian et al., 2019). Indeed, in this study, many antioxidant-related genes, including superoxide dismutase (SOD), peroxidase (POD), catalase (CAT) and thioredoxin were significantly down-regulated in HTAS (Table 2), which suggests a conservative progress of energy-saving response to deal with the toxic substance accumulation in HTAS and the difference in fitness cost between the two strains. Similar results were observed in detoxication-related genes including glutathione S-transferase (GST) and cytochrome P450. The detoxification pathways identified in our study will provide new insights for the investigation of the molecular mechanisms of fitness cost under environmental stresses. FIGURE 8 | qRT-PCR analysis of the expression levels of 10 Unigenes. X axis represents different Unigenes. CL333.Contig6_All, fructose-2,6-bisphosphatase; CL2646.Contig2_All, calcium-binding protein; CL2903.Contig4_All, mucolipin; CL180.Contig20_All, single-stranded DNA-binding protein; CL606.Contig8_All, stromal interaction molecule; CL859.Contig7_All, DNA polymerase; CL517.Contig1_All, RING finger and SPRY domain-containing protein; CL1204.Contig2_All, long-chain fatty acid; CL1418.Contig2_All, HSP70 and CL2746.Contig2_All, HSP20. Y axis represents the relative expression levels of genes.
Using high-throughput MS-based proteomics, 500 DEPs were determined including 225 up-regulated DEPs and 275 downregulated DEPs. These 348 differentially expressed proteins were allocated to the reference pathways in KEGG, where 23 pathways were enriched (Table 3). Sphingolipid signaling pathway has been proved to be involved in defensive mechanisms (Thuleau et al., 2013) which may be contributed to activating the defensive response of N. barkeri by heat stress. The activities of antigen process and presentation, related to immune response (Kelly and Trowsdale, 2018) were modulated by heat acclimation. The upregulation of immune-related genes including serpin B10 (1.71-fold) under heat acclimation further suggested that heat stress facilitated the components of immune defense.
Correlation analysis of differentially expressed proteins and genes displayed that a subset of genes and proteins were expressed with the same trend related to heat acclimation (Supplementary Table S7). These genes including HSP70 (CL1418.Contig6), HSP90 (Unigene15778), isocitrate dehydrogenase (CL4252.Contig1), vitellogenin (CL125.Contig1) and ATP-dependent DNA helicase (Unigene12306) were up-regulated, while cytochrome P450 (CL1691.Contig3), cuticle protein (CL1299.Contig2), choline dehydrogenase (CL854.Contig1), endochitinase (CL4004.Contig2) and many unknown and uncharacterized protein were among the genes down-regulated at both the transcriptional and translational levels. Cuticular lipids are important to the fitness of the insects and a major barrier to water loss. The cuticle protein was downregulated in HTAS which may imply a fitness cost of heat acclimation in desiccation situation (Huang et al., 2019). Physiologically similar phenotypes were identified in D. melanogaster which the proteomic analysis (2D-DIGE) showed the physiological changes of the acclimated D. melanogaster adults could drastically alter expressed proteins especially HSP70 (Colinet et al., 2013). The up-regulated gene isocitrate dehydrogenase as an essential enzyme catalyzes oxidative decarboxylation of isocitrate to α-ketoglutarate and CO 2 with a concomitant reduction of NADP + to NADPH (Wise et al., 2011). Studies have been confirmed that isocitrate dehydrogenases enzymes are involved in the thermal properties in the bacterium, Colwellia maris (Mouri and Takada, 2018). However, there are some genes, including lipoprotein receptor (Unigene8357), glutamine synthetase (CL2142.Contig1) and serine/threonine-protein kinase (CL3245.Contig2), were upregulated at the translational level but down-regulated at the transcriptional level (Supplementary Table S8). This effect may contribute to differences between genes and proteins in the posttranscriptional regulatory mechanism (He and Klionsky, 2009). Four pathways including amoebiasis, glyoxylate and dicarboxylate metabolism, citrate cycle (TCA cycle) and protein digestion and absorption were observed enriched in both DEGs and DEPs in HTAS. Citrate cycle (TCA cycle) metabolism was enriched in plant (Ren et al., 2018; and thermophilic bacteria (Cordova et al., 2017) and yeast  to obtain more ATP production and more NADH, since detoxification of furfural or phenolic compounds is an energy-consuming process under heat stress which indicated resistance to high temperature of N. barkeri was a process of accelerating energy consumption.
In this study, comparative analysis between the CS and the HTAS N. barkeri with 6-years heat acclimation of N. barkeri was firstly performed at a transcriptomic and proteomic level. The results of the GO and KEGG pathway enrichment analyses of 5,374 identified DEGs indicated the DEGs related to metabolism and immunity were enriched after heat acclimation. In addition, the transcriptome results revealed a number of genes that were potentially relevant to HSPs, antioxidant, and detoxication in N. barkeri. A majority of antioxidant-related genes and detoxication-related genes were down-regulated in HTAS, which suggested a conservative progress of energy-saving response and the balance of fitness cost. Correlation of DEGs and DEPs showed four pathways related to metabolism and ATP production were enriched, indicating the importance of energy and metabolism at high temperature. Genes relevant to HSPs and ATP production were up-regulated at both the transcriptional and translational levels. Excellently thermal tolerance in HTAS N. barkeri was a result of long-term heat acclimations and heat hardenings which is different from the response to a sudden heat stimulus. Six years of laboratory heat acclimation is likely to have had a profound impact on the developmental physiology of N. barkeri. Our results contributed to the better knowledge of the complex physiological response mechanism of heat acclimation and provide insight into the heat tolerance of insects.

DATA AVAILABILITY STATEMENT
The datasets generated for this study can be found in the transcriptome raw reads were deposited in the Sequence Read Archive (SRA) bioproject number PRJNA492849. The mass spectrometry proteomics data have been deposited with the ProteomeXchange Consortium via the PRIDE partner repository with the data set identifier PXD011379.

ETHICS STATEMENT
The research project was conducted on mite species that are not subjected to any specific ethical issue and legislation.

AUTHOR CONTRIBUTIONS
CT, YL, and HL conceived and designed the experiments. CT, YL, JH, and WC performed the experiments, analyzed the data, and drafted the manuscript. ZW and HL participated in manuscript drafted and modification. All authors read and approved the final manuscript.

ACKNOWLEDGMENTS
We wish to thank Dr. Guangyun Li (School of Biological Sciences, University of Auckland, Auckland, New Zealand) and Dr. Yan Shi for assistance with improving the manuscript.