Transcriptome Analysis of the Japanese Pine Sawyer Beetle, Monochamus alternatus, Infected with the Entomopathogenic Fungus Metarhizium anisopliae JEF-197

The Japanese pine sawyer (JPS) beetle, Monochamus alternatus Hope (Coleoptera: Cerambycidae), damages pine trees and transmits the pine wilt nematode, Bursaphelenchus xylophilus Nickle. Chemical agents have been used to control JPS beetle, but due to various issues, efforts are being made to replace these chemical agents with entomopathogenic fungi. We investigated the expression of immune-related genes in JPS beetle in response to infection with JEF-197, a Metarhizium anisopliae isolate, using RNA-seq. RNA samples were obtained from JEF-197, JPS adults treated with JEF-197, and non-treated JPS adults on the 8th day after fungal treatment, and RNA-seq was performed using Illumina sequencing. JPS beetle transcriptome was assembled de novo and differentially expressed gene (DEG) analysis was performed. There were 719 and 1953 up- and downregulated unigenes upon JEF-197 infection, respectively. Upregulated contigs included genes involved in RNA transport, ribosome biogenesis in eukaryotes, spliceosome-related genes, and genes involved in immune-related signaling pathways such as the Toll and Imd pathways. Forty-two fungal DEGs related to energy and protein metabolism were upregulated, and genes involved in the stress response were also upregulated in the infected JPS beetles. Together, our results indicate that infection of JPS beetles by JEF-197 induces the expression of immune-related genes.


Introduction
The Japanese pine sawyer (JPS) beetle, Monochamus alternatus Hope (Coleoptera: Cerambycidae), is a major pest that causes serious damage to pine trees by mediating transmission of the pine wilt nematode, Bursaphelenchus xylophilus Nickle (Aphelenchida: Aphelenchoididae) that causes pine wilt disease [1]. This beetle has spread to East Asian countries including Korea, Japan, and China, in addition to Portugal in Europe [2,3]. Disease-causing nematodes or insects that mediate the transmission of these nematodes need to be controlled to prevent the spread of pine wilt disease. The main means of controlling JPS beetles are chemical agents such as neonicotinoids, fenitrothions, thiacloprids, and thiamethoxams [4,5]. However, these chemical agents are toxic to the environment and can induce insect resistance. Therefore, natural microbial materials have been investigated as alternatives to these chemical agents [6,7]. Among natural microbial agents, entomopathogenic fungi have attracted great interest as environmentally sound biopesticides [8]. The entomopathogenic fungus Beauveria bassiana was developed as a protein coding in silico cDNA library of JPS adults was constructed by de novo assembly, and DEGs were analyzed by mapping short reads to the in silico cDNA library. Changes in metabolism were confirmed by analyzing pathways and performing gene ontology (GO) analysis. In addition, changes in the expression levels of JPS adult immune genes upon fungal infection were assessed. Secondly, we confirmed that the expression level of specific fungal genes changed significantly during fungal infection. Candidate genes of M. anisopliae that potentially played a role in killing JPS adults were identified through bioinformatic analyses.

Japanese Pine Sawyer Beetles and the Entomopathogenic Fungus JEF-197
Adult Japanese pine sawyer (JPS) beetles, M. alternatus Hope, were obtained from the insect rearing company Osangkinsect Co. (http://www.osang.com, accessed on 16 April 2021, Korea). Final stage larvae were stored at 10 • C to induce dormancy. Larvae were maintained at 20 • C for 7 days for dormancy awakening and then moved to 25 • C for pupation and emergence. Rearing conditions were 60% relative humidity (RH) and a photoperiod (light:dark) of 16:8. All experiments involved JPS adults within 10 days of emergence.

Treatment of JPS Adults with M. anisopliae JEF-197
One JPS adult was placed in a plastic cup (lid 100 mm diameter × bottom 60 mm diameter × height 150 mm), and then sprayed either with 1 mL of the conidial suspension (1.0 × 10 7 conidia/mL using 0.03% siloxane solution) or 1 mL 0.03% siloxane solution (nontreated control). Plastic cups with JPS adults were maintained at 25 ± 2 • C at 95% RH, and an approximately 10 cm long pine stem was placed in the dish as a food source. Mortalities of JPS adults were assessed on a daily basis. Fifteen JPS adults were included in each treatment group, and experiments were repeated in triplicate. Differences in mortality of JPS adults among groups were analyzed using t-tests in SPSS ver. 19.0 (SPSS Inc., Chicago, IL, USA) using a significance level of 0.05 (α), and LT 50 was calculated by probit analysis.

RNA Extraction and Construction of RNA-Seq Libraries
JEF-197 fungal samples were cultured in 1/4 SDA medium and harvested on day 8 after inoculation of the medium. JPS samples were treated with the fungus as described above, and then three live JPS adults from fungal-treated and non-treated control groups were obtained on the 8th day after treatment, respectively. To analyze changes in the expression of specific immune-related genes over time, additional JPS samples were collected from fungus-treated and non-treated control groups at 6, 24, 48, 72, 96, 120, and 192 h after treatment as described above.
RNA was extracted from all samples on the 8th day after fungal treatment. Individual adult JPS beetles were placed in 15 mL conical tubes (SPL life sciences, Pocheon, Korea) with 5 mL Trizol TM reagent (Molecular Research Center Inc., Cincinnati, OH, USA). Adult JPS beetles were then ground using an iron pestle, and 400 µL of ground insect sample in Trizol reagent was transferred to a 1.5 mL microtube and ground again using an Ultra grinder BTM (Taeshin Bio Science, Namyangju, Korea). Total RNA from fungus-treated and non-treated control JPS adults was extracted using Trizol reagent according to the manufacturer's instructions. The integrity of the extracted RNA was examined using an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). Sequencing libraries of the samples were made using the Truseq RNA kit (Illumina, San Diego, CA, USA) according to the manufacturer's protocol. cDNA library construction and Illumina sequencing of the samples were performed at Macrogen Corporation (Seoul, Korea). First, poly-A containing mRNA molecules were purified using poly-T oligo-attached magnetic beads followed by fragmentation into small pieces under elevated temperature with divalent cations. Cleaved RNA fragments were reverse-transcribed into first strand cDNA with random primers. This was followed by second strand cDNA synthesis using DNA polymerase I and RNase H. These cDNA fragments were subjected to an end repair process, addition of a single dATP, and ligation of sequencing adapters. Products were then purified and enriched by PCR to create the final cDNA library. All samples were sequenced on an Illumina HiSeq2000 sequencer (Illumina, San Diego, CA, USA) to generate high-throughput transcriptome sequence data with an average read length of 101 bp.

De Novo Transcriptome Assembly and Differentially Expressed Gene (DEG) Analysis
Three biological replicates from each group were sequenced, and the qualities of the Illumina short reads were checked using fastQC v.0.11.8 [40]. Short reads were filtered to remove low-quality sequences using the trimmomatic program [41], and JEF-197 and JPS beetle contigs were assembled using the Trinity de novo assembler v.2.8.5 [42] with default options. Assembled contigs with protein coding capacity were identified by TransDecoder v5.5.0 (https://github.com/TransDecoder, accessed on 16 April 2021). The numbers of short reads that mapped to contigs were quantified by Kallisto v.0.45.0 [43] to calculate transcript per million (TPM) values. Among the isoforms of contigs, those with the highest TPM values were used, and those with TPM values of 0 were removed. Filtered contigs were subjected to DEG analysis under the condition of FDR < 0.05 using edgeR [44].

Functional Annotation and Gene Set Enrichment Analysis
JEF-197 and JPS beetle contigs were identified using an E-value threshold of 1 × 10 −5 using the blastx function of Blast2GO v.5.2.5 [45] and the NCBI insect nr database. Immunerelated genes in JPS beetle contigs were identified as described above using a database of immune-related genes of arthropods (immunoDB; http://cegg.unige.ch/Insecta/immunodb, accessed on 16 April 2021). Among the contigs identified in the immunoDB, contigs with e-values lower than 1 × 10 −100 were used to analyze the expression level of immunerelated genes. In addition, genes related to Toll, IMD, and JAK/STAT signaling, which are the major immune signaling pathways in insects, were analyzed, including contigs with E-values higher than 1 × 10 −100 . Annotated contigs encoding immune-related genes were classified by immunological class and log 2 FC value (fungal-treated/non-treated control JPS adults).
Analysis of gene ontology (GO) and KEGG (Kyoto Encyclopedia of Genes and Genomes) pathways was performed by grouping upregulated and downregulated genes. GO analysis of DEGs was performed using the InterProScan function of Blast2GO, which interrogates the CDD, HAMAP, HMMPanther, HMMPfam, HMMPIR, FPrintScan, and BlastProDom databases (https://www.ebi.ac.uk/interpro/, accessed on 16 April 2021). After annotation, we assessed whether upregulated and downregulated contigs belonged to one or more of the following three GO groups: biological process, cellular component, and molecular function at GO level 3. KEGG pathway analysis was performed using the BBH method and insect databases present in the KEGG Automatic Annotation Server (KAAS, http://www.genom.jp/tools/kaas/, accessed on 16 April 2021).
DEGs in fungus-infected JPS beetles were identified based on the Tribolium castaneum protein database (Tcas5.2, Ensembl database) and E-value threshold of 1 × 10 −10 using the blastx function of Blast2GO. DEGs of fungus-infected JPS beetles were classified into eight groups based on up-and downregulated genes (|log 2 FC| < 1, 1-2, 2-3, and >3) and were then subjected to GO enrichment analysis by g:profiler [46] using a Benjamini-Hochberg false discovery rate (FDR) < 0.05 to identify GO terms and KEGG pathways associated with the various groups of genes.

Gene Expression Validation by qRT-PCR
Ten qPCR primers were used in qRT-PCR analyses to validate RNA-seq results. Primers were designed based on the coding regions of contigs using the online program SnapDragon (https://www.flyrnai.org/snapdragon, accessed on 16 April 2021). Sequences of all primers used in this study are listed in Table S1.
The actin gene (JPS_TRINITY_DN629_c0_g1) was used as an amplification control.
Extracted RNA from JPS beetles was used as a template to synthesize cDNA for RT-PCR and qRT-PCR. One microgram of total RNA from each sample was subjected to reverse transcription with an oligo (dT) 15 primer (Promega, MI, USA) using Accupower ® RT PreMix (Bioneer, Daejeon, Korea). RT-PCR was conducted using the following conditions: 94 • C for 3 min, 30 cycles of 94 • C for 30 s, 60 • C for 30 s, and 74 • C for 30 s, and a final extension step at 74 • C for 3 min using Accupower ® PCR PreMix (Bioneer, Daejeon, Korea). PCR products were identified by 1.5% agarose gel electrophoresis (data not shown). Realtime PCR was performed using the Thermo Scientific Verso SYBR Green 1-step qRT-PCR ROX Mix kit (Thermo-Fisher Scientific, Carlsbad, CA, USA) and the 96-well Bio-Rad CFX96 Real-Time PCR System (Bio-Rad, Hercules, CA, USA). Total RNA samples that were not reverse transcribed were used as additional negative controls for PCR. PCR conditions were as follows: 95 • C for 2 min followed by 40 cycles of 95 • C for 5 s and 60 • C for 15 s. The actin gene of JPS beetle was used to normalize the expression level of target genes. Melting curve analysis was performed to assess non-specific amplification. Relative gene expression (fold change) was calculated using the 2 −∆∆Ct method. All experiments were performed in triplicate. Statistical analyses were performed using Student's t-test, and a p-value < 0.05 was considered to indicate a significant difference.

Time Course of Fungal Virulence
JPS adults sprayed with a conidial suspension (1.0 × 10 7 conidia/mL) of JEF-197 had a mortality of 53.3% on the 8th day after fungal inoculation (Figure 1a), and the LT 50 of adults was 7.86 days, which was confirmed by probit analysis (χ 2 = 8.774, df = 10, p = 0.554). On the 10th day after fungal treatment, the mortality of non-treated control and fungus-treated JPS adults was 11.1% and 66.7%, respectively, which was significantly different (t = 9.423, p < 0.001). On the 4th day after death, white mycelia and green conidia were observed on the surface of JPS cadavers (Figure 1b). After the fungal treatment, the time point of LT 50 in the infected JPS beetles was determined to extract RNAs and analyze gene expression to figure out the response of beetle when the fungus was actively infecting.
In this work, we analyzed gene expression in JPS beetles and the fungus at the same time points of LT 50 because the fungus in the infected insect needed to grow enough for analysis. In the results, significant responses of JPS genes against fungal infection were identified; however, the fungal gene expression profile was relatively simple, probably due to the small amounts of fungal transcripts obtained from the infected JPS. Because the immune response is induced at the onset of infection, the immunity of insects at the LT 50 reflects the late stage of infection. However, it was identified for the persistent immune response of insects against fungal infection. These are described in Sections 3.4 and 3.6. In order to analyze the reaction to the fungi more clearly, an improvement in the method of the experiment was needed. Previous studies have used transcriptome analysis to investigate interactions based on sequencing of RNA extracted from the pathogen and host at different time points [47]. However, using dual RNA-seq, it is possible to simultaneously measure multiple transcripts without physically separating cells [48]. Individual studies are needed in consideration of the timing of analysis of fungi and insects. the experiment was needed. Previous studies have used transcriptome analysis to investigate interactions based on sequencing of RNA extracted from the pathogen and host at different time points [47]. However, using dual RNA-seq, it is possible to simultaneously measure multiple transcripts without physically separating cells [48]. Individual studies are needed in consideration of the timing of analysis of fungi and insects.

Construction of in Silico cDNA Libraries of M. anisopliae JEF-197 and JAPANESE Pine Sawyer Beetle
A total of 55,816,117 raw RNA-seq reads were obtained from JEF-197 (NCBI SRA accession# PRJNA691966) while 51,085,191 and 50,416,917 raw reads were obtained from negative control and fungus-inoculated JPS beetles (NCBI SRA accession# PRJNA691967), respectively (Table S2). The quality of raw reads was assessed using the fastQC program. The de novo assembled transcript sequences of JEF-197 and JPS adults were filtered to select contigs with protein coding capacity; 35,334 and 52,306 contigs including 9286 and 19,046 unigenes were obtained from JEF-197 and JPS beetles, respectively. The N50 of JEF-197 was 1500 bp while that of JPS beetles was 2007 bp (Table S3). JEF-197 and JPS beetles had 59.1% and 65.9% contigs with a contig length less than 1 kb, respectively ( Figure S1).

Differentially Expressed Genes in Japanese Pine Sawyer Beetle as a Result of Fungal Infection
To accurately measure gene expression levels, only contigs with the highest TPM were retained, and genes with a TPM value of 0 were removed. Finally, 14,156 contigs were filtered for DEG analysis. Illumina short reads of negative control and fungustreated JPS beetles were mapped to the 14,156 JPS beetle unigenes. Gene expression levels of 2672 JPS beetle unigenes showed significant differences, and there were 719 and 1953

Construction of in Silico cDNA Libraries of M. anisopliae JEF-197 and JAPANESE Pine Sawyer Beetle
A total of 55,816,117 raw RNA-seq reads were obtained from JEF-197 (NCBI SRA accession# PRJNA691966) while 51,085,191 and 50,416,917 raw reads were obtained from negative control and fungus-inoculated JPS beetles (NCBI SRA accession# PRJNA691967), respectively (Table S2). The quality of raw reads was assessed using the fastQC program. The de novo assembled transcript sequences of JEF-197 and JPS adults were filtered to select contigs with protein coding capacity; 35,334 and 52,306 contigs including 9286 and 19,046 unigenes were obtained from JEF-197 and JPS beetles, respectively. The N50 of JEF-197 was 1500 bp while that of JPS beetles was 2007 bp (Table S3). JEF-197 and JPS beetles had 59.1% and 65.9% contigs with a contig length less than 1 kb, respectively ( Figure S1).

Differentially Expressed Genes in Japanese Pine Sawyer Beetle as a Result of Fungal Infection
To accurately measure gene expression levels, only contigs with the highest TPM were retained, and genes with a TPM value of 0 were removed. Finally, 14,156 contigs were filtered for DEG analysis. Illumina short reads of negative control and fungus-treated JPS beetles were mapped to the 14,156 JPS beetle unigenes. Gene expression levels of 2672 JPS beetle unigenes showed significant differences, and there were 719 and 1953 unigenes with increased and decreased expression, respectively. The number of unigenes that were between 2-and 8-fold upregulated or downregulated genes was 614 (85.4%) and 1500 (76.8%), respectively (Figure 2a). Gene expression levels were compared using the condition of FDR < 0.05 in edgeR. There were more unigenes with decreased expression levels than with increased expression levels (Figure 2b,c). To validate the RNA-seq results, the expression levels of seven JPS beetle genes were assessed using qRT-PCR. Gene expression levels measured by qRT-PCR were in a good accordance with those inferred based on RNA-seq ( Figure S2). unigenes with increased and decreased expression, respectively. The number of unigenes that were between 2-and 8-fold upregulated or downregulated genes was 614 (85.4%) and 1500 (76.8%), respectively (Figure 2a). Gene expression levels were compared using the condition of FDR < 0.05 in edgeR. There were more unigenes with decreased expression levels than with increased expression levels (Figure 2b,c). To validate the RNA-seq results, the expression levels of seven JPS beetle genes were assessed using qRT-PCR. Gene expression levels measured by qRT-PCR were in a good accordance with those inferred based on RNA-seq ( Figure S2).

Changes in Japanese Pine Sawyer Beetle Gene Expression after Fungal Treatment
Functions of the 2672 DEGs were annotated by InterProScan based on the EMBL database. A total of 578 (48.3%) upregulated and 205 (24.5%) downregulated JPS beetle contigs were annotated with GO functions. As a result, a total 246 and 239 genes were annotated with the GO term of biological process, while 145 and 90 genes were annotated with the GO term cellular component function and 327 and 363 genes with the GO term molecular function, respectively ( Figure 3). There were more upregulated unigenes than downregulated unigenes annotated as being involved in biological processes. In particular, many upregulated genes were annotated as being involved in cellular component organ-

Changes in Japanese Pine Sawyer Beetle Gene Expression after Fungal Treatment
Functions of the 2672 DEGs were annotated by InterProScan based on the EMBL database. A total of 578 (48.3%) upregulated and 205 (24.5%) downregulated JPS beetle contigs were annotated with GO functions. As a result, a total 246 and 239 genes were annotated with the GO term of biological process, while 145 and 90 genes were annotated with the GO term cellular component function and 327 and 363 genes with the GO term molecular function, respectively (Figure 3). There were more upregulated unigenes than downregulated unigenes annotated as being involved in biological processes. In particular, many upregulated genes were annotated as being involved in cellular component organization or biogenesis ( degeneration-multiple diseases pathways (KEGG:05022, 12), amyotrophic lateral sclerosis (KEGG:05014, 11), Toll and Imd signaling pathways (KEGG:04624, 10), and Shigellosis (KEGG:05131, 10). The 156 downregulated unigenes were involved in 246 pathways including axon regeneration (KEGG:04361, 13), human papillomavirus infection (KEGG:05165, 13), pathways in cancer (KEGG:05200, 12), Wnt signaling pathway (KEGG:04310, 9), axon guidance (KEGG:04360, 9), and the MAPK signaling pathway (KEGG:04010, 8) (Table S4). The 652 and 1218 upregulated and downregulated unigenes were annotated using the genome of T. castaneum, respectively, and GO enrichment analysis showed that unigenes that were more than 2-fold upregulated were involved in 28 biological processes, 46 cellular components, and 8 molecular functions while unigenes downregulated less than 0.5-fold were involved in 43 biological processes and 10 molecular functions. In the GO enrichment analysis, both up-and downregulated unigenes were annotated with the GO terms of biological processes and molecular function, while only upregulated unigenes received the GO annotation of cellular components. As shown in the GO analysis results, JPS adults infected by JEF-197 showed a decrease in the expression of genes involved in signaling and transmembrane transport. Signaling is critical for information  (Table S4).
The 652 and 1218 upregulated and downregulated unigenes were annotated using the genome of T. castaneum, respectively, and GO enrichment analysis showed that unigenes that were more than 2-fold upregulated were involved in 28 biological processes, 46 cellular components, and 8 molecular functions while unigenes downregulated less than 0.5-fold were involved in 43 biological processes and 10 molecular functions. In the GO enrichment analysis, both up-and downregulated unigenes were annotated with the GO terms of biological processes and molecular function, while only upregulated unigenes received the GO annotation of cellular components. As shown in the GO analysis results, JPS adults infected by JEF-197 showed a decrease in the expression of genes involved in signaling and transmembrane transport. Signaling is critical for information transfer in vivo, while transmembrane transport is important for sequestration of pathogens and manipulation of the transmembrane transport mechanisms of the host [49][50][51]. We predict that signal transmission and transport in JPS beetles are controlled by the fungus upon infection. In particular, pathways affected by the expression of eight gene expression level groups were identified. Genes involved in RNA transport (KEGG:03013) were 2-to 4-fold upregulated. Genes involved in the Toll and Imd signaling pathways (KEGG:04624), ribosome biogenesis in eukaryotes (KEGG:03008), and the spliceosome (KEGG:03040) were 4-to 8-fold upregulated. Genes involved in metabolic pathways (KEGG:01100), cysteine and methionine metabolism (KEGG:00270), glycine, serine, and threonine metabolism (KEGG:00260), purine metabolism (KEGG:00230), one carbon pool by folate (KEGG:00670), and ubiquinone and other terpenoid-quinone biosynthesis (KEGG:00130) pathways were 8-fold upregulated. Genes involved in the hedgehog signaling pathway-fly (KEGG:00670), Wnt signaling pathway (KEGG:00670), and Notch signaling pathway (KEGG:00670) were 0.5-to 0.25-fold downregulated (Table S5).
As mentioned above, certain genes involved in Toll and Imd signaling pathways, which are the main immune signaling pathways in insects, were 4-to 8-fold upregulated. Activation of the Toll and Imd signaling pathways is a typical immune response in fungusinfected insects [52]. Expression levels of the Gram-negative bacteria binding-protein 3 (GNBP3), modular serine protease (ModSP), and Spätzle (Spz) that recognizes β-1,3glucans of fungi were increased, and the expression levels of Tube, Pelle, and Cactus genes, which encode cytosolic components, were also increased. Expression of IMD pathway gene in JPS adults was increased by JEF-197 treatment. In particular, levels of peptidoglycan recognition protein LC (PGRP-LC), immune deficiency (IMD) protein, and a caspase (CASP8), which are involved in fungal recognition, were upregulated. These two pathways were activated in defensive responses, but genes related to antimicrobial peptides (AMPs) were not annotated ( Figure 4a, Table S6). Among DEGs, no genes related to the production of AMPs were annotated. Cecropin and defensin genes, which are AMPs, were identified, but were less than 2-fold upregulated (Table S7). The reason why AMPs were not identified may be because assembly was performed with default settings with the minimum contig length set to 200 bp, and several types of AMPs are short-length sequences. However, expression levels of genes in the Toll and IMD signaling pathway should increase prior to AMP production. Our results indicate that the immune response of JPS beetles was activated in response to fungal infection, and AMP production was likely induced. In previous studies, AMPs such as attacins, cecropin, coleoptericin, and lysozymes were identified in JPS beetles [36]; further studies are needed to determine the repertoire of AMPs produced by JPS beetles in response to fungal infection.
Upregulated genes involved in JPS beetle defense against fungi were assessed by KEGG analysis, and Toll and IMD signaling pathways were identified as playing major roles (Figure 4b, Table S6). The Toll signaling pathway is recognized and initiated by GNBP3, which binds to β-1,3-glucan, a major component of the fungal cell wall [53][54][55]. ModSP activation induces sequential activation of Clip-sectine proteases and Spätzle-processing enzymes [56,57]. Toll is activated by binding of Spätzle cleaved by proteolysis [58]. When Toll is activated, the Toll-interleukin1-resistance (TIR) domain of the Toll receptor recruits a cascade of signaling adapters and kinases to induce phosphorylation of cactus [59]. Toll signaling ultimately induces activation of AMP genes and dissociation of NF-κB protein from cactus [60]. The IMD pathway can be activated as part of the antifungal response, and the Rel/NF-κB transcription factor can be induced in later stages of fungal infection [61,62]. The IMD pathway is activated by recognition of microbial-derived molecules by the peptidoglycan recognition protein LC (PGRP-LC) [63]. Components of the IMD pathway such as Imd, dTAK1, Ird5, Kenny, CASP8, and Relish, which were genes identified as upregulated in our study, regulate the expression of antimicrobial peptide genes by activating Rel/NF-κB-like transcription factors [64]. In addition, the expression of genes involved in the DUOX pathway was increased, suggesting that the bactericidal ROS production pathway was activated in addition to the IMD pathway. Microorganisms capable of destroying DUOX-dependent ROS are regulated by IMD-dependent AMP and may play a complementary role in the ROS system, and ROS plays an important role in the control of intestinal bacteria [65,66]. capable of destroying DUOX-dependent ROS are regulated by IMD-dependent AMP and may play a complementary role in the ROS system, and ROS plays an important role in the control of intestinal bacteria [65,66]. The Wnt signaling pathway controls the balance between differentiation and proliferation and is involved in the cell cycle and gene transcription [67]. In KEGG pathway analysis, o-palmitoleoyl transferase, which affects both Wnt and Wnt by binding to the frizzled protein to transmit a signal to transcription factor 7 in the cells, was upregulated. In addition, frizzled expression results in gene transcription via MAPK signaling pathway activation. However, genes in this pathway were downregulated, suggesting that Wnt was not recognized by the frizzled receptor. This indicates that the Wnt signaling pathway was not fully functional and that cell cycle and gene transcription activity in JPS beetles was reduced as a result of fungal infection. This is consistent with the decrease in expression of genes with the GO annotation of 'Signaling'. Downregulated contigs were involved in the Notch signaling pathway, and the expression level of notch was also decreased, suggested decreased activation of the MAPK signaling pathway. Certain genes involved in the MAPK signaling pathway were found to be downregulated, suggesting reduced activation of Wnt and Notch signaling pathways and a decrease in cell cycling and gene transcription [68,69]. Fungal infection appears to reduce the expression of genes in the MAPK signaling pathway, which controls the cell cycle and gene transcription via Wnt and Notch signaling pathways. In our study, active immunity of JPS adults to fungal infection decreased gradually after the eighth day of infection, and infected insects ultimately died. The Wnt signaling pathway controls the balance between differentiation and proliferation and is involved in the cell cycle and gene transcription [67]. In KEGG pathway analysis, o-palmitoleoyl transferase, which affects both Wnt and Wnt by binding to the frizzled protein to transmit a signal to transcription factor 7 in the cells, was upregulated. In addition, frizzled expression results in gene transcription via MAPK signaling pathway activation. However, genes in this pathway were downregulated, suggesting that Wnt was not recognized by the frizzled receptor. This indicates that the Wnt signaling pathway was not fully functional and that cell cycle and gene transcription activity in JPS beetles was reduced as a result of fungal infection. This is consistent with the decrease in expression of genes with the GO annotation of 'Signaling'. Downregulated contigs were involved in the Notch signaling pathway, and the expression level of notch was also decreased, suggested decreased activation of the MAPK signaling pathway. Certain genes involved in the MAPK signaling pathway were found to be downregulated, suggesting reduced activation of Wnt and Notch signaling pathways and a decrease in cell cycling and gene transcription [68,69]. Fungal infection appears to reduce the expression of genes in the MAPK signaling pathway, which controls the cell cycle and gene transcription via Wnt and Notch signaling pathways. In our study, active immunity of JPS adults to fungal infection decreased gradually after the eighth day of infection, and infected insects ultimately died.

Changes in the Expression Level of Japanese Pine Sawyer Beetle Immune-Related Genes
Among the differentially expressed JPS beetle unigenes, 172 unigenes were annotated as immune-related genes base on immunoDB (E-value < 1 × 10 −100 ) ( Table S7). The immune-related genes belonged to 27 immune groups and included genes involved in autophagy, beta glucan binding, as well as genes encoding caspases, clip-domain serine proteases, c-type lectins, IMD pathway members, inhibitors of apoptosis, JAK/STAT pathway members, scavenger receptors class B proteins, serine protease inhibitors, thioredoxin peroxidases, and Toll pathway members ( Figure 5). The group of gene expressions in all immune classes, except catalases, MD2-like receptors, and superoxide dissidents, showed significant differences.

Changes in the Expression Level of Japanese Pine Sawyer Beetle Immune-Related Genes
Among the differentially expressed JPS beetle unigenes, 172 unigenes were annotated as immune-related genes base on immunoDB (E-value < 1 × 10 −100 ) ( Table S7). The immune-related genes belonged to 27 immune groups and included genes involved in autophagy, beta glucan binding, as well as genes encoding caspases, clip-domain serine proteases, c-type lectins, IMD pathway members, inhibitors of apoptosis, JAK/STAT pathway members, scavenger receptors class B proteins, serine protease inhibitors, thioredoxin peroxidases, and Toll pathway members ( Figure 5). The group of gene expressions in all immune classes, except catalases, MD2-like receptors, and superoxide dissidents, showed significant differences. Figure 5. Expression of immune-related genes in JPS beetles in response to M. anisopliae JEF-197 fungal treatment. Immunerelated genes in JPS beetles whose expression was altered by JEF-197 fungal treatment were determined by Blast2GO based on an arthropod immune-related gene database in immunoDB (E-value < 1.0 × 10 −100 ). Genes that showed over 2fold upregulation or 0.5-fold downregulation were considered upregulated and downregulated, respectively. Differences between groups of gene expression levels in adults with JPS within the immune class were analyzed using a one-way ANOVA in SPSS using a significance level of 0.05 (α). * represents a significant difference in the level of gene expression in the two groups, and ** represents a significant difference in the level of gene expression in the three groups.
Gene expression levels of JPS beetle immune-related unigenes were measured at 6, 24, 48, 72, and 96 h after treatment with JEF-197. Among IMD pathway members, inhibitor of nuclear factor kappa-B kinase gene expression was not significantly altered compared with baseline 6 and 24 h after fungal treatment; however, over 4-fold of upregulation was observed by qRT-PCR at 48 h after treatment, followed by attenuation of the expression level of this gene. Mitogen-activated protein kinase gene expression decreased after fungal infection ( Figure S3a,b). Among JAK/STAT pathway members, gene expression of the tyrosine-protein kinase hopscotch decreased slightly up to 24 h after fungal treatment; however, over 1000-fold of upregulation was observed by pRT-PCR at 48 h after treatment, followed by attenuation over time. Cytokine receptor expression continued to decrease after infection, while expression of the signal transducer and activator of transcription (STAT) increased ( Figure S3c-e). The JAK/STAT pathway is known to be activated in Figure 5. Expression of immune-related genes in JPS beetles in response to M. anisopliae JEF-197 fungal treatment. Immunerelated genes in JPS beetles whose expression was altered by JEF-197 fungal treatment were determined by Blast2GO based on an arthropod immune-related gene database in immunoDB (E-value < 1.0 × 10 −100 ). Genes that showed over 2-fold upregulation or 0.5-fold downregulation were considered upregulated and downregulated, respectively. Differences between groups of gene expression levels in adults with JPS within the immune class were analyzed using a one-way ANOVA in SPSS using a significance level of 0.05 (α). * represents a significant difference in the level of gene expression in the two groups, and ** represents a significant difference in the level of gene expression in the three groups.
Gene expression levels of JPS beetle immune-related unigenes were measured at 6, 24, 48, 72, and 96 h after treatment with JEF-197. Among IMD pathway members, inhibitor of nuclear factor kappa-B kinase gene expression was not significantly altered compared with baseline 6 and 24 h after fungal treatment; however, over 4-fold of upregulation was observed by qRT-PCR at 48 h after treatment, followed by attenuation of the expression level of this gene. Mitogen-activated protein kinase gene expression decreased after fungal infection ( Figure S3a,b). Among JAK/STAT pathway members, gene expression of the tyrosine-protein kinase hopscotch decreased slightly up to 24 h after fungal treatment; however, over 1000-fold of upregulation was observed by pRT-PCR at 48 h after treatment, followed by attenuation over time. Cytokine receptor expression continued to decrease after infection, while expression of the signal transducer and activator of transcription (STAT) increased ( Figure S3c-e). The JAK/STAT pathway is known to be activated in response to viruses; however, the expression level of STAT increased steadily after treatment. Although this gene is only one of the JAK/STAT genes, a previous study suggested that fungi cannot activate the JAK/STAT pathway. Our result as one of the highlights suggests that fungal infection may induce the expression of JAK/STAT pathway-related genes. Gene expression of Toll-like receptor 7 increased slightly 48 h after fungal treatment and decreased thereafter, while that of the Toll-like receptor Tollo decreased continuously ( Figure S3f,g). Among the beta-glucan binding proteins, the expression of beta-1,3-glucan-binding protein (GNBP1) increased after 72 h after treatment as did that of GNBP3 after 24 h after treatment, while the expression of beta-1,3-glucan-binding protein 1 (GNBP2) did not change significantly after 72 h after treatment ( Figure S3h,j). Although the change in expression of the Tolllike receptor appears to be insignificant, KEGG enrichment analysis result showed that expression of genes involved in the Toll pathway was increased after fungal treatment; increased expression of GNBP3, which recognizes fungi, may be the trigger that induces activation of the Toll pathway.

Expression of M. anisopliae JEF-197 in Japanese Pine Sawyer Beetles
Three replicate RNA-seq short reads of JEF-197 and fungal-treated JPS beetles were mapped to the contigs of JEF-197. Contigs were filtered, and the filtered contigs were analyzed for DEGs as described above.
A total of 73 DEGs were identified in JEF-197-infected JPS adults, and 42 of these DEGs were upregulated (Table S8). Of the upregulated contigs, 17 contigs were annotated as functional genes, while 17 hypothetical proteins and 8 no-hit contigs were identified, respectively. Functional genes were involved in energy metabolism, gene transcription regulation, material transport and degradation, and the stress response.
Energy metabolism-related genes such as oxidoreductase FAD-binding domain protein, ubiquinol-cytochrome c reductase complex subunit, and ATP synthase beta chain precursor were upregulated. Fungi use ATP synthesis to gain energy for invasion of JPS beetles. DNA-directed RNA polymerase II, Spt20 family protein, and methyltransferase sirN-like protein, which are genes related to gene transcription regulation, were also upregulated. Recent studies have shown the modulation of pathogenesis and virulence by transcriptional factors of entomopathogenic fungi [70][71][72]. Regulation of activation of gene transcription can contribute to toxicity in the host [73,74], while the recQ family helicase is involved in DNA repair and replication [75]. The increase in expression level of these genes suggests that transcription-, repair-, and replication-related genes play important roles in the pathogenesis and proliferation of the fungus JEF-197 in JPS beetles. The general substrate transporter, which transports small solutes in response to a chemical osmotic ion gradient [76], and the protein BFR2, which regulates protein transport or induces mass cell proliferation [77], were also upregulated in JEF-197.
The response of fungi to stress is essential for fungal growth in an insect host. Superoxide dismutase (SOD) and ubiquitin were identified as upregulated genes. SOD, an enzyme found in all living cells, breaks down superoxide and prevents tissue damage. It also acts as an enzyme that protects against reactive oxygen species (ROS) formed by solar ultraviolet radiation (UA-A and UA-B) and protects conidia against adverse environmental conditions [30,78]. In addition, catalase activity protects against host-derived H 2 O 2 [79]. In this study, we predicted that the expression of SOD would increase to protect the fungi from H 2 O 2 produced by JPS beetles. Ubiquitin binds to other proteins and promotes the breakdown of proteins, and its expression is increased in response to stress. In particular, damaged proteins that cannot be recovered due to heat shock are degraded through the ubiquitin-proteasome pathway [80]. We assumed that JEF-197 would have evolved strategies to overcome host-cell-induced stress [81]. In addition, we predicted that this fungus would overcome host-induced stress by increasing energy and protein metabolism and by proliferating in the insect host.
The 42 upregulated DEGs of JEF-197 are candidate genes that likely play crucial roles in invasion of the insect host. Upregulation of genes involved in energy and protein metabolism as well as responses to stressors suggests that these genes are involved in invasion. Defense and proliferation likely occur simultaneously to facilitate the growth of the fungus, and changes in expression patterns of certain fungal genes brought about by methyltransferases likely play a role in toxogenesis in insects. Proteins with unknown functions, such as domain of unknown function (DUF) and hypothetical proteins, showed elevated expression in JPS beetles. The functions of these genes should be clarified in future studies.
In the future, there is a need for practical verification through proteomics and antibody analysis of similar species based on the result of transcriptome analysis. Protein profiles have been proposed as more reliable predictors than transcription profiles [82]. Recently, gene and protein expression studies had mutual connectivity, and transcriptomics data and proteomics data could be complemented by providing protein database and protein level verification, respectively [83,84]. The results of present study, transcriptomic analysis of JPS and M. anisopliae that interacted after fungal treatment, will be used as a basis for further study of genomics-inspired proteomics. In the future, we would like to study the correlation between the expression level of the major genes identified in this study and the actual protein production.

Conclusions
To investigate the interaction between M. anisopliae JEF-197 and JPS beetles, we performed RNA-seq analysis of RNA samples obtained from JPS beetles on day 8 after fungal treatment. JPS genes related to RNA transport, ribosome biogenesis, and pathways important for protein production of spliceosomes were upregulated in the JEF-197-infected JPS adults, while JPS genes involved in pathways related to the cell cycle and gene transduction were downregulated. Immune-related genes in the Toll and Imd pathway were upregulated in response to fungal infection, consistent with previous findings. In JEF-197, genes involved in energy, protein metabolism, and stress were upregulated. Although the 42 genes of M. anisopliae identified in this study cannot fully explain the pathogenetic mechanisms of this fungus, these 42 genes are candidate genes that play an important role in the fungal invasion to JPS beetles and can provide insight into fungal pathogenesis in future studies.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/jof7050373/s1, Figure S1: Distribution of the assembled JPS transcripts according to length of contig, Figure S2: Validation using qRT-PCR for the expression level of the selected seven genes by RNA-seq analysis, Figure S3: The expression levels of immune-related genes in Metarhizium anisopliae JEF-197-treated JPS, Table S1: A list of primers used in qRT-PCR to validate the data of RNA-seq analysis and immune-related gene expression in Metarhizium anisopliae JEF-197-treated Japanese pine sawyer adult, Table S2: Summary of obtained read of non-treated control and fungus-treated Japanese pine sawyer adults by Illumina sequencing, Table S3: In silico cDNA library of assembled M. anisopliae JEF-197 and Japanese pine sawyer, Table S4: KEGG pathway annotation by up-and downregulated unigenes of Metarhizium anisopliae JEF-197-treated Japanese pine sawyer, Table S5: Enriched GO term and KEGG of Metarhizium anisopliae JEF-197-treated Japanese pine sawyer using a g:profiler, Table S6: DEGs of Japanese pine sawyer treated with Metarhizium anisopliae JEF-197 annotated in the IMD and Toll pathway, Table S7: List of immune-related genes in Japanese pine sawyer, Table S8: The DEGs of Metarhizium anisopliae JEF-197 in Japanese pine sawyer.