The S. carpocapsae (strain All) was supplied by the Guangdong Institute of Applied Biological Resources, Guangdong Province, China (Han et al., 1997). We used Galleria mellonella larvae as the hosts of S. carpocapsae (Kaya and Stock, 1997). The infected hosts were cultured in a 25°C incubator. These nematodes were selected as the controls and original nematodes for the heat shock tests. After the culture period, the infective juveniles (IJs) were introduced into sterile distilled water in a culture dish in preparation for the next step.

We measured the nematode heat resistance based on procedures used in a previous study (Yaari et al., 2016). We exposed the nematodes to 32°C for 3 h. They were then transferred to 23°C for 1 h to recover and finally exposed to 35°C for 5 h. We then selected 10 drops with 10-μl suspensions in each drop and calculated the average nematode mortality. Each drop contained 80–100 nematodes. The surviving nematodes were cultured in the host in normal condition, and the offspring were tested in the same manner. As the survival rate exceeded 50%, the high temperature used in the experiment was increased from 35 to 38°C with successive generations of nematodes. The nematodes whose survival rate exceeded 50% at 38°C were collected as the treatment group. Samples for RNA extraction were taken from all instars. Both the heat shock treatment and the control group were prepared in three replicates. Approximately 2,000–2,200 nematodes, representing all growth stages, were collected for each replicate.

Total RNA for sequencing was isolated using the RNeasy^® plus Micro Kit according to manufacturer’s instructions (Qiagen, Hilden, Germany). A NanoVue spectrophotometer (GE Healthcare Bio-Science, Uppsala, Sweden) was used to assess the quality of total RNA. The degradation and quality of the RNA were tested by electrophoresis on a 1.0% agarose gel. To isolate poly(A) mRNA, magnetic beads with oligo (dT) were used. After the fragmentation of the mRNA under the effect of the fragmentation buffer, the cDNA was synthesized using the template mRNA fragments. These short cDNA fragments were resolved and purified by ethidium bromide buffer for end reparation and connected with adapters. After the cDNA samples were selected for PCR amplification and quantified and qualified using an Agilent 2100 Bioanaylzer and ABI StepOnePlus Real-Time PCR System (Agilent Technologies, Santa Clara, CA, United States). The library was sequenced using an BGIseq PE150 system (Shenzhen, China). The cDNA library construction and sequencing were performed by Beijing Genomics Institution (Shenzhen, China).

The adaptors and low-quality sequences from the raw reads were removed using FASTX¹. SOAP2 (version 2.2.1)² was used to filter out potential rRNA reads, and the data were filtered to be clean reads. The filtered clean reads were mapped to the reference genome, S. carpocapsae, using HISAT2 (Hierarchical Indexing for Spliced Alignment of Transcripts)³ (Rougon-Cardoso et al., 2016). The assembled de novo transcriptome was implemented using Trinity (Grabherr et al., 2011). The reads with overlap length were combined to structure the contigs and then mapped back to contigs. Trinity connected the contigs from the same transcript that could not be extended on either end and formed the defined unigenes. The unigenes obtained were distributed into two groups: singletons and clusters. With an E-value cutoff of 10^–5, the sequences were compared based on the databases of Nr and Kyoto Encyclopedia of Genes and Genomes (KEGG) database. The software Blast2GO was used to analyze the gene annotation with Gene Ontology (GO) terms⁴ with default parameters. The KEGG database was used to perform pathway annotation analysis⁵. Selected from the Nr annotation, the interesting genes were analyzed using NCBI BLASTx⁶.

After the clean reads were aligned to the genome sequence, the expression level of each sample was calculated using RSEM (Li and Dewey, 2011) based on the RPKM (Reads Per Kilobase per Million mapped reads) value. Next, significance levels (P-values) of DEGs between treatment and control groups were calculated using edgeR (Robinson et al., 2010), an R package⁷ recommended for use in cases without biological or technical replicates. According to the recommendation of the edgeR manual⁸, all of the parameters were set. By controlling the false discovery rate (FDR), two methods were used to correct significant levels (P-values): the Benjamini and Hochberg and Storey and Tibshirani methods (Benjamini and Hochberg, 1995; Storey and Tibshirani, 2003). For the accuracy of the improved DEGs, we defined the fold change of gene expression ≥ 2 and adjusted the P ≤ 0.001 as the criteria for significant differences in gene expression levels. To deeply extract the DEGs and pathway data that were related to the treatment and control, GO and KEGG pathway cluster analyses were conducted using Blast2GO and KOBAS, respectively. After revision by FDR, the notable levels were defined as those with a P ≤ 0.001.

To convey the results of the transcriptome, quantitative real-time PCR (qPCR) was conducted for the key genes in the heat reaction, including SOD and the heat shock proteins (HSPs). Following alignment using NCBI BLASTx, the isolated sequences were preliminarily confirmed, and primers of these genes are designed using Primer 3.0⁹ and delineated in Supplementary Table S2. The expression abundances of mRNA were detected by qPCR using ABI StepOne Plus^TM Real-Time PCR System (Life Technologies, Woodlands, Singapore) with GoTaq^® qPCR master mix (Promega, Madison, WI, United States). The templates for qPCR were the same as those for the RNA-Seq. qPCR detection was based on a 20-μl reaction, composed of 10 μl master mix, 7 μl ddH₂O, 1 μl template cDNA, and 1 μl of each primer. The cycling program is as follows: 95°C for 2 min, 40 cycles of 95°C for 15 s, and 60°C for 30 s, 60°C for 15 s, and 95°C for 30 s. The relative expression was calculated using the 2^{–Δ Δ Ct} method (Pfaffl, 2001) with Actin set as the reference gene (Yan et al., 2011). Fold changes were clarified after the relative expression values were standardized by the lowest value.

The raw data were filtered after sequencing, and the low-quality reads were discarded. We obtained 43.9 and 42.9 million clean reads with an average of approximately 96.6 and 96.1% Q20 bases for the treatment and control groups, respectively. The filtered clean data were mapped to the S. carpocapsae genome. The percentage of mapping was 82.2% in the control and 88.9% in the treatment group. A total of 72.3% (15,850) and 65.5% (14,354) of all the genes (21,917) in the treatment and control groups were obtained as expressed genes (FPKM > 1). For functional annotation, the sequences were aligned to the database, including Nr, AnimalTFDB 2.0, GO, and KEGG. The number of annotated genes in the TFDB database was 2,206, as well as 13,319, 12,525, and 6,902 genes in the KEGG, Nr, and GO databases, respectively.

A total of 7,891 determined genes have been significantly expressed between the treatment and control groups, including 4,902 upregulated and 2,989 downregulated DEGs (Supplementary Figure S1). The top 10 most upregulated and downregulated genes are shown in Supplementary Table S1. Among them, the top three upregulated DEGs were ladderlectin-like protein, proline-rich extension-like protein, and synapse-associated protein. The three most downregulated DEGs were WD repeat domain 39, followed by decaprenyl-diphosphate synthase and lipase-like protein. In addition, some longevity regulating pathway and phosphoinositide 3-kinase (PI3K)-Akt signaling pathway genes were conspicuously expressed, such as insulin-like receptor and nematode cuticle collagen domain protein. Finally, several DEGs were detected as hypothetical or predicted proteins, indicating that their functions are not known and merit investigation.

We classified significant enriched GO terms involved in the heat shock reaction by their determination as DEGs of the treatment and control groups. All of the DEGs were categorized into three GO types: biological process, cellular component, and molecular function, and these contained 19, 15, and 10 terms, respectively (Figure 1). In the biological process category, the top three terms with the largest number of DEGs included cellular process (874), metabolic process (723), and biological regulation (488). Among the cellular component category, the membrane (1,214), membrane part (1,143), and cell (988) terms contained the greatest number of DEGs. Binding, catalytic activity, and transporter activity were the top three terms in the molecular function category. Among all of these terms, the three least common terms were pigmentation (biological process), virion (cellular component), and virion part (cellular component), which all included four DEGs.

In the KEGG pathway enrichment analysis, 44 terms were classified to six categories: cellular processes, environmental information processing, genetic information processing, human diseases, metabolism, and organismal systems (Figure 2). Among them, signal transduction (791), global and overview maps (637), and cancers: specific types (620) were top three DEG terms. These terms were categorized to environmental information processing, metabolism, and human diseases, respectively. The term “biosynthesis of other secondary metabolites” included four DEGs, while the lipid metabolism term included 229 DEGs.

Annotation of the DEGs showed that some genes were clearly associated with the heat shock response. After isolation, some DEGs were identified as heat shock-related DEGs (HRDs) to comprehensively understand the stress response of the HRDs to heat shock. Based on the preliminary categories of heat shock genes in nematodes in previous studies (Yaari et al., 2016), we isolated and identified four pathways: Longevity regulating pathway, IIS (insulin/insulin-like signaling) pathway, Drug metabolism-CYP450 pathway, and Fatty acid metabolism pathway (Table 1). Among the four pathways, 25 HRDs were selected, including 17 upregulated and eight downregulated. The top three most upregulated HRDs were Fatty acids protein 6 (8.3), Acyl-CoA synthetase (7.4), and LIPS17 (7.1). DAF-2 (-11.2) was the most downregulated HRD.

Ten HRDs, which were clustered to different groups and key pathways, were isolated to validate the RNA-Seq results using qPCR. The expression trends of all 10 selected HRDs determined by qPCR were consistent with those detected by RNA-Seq, indicating that the transcriptome analysis was reliable (Figure 3).