Mitochondrial DNA analysis reveals spatial genetic structure and high genetic diversity of Massicus raddei (Blessig) (Coleoptera: Cerambycidae) in China

Abstract The oak longhorned beetle (OLB), Massicus raddei (Blessig, 1872) (Coleoptera: Cerambycidae), is widely distributed in Asia (China, the Korean Peninsula, Japan, Vietnam and the Russian Far‐East), but pest outbreaks have occurred only in Liaoning Province and Jilin Province of China. In order to explore possible mechanisms of local population outbreaks and characterize the genetic diversity and genetic structure of M. raddei across its range in China, three mitochondrial genes (COI, Cytb, and COII) were sequenced and analyzed for seven M. raddei populations collected from six provinces in China. From these different populations, we found a high haplotype and nucleotide diversity. Haplotype networks and phylogenetic analyses both demonstrate apparent genetic diversification between SC (southern China) and NC (northern China) population groups. A set of 21 pairwise comparisons for Fst (pairwise fixation indices) and Nm (genetic flow index) showed significant genetic differentiation and limited gene flow except for two pairs, Shandong (SD) and Liaoning (LN), and Anhui (AH) and Henan (HN). This pattern suggested that the periodic outbreak of the LN population could not be attributed to the absence of genetic flow with other spatial populations and that regional environmental factors might be responsible. AMOVA (Analysis of molecular variance) showed that the greater molecular genetic variation was among populations. Based on Tajima's D statistic, Fu's Fs, and the mismatch distribution test, we determined that the seven populations sampled were stable and had not experienced any recent population expansion. The fact that all the sampled populations showed only unique haplotypes and lacked shared or ancestral haplotypes, as well as the nonstar‐like distribution of haplotype network for concatenated genes, collectively provided powerful evidence of the stable and isolated nature of most populations. The high genetic differentiation and spatial genetic structuring among populations are both likely related to the beetle's moderate flight capacity, regional variation in host tree species and microclimate, as well as the geographic distance between sampling sites.

population outbreaks and characterize the genetic diversity and genetic structure of M. raddei across its range in China, three mitochondrial genes (COI, Cytb, and COII) were sequenced and analyzed for seven M. raddei populations collected from six provinces in China. From these different populations, we found a high haplotype and nucleotide diversity. Haplotype networks and phylogenetic analyses both demonstrate apparent genetic diversification between SC (southern China) and NC (northern China) population groups. A set of 21 pairwise comparisons for Fst (pairwise fixation indices) and Nm (genetic flow index) showed significant genetic differentiation and limited gene flow except for two pairs, Shandong (SD) and Liaoning (LN), and Anhui (AH) and Henan (HN). This pattern suggested that the periodic outbreak of the LN population could not be attributed to the absence of genetic flow with other spatial populations and that regional environmental factors might be responsible. AMOVA (Analysis of molecular variance) showed that the greater molecular genetic variation was among populations. Based on Tajima's D statistic, Fu's Fs, and the mismatch distribution test, we determined that the seven populations sampled were stable and had not experienced any recent population expansion. The fact that all the sampled populations showed only unique haplotypes and lacked shared or ancestral haplotypes, as well as the nonstar-like distribution of haplotype network for concatenated genes, collectively provided powerful evidence of the stable and isolated nature of most populations. The high genetic differentiation and spatial genetic structuring among populations are both likely related to the beetle's moderate flight capacity, regional variation in host tree species and microclimate, as well as the geographic distance between sampling sites.

| INTRODUC TI ON
Massicus raddei (Blessig) (Coleoptera: Cerambycidae) is known as a serious trunk borer that mainly infests species of oaks (Chen et al., 1959). The larvae of M. raddei bore into the xylem of host trees, creating large galleries that reduce the transportation of water and nutrition, causing crown dieback and ultimate death (Sun, 2001). M. raddei has a wide distribution in the oriental region, including most provinces in China (Sun, 2001), the Korean Peninsula (Cho et al., 2010;Kim et al., 2014;Lim et al., 2014), Japan (Kojima, 1931;Leksono et al., 2006), Vietnam (Nga et al., 2014), and the Russian Far-East (Anisimov & Bezborodov, 2017;Mokuroku Database, 2017). Although this borer species has not yet been reported from other continents, there still exists a risk of invasion vectored by transplanted host plants or timbers infested with its larvae or pupae, this species has recently been added to the "alert list" of the European and Mediterranean Plant Protection Organization and the United States Department of Agriculture (EPPO, 2015;USDA, 2015). This borer shows ecological differences in its life history, host trees, and timing of emergence among populations in various regions (EPPO, 2018). In heavily infested stands of Quercus mongolica Fisch.ex Leddb and Quercus liaotungensis Mary in northeastern China (Liaoning Province), one generation is completed in 3 years with six larval instars (Wang et al., 2012). M. raddei was first reported as damaging oak forests in Liaoning Province, with nearly perfectly synchronized mass emergences noted in every 3 years from 1993 to 2017. However, in interval 2 years, there were relatively scarce adults to emerge Wang & Shao, 2013;Wang et al., 2012;Yang et al., 2011). In Shandong Province, Anhui Province, and Henan Province, M. raddei attacks Quercus acutissima Carruth, while in more southern places such as Fujian and Yunnan provinces, the common host is Quercus glauca Thunb. In addition, fewer than 3 years might be needed in these sites (EPPO, 2018). In other distributional areas apart from Liaoning and Jilin provinces, M. raddei does not display periodic emergence, and adults emerge every year without outbreaks and do not pose great threat to the growth of oak forests (Luo et al., 2005;Wei et al., 2009;State Forestry Administration, 2016;EPPO, 2018). Therefore, there are reasons to suspect that different ecological patterns of phenology and density of M. raddei might be influenced by limited genetic flow among spatial populations, variation in host tree adaptation, and regional climatic conditions, as well as the regulatory effect of natural enemies (such as Dastarcus helophoroides Fairmaire or Cerchysiella raddeii Yang found in oak forests of Liaoning Province) (Gao et al., 2003;Yang et al., 2013).
The periodicity, involving synchronized adult emergences, fixed life-cycle length, and intervals between emergences with no or scarce adults present, is rare in insects (Heliövaara et al., 1994).
Of these, periodical cicadas Magicicada spp. and some periodical moths are typical study objects for deep comprehension of periodicity development (Ito et al., 2015;Martin & Simon, 1990;Várkonyi et al., 2002;Williams & Simon, 1995). Previously, Yoshimura (1997) proposed that limited gene flow restricted by geographical barriers and extreme low temperature events during the Pleistocene might be driving factors for the development of periodical cicadas from nonperiodical species Okanagana spp. In addition, the adults of the noctuid genus Xestia exhibit periodic dynamics in eastern Lapland but not in western Lapland, Kankare et al. (2002) discovered that substantial genetic differentiation and isolation did not exist between nonperiodical and periodic cohorts, which might not be a plausible explanation for its periodicity evolution. As for M. raddei, whether there exist in subtle genetic differentiation and limited gene flow between population of Liaoning Province and those of other distributional regions is still a question, which is speculated to reinforce its development of their particular periodic adult emergence pattern in Liaoning Province.
Recently, mitochondrial DNA markers COI (cytochrome oxidase subunit I), Cytb (cytochrome b), and COII (cytochrome oxidase subunit II) genes have been developed into the common molecular tool for many phylogenetic studies of beetles. Eight microsatellite loci and a partial mtDNA COI gene were used as molecular markers to elucidate patterns of genetic structure among 32 populations of the endangered beetle Rosalia alpina Linnaeus (Coleoptera: Cerambycidae) across central and south-east Europe, to clarify its biology and recent population expansion, and to devise more effective conservation measures (Drag et al., 2015). Tang et al. (2016) used the combination of three mitochondrial genes to study the genetic structure of Galerucella birmanica Jacoby (Coleoptera: Chrysomelidae) in its main distributional habitats. Moreover, Javal et al. (2017)

investigated invasion history of
Anoplophora glabripennis Motschulsky (Coleoptera: Cerambycidae) across Europe based on COI gene. Lee et al. (2020) presented the entire genetic structure among ALB (A. glabripennis) populations of South Korea and confirmed that ALB invasion has occurred even within the species' native ranges based on COI mitochondrial sequences. Studies on M. raddei have focused mostly on its biology, damage pattern, field surveys for natural enemies, and integrated control measures (Li et al., 2017;Tang, 2011;Wang et al., 2010;Yang et al., 2013Yang et al., , 2014. However, as yet there has been no further researches into genetic aspects of M. raddei populations in China. In this paper, we used three mitochondrial COI, Cytb, and

| Specimen collection
A total of 40 adults of M. raddei were collected by a specially designed black light  at night from early May to late August in 2019. The seven locations mostly cover the range of this pest which differs in host species in China ( Figure 1; Table 1). The observation of sampled adults was conducted with a Nikon SMZ1500 stereomicroscope to distinguish main morphological characteristics. Then, specimens were identified according to the morphological description and identification key of M. raddei adults in the book named "Economic Insect Fauna of China (Volume 1) Cerambycidae" (Chen et al., 1959), and the original description of M. raddei (Blessig, 1872 (Figure 1). ArcGIS platform version 10.2 was used to produce a distribution map based on the geographical coordinates of the localities. All samples were stored at −80°C before DNA extraction.

| DNA extraction, amplification, and sequencing
Genomic DNA was extracted from mesothorax muscles (approximately 0.1 g) after the epicuticle was removed, using the DNeasy Tissue & Blood Kit (QIAGEN) in accordance with the manufacturer's instructions. Extraction yield was measured using a Nanodrop-1000 Spectrophotometer . Subsequently, products were diluted to 30-40 ng/μL and stored at −20°C before polymerase chain reactions (PCR).
Based on the mitochondrial genome sequence of M. raddei (Wang et al., 2016), we designed specific primers for the partial sequences of the COI, Cytb, and COII using Primer version 3.0 (http://prime r3.ut. ee/) (Kõressaar et al., 2018), and these primers were synthesized by F I G U R E 1 Map illustrating the distribution and sampled localities of Massicus raddei (Blessig) in China. NC and SC population groups are labeled by black and red points, respectively. The sampling provinces and outbreak regions are labeled, and the population codes are shown in yellow words Invitrogen Biotechnology Co., Ltd. Details about primer sequences and product size are given in Table S1.
Polymerase chain reactions (PCR) were performed using a S1000 Thermal Cycler (BIO-RAD) in a total volume of 25 μl containing 12.5 μl 2× Green-Mix, 1 μl forward and reverse primer (10 μmol/L), 1 μl genomic DNA, and 9.5 μl ddH 2 O. PCR amplification was employed with an initial denaturation at 94°C for 5 min, followed by 35 amplification cycles consisting of 95°C for 30, 30 s at a primer-specific annealing temperature (52°C for COI gene, 53°C for Cytb and COII genes) and an extension at 72°C for 40 s, and then a final elongation at 72°C for 10 min. The products were visualized by 1% agarose gel electrophoresis and gel imaging system (UVP, Biolmaging Systems), and sequencing was implemented by Invitrogen Biotechnology Co., Ltd.

| Mitochondrial DNA statistical analysis
Sequencing data were analyzed and edited using Geneious R11.1.5 (https://www.genei ous.com; Kearse et al., 2012) and Bioedit version 7.2.5 (Clustal W method, default settings) (Hall, 1999). In order to corroborate that the correct gene fragments were obtained, sequence homology analyses were done using the BLASTN tool, which is available online from the NCBI (http://blast.ncbi.nlm.nih.gov.). Sequences of COI (679 bp), COII (369 bp), and Cytb (462 bp) of M. raddei were deposited in NCBI Genbank (Genbank accession numbers: COI, from MN796130 to MN796169; COII, from MN796170 to MN796209; Cytb, from MN796210 to MN796249). Subsequently, these three partial genes from the same individuals were concatenated to yield a combined sequence of 1,512 bp as a single locus because of the lack of recombination of mtDNA (Moraes et al., 2002). The nucleotide composition, number of conserved sites, variable sites, and parsimony informative sites of the three single partial gene sequences were then analyzed with Mega version 7.0 (https:// www.megas oftwa re.net/citat ions).

| Intrapopulation genetic diversity
Haplotypes were identified based on single and concatenated genes, respectively. Genetic diversity was estimated by number of haplotypes (n), haplotype diversity (Hd), nucleotide diversity (Pi), and the nucleotide average difference (K) for concatenated genes with the DnaSP version 5.0 (Librado & Rozas, 2009). Fst (using haplotype frequencies) was conducted with 10,000 permutations to estimate the hierarchical genetic structure using the Arlequin version 3.5 (Excoffier & Lischer, 2010). Only two hierarchical levels were defined: (a) among populations and (b) within populations. Following the criterion for genetic differentiation of Wright (1978), we defined genetic differentiation as low for Fst < 0.05, moderate for 0.05 < Fst < 0.15, high for 0.15 < Fst < 0.25, and very high for Fst > 0.25 (Govindaraju, 1989). With reference to the criterion for gene flow values, we defined genetic flow as low for Nm < 1, high for 1 < Nm < 4, and very high for Nm > 4 ( Boivin et al., 2004). The geographical distances were calculated according to longitude and latitude data using Geographic Distance Matrix Generator version 1.2.3 (Ersts, 2018). To examine whether any isolation-by-distance (IBD) effect occurred, matrices of pairwise genetic distances and geographical distances (km), matrices of genetic differentiation estimated by [Fst/(1 − Fst)], and the natural logarithm of geographical distance data (In Km) between all the sampling sites were analyzed for their degree of correlation by two mantel tests, with significance tests performed with 9,999 permutations. The analysis was implemented in GenAIEX 6.41 (Peakall & Smouse, 2012). TA B L E 1 Sampling information of seven geographical populations of Massicus raddei (Blessig)

| Phylogenetic relationships among populations and haplotypes
Applying Mega version 7.0 (https://www.megas oftwa re.net/citat ions.), ML phylogenetic analyses were used to identify major clades and evaluate the relationships among the seven populations by concatenated genetic sequences. An ML phylogenetic tree was constructed using General Time Reversible model with 1,000 bootstrap replications (Kumar et al., 2016). The Bayesian tree was obtained by running Mrbayes v3.1.2 (Ronquist & Huelsenbeck, 2003). The MCMC method was used to calculate 3 million generations, sampling once every 100 generations to ensure independence of sampling.
The original 3,000 trees were discarded as burn in.
Haplotype networks of M. raddei were generated using a median-joining algorithm in Network 5.0 (Bandelt et al., 1999), and the geographical distribution maps of haplotypes were visualized

| Demographic history and neutrality test
With the DnaSP version 5.0 software (Librado & Rozas, 2009), the Tajima's D and Fu's Fs values were estimated to test demographic expansion (Kimura, 1983;Tajima, 1989), and we assessed their significance with 1,000 permutations at levels of each geographical populations and overall population. The demographic history of whether M. raddei experienced a range expansion was also investigated by mismatch distributions of pairwise differences among sequences calculated with parametric bootstrapping (1,000 replications). The expected values were calculated assuming that sudden population growth model through coalescent simulations. According to simulations, demographically stable or admixed populations must present a multimodal distribution, whereas populations that have undergone a recent expansion generally exhibit a unimodal distribution (Harpending et al., 1998).

| Sequence analyses
Based on the specimen identification, all sampling adults were determined as the adults of M. raddei ( Figure S1). Forty partial COI, Cytb, and COII gene sequences of M. raddei with a length of 679, 462, and 369 bp were obtained, respectively. Nucleotide compositions among the seven geographical populations were basically consistent. The total content of A + T was much higher than that of G + C, which was in accordance with the feature of mitochondrial nucleotide composition (Jermiin & Crozier, 1994) (Table 2).   (Table 4). Moreover, AMOVA results indicated that the larger proportion of molecular genetic variation was found among populations (62.97%), and the remaining variation came within populations (37.03%). Exact tests showed significant variation on two levels (p < .001) ( Table 5).

| Genetic differentiation and genetic structure
Through two mantel tests, no significantly positive relationships between geographical distances and genetic diversification were found over seven geographical populations of China based on concatenated genes. R = .289, p = .138 and R = .257, p = .100 were detected between average pairwise genetic distances and geographic distances and between lnKm and Fst/(1 − Fst)) ( Figure 3).    (Figures 2 and 4).

| Demographic history
All Tajima Table 3). Distributions of pairwise differences (mismatch distributions) of overall populations through concatenated and single gene sequences presented a multi-peak distribution ( Figure 5), suggesting that the seven geographical populations of M. raddei did not undergo any demographic expansion recently (Harpending et al., 1998;Rogers & Harpending, 1992). Insignificant values on the basis of Tajima'D and Fu's Fs also supported this interpretation.

| D ISCUSS I ON
The present study used three genes of mitochondrial DNA to provide basic molecular information to allow better insight into factors contributing to genetic diversification among populations and possible mechanisms triggering the evolution of different patterns of adult emergence of M. raddei. High genetic diversity and distinct genetic differentiation among sampled populations were identified through the mitochondrial analyses, and the assumption that particular periodic behaviors and outbreaks of M. raddei in Liaoning Province mainly result from the absence of genetic flow with other populations and local genetic isolation was preliminarily denied.

| Genetic diversity
The largest nucleotide diversity and differences are detected in HN population, and four beetle specimens occupy one unique haplotype, respectively. All results provide a possibility that Xixia County might constitute the ancestral range of this species in Henan Province.
The homogeneity is achieved at small spatial scales, but limited dispersal ability reinforces genetic divergence over long distances (Peterson & Denno, 1998). We have reason to assume that there are positive correlations between genetic and geographic distances among populations of M. raddei. However, our current data were unable to obtain statistically significant results in mantel tests, which mainly resulted from our small sample size.
Furthermore, high genetic differentiation and obvious genetic diversity among the seven populations might conversely determine great variance in fecundity, larval capacity of degradation and nutrient acquisition along with adaptability to extreme environmental conditions, which further leads to the differences of regional population density and destructive levels to oak forest ecosystems. inces. By coincidence, among those regional populations, a periodic (once per 3 years) adult emergence phenology was also observed (EPPO, 2018). It is therefore assumed that there exists a causal relationship between the particular pattern of adult emergence and serious infestation. Previous studies have proposed that, within an insect species, the populations with annual emergence are considerably less dense than periodical populations (Bulmer, 1977;Lloyd & Dybas, 1966;Martin & Simon, 1990), and extraordinarily high densities and epidemic outbreaks in M. raddei populations of northeastern China are likely to be attributed to predation avoidance or predator satiation achieved by periodic behavior (Williams & Simon, 1995). In the present study, the assumption that serious periodic outbreak of population in Liaoning Province might be a result of local evolution in the absence of gene flow was preliminarily ruled out. Future researches should pay more attention on its biology, host-parasitoid interaction, and climatic factors, which may trigger proto-periodicity and favor the development and perfection of periodicity, in order to clarify the origin and maintenance mechanisms behind this interesting ecological phenomenon.
In conclusion, the present work, a preliminary attempt to investigate population structure across distributional areas of M. raddei in China. These results not only expand our genetic knowledge of this borer species, but also provide some new perspectives and theoretical basis for subsequent specialized studies in some regional populations. For instance, all spatial populations of Henan Province can be regarded as one research focus to evaluate the possibility of Xixia County as the ancestral center. More local populations of Shandong and Liaoning provinces should be collected to verify frequent gene flow and the lack of local genetic isolation between periodical and non-periodical populations once again, and then provide more detailed information to obtain a comprehensive understanding of contributing ecological and artificial factors shaping the phenomenon.
Moreover, future research employing more extensive samples and covering wider distributional sites, along with additional mitochondrial and nuclear molecular markers are needed to allow a better understanding of genetic structure among populations all over China.

ACK N OWLED G M ENTS
This work was supported by the Fundamental Research Funds of Chinese Academy of Forestry (CAFYBB2018ZB001). We thank Hao Liu and Xuewei Li for their assistance on specimen collection. We also thank Yingqiao Dang for her review and valuable comments on the early version of the manuscript.

CO N FLI C T O F I NTE R E S T
None declared.  Table S1.