When sampling each population, we chose three to five locations with an interval of more than 5 m and collected five individuals at least for each location. Collected samples were placed in absolute alcohol and stored at −80°C until DNA extraction. A total of 167 SAA adults were selected from 23 natural populations across China (Table 1 and Figure 1). Voucher specimens (NO. VAphA001–VAphA167) were deposited in the Entomological Museum of China Agricultural University. We obtained whole genomic DNA from individual specimen using the DNeasy Blood and Tissue Kit (Qiagen, Hilden, Germany) according to the manufacturer’s protocol. The quality and concentration of extracted DNA were measured using Qubit 3.0 (Invitrogen, United States). Each specimen was identified morphologically (Zhang and Zhong, 1983) and checked by BLAST search in the GenBank database using their COI gene sequence (with percent identity ≥ 99.00%), which was acquired through PCR sequencing under the following conditions: initial denaturation at 95°C for 1 min; 40 cycles at 98°C for 10 s; 45°C for 50 s; 68°C for 1 min; and a final single extension at 72°C for 2 min (with primers, reagents, and reaction system shown in Supplementary Tables 1, 2).

A 350 bp insert size library was constructed from each specimen and sequenced in 150 bp pair-end mode to obtain ∼3 Gb of data using an Illumina NovaSeq 6,000 platform (San Diego, CA, United States) at Berry Genomics (Beijing, China). Raw reads were trimmed off adapters using Trimmomatic (Bolger et al., 2014), and low-quality or short reads were discarded with Prinseq (Edwards, 2011). The remaining high-quality reads were mapped onto the published reference mitogenome (GenBank accession: MK766411; Liu et al., 2019) using Geneious 10.1.3,¹ allowing for mismatches of up to 2% with a maximum gap size of 5 bp and a minimum overlap of 30 bp. The consensus sequences were generated as novel individual mitogenomes under the threshold of 75% uniformity at each site. Sequences obtained from mapping were annotated in Geneious according to the reference sequence. Annotation of protein-coding genes (PCGs) and ribosomal RNA genes (rRNAs) was aligned with homologous genes of the reference mitogenome. Transfer RNA genes (tRNAs) were rechecked using the MITOS Web Server (Bernt et al., 2013). Each PCG was independently aligned using TranslatorX online platform with the MAFFT method, with all stop codons removed (Abascal et al., 2010), and the PCG dataset was finally obtained. Besides, we also performed ABGD (Automatic Barcode Gap Definition) analysis with the aligned PCG dataset (Pmin = 0.001, Pmax = 0.02) and confirmed that all sequences belong to the same species T. trifolii. Analyses in this study were all based on the PCG dataset.

The number of polymorphic sites (S), the number of haplotypes (Nh), nucleotide diversity (π), and haplotype diversity (Hd) were calculated using DnaSP 6.0 (Rozas et al., 2017). To check the differentiation levels in different populations and lineages, we calculated pairwise genetic differentiation (F_ST) and obtained p-values for F_ST values with 3,000 permutations using Arlequin 3.5 (Excoffier and Lischer, 2010), and calculated genetic distance (uncorrected p-distance) using MEGA7 (Kumar et al., 2016). Mantel tests were used to test the isolation-by-distance model for all geographic populations except Changde, Hunan population (HNCD) including only two samples. Correlations between F_ST and linear geographic distance (in km) were plotted with 999 replicates in the ade4 v. 1.7 module of R.

Phylogenetic trees were constructed by the Bayesian inference (BI) method using PhyloBayes MPI 1.5a (Lartillot et al., 2013) and the maximum-likelihood (ML) method using IQ-TREE 1.6.5 (Trifinopoulos et al., 2016). PartitionFinder2 (Lanfear et al., 2016) on the CIPRES Science Gateway (Miller et al., 2010) was used to select the optimal partition schemes and substitution models. Thirty-nine (13 × 3) codon partitions were predefined for this dataset, and we used the “greedy” algorithm with branch lengths estimated as “unlinked” and the Akaike information criterion to select the partition schemes and models (summarized in Supplementary Table 3). For the ML analysis, the best evolutionary models predicted before were selected under the corrected Bayesian information criterion (BIC) in IQ-TREE, and the ML tree was obtained with IQ-TREE using an ultrafast bootstrap approximation approach with 10,000 replicates. For BI analysis, we chose the site-heterogeneous mixture model (CAT + GTR), with four independent chains starting from a random tree and running for 30,000 cycles. Trees being sampled every 10 cycles and the initial 25% trees of each MCMC run were discarded as burn-in. A consensus tree was computed from the remaining trees combined from two runs which converged at a MaxDiff < 0.1.

Clustering of individuals was performed with BAPS v. 6.0 (Cheng et al., 2013) under the spatial clustering of groups of individuals. Based on the same dataset, the median-joining network of haplotypes was constructed in PopART (Leigh and Bryant, 2015) to further analyze the relationship among different haplotypes. The frequency of haplotypes was obtained using Arlequin software. We also carried out the hierarchical analysis of molecular variance (AMOVA) among populations and between defined lineages (the result of BAPS) with the Arlequin program to gain a better insight of the biogeographical processes shaping the genetic structure. The significance of differentiation between pairs of populations/groups was evaluated by 1,000 bootstrap replications (Excoffier et al., 1992).

We estimated divergence time using BEAST 2.4.8 (Bouckaert et al., 2014). Given that there were no fossils available to calibrate the nodes, we used a molecular clock for divergence time estimation. COI gene was evaluated as a separate partition with the substitution rate fixed at 1.77% per site per million years (Mya) (Papadopoulou et al., 2010). The other 12 PCGs were applied as another partition with the evolutionary rate estimated accordingly. An uncorrelated lognormal relaxed clock model was applied with the birth-death model for the tree prior. Two independent MCMC runs per 500 million generations were performed with tree sampling every 50,000 generations. Tracer v. 1.7 (Rambaut et al., 2018) was used to verify whether the MCMC runs reached a stationary distribution based on the effective sample sizes (ESS) of each estimated parameter. We required ESS for the posterior, prior, and tree likelihood to be at least 200. The resulting trees for each replicate were combined using LogCombiner program in BEAST, and TreeAnnotator program in BEAST was used to calculate the consensus tree and annotate the divergence times with “mean height” after discarding the initial 25% trees as burn-in.

To test whether expansion occurred, we performed mismatch distribution analyses in Arlequin with 1,000 bootstrap replicates to detect expansion through the linear fitting of observed and simulated curves with the significance of two parameters, the sum of squares deviations (SSD) and Harpending’s raggedness index (r) (Excoffier and Lischer, 2010). Bayesian skyline plots (BSPs) were obtained for different lineages in BEAST under Coalescent Bayesian Skyline model for the tree prior. The same partitions and substitution rates were applied consistent with these in the estimation of divergence time. BEAST analysis was conducted for 300 million generations with sampling trees every 30,000 generations. We determined the effective population size through time using Tracer, discarding the initial 10% of generations as burn-in.

We obtained the whole mitogenome sequences for 167 individuals from 23 SAA populations across China (Table 1), with the size of mitogenome ranging from 16,068 to 16,083 bp. We detected 87 haplotypes in total (Table 1 and Supplementary Table 4), with 202 polymorphic sites (S) in average. Overall, populations showed a high haplotype diversity (Hd = 0.976) and a low nucleotide diversity (π = 0.00197). Among all the different geographic populations, GSQY showed the highest nucleotide diversity (π = 0.00253) and three populations showed lowest (π = 0.00000), including CQBB, HNCD, and YNBC.

We estimated pairwise F_ST and uncorrected p-distance between 23 different geographical populations (Table 2). The F_ST value between CQBB, HNCD, and YNBC was the highest (1.0000), and the lowest was between HNCD and NGNZ (−0.2310). For the p-distance, the value between HNZZ and YNZK was the highest (0.0038), and the lowest was between SCLS and YNBC (0.0000). Populations from northern China showed low p-distance values from each other. HBSY, HNZZ, and XJCJ populations showed higher F_ST and p-distance values than those from other populations. We also calculated the p-distance between 87 haplotypes, and the results ranged from 0.0001 to 0.0045 (Supplementary Table 5). The Mantel test indicated that the correlation between geographic distance and genetic differentiation was not significant (r = −0.0374, p = 0.7120, Supplementary Figure 1).

Phylogenetic relationships of major nodes in the ML, BI, and BEAST trees showed similar topologies, even though the nodal supports of some pivotal nodes were not high (Figure 2). The monophyly of the three lineages was well supported: the first lineage included all samples from XJCJ and one individual from another lineage, the second lineage included samples mainly from middle China, and the last lineage included samples from various populations, ranging from northeast to southwest of China. Therefore, we assigned the name of lineages according to their main geographic region, i.e., Xinjiang (XJ), middle China (MC), and other populations within China (OC) lineages, respectively (Figure 2). In MC lineage, most individuals from Ningxia area diverged first followed by samples from nearby cities. Individuals of OC lineage that diverged first were all from north and northeast of China, and individuals that diverged later were mostly from northwest and southwest of China. BAPS result also showed that all samples were divided into three clusters (Figure 1); XJ cluster (blue color in Figure 1) was only distributed in Xinjiang province, the MC cluster (red color in Figure 1) was mainly distributed in the north and middle of China, and the OC cluster (green color in Figure 1) was distributed as a majority (the haplotype from this cluster can be found in every geographical population).

Among the lineages revealed by the phylogenetic analyses, OC lineage with the largest sample size showed the lowest diversity, followed by MC lineage (Table 3), and XJ lineage exhibited the highest values even though it comprised the fewest individuals (Table 3). The p-distance among three lineages varied from 0.0031 to 0.0041 (Table 4). AMOVA analysis showed that the greatest proportion of variation was among lineages (81.79%) and the genetic differentiation among all of the populations was F_ST = 0.877, p < 0.0001 (Table 5).

The median-joining network of haplotypes revealed three haplotype groups (Figure 3), consistent with the results of BAPS and phylogenetic trees. For the OC lineage, there were several star-like topologies observed, indicating the recent expansion in this lineage. Haplotype H6 was the most frequently shared haplotype, which included 18 individuals from 10 different populations, followed by haplotype H9, which included 10 individuals from three populations. There were many mutation steps observed between the three lineages in median-joining network.

Our divergence time estimation suggested that the most recent common ancestor (MRCA) of SAA was 83.2 Ka (kilo-annum) and diverged into the MRCA of XJ and MC, and OC lineages with the ages of 75.3 and 27.7 Ka (Figure 1). The MRCA ages of XJ and MC were 63.5 and 25.1 Ka, respectively. Thus, all three lineages diverged during the late Pleistocene Quaternary.

In the mismatch distribution analyses (Figure 4), the observed and simulated curves did not significantly differ from each other in OC lineage with the parameters of SSD and r being small and insignificant (P > 0.05). Hence, the population expansion process was supported in this lineage. In contrast, XJ and MC lineages exhibited relatively large and significant parameters of SSD and r. Thus, the population expansion event was not supported in these two lineages. BSP analysis (Figure 4) revealed that OC lineage and total samples experienced a significant increase in effective population size about 1500–2000 years ago, which was a very recent expansion event. However, the population size of MC and XJ groups showed a slight decrease after a long steady period.