Selection and Validation of Suitable Reference Genes for RT-qPCR Analysis in Bursaphelenchus Xylophilus

Quantitative reverse transcription polymerase chain reaction (RT-qPCR) is a powerful technique for studying gene expression, and it is widely used in molecular biology and genetic research. The key to quantitative accuracy depends on the stability of the reference genes used for data normalization under different experimental conditions. The pinewood nematode (PWN; Bursaphelenchus xylophilus) is the causal agent of a devastating disease, the pine wood disease (PWD), demanding extensive and prompt research to understand the molecular mechanism of PWD, but the identication of reference PWN genes for standardized RT-qPCR has not been reported yet. In this study, we have analyzed eight candidate reference genes of PWN under different temperature conditions and at different developmental stages. Delta Ct method, GeNorm, NormFinder, BestKeeper and RefFinder algorithms were used to evaluate the expression stability of these genes. 18SrRNA and HIS were selected for gene expression research under temperature treatments, while EF1γ and 18SrRNA for gene expression studies across different developmental stages. In general, the results indicate that 18SrRNA is the most stable gene for RT-qPCR under different conditions. Finally, we use arginine kinase gene (AK) in different temperature and heat shock protein 90 (HSP90) in different developmental stages to conrm the stability of expression. The systematic analysis of qRT-PCR reference gene selection of B. xylophilus will be helpful for future functional analysis and exploration of genetic resources of B. xylophilus.


Introduction
Pine wilt disease (PWD) has devastated a serious threat to forest safety and ecosystem stability in parts of Asia and Europe and has been quarantined globally 1 #64}. In China, the disease is caused by an invasive nematode species, the pine wood nematode (PWN), Bursaphelenchus xylophilus.
The destructiveness of PWN is usually closely related to its native vectors beetles and ophiostomatoid fungal species 3 . PWN feeding on the epithelial tissue of the tree host and ophiostomatoid (blue -stain) fungi. The fungi reproduce on the damaged tree host 4,5 . PWN undergoes four propagative larval stages (L1-L4) from development to reproduction 6 .
However, under adverse conditions, such as scarce food and extreme temperatures , PWN will molt from L 2 into dispersed third juveniles (L III ) and enter a dispersal stage of its life cycle 7 .
When the environmental conditions are suitable, L III are attracted by the insect vector, Monochamus alternatus, and gather around their pupal chambers 2,5 . Then PWN makes the development of the fourth juveniles (L IV ) dispersal stages synchronized with the late pupal or early adult phase of vector beetles, in order to these L IV can enter the tracheal system of M. alternatus and spread to new pine hosts 2,5 . Since the pathogenesis of B. xylophilus is complex and there are strong interactions between vector beetles, For the control of this destructive disease, a lot of effort has been spent on the control of pine wood nematodes, including the chemical control being one of the most common approach. However, the indiscriminate in addition to posing environmental hazards, have led to the emergence of pesticide resistance in PWN 8 . Developing environmentally friendly pest control methods, such as RNAi or nematophagous fungal formulations, is a necessary alternative solution. Considering the economic importance and the urgency of developing effective and e cient control strategies, further ecological, physiological, and molecular biological explorations on B. xylophilus has been extensively conducted over the past decades [9][10][11][12][13][14] . In the continual research process, further studies involving various aspects of the molecular mechanism of pine wood nematodes, including development, immunity, energy metabolism, reproduction, chemoreception, and other biological process are required [15][16][17][18] . Analysis of the regulation and expression patterns of functional genes is one of the major approaches in such explorations.
Quantitative reverse transcription polymerase chain reaction (RT-qPCR) has been widely used for gene expression analysis owing to its high speci city, rapidity, reproducibility and low costs [19][20][21] . An indispensable factor for the RT-qPCR method is to select a reference gene with stable expression and use it as a reference to reduce the differences between samples through analysis and calculation to ensure the con dence of the nal result 22 .  35 . The washed isolates of nematodes were poured into a glass dish and allowed to stand for half an hour. When the nematode eggs got settled at the base, the excess liquid was drained off and PBST containing streptomycin sulfate and chloramphenicol were added and the glass dishes were incubated for 24h at 25℃ to hatch the eggs and get L 2 larvae.
Collection of third-stages reproductive juveniles (L3). The L 2 nematodes were cultured on B. cinerea -PDA culture plates and incubated for 18 hours at 25℃. The L 3 nematodes were isolated from the cultures using the Baermann funnel technique and cleaned as described earlier followed by freezing in liquid nitrogen and storage at -80℃ for later use.
Collection of fourth-stages reproductive juveniles (L4). The earlier prepared L 2 nematode cultures were incubated for 42 hours to get them developed into fourth-stage reproductive juveniles (L 4 ). Once isolated and cleaned as described earlier, the L 4 nematodes were immediately frozen in liquid nitrogen and stored at -80℃.
Collection of third-stages dispersal juveniles (Lш). The wild wilted pine wood was cut into 35 cm thin strips, and the nematodes were separated by the Baermann funnel technique 36 , whereas any of the non-Lш nematode isolated, were removed during microscopic observations. To further purify the isolates, sucrose otation centrifugation method was used to remove fungal hyphae and impurities 37 . These puri ed supernatants containing Lш nematodes were frozen in liquid nitrogen and stored at -80℃ for later use.
Collection of fourth-stages dispersal juveniles (LIV). The newly eclosing M. alternatus beetles were dissected and placed in a 35 mm Petri dish with sterilized ddH 2 O for 30 min. After half an hour, a large number of L IV nematode were clearly visible oating in the solution. these were isolated and cleaned before freezing in liquid nitrogen and storage at -80℃.
Collection of reproductive adult. L 2 nematodes as collected earlier, were inoculated to B. cinerea cultures and kept for 72 hours. Once the larval stages developed in to adults, nematodes were isolated by Baermann funnel method and cleaned as previously described. The adult nematodes were stored in 1.5 mL centrifuge tubes and stored at -80℃.
Collection of nematodes treated with different temperatures. To study the effect of different temperatures, randomly collected nematodes from B. cinerea culture plates were isolated and cleaned as described earlier. These cultures were exposed to 4℃, 10℃, 15℃, 20℃, 25℃, 30℃ and 35℃ for 48 hours, respectively. Another dish of random culture of nematodes was incubated at 50℃ for 20 min and then shifted at 25℃ for next 24 hours. Each temperature treatment was replicated thrice and after incubation all samples were stored in liquid nitrogen at -80℃. Assessing the expression stability of reference genes. The expression stabilities of the candidate reference genes were estimated and ranked using four different statistical algorithms, geNorm, NormFinder, BestKeeper, and the delta Ct method, and a web-based analysis tool, RefFinder. The geNorm evaluates the most stable reference genes and determines the optimal number of reference genes required for normalization 31 . When the expression stability value (MV) is less than 1.5, the expression of this gene is stable 40 . The NormFinder ranks candidate reference genes by calculating their stability value (SV) 32 . The BestKeeper calculates the stability of the candidate genes based on the standard deviation (SD) and coe cient of variation (CV), the smaller SD and CV, the more suitable they are as reference genes 30 . The comparative delta Ct method calculates the delta Ct of the candidate reference genes, and smaller ΔCt means more stable gene expression 29 . Finally, we comprehensively ranked the candidate reference genes based on the geomean (GM) values of the above results from the four different statistical algorithms, using the web-based analysis tool RefFinder 33 .
Validation of reference genes. To con rm the stability of expression of the selected reference genes, veri cation experiments were carried out in samples from different temperature and different developmental stages. For the developmental stage studies, we quanti ed the relative expression of arginine kinase genes (AK), which have been shown to be differentially expressed at different developmental stages 41,42 . For the different temperatures studies, we quanti ed the relative expression of the thermotolerance-related gene heat shock protein 90 (HSP90), which is essential for the survival of B. xylophilus 43

Results
Performance of RT-qPCR primers. In order to determine the availability of reference genes, rst verify the speci city and e ciency of its ampli cation primers. The 2.0% agarose gel image showed that the PCR ampli cations of each primer pair had a single expected size (Fig. 1A), and the melting curves showed a single peak ( Fig. 1B; Fig.S1). The primer e ciency (E) was between 106.43% and 219.36%, and the linear regression coe cient (R 2 ) values were between 0.983 and 0.998 (Table 1). The comprehensive results showed that the designed primers could accurately amplify the candidate reference genes.
Expression pro les of candidate reference genes. We examined the expression pro les of the candidate reference genes to have an overall representation of primer variability under different experimental conditions (Fig. 2). All the Ct values of the candidate reference genes were between 13 and 26, which indicated that the eight reference genes expressed high level under different experimental conditions. PMP-2 and HIS had average Ct values >23 cycles, followed by UBCE, TUB, EF1γ, UBQ, Actin and 18SrRNA whose average Ct values ranged between 13-23 cycles (Fig. 2). Meanwhile, it was obviously that the degree of variation in the expression of some reference genes varies with experimental conditions. For example, UBQ varied less (~1 cycle) between temperature treatments than across developmental stages (>4 cycles).
Stability of candidate reference genes. To further evaluate the expression stability of the eight candidate reference genes under two experimental conditions (temperature and developmental stage), ΔCt method, BestKeeper, NormFinder, and geNorm were employed, and then RefFinder was used to calculate an overall stability ranking.
For different temperature treatments, the analysis of the result of ΔCt method, BestKeeper, and NormFinder obviously showed that HIS, 18SrRNA and TUB were the most stable reference genes ( Table  2). The geNorm showed that all eight candidate genes were suitable for reference genes (MV < 1.5), while 18SrRNA and UBQ were the most stable genes (Fig. 3A). EF1γ and PMP-2 were the most unstable reference genes in the results of the four statistical algorithms. The stability of genes from most to least stable ranked by RefFinder was HIS>18SrRNA>TUB>UBQ>UBCE>Actin>PMP-2>EF1γ (Fig. 4A). Pairwise variations values of the expression stability of candidate reference genes calculated by geNorm was Vvalue < 0.15 at V2/3, which indicated that only two references genes were required to normalize the target gene data (Fig. 5). Thus, the best combination of reference genes for temperature treatments of B. xylophilus was HIS and 18SrRNA.
For the different developmental stages, the ΔCt method and NormFinder represented that 18SrRNA and EF1γ were the most stable genes, whereas TUB and HIS exhibited the greatest variation ( Table 2). Based on BestKeeper, 18SrRNA and Actin were the most stable reference genes (Table 2). Whereas, the analysis result of BestKeeper indicated that HIS and UBQ were relatively unstable (SD > 1) ( Table 2). Additionally, GeNorm calculated the lowest M value for the 18SrRNA/EF1γ pair, suggesting that they are the most stable reference genes (Fig. 3B). RefFinder ranked the genes from most to least stable as: EF1γ>18SrRNA>UBQ>actin>UBCE>PMP-2>TUB>HIS (Fig. 4B). The geNorm pairwise variation analysis also revealed that the rst V-value <0.15 emerged at V2/3 (Fig. 5). Thus, the best combination of reference genes for developmental stages samples of B.xylophilus was EF1γ and 18SrRNA.
Validation of selected reference genes in B. xylophilus.To validate the reference genes selected in B. xylophilus under different experimental conditions, the expression of HSP90 and AK were checked under different temperature treatments and developmental stages. The proposed stable gene combinations (HIS/18SrRNA or EF1γ/18SrRNA, respectively) (Fig. 6) and the least stable combinations (EF1γ/PMP-2 or TUB/HIS, respectively) ( Fig. 7) were used to normalize the expression levels of those two target genes.
The expression patterns of HSP90 under different temperature treatments were inconsistent when normalized with the two most stable or two least stable reference genes. HSP90 expression was found to be 2.6 fold (p < 0.05) higher at 50℃ than 25℃ when the proposed gene combination (HIS/18SrRNA) was used (Fig. 6A). When the least stable combination (EF1γ/PMP-2) was used, HSP90 expression was only 1.5 fold (p > 0.05) higher at 50℃ than in the controls under the same conditions (Fig. 6B).
At different developmental stages, AK expression was sharply up regulated at L 4 stage and adult stage.
However, when AK expression at L 4 stage was 214 fold higher than at L III stage when it was normalized to the two most stable reference genes, and only 80 fold higher when normalized to the two least stable reference genes (Fig. 7).

Discussion
B. xylophilus a serious pest of pine trees causing Pine wilt disease (PWD), which has been considered as the most serious forestry disaster in China in recent years. Due to the signi cance of PWD and the economic losses occurring due to this pest, many in depth studies are being initiated to understand the disease cycle and most of them involve the qRT-PCR technique to quantify the expression of target genes e ciently and accurately. This technique relies on stable expression of reference genes to eliminate the impact of differences in RNA quality, reverse transcription e ciency, and human factors in different samples on quantitative results. Hence, selection of the most appropriate and stable reference genes is the basis of qRT-PCR success 44  Regarding the stability of eight reference genes selected reported by different analyses, different results were demonstrated (Table 2). It is di cult to identify reference genes which are generally stable across these four algorithms, as each program has its own strengths and appropriate application conditions 29,54,55 . When compared across different temperature treatments, the results of geNorm is quite distinct from those of delta Ct, BestKeeper, and NormFinder (Table 2), nonetheless, there also a partial consistency was noted. As shown in Table 2, the stable combination of reference genes for temperature treatments was HIS and 18SrRNA, whereas EF1γ and PMP-2 were the most unstable reference genes in the results of the four statistical algorithms. In the analysis of the developmental stages, the ΔCt method and NormFinder represented that 18SrRNA and EF1γ were the most stable genes. Based on BestKeeper, 18SrRNA and Actin were the appropriate reference genes (Table 2).
Therefore, when diverse methods are used, a comprehensive ranking must be made based on the results of different methods. In the latest research, in order to overcome the differences in the results of different algorithms and obtain the nal ranking, the web-based analysis tool RefFinder was used. RefFinder obtained the comprehensive overall ranking based on the Geomean of the ranking values from abovementioned algorithms 56 . However, one disadvantage of RefFinder is that the results of these algorithms are not weighted according to the unavailability of their cut-offs and appropriate weights 57 .
Based on the results from the comprehensive analysis using the ve algorithms, HIS and 18SrRNA, EF1γ and 18SrRNA are pairs of the reliable combination, which will facilitate the molecular mechanisms conferring functional and development genomics studies of the nematode in future studies.
Further, according to the results of gene expression studies, we found that traditional reference genes are not always stable and consistent with earlier studies. Many studies have shown that some common reference genes have different expression stability under different experimental conditions 52,[58][59][60][61] . This study further emphasizes the shortcomings of simple use of commonly reference genes without a comprehensive evaluation (Table 2; Fig. 4). For example, the commonly used internal reference gene ACT has unstable expression under different temperature treatments and different developmental stages, especially different temperature treatments. Similarly, HIS gene expression is stable under different temperature treatments, but extremely unstable at different developmental stages ( Table 2; Fig. 2). These data indicate that under different conditions, the expression levels of traditional reference genes may vary signi cantly, which requires that all reference genes should be validated.

Conclusions
Under a series of experimental conditions, we identi ed suitable and reliable reference genes for the normalization of gene expression in B. xylophilus in conclusion, 18SrRNA and HIS can give stable gene expressions under temperatures while, EF1γ and 18SrRNA were found a good selection as reference gene in expression studies across different developmental stages. This study identi ed several suitable internal reference genes for the standardized RT-qPCR analysis of pine wood nematode, which will contribute to further explorations of molecular mechanism and functional genomics of pine wood nematode with pest control implications.