Transcriptome-Based Analysis Reveals a Crucial Role of BxGPCR17454 in Low Temperature Response of Pine Wood Nematode (Bursaphelenchus xylophilus)

Background: The causal agent of pine wilt disease is the pine wood nematode (PWN) (Bursaphelenchus xylophilus), whose ability to adapt different ecological niches is a crucial determinant of their invasion to colder regions. To discover the molecular mechanism of low temperature response mechanism, we attempted to study the molecular response patterns under low temperature from B. xylophilus with a comprehensive RNA sequencing analysis and validated the differentially expressed genes (DEGs) with quantitative real-time polymerase chain reaction (qRT-PCR). Bioinformatic software was utilized to isolate and identify the low-temperature-related BxGPCR genes. Transcript abundance of six low-temperature-related BxGPCR genes and function of one of the BxGPCR genes are studied by qRT-PCR and RNA interference. Results: The results showed that we detected 432 DEGs through RNA sequencing between low-temperature-treated and ambient-temperature-treated groups nematodes. The transcript level of 6 low-temperature-related BxGPCR genes increased at low temperature. And, the survival rates of BxGPCR17454 silenced B. xylophilus revealed a significant decrease at low temperature. Conclusion: in conclusion, this transcriptome-based study revealed a crucial role of BxGPCR17454 in low temperature response process of pine wood nematode. These discoveries would assist the development of management and methods for efficient control of this devastating pine tree pest.


Introduction
As one of the most dangerous plant pests in the world, pine wood nematode (Bursaphelenchus xylophilus, PWN) causes devastating pine wilt diseases to the pine trees in Asia, Europe, and North America [1]. Unfortunately, the infestation area of PWN, in all probability, will continue expanding to colder regions of Asia and Europe [2][3][4][5]. For many parasites, PWN included, their ability to adapt different ecological niches is a crucial determinant of their invasion to new regions [6]. The formation of low-temperature-induced diapause stage of PWN is one of the most pivotal factors for PWN's survival and expansion. In the autumn, PWN will gradually stop their growth and gradually enter winter diapause in response to the decreased environmental temperature. Additionally, the second-stage propagative juveniles (J 2 ) gradually turn into specialized third-stage dauer larva (DL 3 ), so that they can enter winter diapause together with the adults [7][8][9]. All the low-temperature-induced nematodes revealed similar characteristics: extended lifespan and accumulated fat content [8].
Although very little is known about the low temperature-sensing mechanism in PWN, much effort has been invested in probing the functions of some genes that respond to the low temperature 2 of 13 in the model nematode Caenorhabditis elegans. The low temperature response process of C. elegans is regulated by many genes and groups including Cyclic guanosine monophosphate (cGMP) pathway, insulin-signaling pathway [10], G protein signaling [11], and phospholipid saturation [12], among others. The analysis of this widely studied model nematode inspired our analysis for the mechanism study of low temperature sensation in PWN. Some efforts have been made in the study of molecular mechanism of PWN's response to low temperature in our previous research: we proved the key role of patched-related protein gene Bx-DAF6 [13], cGMP pathway [14], and stearoyl-coA desaturase gene Bx-SCD [15]. However, the mechanism of low-temperature-induced lifespan extension of PWN remains elusive.
In this study, we performed RNA sequencing on mixed-stage PWN between ambient temperature treated and low temperature treated groups. The results of RNA sequencing illustrated an essential role of G protein signaling during PWN's early response to the low temperature. G protein coupled receptors (GPCRs) are a functionally diverse superfamily which are known as widely related to the environmental signal detection process [16][17][18]. Thus, we identified six GPCRs genes in G protein signaling from PWN and studied their response pattern to low temperature. We also validated the function of one GPCRs gene in low temperature-sensing with RNAi method. The results reveal that BxGPCR17454 plays a key role in PWN's response to low temperature. This means the BxGPCR17454 in PWN may have the potential to be treated as a target for the management of this dangerous pest.

RNA Sequencing and General Transcription Patterns
We performed RNA sequencing and obtained 312,138,316 high quality clean reads in all six samples after removal of rRNA and low-quality reads ( Table 1). The read counts were converted to FPKM. Dendrogram clustering of the biological duplicates for each RNA sequencing condition reveals conserved alignment between replicates. ( Figure 1A). The proportion of total reads that mapped to the reference genome ranged from 70.47% to 75.61%, and correlation values were significantly higher between the duplicated samples than among the treatments ( Figure 1B). We performed RNA sequencing and obtained 312,138,316 high quality clean reads in all six samples after removal of rRNA and low-quality reads ( Table 1). The read counts were converted to FPKM. Dendrogram clustering of the biological duplicates for each RNA sequencing condition reveals conserved alignment between replicates. ( Figure 1A). The proportion of total reads that mapped to the reference genome ranged from 70.47% to 75.61%, and correlation values were significantly higher between the duplicated samples than among the treatments ( Figure 1B).

B. xylophilus Genes Differentially Expressed in Response to Low Temperature
Between low-temperature-treated and CK nematodes, we identified 432 DEGs of which 314 were up-regulated, while 118 DEGs were down-regulated ( Figure 2A, Table S3). Among all the DEGs, the gene expression patterns were similar between the biological duplicates but significantly different between the low-temperature-treated and CK group. ( Figure 2B). GO enrichment showed that the DEGs mainly distributed in the single-organism process, cellular process, metabolic process, developmental process, and response to stimulus in biological process terms and binding and catalytic activity in molecular function terms ( Figure 2C). Our results also indicated that low temperature influenced the expression of several gene groups with the putative function of environmental signal detection and stress adaptation (Table S3). These genes groups include four cytochrome P450s (BXY_0803400, BXY_0817800, BXY_1185500 and BXY_1697600), six heat shock proteins (BXY_0165400, BXY_0586000, BXY_0640100, BXY_0768000, BXY_1274600 and BXY_1563600) and 16 GPCRs (Table S4), among others.

Validation of DEGs by Quantitative Real-Time Polymerase Chain Reaction (qRT-PCR)
12 DEGs were selected randomly for qRT-PCR validation. In general, the patterns of up-regulation and down-regulation were consistent with those obtained from the RNA sequencing analysis. Thus, qRT-PCR analysis confirmed that the changes detected in the RNA sequencing analysis were reliable ( Figure 3, Table S2).

Identification and Transcript Abundance Analysis of 6 Low-temperature-related BxGPCRs
GPCRs are a group of genes which has been reported as closely related to the environmental signal detection process [16][17][18]. For this reason, we studied the relationship of low temperature stress and BxGPCRs expression. Structurally, GPCRs are also known as 7TM receptors because they have seven transmembrane (TM) domains [18,19]. Thus, we then used conserved domain analysis and transmembrane region prediction to filter the 16 differently expressed genes with GPCRs putative function description. BxGPCRs with incomplete conserved domain (Table S5) or with numbers of TM not equal to seven ( Figure S1) are not further considered. There are six low-temperature-related BxGPCRs left after this filtering. According to these six genes' ID BXY_1375800, BXY_1672700, BXY_0593200, BXY_0521000, BXY_1391500, and BXY_1745400, we named them BxGPCR13758, BxGPCR16727, BxGPCR05932, BxGPCR05210, BxGPCR13915, and BxGPCR17454, respectively. Blastp results showed that the deduced amino acid sequence of these six low-temperature-related BxGPCRs have a relatively high level of identity with other GPCR proteins of several nematodes. On this basis, we downloaded selected homologous amino acid sequences from NCBI (Table S6)

Validation of DEGs by Quantitative Real-Time Polymerase Chain Reaction (qRT-PCR)
12 DEGs were selected randomly for qRT-PCR validation. In general, the patterns of up-regulation and down-regulation were consistent with those obtained from the RNA sequencing analysis. Thus, qRT-PCR analysis confirmed that the changes detected in the RNA sequencing analysis were reliable ( Figure 3, Table S2).

Identification and Transcript Abundance Analysis of 6 Low-temperature-related BxGPCRs
GPCRs are a group of genes which has been reported as closely related to the environmental signal detection process [16][17][18]. For this reason, we studied the relationship of low temperature stress and BxGPCRs expression. Structurally, GPCRs are also known as 7TM receptors because they have seven transmembrane (TM) domains [18,19]. Thus, we then used conserved domain analysis and transmembrane region prediction to filter the 16 differently expressed genes with GPCRs putative function description. BxGPCRs with incomplete conserved domain (Table S5) or with numbers of TM not equal to seven ( Figure S1) are not further considered. There are six low-temperature-related BxGPCRs left after this filtering. According to these six genes' ID BXY_1375800, BXY_1672700, BXY_0593200, BXY_0521000, BXY_1391500, and BXY_1745400, we named them BxGPCR13758, BxGPCR16727, BxGPCR05932, BxGPCR05210, BxGPCR13915, and BxGPCR17454, respectively. Blastp results showed that the deduced amino acid sequence of these six low-temperature-related BxGPCRs have a relatively high level of identity with other GPCR proteins of several nematodes. On this basis, we downloaded selected homologous amino acid sequences from NCBI (Table S6), performed alignment of the deduced amino acid sequences ( Figure  S2) and constructed phylogenetic tree ( Figure 4A).
To validate the transcript pattern of BxGPCRs under low temperature, we measured transcript levels of six BxGPCRs at 5 °C for 1, 3, 5, and 7 days, and 25 °C for 1, 3, 5, and 7 days, respectively, with qRT-PCR ( Figure 4B). The result showed that all of 6 BxGPCRs genes revealed higher transcript levels at low temperature than regular temperature over 7 days. The gene expression culminates at 3 days.

RNAi Validation of Low-temperature-related BxGPCR
The FITC treated nematodes reflected green fluorescence after exposed under ultraviolet light ( Figure 5A), which indicated that the dsRNA can be fully absorbed by nematodes with soaking method. According to the result of expression changes between BxGPCR17454 dsRNA-treated and dsRNA-free nematodes, BxGPCR17454 gene was significantly silenced after dsRNA treatment. The mean expression level of BxGPCR17454 was decreased to 11.85% compared to the dsRNA-free nematodes. However, the dsRNA had no obvious effect on the transcript level of another internal To validate the transcript pattern of BxGPCRs under low temperature, we measured transcript levels of six BxGPCRs at 5 • C for 1, 3, 5, and 7 days, and 25 • C for 1, 3, 5, and 7 days, respectively, with qRT-PCR ( Figure 4B). The result showed that all of 6 BxGPCRs genes revealed higher transcript levels at low temperature than regular temperature over 7 days. The gene expression culminates at 3 days.

RNAi Validation of Low-temperature-related BxGPCR
The FITC treated nematodes reflected green fluorescence after exposed under ultraviolet light ( Figure 5A), which indicated that the dsRNA can be fully absorbed by nematodes with soaking method. According to the result of expression changes between BxGPCR17454 dsRNA-treated and dsRNA-free nematodes, BxGPCR17454 gene was significantly silenced after dsRNA treatment. The mean expression level of BxGPCR17454 was decreased to 11.85% compared to the dsRNA-free nematodes. However, the dsRNA had no obvious effect on the transcript level of another internal control gene β-actin ( Figure 5B). The alignment of six BxGPCRs nucleic sequences also illustrated that RNAi primer sequences is unique for BxGPCR17454 ( Figure S3). These results indicated that the BxGPCR17454 RNAi was potent and specific in this study. control gene β-actin ( Figure 5B). The alignment of six BxGPCRs nucleic sequences also illustrated that RNAi primer sequences is unique for BxGPCR17454 ( Figure S3). These results indicated that the BxGPCR17454 RNAi was potent and specific in this study.
The results of the survival rates calculation showed that B. xylophilus revealed a significantly decreased survival rate after RNAi of BxGPCR17454 under low temperature compared to dsRNA-free group over 30 days ( Figure 5C). This indicated that gene silencing of BxGPCR17454 can significantly decrease the survival rate of B. xylophilus at low temperature. BxGPCRs plays a key role in the low temperature response process of B. xylophilus. Although the transcript level of BxGPCR17454 was evaluated in the first 7 days according to the qRT-PCR results, there are no differences between the survival rate of dsRNA-free group and dsRNA-treated groups in the first 15 days. This may be because there are other downstream genes are influenced by the RNAi of BxGPCR17454. All these genes together with BxGPCR17454 influenced the survival rate of B. xylophilus under low temperature. Further studies need to be carried out to detect the downstream genes of BxGPCR17454 in the future.  . (B). BxGPCR17454 RNAi efficiency was validated by expression changes between BxGPCR17454 dsRNA-treated and dsRNA-free nematodes. BxGPCR17454 gene was significant silenced after dsRNA treatment, while dsRNA had no obvious effect on the transcript level of internal control gene β-actin. (C). BxGPCR17454 dsRNA-treated nematodes showed significate decreased survival rate under low temperature. Data represent mean values ± SD from different repetitions. Asterisks indicates statistically significant differences (* p < 0.05, ** p < 0.01, Student's t-test) was found between the dsRNA-treated and dsRNA-free groups. The results of the survival rates calculation showed that B. xylophilus revealed a significantly decreased survival rate after RNAi of BxGPCR17454 under low temperature compared to dsRNA-free group over 30 days ( Figure 5C). This indicated that gene silencing of BxGPCR17454 can significantly decrease the survival rate of B. xylophilus at low temperature. BxGPCRs plays a key role in the low temperature response process of B. xylophilus. Although the transcript level of BxGPCR17454 was evaluated in the first 7 days according to the qRT-PCR results, there are no differences between the survival rate of dsRNA-free group and dsRNA-treated groups in the first 15 days. This may be because there are other downstream genes are influenced by the RNAi of BxGPCR17454. All these genes together with BxGPCR17454 influenced the survival rate of B. xylophilus under low temperature. Further studies need to be carried out to detect the downstream genes of BxGPCR17454 in the future.

Discussion
It is commonly known that low temperature response processes of both homeotherms and poikilotherms are not passive thermodynamic process but active genetic-promoted processes [20]. In this study, our results coincided with this theory. Transcriptome analysis revealed that low temperature response process of this plant parasite nematode is also a genetic process regulated by many genes. The past 15 years has witnessed a rapid progress in the understanding of how low temperature influenced the genetic response patterns of many plant pests with RNA sequencing method [21][22][23][24][25][26]. However, little is known about how B. xylophilus, which is also one of the most destructive pine tree pests, genetically responded to low temperature.
The identification and characterization of genes involved in low temperature responses is essential to elucidate low temperature defense mechanisms and develop effective control strategies. Changes in the transcriptome of B. xylophilus during low temperature directly reflect the impact of low temperature on the genetic activities. In this article, we conducted a comprehensive transcriptome analysis and characterized the gene expression profiles of B. xylophilus under low temperatures. Through the analysis of DEGs sets, we identified transcriptome changes in B. xylophilus in response to low temperature and validated the transcriptome result with qRT-PCR. The DEG analysis revealed that many gene groups including heat shock proteins, cytochrome P450s and GPCRs, among others, may potentially have played a crucial role in low temperature response process of B. xylophilus. On this basis, six low-temperature-related BxGPCRs genes were identified by phylogenetic analysis. In addition, we found these six low-temperature-related BxGPCRs genes revealed higher transcript abundance under low temperature than ambient temperature. Furthermore, the RNAi method was utilized to study the functions of one of the low-temperature-related BxGPCRs which we named BxGPCR17454. The results indicated that interference of BxGPCR17454 can significantly reduce the survival rate of B. xylophilus under low temperature. Taking all these factors into consideration, we hypothesized that G protein signaling may play a crucial role for the low temperature response process in B. xylophilus. This is widely consistent with previous research results in other organisms such as nematodes [10,11,27], bacterium [28], plants [29,30], fish [31], and mammals [32][33][34]. Interference of BxGPCR17451 may represent a novel approach for the prevention of this dangerous plant pest. Apart from G protein signaling, others DEGs such as heat shock proteins, are widely investigated as a big gene group regarding their role in thermal stress response, development and lifespan regulation [35][36][37][38]. Cytochrome P450s are not only reported as an import detoxification gene group [39], but also thermal-stress-related genes [40,41]. Further validations are expected to be made in these gene groups in the future. Other gene groups and pathways that have been reported to have a low temperature response function in model organism C. elegans include TRP channels [20,42,43], the cAMP pathway [44], and fatty acid metabolism pathway [45,46]; their orthologous genes in B. xylophilus may also share similar low temperature response functions. Further efforts are also expected to be made in the relevant fields in the future.
This study focused on the molecular response pattern of B. xylophilus under low temperature, identification and response pattern validation of six low-temperature-related BxGPCRs genes as well as functional validation of BxGPCR17454. The results shaped an important role of a G protein signaling gene BxGPCR17454 and potential value of other low-temperature-related genes in the low temperature response process of B. xylophilus. These discoveries would assist the development of management and methods for efficient control of this pine tree pest.

Sample Preparation
B. xylophilus, maintained in the Forestry Protection Laboratory of Northeast Forestry University, Harbin, China, were kindly provided by the Chinese Academy of Forestry, Beijing, China. The nematodes were cultured on Botrytis cinerea for 5-7 days at 25 • C. The Baermann funnel technique was used to isolate B.xylophilus from potato dextrose agar medium plates. About 60,000 nematodes were mixed evenly in 6 mL distilled water. Then, suspensions were separated into six 1.5 mL centrifuge tubes equally. Three tubes were cultured at 5 • C for 24 h as treated group. The other three tubes were cultured at 25 • C for 24 h as CK group. After the incubation was finished, all samples were collected in the bottom of the tubes with centrifuge. Then, the tubes were all transferred in liquid nitrogen immediately for the next analysis.

RNA Sequencing
Total RNA of above six samples was extracted with Trizol reagent (Thermo Fisher Scientific, Shanghai, China). RNA-seq library was then constructed with NEB #7530 Kit (New England Biolabs, Ipswich, MA, USA) as following method: mRNA was enriched by Oligo(dT) beads. First, the enriched mRNA was fragmented into short fragments using NEBNext First Strand Synthesis Reaction Buffer. Random primers together with ProtoScript II Reverse Transcriptase and Marine RNase Inhibitor were then used to reverse transcript mRNA into cDNA with following procedure: 25

Sequencing Data Analysis
Raw reads containing adapters, containing more than 10% of unknown nucleotides or containing more than 50% of low quality (Q-value ≤ 20) bases were removed by fastp [47] (Plant Pathology, SCRI, Invergowrie, Dundee DD2 5DA, UK. version 0.18.0) in order to get high quality clean reads. Short reads alignment tool Bowtie2 [48] was used for mapping reads to ribosome RNA (rRNA) database with default parameters. The rRNA mapped reads were removed. The remaining rRNA removed reads of each sample were then mapped to reference genome by TopHat2 [49] (Center for Bioinformatics and Computational Biology, University of Maryland, College Park, MD, USA. version 2.0.3.12). The mapping options used with TopHat are as follows: maximum read mismatch is 2, the distance between mate-pair reads is 50bp and the error of distance between mate-pair reads is ±80bp. The reference genome of B. xylophilus which was submitted directly to WormBase database by Taisei Kikuchi in 2011 [7] was downloaded from WormBase Parasite (BioProject PRJEA64437). The reference genome has 17,704 transcripts in total. We used RSEM software [50] (Department of Computer Sciences, University of Wisconsin-Madison, Madison, WI, USA) to calculate the gene expression level which was normalized by using FPKM (Fragments Per Kilobase of transcript per Million mapped reads) method. The edgeR package (Bioconductor, Roswell Park Cancer Institute, Buffalo NY, USA. http://www.r-project.org/) was used to identify differentially expressed genes with raw counts from RSEM across two groups. We identified genes with a fold change ≥2 and a false discovery rate (FDR) <0.05 in a comparison as significant differentially expressed genes (DEGs). BLASTX against the NCBInr database and Gene Ontology database was used to analysis gene description and annotations. Clustering, Pearson correlation analysis, volcano plot, heatmap, and GO categories are performed by R (Bioconductor, Roswell Park Cancer Institute, Buffalo NY, USA. version 3.2.1).

Validation of DEGs by Quantitative Real-Time Polymerase Chain Reaction (qRT-PCR)
Quantitative real-time Polymerase Chain Reaction (qRT-PCR) was performed to validate the expression of DEGs. Gene-specific primers (Table S1) for 12 randomly selected DEGs were designed by Sangon Online Primer Designer (available online: https://www.sangon.com/newPrimerDesign) (Sangon Biotech Co. Shanghai, China). The Bx28S gene was used as the internal control gene (Table S1). Double-stranded cDNA was obtained using GoTaq 2-Step RT-qPCR System (Promega, Madison, WI, USA) according to instructions of the manufacturer. We used the cycle threshold data to calculate the fold of relative transcript abundance (treated group/CK group). The PCR program used in this analysis was first step: 95 for 2 min, second step: 95 for 15 s and 60 for 1 min, in 40 cycles.

Identification of Low-temperature-related BxGPCRs
According to the analyzed sequencing data, DEGs with a NR description of GPCRs were selected as candidate low-temperature-related BxGPCRs. Candidate genes with incomplete conserve domain were removed. The candidate genes whose predicted transmembrane domains unequal to seven were also removed, considering seven transmembrane domains is a typical feature of GPCRs [51,52]. Conserved domains were analyzed by NCBI Batch CD (available online: https: //www.ncbi.nlm.nih.gov/Structure/bwrpsb/bwrpsb.cgi) and transmembrane domains were predicted by the online tool TMHMM from the Center for Biological Sequence website (available online: http://www.cbs.dtu.dk/services/TMHMM/). BxGPCRs homologous sequences of the other organisms were obtained from BLASTP (available online: https://blast.ncbi.nlm.nih.gov/Blast.cgi). Alignments of deduced amino acid sequences were performed by Muscle (EBI, UK) [53] and visualized by Jalview (version 2. 10. 4b 1). Phylogenetic analyses were performed with Mega 7.0 software using a maximum likelihood tree.

Transcript Abundance Analysis of 6 BxGPCRs under Low Temperature
The transcript abundance was analyzed as the method described before [14]. About 80,000 nematodes were mixed evenly in 8 mL distilled water. Then we separated the suspensions into eight 1.5 mL centrifuge tubes equally. Four tubes of nematodes were cultured at 5 • C for 1, 3, 5, and 7 days. The other four tubes of nematodes were cultured at 25 • C for 1, 3, 5, and 7 days. Total RNA of the nematodes cultured in these eight tubes were extracted and synthesized into cDNA with the method as described above after the incubation was finished. The transcript abundance and qRT-PCR Primers (Table S2) were also calculated and designed as described above. The whole transcript abundance analysis was analyzed with three replicates as three independent trials.

RNAi Validation of Low-Temperature-Related BxGPCR
We utilized RNAi method to study the function of one of the low-temperature-related BxGPCRs BxGPCR17454. The nematode RNAi was performed with soaking method as outlined in previous papers [14,39,[54][55][56]. Mixed-stage B. xylophilus (male: female: juvenile ratio is approximately 1:1:2) was used in this research. MAXIscript T7/T3 RNA Synthesis Kit (Ambion, Tokyo, Japan) was used to obtain dsRNA with the following primers: i-BxGPCR17454-F (GCT AAT ACG ACT CAC TAT AGG GGG AAT CAT TGT ATA CCA GCC AAT TGC) and i-BxGPCR17454-R (AGT AAT ACG ACT CAC TAT AGG GTG GCC CGA TTG AAA TCC TGC). The size of BxGPCR17454 dsRNA for RNAi is 350 bp. The dsRNA-treated nematodes were soaked in ddH 2 O with final concentration of 2 mg/mL dsRNA corresponding to BxGPCR17454 sequence for 24 h at 25 • C to uptake the dsRNA. The uptake of the dsRNA was monitored by a final nematode treatment of ddH 2 O containing 1 mg/mL FITC (fluorescein isothiocyanate) for 24 h. The CK nematodes were soaked in ddH 2 O only. After intermittent stirring for 24 h at 25 • C, the nematodes were washed with ddH 2 O to remove the external FITC and dsRNA. A fluorescence microscope was then used to take pictures of the nematode (ZEISS, Germany) to detect the uptake of FITC. Nematodes were divided into two groups after dsRNA uptake: the first group was used to examine the efficiency of RNAi with qPCR, and the second one was used to calculate survival rate under 5 • C environments. Total RNA was extracted from the CK nematodes and the dsRNA-treated nematodes after dsRNA. qRT-PCR was then performed as mentioned above with qRT-PCR primers of BxGPCR17454 and internal control primers Bx28S and β-actin [57] (Table S2). Survival rates of nematodes were calculated with 100 nematodes in 1 mL of ddH 2 O within each 1.5 mL centrifuge tube performed in triplicate. Then the survival rate of nematodes in each centrifuge tube was monitored every 3 days over 30 days. The survival rates of each tube were calculated as follows: number of living nematodes in each tube / number of total nematodes in each tube * 100%.