Mitogenome Analysis of Four Lamiinae Species (Coleoptera: Cerambycidae) and Gene Expression Responses by Monochamus alternatus When Infected with the Parasitic Nematode, Bursaphelenchus mucronatus

Simple Summary The longicorn beetle, Monochamus alternatus, is a major vector for the transmission of pine wilt disease, which is caused by a nematode pathogen, Bursaphelenchus xylophilus (or also possibly by B. mucronatus) that is spread by the beetle as it feeds on pine trees. In this study, the mitochondrial genome sequences of four longicorn species (Coleoptera: Cerambycidae: Lamiinae) were determined to further elaborate the phylogenetic relationships of Lamiinae. RT-qPCR was also used to assess the expression of eight mitochondrial protein-coding genes in M. alternatus when carrying B. mucronatus or not, so as to explore the relationship between these two species. The results showed that expression of mitochondria-encoded genes was elevated in M. alternatus beetles that were infected with B. mucronatus, suggesting that B. mucronatus putatively activates an immune response, which significantly affects the metabolic processes of M. alternatus. These results are of significance for further understanding the phylogenetic relationships of longicorn beetles and controlling the spread of pine wilt disease. Abstract We determined the mitochondrial gene sequence of Monochamus alternatus and three other mitogenomes of Lamiinae (Insect: Coleoptera: Cerambycidae) belonging to three genera (Aulaconotus, Apriona and Paraglenea) to enrich the mitochondrial genome database of Lamiinae and further explore the phylogenetic relationships within the subfamily. Phylogenetic trees of the Lamiinae were built using the Bayesian inference (BI) and maximum likelihood (ML) methods and the monophyly of Monochamus, Anoplophora, and Batocera genera was supported. Anoplophora chinensis, An. glabripennis and Aristobia reticulator were closely related, suggesting they may also be potential vectors for the transmission of the pine wood pathogenic nematode (Bursaphelenchus xylophilus) in addition to M. alternatus, a well-known vector of pine wilt disease. There is a special symbiotic relationship between M. alternatus and Bursaphelenchus xylophilus. As the native sympatric sibling species of B. xylophilus, B. mucronatus also has a specific relationship that is often overlooked. The analysis of mitochondrial gene expression aimed to explore the effect of B. mucronatus on the energy metabolism of the respiratory chain of M. alternatus adults. Using RT-qPCR, we determined and analyzed the expression of eight mitochondrial protein-coding genes (COI, COII, COIII, ND1, ND4, ND5, ATP6, and Cty b) between M. alternatus infected by B. mucronatus and M. alternatus without the nematode. Expression of all the eight mitochondrial genes were up-regulated, particularly the ND4 and ND5 gene, which were up-regulated by 4–5-fold (p < 0.01). Since longicorn beetles have immune responses to nematodes, we believe that their relationship should not be viewed as symbiotic, but classed as parasitic.


Introduction
Cerambycidae, also known as longicorn beetles, are a very large group within the Coleoptera with over 4000 genera and 35,000 species world-wide and are important components of forest ecosystems [1,2]. Lamiinae is considered to be the most diverse subfamily in Cerambycidae, with about 20,000 described species classified into 80 tribes [3]. The longicorn beetles of Lamiinae are widely distributed and most of the species are major pests in agriculture and forestry. Scholars have formerly relied heavily on morphology [4,5], anatomy and physiology [6], and other methods to establish relationships among species and create classification systems [7]. However, the evolutionary relationships among Lamiinae are hard to describe based on traditional morphological data, but with the development of molecular markers, these relationships can now be assessed with much greater precision [8]. Previous studies have used phylogenetic methods to assess the relationships within the subfamilies of Cerambycidae, and the monophyly of Lamiinae has been corroborated [3,[9][10][11]. However, very little other information on mitochondrial genomes has been reported for this subfamily, thus, we aim to enrich the database of Lamiinae and provide useful clues for the internal phylogenetic relationship of Lamiinae.
Monochamus alternatus, belongs to the subfamily Lamiinae. In East Asia, M. alternatus is believed to be the main vector for the pathogenic nematode Bursaphelenchus xylophilus [12] that infects pine trees and causes pine wilt disease [13]. B. xylophilus is transported to host pine trees by Monochamus beetles and after nematodes develop rapidly and form dispersal juveniles, these juveniles enter the tracheal system of the beetle for transport to new pine trees hosts, thus spreading the disease [14]. Bursaphelenchus mucronatus, a sister species to B. xylophilus, is a migratory endoparasitic nematode [14][15][16]. However, as a native sympatric sibling nematode to B. xylophilus, in both its morphological and biological characteristics [17], B. mucronatus shows weaker pathogenicity to forests [18]. However, various scholars believe that B. mucronatus can also cause pine wilt disease, and that this species jointly forms a complex with B. xylophilus [19,20]. Previous studies have revealed the physiological changes and symptom development of pine wilt disease, such as decreased photosynthesis, cambium destruction, traumatic resin canal formation, etc. [21]. Much attention has been given to the effects on the health of pine trees, but few studies have focused on the relationships between longicorn beetles and nematodes.
Mitochondria are mainly known for their crucial roles in aerobic ATP production involving the citric acid cycle and oxidative phosphorylation [22]. However, these organelles are involved in a variety of other processes including apoptosis and immunity [23]. In addition, the mitochondrial genome is widely considered to be a highly informative molecular marker for comparative genomic research [24,25], and is widely used to study phylogenetic relationships due to its maternal inheritance and high evolutionary rate [26]. There have already been several studies related to mitochondrial gene expression [27][28][29][30]. Previous studies have already demonstrated that gene expression in M. alternatus was different when infected with B. xylophilus [31,32]. Zhou et al. [31] proved that M. alternatus infected with B. xylophilus showed an increase in the expression of antioxidant genes in order to obtain immune tolerance. In addition, Ning et al. [32] found that there were many differentially expressed microRNAs (miRNAs) in M. alternatus with nematode infections and the expression of neuron-related miRNAs increased significantly. However, current research on mitochondrial genes related to M. alternatus and pine wood nematodes has focused only on genome sequencing and the regulation of nuclear genes, and there has not been much work on mitochondrial gene expression. Li et al. [33] indicated that there was no difference in the transcript level of M. alternatus infected with B. xylophilus or not. Considering that mitochondria are essential for the production of adenosine triphosphate (ATP) through oxidative phosphorylation [34], the energy metabolism of the longicorn beetles may change and manifest as a change in the transcript level of mitochondrial protein coding genes when carrying nematode. In the present study, the mitochondrial genome was not only used to study the phylogenetic relationships of Lamiinae, but also to figure out the relative expression of mitochondrial protein coding genes in response to B. mucronatus.
The present study aims to determine if there is a difference in mitochondrial gene expression in M. alternatus when infected with B. mucronatus. We examine how B. mucronatus changes the energy metabolism of M. alternatus, and further explore the parasitic relationship. We also provide new ideas about pine wilt disease by studying mitochondrial gene expression responses when M. alternatus beetles are infected by pine wood nematodes. Because molecular data on the species within the subfamily Lamiinae is still scarce, we also sequenced four mitochondrial genomes from Lamiinae species in order to enrich the mitochondrial genome sequence database and we built phylogenetic trees to explore the phylogenetic relationships within the Lamiinae subfamily.

Sampling Collection and Microscopic Examination
Samples of live M. alternatus were captured by traps set in the Taihuyuan pine wood nematode epidemic area of Lin'an (30 •  After morphological identification, the specimens of M. alternatus were then cultivated in the laboratory under artificial feeding conditions. Cultivation was in a plastic incubator that was placed in an artificial climate box; fresh Masson pine was provided every 3 days, the temperature in the box was maintained at 20-25 • C and relative humidity was 60-70%. Specimens of the other three longicorn beetles were deposited at −40 • C in the Animal Specimen Museum, College of Life Sciences and Chemistry, Zhejiang Normal University, China. After experimentation, all the pine wood was incinerated.
All M. alternatus were quickly dissected after a week and small pieces of abdominal tissue, thorax flight muscle, and intestinal tract were taken for microscopical examination. We used an OLYMPUS-SZX16 stereo microscope to observe abdominal tissue samples to directly determine the presence/absence of nematodes. If nematodes were found in M. alternatus tissue, then tissues from that individual were assigned to the experimental group whereas beetles without nematodes were placed into the control group. The excised digestive tract and flight muscle tissue samples were quickly frozen in liquid nitrogen and stored in an ultra-low temperature refrigerator at −80 • C until use.

DNA Extraction, Sequencing, and Molecular Identification
Total genomic DNA was extracted from thorax flight muscle tissues of each specimen by an Ezup Column Animal Genomic DNA purification Kit (Sangon Biotech Company, Shanghai, China). After bright and clear bands were detected with 1% agarose gel electrophoresis, the DNA was then stored at −20 • C. Cycling conditions and reaction volume of PCR amplifications were as in Cheng et al. and Gao et al. [35,36]. All PCR products were obtained using an automated DNA sequencer (ABI 3730) from Sangon Biotech Company (Shanghai, China) for both mitochondrial strands.
Referring to the method of amplifying insect genomes by Simon et al. [37] and combined with the M. alternatus mitochondrial genome published in the GenBank database, Clustal W software in Mega 7.0 [38] was used to align the primer sequences. Eleven pairs of universal primers were modified and 11 adjacent or overlapping fragments were amplified [38][39][40][41]. Using qualified DNA samples as templates, the mitochondrial genes of M. alternatus obtained from the Lin'an District of Hangzhou were sequenced via PCR. Primers were also used to amplify the mitochondrial genome fragments of the nematodes collected and were then further modified to distinguish between B. xylophilus and B. mucronatus [42]. After using these detection primers in PCR experiments, the appearance of a bright and clear band after electrophoresis meant that sample contained the corresponding nematode species. This step also eliminated the possibility that the control samples were carrying nematodes. Through this step, we ensured that the longicorn beetles in the control group do not contain nematodes, and that the beetles in the experimental group can be detected as B. mucronatus by both microscope examination and PCR.

Sequence Annotation and Analyses
All mitochondrial genome sequences were manually proofread and spliced using DNASTAR Package v.6.0 software [43] and tRNA gene annotation was completed through the online website MITOS (http://mitos.bioinf.uni-leipzig.de/index.py, accessed on 15 January 2021) [44]. The complete sequences of some Cerambycidae species were obtained from the GenBank database [1,2,11,[45][46][47][48][49][50][51][52][53], and homologous sequence alignments were performed by the Clustal W program of Mega 7.0 [38] to determine the positions of the two rRNA genes and the AT-rich region. The four mitogenomes were deposited in GenBank with the accession numbers MW858148, MW858150-MW858152. Also, the nucleic acid sequences of the 13 protein-coding genes (PCGs) were translated into amino acid sequences on the basis of the invertebrate genetic codes [54]. Mega 7.0 [38] software was used to find the open reading frame of each of the 13 PCGs, and to calculate the AT content and codon usage rate of the newly determined mitochondrial sequences. The formula used for skew analysis of gene composition was as follow: AT-skew = (A − T)/(A + T) and GC-skew = (G − C)/(G + C) [55]. The relative synonymous codon usage (RSCU) of four longicorn beetle mitochondrial genomes was calculated using PhyloSuite [56].

Comparison of Protein Coding Genes and Phylogenetic Analyses
To study the phylogenetic relationships of Cerambycidae, a nucleotide dataset (13P33) of the 13 protein-coding genes of 33 mitogenomes was used according to Zhang et al. [57] ( Table 1). Newly sequenced mitogenomes of Au. atronotatus, Ap. germari, P. fortunei and M. alternatus were included, with Sympiezomias velatus and Naupactus xanthographus selected as outgroups. The Clustal W program of Mega 7.0 software was used to align the amino acid sequences of the protein-coding genes [38], the conserved regions were selected by Gblock 0.91b [58], and PartionFinder 1.1.1 [59] was used to select partitions based on Bayesian information criterion (BIC). The phylogenetic relationships were analyzed by the Bayesian inference (BI) and maximum likelihood (ML) methods. The partition schemes and best-fit models selected for each data set are provided in Table 2. We used a GTR + I + G model to do ML and BI analyses. ML analysis was performed by RAxML 8.2.0 software based on the GTRGAMMAI best model, and branch confidence was obtained by rapid inference using 1000 repetitions [60]. The BI analysis was run by MrBayes 3.2 [61] based on the GTR + I + G model through 10 million generations, with the first 25% of generations removed as burn in, and the mean standard deviation of Bayesian split frequency was less than 0.01.  ND4L_codon1, ND1_codon1, ND5_codon1, ND4_codon1 TVM + I + G Partition 7 ND4L_codon2, ND1_codon2, ND5_codon2, ND4_codon2 GTR + I + G Partition 8 ND1_codon3 HKY + G Partition 9 ND6_codon2, ND2_codon2, ND3_codon2 TVM + I + G Partition 10 ND2_codon3 TRN + G Partition 11 ND5_codon3, ND4_codon3, ND4L_codon3 HKY + G

RNA Isolation and cDNA Synthesis
Samples (50-100 mg) of M. alternatus intestines frozen at ultra-low temperature were used for RNA extraction according to the instructions of the TaKaRa MiniBEST Universal RNA Extraction Kit (Takara, Japan). The integrity of RNA was verified by the presence of sharp bands for 28S and 18S rRNA using a 1% agarose gel stained with Goldview and run at 130 V, 120 mA for 25 min [62]. Then, the samples were frozen at −80 • C until use. We performed reverse transcription according to the PrimeScript™ RT Master Mix kit instructions (Takara, Japan). Reverse transcription reaction conditions were: 37 • C, reverse transcription for 15 min; 85 • C, reverse transcriptase inactivation for 5 s.

RT-qPCR
We used DNASTAR [43] and Primer Premier 5 software [63] to design primers according to the M. alternatus sequence downloaded from the NCBI website and the newly obtained sequence of M. alternatus in Lin'an, Hangzhou. The β-actin gene was selected as the reference gene [62,64] to design the internal reference primer of M. alternatus. The amplicon was synthesized by Sangon Biotech Company (Shanghai, China), followed by screening of the quantitative primers for 13 protein-coding genes. Due to the strong AT skew and the length of gene fragments in 5 protein-coding genes, their primers failed to amplify. However, the remaining eight protein-coding genes could be used for RT-qPCR experiments ( Table 3). The primers shown in Table 3 were screened by RT-qPCR reactions, and then gene expression in the intestine of M. alternatus was measured in both experimental and control groups. The genes corresponding to each primer pair were assessed in three technical replicates. Reactions were carried out in Primer Script TM RT Master Mix (TaKaRa, Tokyo, Japan), with conditions as follows: initial denaturation at 95 • C for 30 s, and 40 cycles (95 • C for 5 s, 60 • C for 30 s) to generate a melt curve. Table 3. Primer sequences of eight mitochondrial protein-coding genes and β-actin gene designed for RT-qPCR experiment in this study.

Gene
Forward Primer Forward Primer The relative expression of mRNA in the intestinal tract of M. alternatus at a specific stage was analyzed by the 2 −∆∆Ct method and normalized to the β-actin gene. We used SPSS 21.0 software to analyze all data, and data was expressed as mean ± standard deviation (mean ± SE). During the data analysis process, the average value of the control group was set to 1. Gene expression of the experimental and control groups was pre-compared using the t-TEST function, with p < 0.05 considered statistically significant. Finally, Origin 8.0 software was used to plot results.

Nematode Identification
After the microscopic examination, we found that some abdominal tissues contained nematodes, so we proceeded to molecular identification. We used specific primers for detection of B. xylophilus and B. mucronatus genes and performed PCR analysis on the DNA template of M. alternatus. The results showed that only the B. mucronatus primers produced bright and clear bands after gel electrophoresis. Therefore, it was determined that the M. alternatus specimens that were collected contained only B. mucronatus nematodes. M. alternatus tissue samples from individuals without nematodes were separated into the control group, and those containing B. mucronatus were placed in the experimental group.

Components of the Mitochondrial Genomes
We obtained one complete and three almost complete mitogenomes of Lamiinae species representing four genera (Aulaconotus, Apriona, Paraglenea, and Monochamus) (Supplementary Materials, Tables S1-S4). The complete mitochondrial genome of P. fortunei was 15,496 bp in length, and the mitogenomes of Au. atronotatus, Ap. germari, and M. alternatus were 14,491 bp, 14,858 bp, and 14,189 bp in length, respectively.
The mitogenome of P. fortunei contained 13 protein-coding genes, 22 tRNAs, two rRNAs and an A + T rich region, the same as the ancestral mitochondrial genome (Figure 1). The A + T content, AT-skew and GC-skew of corresponding regions (mitogenome, PCGs and rRNAs) of each longicorn beetle species were calculated and are shown in Table 4. The nucleotide composition of the four longicorn beetles mitogenomes was strongly biased towards A and T, which made up 74.2% (Au. atronotatus) to 78.4% (M. alternatus) of the base pairs. All tRNAs of these longicorn beetles showed a typical cloverleaf secondary structure except for trnS1 (AGA) in Au. atronotatus and Ap. Germari, which formed a simple loop instead of the dihydrouridine (DHU) arm ( Figures S1-S4).

Components of the Mitochondrial Genomes
We obtained one complete and three almost complete mitogenomes of Lamiinae species representing four genera (Aulaconotus, Apriona, Paraglenea, and Monochamus) (Supplementary Materials, Tables S1-S4). The complete mitochondrial genome of P. fortunei was 15,496 bp in length, and the mitogenomes of Au. atronotatus, Ap. germari, and M. alternatus were 14,491 bp, 14,858 bp, and 14,189 bp in length, respectively.
The mitogenome of P. fortunei contained 13 protein-coding genes, 22 tRNAs, two rRNAs and an A + T rich region, the same as the ancestral mitochondrial genome ( Figure  1). The A + T content, AT-skew and GC-skew of corresponding regions (mitogenome, PCGs and rRNAs) of each longicorn beetle species were calculated and are shown in Table  4. The nucleotide composition of the four longicorn beetles mitogenomes was strongly biased towards A and T, which made up 74.2% (Au. atronotatus) to 78.4% (M. alternatus) of the base pairs. All tRNAs of these longicorn beetles showed a typical cloverleaf secondary structure except for trnS1 (AGA) in Au. Atronotatus and Ap. Germari, which formed a simple loop instead of the dihydrouridine (DHU) arm ( Figures S1-S4).  Among the four sequenced mitochondrial genomes, 11 protein-coding genes used the typical invertebrate start codon ATN, whereas TTG and GTA were used as the initiation codon in ND1 and ND3. Typical stop codons such as TAA and TAG were widely used. However, we also observed an incomplete stop codon T in five genes: COI (Au.   Among the four sequenced mitochondrial genomes, 11 protein-coding genes used the typical invertebrate start codon ATN, whereas TTG and GTA were used as the initiation codon in ND1 and ND3. Typical stop codons such as TAA and TAG were widely used. However, we also observed an incomplete stop codon T in five genes: COI (Au. atronotatus and M. alternatus), COII (Au. atronotatus, Ap. germari, P. fortunei and M. alternatus), COIII (Au. atronotatus and M. alternatus), ND4 (Ap. germari, P. fortunei and M. alternatus) and ND5 (Au. atronotatus, Ap. germari, P. fortunei and M. alternatus). The RSCU of each of the mitogenomes was calculated and the results indicated that the main codons used were highly similar in the four longicorn beetles (Figure 2). The most used codons (>284) in the PCGs of these mitogenomes were UUA (Leu), AUU (Ile) and UUU (Phe).

Phylogenetic Analyses
Based on the 13 mitochondrial protein-coding genes, the phylogenetic relationships within the subfamily Lamiinae were inferred from BI and ML analyses (Figure 3). The results showed that both ML and BI trees within Lamiinae split into two branches: (1) the clade consisting of Lamiinae sp. 1 ACP-2013, Agapanthia daurica and Au. atronotatus was the earliest diverging clade, and (2) the other 28 species gathered together in another branch. Lamiinae and the outgroup were well separated. In addition, the monophyly of Monochamus, Anoplophora, Batocera genera was supported. results showed that both ML and BI trees within Lamiinae split into two branches: (1) the clade consisting of Lamiinae sp. 1 ACP-2013, Agapanthia daurica and Au. atronotatus was the earliest diverging clade, and (2) the other 28 species gathered together in another branch. Lamiinae and the outgroup were well separated. In addition, the monophyly of Monochamus, Anoplophora, Batocera genera was supported.

Transcriptions of Mitochondrial Protein-Coding Genes
The relative transcription of eight mitochondrial protein-coding genes from the intestines of control and experimental groups of M. alternatus were quantitatively determined by RT-qPCR to assess the effects of the presence versus absence of B. mucronatus nematodes on the transcription of the beetle genes ( Figure 4). The results indicated that in the presence of B. mucronatus, transcript levels of COI, COII, COIII, ND1, ND4, ND5, ATP6 and Cyt b were up-regulated in M. alternatus to varying degrees. The transcription level of ND4 and ND5 genes increased by 4.25 ± 1.01, and 3.35 ± 0.78-fold, respectively (p < 0.01). Transcript levels of Cyt b, COI, and ATP6 genes in beetle intestine also increased significantly by 5.11 ± 2.12, 2.64 ± 0.51 and 2.64 ± 0.71-fold, as compared to the control group, respectively (p < 0.05). The remaining three genes (COII, COIII, ND1) showed no significant differences in their transcription levels between control and nematode-infected groups.
Insects 2021, 12, x FOR PEER REVIEW 10 of 16 determined by RT-qPCR to assess the effects of the presence versus absence of B. mucronatus nematodes on the transcription of the beetle genes ( Figure 4). The results indicated that in the presence of B. mucronatus, transcript levels of COI, COII, COIII, ND1, ND4, ND5, ATP6 and Cyt b were up-regulated in M. alternatus to varying degrees. The transcription level of ND4 and ND5 genes increased by 4.25 ± 1.01, and 3.35 ± 0.78-fold, respectively (p < 0.01). Transcript levels of Cyt b, COI, and ATP6 genes in beetle intestine also increased significantly by 5.11 ± 2.12, 2.64 ± 0.51 and 2.64 ± 0.71-fold, as compared to the control group, respectively (p < 0.05). The remaining three genes (COII, COIII, ND1) showed no significant differences in their transcription levels between control and nematode-infected groups.

Phylogenetic Analyses
The Lamiinae subfamily is widely distributed and contains the largest number of species among the Cerambycidae, but the phylogenetic relationships within Lamiinae are still controversial [2]. However, the subfamily is considered to be a monophyletic group, as has been shown in previous studies [2,[9][10][11]. In our study, two species within the subfamily Entiminae were selected as the outgroups, and the results showed that all the Lamiinae species were well separated with a high degree of confidence. Both ML and BI phylogenetic trees showed the formation of two large branches within the Lamiinae subfamily ( Figure 3). Batocera lineolata (JN986793 and MF521888) was closely related to Apriona swainsoni, which was consistent with the results of Wang et al. [2]. We also found that the clade of ((An. chinensis + An. glabripennis) + Aristobia reticulator) was a sister clade to the genus Monochamus with 100% confidence both in ML and BI analyses. These longicorn beetles showed a close phylogenetic relationship, suggesting that An. chinensis, An. glabripennis and Ar. reticulator may all be potential vectors for the transmission of pine wood nematodes. Our study increases the richness of the Cerambycidae genome

Phylogenetic Analyses
The Lamiinae subfamily is widely distributed and contains the largest number of species among the Cerambycidae, but the phylogenetic relationships within Lamiinae are still controversial [2]. However, the subfamily is considered to be a monophyletic group, as has been shown in previous studies [2,[9][10][11]. In our study, two species within the subfamily Entiminae were selected as the outgroups, and the results showed that all the Lamiinae species were well separated with a high degree of confidence. Both ML and BI phylogenetic trees showed the formation of two large branches within the Lamiinae subfamily ( Figure 3). Batocera lineolata (JN986793 and MF521888) was closely related to Apriona swainsoni, which was consistent with the results of Wang et al. [2]. We also found that the clade of ((An. chinensis + An. glabripennis) + Aristobia reticulator) was a sister clade to the genus Monochamus with 100% confidence both in ML and BI analyses. These longicorn beetles showed a close phylogenetic relationship, suggesting that An. chinensis, An. glabripennis and Ar. reticulator may all be potential vectors for the transmission of pine wood nematodes. Our study increases the richness of the Cerambycidae genome information and can assist in phylogenetic, molecular systematics and evolutionary studies of Cerambycidae. However, molecular data on species of Lamiinae are still scarce and more information is needed to fully explore the phylogeny within Lamiinae.

Mitochondrial Gene Expression
The relative transcription of eight mitochondrial protein-coding genes in the intestines of the control and experimental groups of M. alternatus was quantitatively determined by RT-qPCR to compare the effects of the presence versus absence of B. mucronatus nematodes on mitochondrial gene expression in M. alternatus. The results showed that transcription of all eight genes in the experimental group was up-regulated to varying degrees.
Mitochondria are essential to the process of generating adenosine triphosphate (ATP) through oxidative phosphorylation [34]. The mitochondrial respiratory chain responsible for generation of cellular ATP consists of five multi-subunit enzyme complexes (Complex I-V) located at the inner mitochondrial membrane. Four of these, Complexes I, III, IV and V, are bipartite and consist of subunits derived from both mitochondrial DNA (mtDNA) and nuclear DNA [65][66][67]. Thirteen protein-coding genes are encoded by the mitogenome and whether they are correctly expressed is very important for maintaining mitochondrial respiratory chain function [68]. Changes in the levels of mtDNA-encoded mRNA occur rapidly in response to changes in energy requirements of the cell [66]. When infected with B. mucronatus, the expression of COI, ND4, ND5, ATP6, and Cyt b genes in M. alternatus were all significantly up-regulated, and ND4 and ND5 genes increased by 4.25 ± 1.01 and 3.35 ± 0.78-fold compared to control groups (p < 0.01). These results support the viewpoint proposed by Li et al. [33] that a large number of stress-responsive genes were up-regulated in M. alternatus when infected with pine wood nematodes. These five genes are important subunits of Complex I, III, IV, and V and can significantly affect metabolic energy processes. ND4 and ND5 genes are believed to have originated from a common ancestor and have similar functions in the synthesis of Complex I [69]. Complex I is the main point for electrons to enter the electron transport chain, and is key for cells to produce large amounts of ATP through oxidative phosphorylation [34]. Complex I is also one of the main sources of reactive oxygen species (ROS) [67,70]. The Cyt b gene is responsible for electron transfer in Complex III. Cytochrome bc1 encoded by Cyt b is the only subunit of Complex III that is encoded by a mitochondrial gene, and has an important relationship with the composition of the functional Complex III [71]. The COI gene is one of the subunits of Complex IV [72]. Complex IV catalyzes the oxidation of cytochrome C and the reduction of oxygen to water, it is an enzyme at the end of the electron transport chain and it also functions as a proton pump [73]. The ATP6 gene is one of the subunits of Complex V. Under normal physiological conditions, this multi-subunit protein utilizes energy in the form of the proton gradient, generated by the electron transport chain, to phosphorylate ADP to form ATP [74].
B. mucronatus and B. xylophilus belong to the same genus, and in the past, B. mucronatus has attracted little attention as compared with B. xylophilus. Some previous studies have already shown that the expression of several genes in M. alternatus was altered when infected by B. xylophilus. For example, Zhou et al. [31] showed that in order to obtain immune tolerance, M. alternatus infected with B. xylophilus would induce an increase in the expression of antioxidant genes. Through the increased expression of these antioxidant genes, the level of ROS in the trachea of longicorn beetles increases and reaches a state of equilibrium [31]. Coincidentally, as one of the main sources of ROS, the genes ND4 and ND5, which code for two subunits of Complex I, were found to be significantly up-regulated in this study. Besides, mitochondria can simultaneously produce ATP and macromolecule synthesis, which enables mitochondria to meet the metabolic needs of different immune cells [75]. Accumulating evidence has shown that cellular energy metabolism pathways are altered during the differentiation and activation of immune cells, revealing mechanisms of metabolic regulation and roles in immune functions [76]. In our study, analysis of the expression of mitochondrial-encoded genes of the electron transport chain complex showed up-regulation of transcript levels in M. alternatus infected with B. mucronatus, indicating increased synthesis of protein subunits in Complexes I, III, IV, and V. This may indicate that M. alternatus is gaining immune tolerance to the parasitic nematode. Besides, a hypothesis has been proposed that metabolism can be an input to the system and act as a signal to control the immune functions of a cell [77]. Loftus et al. also proposed that cellular metabolism has direct roles in regulating immune cell function. Ning et al. showed that M. alternatus with and without B. xylophilus nematodes had different metabolic characteristics [32]. There were many differentially expressed microRNAs in M. alternatus infected with B. xylophilus and their functions were mostly directed towards regulating genes associated with metabolism and immunity [32], which is also consistent with the significant increase in mitochondrial gene expression described in our research results. The expression of certain miRNAs related to apoptosis, material metabolism, and olfactory neurons in M. alternatus with B. xylophilus increased significantly [32]. Li et al. found that in B. xylophilus infected beetles the differentially expressed gene functions in adult M. alternatus tissues mainly point to pressure regulation in the trachea, hormone response, and tissue and DNA repair [33]. By enhancing the expression of various genes, the tolerance of M. alternatus to B. xylophilus is improved, and a stable state of life is also maintained [22].
Combining these research results with ours, indicated that the up-regulated mitochondrial protein coding genes expressed in long beetles carrying B. mucronatus may play a key role in the process of carrying nematodes. The relationship formed between M. alternatus and B. mucronatus could be related to the expression of mitochondrial protein-coding genes. The results showed that when carrying nematodes, there was an enhancement in mitochondrial metabolism, and the up-regulated genes were all related to immune function, suggesting that the immune function of longicorn beetles is affected and they could have created immune measures against nematodes to protect themselves. Previous studies have shown that in the process of M. alternatus responding to pine wood nematodes, M. alternatus is in a state of suppressed immunity as a whole. Some immune genes in the epidermis and trachea are up-regulated, and the expression level of most immune genes remains unchanged. M. alternatus maintains immune homeostasis between the two through the regulation of ROS [31,78]. Some researchers believe that longicorn beetles and nematodes have a symbiotic relationship [33,79]. However, as the existence of B. mucronatus increased the transcription level of immune-related genes in M. alternatus, which manifested as those M. alternatus maintaining a high metabolic level in the intestine, we prefer to suggest the relationship between them is more like a parasite. We speculate that more efforts were made by M. alternatus to offset the negative effects of nematodes. Werner also suggested through molecular phylogeny that the relationship between nematodes and beetles should be host switching rather than nematode-beetle coevolution [80]. The other mechanisms need to be studied in the future.

Conclusions
The complete mitochondrial genome sequence of P. fortunei and three nearly complete mitogenomes of Au. atronotatus, M. alternatus, Ap. germari were determined. The phylogenetic tree supported the monophyly of Monochamus, Anoplophora, Batocera genera. Since the clade consisting of An. chinensis, An. glabripennis and Ar. reticulator was a sister clade to the one containing Monochamus with 100% confidence, it is suggested that all four of these species may be potential vectors for the transmission of pine wood nematodes. All of the eight mitochondrial genes that were assessed showed an up-regulation trend, with CO I, ATP6 and Cyt b being significantly up-regulated (p < 0.05) with transcript levels increased by 5.11 ± 2.12, 2.64 ± 0.51 and 2.64 ± 0.71-fold, respectively. ND4 and ND5 genes were also significantly up-regulated (p < 0.01) with 4.25 ± 1.01 and 3.35 ± 0.78fold increases, respectively. Since five genes corresponding to important subunits of Complex I, III, IV, and V were affected, it can be proposed that the increase in gene expression is a response to the changing demands of cellular energy requirements. To sum up, the M. alternatus that carries the B. mucronatus nematode still maintains a high metabolic level in the intestine, and this could be related to an autoimmune or symbiotic relationship of longicorn beetles against nematodes or it may indicate a parasitic relationship between longicorn beetles and nematodes.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/insects12050453/s1, Table S1. Location of features in the mitochondrial genome of Au. atronotatus. Table S2. Location of features in the mitochondrial genome of Ap. germari. Table S3. Location of features in the mitochondrial genome of P. fortune. Table S4. Location of features in the mitochondrial genome of M. alternatus. Figure S1. Inferred secondary structures of the tRNA genes in Au. atronotatus. Figure S2. Inferred secondary structures of the tRNA genes in Ap. germari. Figure  S3. Inferred secondary structures of the tRNA genes in P. fortune. Figure S4.  The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Institutional Review Board Statement: Not applicable.
Data Availability Statement: Not applicable.