The damaged branches of host trees were collected in wild fruits forests and orchards from May to August 2018 to 2021. Branches were then cultured in vitro at room temperature. Once adult beetles emerged, they were placed into the centrifugal tubes and stored at −80°C prior to DNA extraction. Genomic DNA was extracted from an individual without the abdomen and elytra using a TIANamp Genomic DNA Kit (TIANGEN, Beijing, China) according to the manufacturer’s instructions. Voucher specimens were deposited at the insect museum of the Chinese Academy of Forestry, Beijing, China. Hundreds of A. mali adult individuals were obtained from 17 locations in China (Table 1; Figure 1).

Three mitochondrial protein-coding genes (COI, COII, and CytB) were utilized as molecular markers in the present study. Based on the published mitochondrial genome sequence of A. mali (Sun et al., 2020), we designed specific primers for the amplification of complete COI, COII, and CytB sequences using the software Primer Premier version 5 (Lalitha, 2000). Details of the primer sequences and product sizes are provided in Supplementary Table S1.

Polymerase Chain Reaction (PCR) amplification was conducted in a total volume of 25 μl containing 1.5 μl genomic DNA, 1 μl each of forward and reverse primers (10 μmol/L), 12.5 μl 2× Flash PCR MasterMix (CWBIO, China), and 9 μl ddH₂O. PCR amplification was performed with an initial denaturation at 98°C for 1 min; 35 cycles of 10 s at 94°C, 15 s at a primer-specific annealing temperature (46°C for COI, 49°C for COII, and 53°C for CytB), and extension at 72°C for 30 s; and a final elongation at 72°C for 1 min. The PCR products were visualized on 1.0% agarose gels under UV light. The sequencing reaction was carried out using an ABI BigDye Terminator version 3.1 cycle sequencing kit on an ABI 3730XL (Tsingke Biotechnology Co., Ltd., Beijing, China) for Sanger sequencing in both directions.

To confirm that the correct sequences were obtained, we compared each sequence with the non-redundant (nr) database of National Center for Biotechnology Information (NCBI) using BLASTN. Three mitochondrial gene sequences were analyzed and edited using Geneious R11.1.5 (https://www.geneious.com) (Kearse et al., 2012) and MEGA version 11 (Tamura et al., 2021). Deletions or insertions in the nucleotide sequences of mitochondrial protein-coding genes (COI, COII, and CytB) were inferred from the amino acid alignment. The sequence of each PCG (excluding stop codons) was then aligned individually as per codon-based multiple alignments using the MAFFT algorithm, which was implemented in the TranslatorX online platform (Abascal et al., 2010). Subsequently, the three genes of an individual were concatenated using MEGA version 11 to generate a combined sequence of 2,913 bp as a single locus because of the lack of recombination of mtDNA (Moraes et al., 2002). The nucleotide composition and the number of variable sites, conserved sites, parsimony informative sites, and singleton variable sites of the three single partial gene sequences, as well as of the combined sequences were analyzed using MEGA version 11.

Genetic diversity was estimated by the number of polymorphic sites (S), the number of haplotypes (h), haplotype diversity (Hd), average number of nucleotide differences (K), and nucleotide diversity (π) for concatenated genes, using DnaSP version 5.0 (Librado and Rozas, 2009) or Arlequin version 3.5 (Excoffier and Lischer, 2010). Besides the analysis of each Tianshan population alone, we combined the nine populations in these regions as a single group (YL) to examine their genetic diversity for comparison with non-Tianshan populations. For the comparisons, we additionally divided the Tianshan populations into natural forest populations (YLNA) and plantation populations (YLPL) (Table 3).

The haplotype sequences for each gene were exported by DnaSP version 5.0. The haplotype sequences were deposited in GenBank under the accession numbers OM177028–OM177057 for COI, and OM212329–OM212346 for COII. All CytB sequences were also deposited in GenBank under the accession numbers OM264916–OM265120, and haplotype sequences are available in the Zenodo database with a DOI of https://doi.org/10.5281/zenodo.6107571.

Based on combined genes sequences, pairwise fixation (F _st ) and gene flow indices (N _m ) were estimated using the DnaSP version 5.0 (Librado and Rozas, 2009). With the criterion for genetic differentiation set by Wright (1984), the genetic differentiation is low for F _st < 0.05, moderate for 0.05 < F _st < 0.15, high for 0.15 < F _st < 0.25, and very high for F _st > 0.25 (Govindaraju, 1989). Similarly, genetic flow was low for N _m < 1, high for 1 < N _m < 4, and very high for N _m > 4 (Boivin et al., 2004).

The analysis of molecular variance (AMOVA) based on F _st (using haplotype frequencies) was performed using 10,000 permutations to estimate the hierarchical genetic structure using Arlequin version 3.5 (Excoffier and Lischer, 2010). Two hierarchical levels, among populations and within populations, were defined.

To test isolation by distance (IBD), the pairwise genetic distance and geographical distance (km) matrices between all 17 sampled populations were compared using the Mantel test with 9,999 permutations (Mantel, 1967). Pairwise average genetic distances between different populations were calculated using the Kimura-2-parameter model in MEGA version 11 (Tamura et al., 2021). Geographical distances were calculated from longitude and latitude data using the Geographic Distance Matrix Generator version 1.2.3 (Ersts, 2018). Mantel analysis was performed using GenAleX version 6.5 (Peakall and Smouse, 2012).

The dataset of mitochondrial haplotypes was analyzed under Bayesian inference (BI) using MrBayes version 3.2.5 (Ronquist and Huelsenbeck, 2003). The congeneric species of Agrilus planipennis were used as an outgroup. For Bayesian analyses, we analyzed the combined genes based on the GTR + I + G model, which was selected using jModelTest version 2.1.7 (Posada, 2008). Then, two simultaneous runs of 10,000,000 generations were conducted for the matrix, each of which was sampled every 1,000 generations with a burn-in of 25%. Trees inferred prior to stationarity were discarded as burn-in, and the remaining trees were used to construct a 50% majority-rule consensus tree. The relationships among haplotypes were further tested by creating a median-joining algorithm using PopART version 1.7 (Leigh and Bryant, 2015).

To explore the demographic history of different lineages, neutral tests using Tajima’s D (Tajima, 1989) and Fu’s Fs (Fu, 1997) for the defined populations were performed using DnaSP version 5.0 (Librado and Rozas, 2009), with significantly negative values indicating population expansion and significantly positive values indicating processes such as population subdivision or recent population bottlenecks (Kimura, 1983; Tajima, 1989).

Additionally, the distribution of pairwise differences (mismatch distributions) between each pair of sequences was analyzed using DnaSP version 5.0 (Librado and Rozas, 2009) with default parameters to detect population expansion through a linear fitting relationship between the observed and simulated curves and infer demographic equilibrium or expansion according to Rogers and Harpending (1992).

The final combined mitochondrial dataset included 2,913 bp of protein-coding regions (COI:1,395 bp, COII:582 bp, CytB:936 bp). No insertions or deletions were detected in these fragments. To ensure the credibility of the data, we deleted a 186 bp sequence containing the poly-T base in the CytB gene. The nucleotide compositions of the 17 geographical populations were consistent. The obtained combined sequences showed base compositions similar to those reported previously in mitochondrial sequences for other insects, with a heavy bias towards A/T (68.7%), which was in accordance with the mitochondrial nucleotide composition (Jermiin and Crozier, 1994).

Every mitochondrial gene displayed several haplotypes, with a total of 30 for COI, 18 for COII, and 37 for CytB, among which 7, 5, and 24 were shared, respectively. Details of the haplotype distribution of mitochondrial COI, COII, CytB, and their combined genes in each A. mali population are shown in Supplementary Tables S3–S6. Among the combined genes, 151 polymorphic sites (including 134 parsimony sites and 17 singleton sites) and 54 haplotypes (including 37 unique haplotypes) were found among all samples (Tables 2, 3). For a single population, the average number of haplotypes (h) ranged from 2 (HM) to 13 (YLK); the average number of nucleotide differences (K) ranged from 0.667 (HM) to 28.974 (CY), with an average of 10.933 (Table 3); haplotype diversity (Hd) within populations ranged from 0.590 (PL) to 0.963 (YLK), with an average of 0.793; and nucleotide diversity (π) ranged from 0.00023 (HM) to 0.00955 (CY). The whole Chinese population exhibited high haplotype diversity and high nucleotide diversity (Hd = 0.924; π = 0.00752).

For the western Tianshan Mountains population, the highest genetic diversity was found in the YLK population (Hd = 0.963, π = 0.00398) from the natural forests, while the lowest (Hd = 0.591, π = 0.00172) was found in the YLN population from plantations. Overall, the genetic diversity of natural forest populations is slightly higher than that of plantation populations (Table 3).

The pairwise F _st values for genetic differentiation varied from −0.0866 to 0.9428, based on the genes combined. The maximum value was detected between the HM and FX populations, and the minimum value was detected between the JZ and YLB populations. Five populations (CF, CY, FX, HM, and PL) exhibited significant genetic differentiation from other populations. The SY population showed higher genetic differences with other populations, whereas there was a certain gene flow between the BY population and the JZ, YLB, YLK, and YLZS populations, respectively. Additionally, no genetic differentiation was observed between the JZ population and the YL group. The N _m values indicated that frequent genetic flow occurred between the YL group and the JZ population (Table 4).

AMOVA showed that overall variation mainly derived from variation among populations (53.25% among 17 populations and 63.09% among nine groups), and genetic variation within populations was slightly lower (Table 5). The Mantel test revealed a significant positive relationship for all populations (R = 0.5590 and p = 0.0001) (Figure 2A) while no significant correlation for the western Tianshan Mountains populations (R = 0.6377, p = 0.952) (Figure 2B) between the average pairwise genetic distance and geographic distance, indicating that genetic differentiation in populations is probably caused by geographical isolation.

Based on the combined genes, the BI phylogenetic trees displayed three parallel clades: Clade 1 containing YL, JZ, SY, BY (these four populations formed the biggest subdivision), and CY populations (another subdivision); Clade 2 containing CF; Clade 3 containing CF, FX, HM, and PL (Figure 3). Within Clade 1, the biggest subdivision was divided into three clusters. Two of them contained pure YL haplotypes despite sharing several haplotypes with other populations, while the other was mixed with SY population.

Median-joining networks inferred from the combined genes also supported similar phylogenetic relationships. Based on combined genes, the median-joining haplotype network for all populations did not exhibit the ancestral haplotype of A. mali across overall populations (Figure 4A), whereas three ancestral haplotypes (H1, H19, and H21) for the network of the YL group were found (Figure 4B). These three ancestral haplotypes accounted for more than half (73/129) of the individuals in the YL group. Furthermore, the ancestral haplotypes of the YL group were more apparent in single genes, namely H1, H16, and H18 in COI; H1, H13, and H14 in COII; and H1, H14, and H15 in CytB. In contrast, the median-joining network of the haplotypes based on a single gene of COII was different from the results obtained using COI, CytB, and combined genes (Figure 4; Supplementary Figure S1), which resulted from a lower nucleotide diversity and less genetic variation caused by the shorter COII sequence.

For the neutral test, despite one significantly positive Tajima’s D value in the YLH population (2.38509, p < 0.05), most Tajima’s D and Fu’s Fs values (Table 3) in the overall population displayed no significance, suggesting that the Chinese population, including the Tianshan population from natural forests, did not display population expansion (Rogers and Harpending, 1992; Harpending et al., 1998).

Distributions of pairwise differences obtained with either single or combined genetic data across populations exhibited declined simulated curves (Figure 5A; Supplementary Figure S2). This suggested that the demographic history of A. mali populations in China fitted equilibrium models and that these insects did not experience recent population expansion. The same result was also obtained among Tianshan populations (Figure 5B).