Human- mediated eco- evolutionary processes of the herbivorous insect Hyalopterus arundiniformis during the Holocene

Aim: Pleistocene climatic fluctuations have been identified as a dominant factor in evolutionary and ecological processes underlying contemporary species distribution. In contrast, anthropogenic activity during the Holocene has been largely neglected in phylogeographic inference. Here, by integrating archaeological evidence, phylogeographic analyses and a climatic niche- based approach, the role of Holocene human activity in determining the evolutionary trajectory of insect pests is investigated using Hyalopterus arundiniformis (HAR) as a model. Eurasia. Methods: We compiled archaeological records and historical data for the host plant (i.e. peach) of HAR. Phylogeographic structure, genetic diversity and gene flow were estimated. Historical demographic dynamics were explored using coalescent-based analyses. The origin and dispersal history were assessed using an approximate Bayesian computation (ABC) approach. The changes in climatic niche breadth during geographic expansion were evaluated via niche comparison analyses. Results: Our results revealed that HAR originated in Southwest China and dispersed to Southeast and Northwest China, and it was then introduced to western Eurasia from Northwest China. The climatic niche shifted significantly with range expansions, and niche expansion explained the largest proportion of niche dynamics. Population size expansion occurred in the Holocene and coincided with the domestication and trans - Eurasian trade of peach.


| INTRODUC TI ON
Contemporary biodiversity on Earth has been profoundly impacted and shaped by the climate oscillations associated with Quaternary glaciations (Gavin et al., 2014;Hewitt, 2000Hewitt, , 2004. Nevertheless, during the Holocene, human activity also impacted diversification and dispersal dynamics for many organisms and triggered both local and global changes in diversity patterns (Bernardo-Madrid et al., 2019;Kemp et al., 2020), impacts that have been largely overlooked in phylogeographic inference. In the Neolithic period [10,000-5000 years before present (BP)], with the transition from hunting and gathering to a farming society, a broad array of animals and crops were domesticated. Subsequently, with the global expansion of agropastoral economies, domestic livestock and crops were spread worldwide along with human activities (Frantz et al., 2020;Purugganan, 2019;Wang et al., 2014). Through domestication and translocations, Holocene human activity changed effective population sizes (Ne) (Purugganan, 2019) and broadened the climatic niche of domestic animals, a process known as the anthropogenic niche construction (Boivin et al., 2016;Kemp et al., 2020;Phelps et al., 2019). Changes in traits under domestication could alter ecological and evolutionary processes within biotic communities, but the effects of domestication on higher trophic levels have not been well documented . In addition, domestic livestock and crops also act as a dispersal vector for other organisms (Auffret et al., 2014;Gippet et al., 2019), such as parasites and herbivorous insects.
Indeed, human-mediated dispersal of crops has also significantly promoted global insect invasions due to the ongoing globalization of trade and travel (Gippet et al., 2019). Therefore, knowledge of human-driven processes, especially for herbivorous insects, will provide a better understanding of the dispersal and species interaction mechanisms.
The Eurasian continent is characterized by heterogeneous geographical features and diverse climatic and environmental conditions. The western and eastern parts of the Eurasian continent are at similar latitudes but have distinct geological histories and climatic conditions, especially for the Pleistocene epoch (Hughes et al., 2013). Previous studies on species widely distributed across Eurasia mainly focussed on the population structure, demographic expansion and variable responses to Pleistocene climatic fluctuations of regional populations (e.g. Song et al., 2016Song et al., , 2020. Nevertheless, this region was also one of the cradles of Neolithic agriculture (Jones & Liu, 2009;Zeder, 2008) and played an essential role in the development and interaction of human civilizations (e.g. Silk Road; Liu, 2010). Based on archaeological and genetic studies, human-mediated movement of farm animals and domesticated plants during the Neolithic era in western or eastern Eurasia has been well documented (e.g. Krak et al., 2019;Ottoni et al., 2013). More recent studies have further proven that human dispersal and Eurasia-wide trade play an essential role in contemporary species distributions (Feuerborn et al., 2021;Kawakubo et al., 2021). Thus, the Eurasian region provides an ideal setting for disentangling how climate-or human-mediated factors have shaped contemporary distributions and past demographic dynamics, which will be fundamental for anticipating their impacts on biological control and future scenarios.
Hyalopterus arundiniformis Ghulamullah (HAR) (Aphididae: Aphidinae) is widely distributed across Eurasia  and has been considered a serious pest that inflicts severe damage on agriculturally important fruit peach (Prunus persica; Blackman & Eastop, 2020;Favret et al., 2017;Lozier et al., 2008). The species is heteroecious and holocyclic, mainly utilizing peach as the primary host for sexual generation and migrating to reeds (Phragmites) as the summer (secondary) hosts for asexual reproduction (Blackman & Eastop, 2020). Based on broad sampling and multispecies coalescent models,  further revealed that HAR contains a small sub-lineage specializing on plum (Prunus salicina) and red plum (Prunus cerasifera f. atropurpurea) and shows morphological differentiation. However, the demographic dynamics and evolutionary history of HAR in Eurasia remain uninvestigated, leaving considerable uncertainty about its migration and expansion.
Host plants are believed to have played a pivotal role in shaping aphid speciation and dispersal (Fang et al., 2018;Liu et al., 2022;Zhang et al., 2012). Evidence from archaeological and genomic studies has suggested that the primary host of HAR (i.e. peach) originated in Southwest China (Su et al., 2015;Wang & Zhuang, 2001;Yu et al., 2018). The domestication of peach dates back as far as 7500 years ago in the lower Yangzi River valley .
Domesticated peach was also spread across Persia (present-day Iran) via the ancient Silk Road in the first centuries before Christ (BC) and further westward into Europe by the Romans (Faust & Timon, 1995).
After several thousand years of domestication and trade, peach is now widespread worldwide and strongly associated with anthropogenic habitats (e.g. roadsides, fields and parks). The cultivation of peach also greatly expanded the distribution of herbivorous insects, such as Myzus persicae (Blackman & Eastop, 2020;Margaritopoulos et al., 2009). Nevertheless, the eco-evolutionary effects of the domestication and Eurasia-wide trade of peach on the demographic history, migration and expansion of herbivorous insects are poorly understood.
Here, based on comprehensive analyses, the eco-evolutionary effect of crop domestication on herbivorous insects is investigated using HAR as a model. We specifically test three hypotheses. First, we test whether human-mediated dispersal of domesticated peach K E Y W O R D S anthropogenic activity, cryptic invasion, herbivorous insect, Holocene, human niche construction, migratory episodes, silk road has significantly expanded the distribution of HAR far beyond the ancestral range. Specifically, we predict that the geographic origin and migratory episodes of HAR were closely associated with peach.
Second, we predict that the Ne of HAR expanded in the Holocene, and the expansion of population size coincided with the domestication and trans-Eurasian trade of peach. Finally, we predict that the climatic niche of HAR broadened during range expansions.

| Sampling and laboratory work
In total, 303 samples of HAR were collected from 2002 to 2021 during field surveys (Table S1). According to the Quick-Start guidelines, total genomic DNA was obtained from one aphid individual per sample using the DNeasy Blood and Tissue Kit (Qiagen, Germany).
Four mitochondrial (mtDNA) fragments, namely cytochrome c oxidase subunit I (COI), tRNA-leucine + cytochrome oxidase c subunit II (tRNA/COII), cytochrome oxidase b (Cytb) and 12S ribosomal RNA (12S rRNA), and one nuclear marker (nuDNA), namely elongation factor-1ɑ (EF-1ɑ), was amplified. The five genes are widely utilized in aphid phylogenetic and phylogeographic studies (e.g. Fang et al., 2017;Fang et al., 2018;Liu et al., 2022;Zhang et al., 2012). All primers and the PCR procedures used in this study are summarized in Table S2. Raw amplified fragments were checked manually for insertions or deletions (indels) and ambiguous positions using SeqMan II (DNAStar, USA). All protein-coding sequences were translated into amino acid sequences to exclude the stop codons by EMBOSS Transeq (Li et al., 2015). Furthermore, the published sequences of HAR deposited in GenBank were retrieved (Table S1). Ultimately, 326 samples of HAR were included, which spanned much of the known geographical ranges across the whole Eurasian continent.
The introns in the EF-1ɑ gene were included in the later analyses.
The DNA sequences of all markers were aligned with MEGA X (Kumar et al., 2018) and subsequently produced the sequences of 658 bp of COI, 739 bp of tRNA/COII, 745 bp of CytB, 339 bp of 12S rRNA and 943 bp of EF-1ɑ. The sequences generated in this study were submitted to the GenBank database, and the accession numbers are provided in Table S1.

| Archaeological records and historical data for peach
It is essential to collect archaeological and historical records to gain insights into the distribution and dispersal processes of species (e.g. Deng et al., 2020;Zhao et al., 2017). However, there is a paucity of archaeological and historical evidence for HAR. As a specialist herbivorous insect, HAR is closely associated with its primary host plant. Therefore, we compiled information from the published literature about the historical distribution of peach. Specifically, we conducted a literature review by searching for publications mentioning peach stones from different archaeological sites. The historical distributions of peach in East Asia during the late Holocene were obtained from Wang and Zhuang (2001) and Huang et al. (2008), who compiled data from historical literature [such as Shi Jing (1100-600 BC) and Han Fei Zi (280-233 BC)] and local gazetteers. Moreover, we compiled historical records of peach in western Eurasia from Faust and Timon (1995). Subsequently, all records were synthesized in terms of archaeological site, time and location. We used the upper limit of the estimated time for archaeological sites and historical records for visualization using ArcGIS 10.2 (ESRI, Redlands, CA, USA). The detailed dispersal routes of peach summarized by Wang and Zhuang (2001) (Figure 1b) were also retrieved. Combing the dispersal routes and archaeological and historical records revealed the temporal-spatial changes in the origin, domestication and trade of peach. The demographic dynamics revealed by the genetic analyses and the historical records of peach were used to delineate the origin and migration scenarios of HAR in Eurasia. In addition, according to the spread of peach and the archaeological evidence, five regions were delimited to facilitate subsequent analysis (detailed information is provided in the Supporting Information).

| Phylogenetic relationships and population structure
Phylogenetic relationships were reconstructed based on concatenated mtDNA and nuDNA datasets with Hyalopterus pruni and Rhopalosiphum padi as outgroups. Phylogenetic analyses were carried out with Bayesian inference (BI) methods implemented in MrBayes 3.2.6 (Ronquist et al., 2012) and maximum likelihood (ML) frameworks using RAxML 7.2.6 (Stamatakis, 2014). We used PartitionFinder 2.1.1 (Lanfear et al., 2017) to identify the best-fitting partitioning strategy and evolution models with the greedy algorithm and the Bayesian information criterion (BIC). The resulting partitioning scheme and substitution models are provided in Table   S3. Two independent runs of the BI analysis with four incrementally heated Metropolis-coupled Monte Carlo Markov chains (MCMCs), including one cold and three heated, were conducted for 10 million generations and sampled every 1000 generations, with the first 25% of the sampled generations discarded as burn-in. The ML analyses were inferred with a GTRCAT evolution model using a rapid 1000 thorough bootstrap replicates. The phylogenetic analyses (including RAxML and MrBayes analyses) were performed on the CIPRES Science Gateway (Miller et al., 2010). For mtDNA, a median-joining network method implemented in Network 10 (Bandelt et al., 1999) was applied to reveal the relationships among haplotypes.
To assess genetic diversity among geographic regions, we calculated the nucleotide diversity (π) and haplotype diversity (Hd) of mtDNA using DnaSP 5.0 (Librado & Rozas, 2009). Furthermore, isolation by distance (IBD) was examined using a Mantel test in the R package "vegan" 2.5-7 (Oksanen et al., 2019) with 999 permutations to evaluate the correlation between pairwise genetic distance and spatial distance. Genetic distance was estimated as the net average

| Gene flow analysis
The migration rates (M) and effective population size (θ) were estimated based on mtDNA datasets using MIGRATE-N 4.4.3 (Beerli & Palczewski, 2010) to evaluate the dispersal pattern among geographic regions. This method used a Bayesian coalescent approach to assess asymmetric gene flow between two populations. The MCMC strategy parameter was setting as following: long-chains = 1, long-inc = 100, long-sample = 1,000,000, burn-in = 100,000, heating = YES: 1: {1.0, 1.5, 3.0, 1,000,000}, and replicate = YES: 5. N e m (the effective number of migrants from each population per generation) was calculated by Mθ. The analyses were run two independent times with different random seeds to ensure the consistency of the results.

| Demographic history analyses
We used the coalescent-based Bayesian skyline plot (BSP) approach implemented in BEAST 2.6.6 (Bouckaert et al., 2014) to infer the changes in Ne since the time to the most recent common ancestor (TMRCA) based on mtDNA datasets. This method considers uncertainty in genealogy using coalescence theory and simultaneously evaluates Ne over time using an MCMC approach. The substitution models were selected based on the result of PartitionFinder (Table   S3). A substitution rate of 0.0177 per site per million years for COI was employed (Papadopoulou et al., 2010) and used to estimate rates for other genes. A relaxed lognormal clock model and a coalescent Bayesian skyline tree prior were selected for each clade. The MCMC analysis was run for 200 million generations, with one tree sampled every 20,000 generations. Effective sample sizes greater than 200 were considered as convergence. A BSP was reconstructed with Tracer 1.7 (Rambaut et al., 2018).

| ABC-based inferences about the origin and migration scenarios
Approximate Bayesian computation (ABC) implemented in DIYABC 2.1.0 (Cornuet et al., 2014) was then applied to infer the geographic origins and migratory episodes of HAR based on mtDNA genes. This method combines simulations and summary statistics to evaluate complex evolutionary scenarios (Cornuet et al., 2014). Different strategies were applied to delimit population units in previous ABC analyses, with some selecting population units identified by Bayesian clustering methods and others using geographic regions (Stone et al., 2017). Here, we aimed to understand whether the migratory history of HAR was associated with the dispersal of its primary host.
Hence, based on information from the archaeological records and migration routes of peach, the populations were delimited into five geographic regions ( Figure S1). For five groups, as many as 120 scenarios could be evaluated. However, it is computationally difficult or impracticable to compare all 120 scenarios in a single set. Thus, a hierarchical approach was implemented by optimizing the process into two successive steps to simplify the competing scenarios (Lombaert et al., 2014). In step 1, three major scenarios for the origin of HAR were assessed: south-western (Models 1-3), south-eastern (Models 4-6) and north-western (Models 7-9). Nine possible models were compared, with each scenario including three sub-models to investigate the relationships among the four groups distributed in East Asia. In step 2, hypotheses about the origin of Central Asia and Europe (AE) group were proposed based on the best-fitting demographic model selected in the last step. Based on genetic diversity, phylogenetic relationships and phylogeographic structure, a total of four biologically plausible scenarios were tested: groups AE and Northwest China (WE) emerged simultaneously from a common ancestor in Southwest China (SW) (Model 10), and group AE emerged from WE and Northeast Asia (NE), respectively (Models 11-13). In addition, we rerun the analysis excluding these individuals to assess the effect of the host-specialized lineage from plum and red plum on the reconstruction of the migration history. The prior distribution of demographic parameters was set as uniform and is listed in Table   S4. We selected HKY as the mutation model and used all possible summary statistics. We performed one million data simulations to build a reference table for each model and compared the models by calculating the relative posterior probabilities based on logistic regression. The model with the highest posterior probability value and nonoverlapping 95% confidence intervals (CIs) was selected as the best model (Lombaert et al., 2014;Stone et al., 2017).

| Niche comparison analyses
We compared the climatic niches in native and nonnative ranges based on the methodological framework of Broennimann et al. (2012) to evaluate changes in climatic niche breadth during the geographic spread of HAR. The results of the ABC-based analysis were utilized as the geographic extent of the native and nonnative ranges of HAR: F I G U R E 1 Integrated diagram of the inferred demographic history of Hyalopterus arundiniformis in Eurasia. (a) Archaeological and historical records of peach across Eurasia. Triangles and circles indicate archaeological records and historical evidence, respectively. (b) Schematic summarizing the routes of the Silk Road (solid grey lines) and the dispersal routes of peach (dashed grey lines, modified from Wang & Zhuang, 2001) and Hyalopterus arundiniformis (orange line). (c) The best-fitting model of the origin and migration scenarios of ABCbased inferences for Hyalopterus arundiniformis. The abbreviations represent different populations from Southwest China (SW), Southeast China (SO), Northwest China (WE), Northeast Asia (NE) and Central Asia and Europe (AE) the native range refers to Southwest China (SW; Figure S2), where HAR originated; in contrast, the nonnative range is the expanded distribution range (i.e. all ranges except for the south-western region).
To better interpret the results in a biological context, eight climate variables that have been widely selected in previous studies were used and divided into three categories: environmental energy, water availability and climate seasonality (Table S5). Using a principal component analysis (PCA), we translated eight climatic variables into a two-dimensional gridded environmental space for the climatic niche analyses using a principal component analysis (PCA). The occurrence densities for both native and nonnative ranges were then estimated using a kernel density function. Niche overlap between occurrence densities in both ranges was calculated using Schoener's D index, which yielded a value ranging from 0 (no overlap) to 1 (complete overlap) Warren et al., 2008). We also applied an approach from niche dynamics to disentangle changes in

| Archaeological evidence
We compiled archaeological and historical evidence for peach from previous literature to evaluate the consilience between the evolutionary dynamics of HAR and its primary host (Figure 1a; Table   S6). The oldest peach stones were found in Kunming, Yunnan, and have been dated to ca. 2.6 million years ago. Peach stones from the Neolithic era were mainly found in Southwest and Southeast China, while fossil data are lacking in Northwest China. Peach stones were also from the Early Jomon period (6700-6400 BP) in Japan. The earliest records of peach in Northwest China were from written records from the Bronze Age (ca. 3000 BP). The archaeological and historical records of peach from Central Asia and Europe did not appear until the first century BC.

| ABC inference of population history
The ABC analyses revealed that Models 3 and 12 were most supported in each step (Figures S2, S3). These models showed the highest posterior probabilities (0.9867 and 0.7032, respectively), and their 95% CIs did not overlap with those of other candidate scenarios (Table S7). In step 1, under the best-supported model, the results suggested that HAR originated in Southwest China and then spread to Southeast and Northwest China ( Figure S2). This model also showed that the NE region was formed by admixture between WE and Southeast China (SO). In the second step, among the four competing scenarios for the source of AE, a more recent divergence from WE was the most likely ( Figure S3). Overall, our results suggested that the dispersal routes of HAR were from Southwest China to Southeast and Northwest China and then further westward to Central Asia and Europe, which are roughly consistent with the dispersal routes of peach (Figure 1b,c). In addition, the dispersal routes of HAR were consistently supported after removing the individuals of the host-specialized lineage.

| Phylogeographic structure and gene flow
The BI and ML analyses based on the concatenated dataset revealed six deep lineage splits, five of which were exclusively presented in East Asia (AL1-AL6, Figure 2a formed the broadest range that spanned the Eurasian continent and included six sub-clades. The most basal subclade of Clade AL6 was distributed from south-western China to western Eurasia (AL6-1; Figure S4). However, the six sub-clades of Clade AL6 showed slight genetic divergence, and all of them spanned a large region and overlapped with each other, suggesting a lack of population structure ( Figure S4). Similarly, the mtDNA haplotype network also revealed six clusters corresponding to six clades (Figure 2c).
The inferred gene flow of HAR between the five regions is shown in Figure S5. SW and WE had relatively higher outward gene flow than in the inward direction, while SO had both substantial outward and inward gene flow. For AE and NE, they all received substantial gene flow from WE. Overall, there are two major paths of gene flow between different regions. One was derived from SW through SO and ended in NE, and the second was mainly from SW to WE and then from WE to AE and NE, respectively.

| Historical demographic dynamics
The BSP analyses were conducted for only three clades (AL4, AL5 and AL6) due to the small sample sizes for the remaining clades. The results revealed significant increases in Ne for all three clades but experienced distinct demographic histories for different lineages

| Climatic niche divergence
In total, 82.13% of the variance in the PCA-env analysis was explained by the first two principal component axes (PC1 = 68.86%, PC2 = 13.27%; Figure S6a). The climatic variables representing environmental energy and water availability were occupied in the first axis, whereas the 2nd axis was related to climate seasonality variables. The climatic niche comparison analysis showed that the niche overlap in both ranges was low (Schoener's D = 0.10). The niche similarity test revealed that the climatic niche was significantly divergent between the two ranges compared with random expectations (p = .27153 is < 0.05; Figure S6b). Concerning climatic niche dynamics, niche expansion explained the largest proportion (78%), and niche stability and unfilling accounted for 22% and 19%, respectively ( Figure S6c). The niche dynamics related to environmental energy and climate seasonality (represented by BIO6 and BIO4, respectively) between the two ranges for HAR are shown in Figure 4.

| Spatial origin and range expansion
Our results revealed that HAR originated in Southwest China, followed by eastward and north-westward dispersal (Figure 1b,c). The south-western mountainous and adjacent regions (i.e. Qinling and Daba Mountains) are characterized by complex topographic conditions and diverse environments (He et al., 2019;Lei et al., 2015), leading to extraordinarily high levels of endemism and species diversity (Lei et al., 2003). These regions also serve as an origin centre and glacial refuge for different species (Ding et al., 2020;Lopez-Pujol et al., 2011;Wang et al., 2020). Previous studies showed that the origin and differentiation of aphids are closely related to their host plants Liu et al., 2022). The oldest peach fossil was reported from the late Pliocene strata of Kunming, Yunnan (south-western China; Su et al., 2015). Previous genomic evidence also indicated that peach originated in Southwest China approximately 2.47 Ma (Yu et al., 2018). Similarly, the inferred origin of the pest Grapholita molesta, which also damages peach, was consistent with that of its original host plant (Wei et al., 2015). Therefore, it is likely that HAR originated in south-western China.
Compared with birds and mammals, the aphid HAR has a poor dispersal ability but is distributed across the entire Eurasian continent, including the Q inghai-Tibet Plateau, southern Asia and western Europe (Figure 2a). Additionally, the lack of phylogeographic F I G U R E 4 Shifts in the climatic niche along BIO6 and BIO4 between native and nonnative ranges of Hyalopterus arundiniformis.
Occurrence densities along the climatic variables are shown for the native (green), nonnative (pink) ranges and overlap between the two ranges (violet). The solid and dashed contour lines show the 100th and 75th quantiles, respectively, of the density in the available climate space. The horizontal bars show niche dynamics present along the x-axis: stability (violet) and expansion (pink). The solid arrows indicate the shift direction of the niche centroid between the two ranges. The shift direction of the average available climatic conditions between the two ranges is shown by the dashed arrows herbivores (Kebe et al., 2017;Silva et al., 2020). In addition, human transportation of host plants usually provides ideal conditions for herbivorous insects to spread over long distances (Gippet et al., 2019), which results in many insects exhibiting a broad geographic range, even across several continents (Kemp et al., 2020). The extensive export of peach to West Eurasia dates back to the time of Emperor Wu of Han (Wang & Zhuang, 2001). In the present study, the increase in the Ne for Clade 6 occurred within this period.

| Demographic dynamics in the Holocene
According to the BSP analysis, the population size expansion for HAR occurred in the Holocene after the Last Glacial Maximum (LGM, 21 kya). This result was inconsistent with those from previous studies of other organisms (Cheng et al., 2016;Wang et al., 2013;Zhao et al., 2012), in which East Asian lineages usually exhibited a "pre-LGM expansion" model (i.e. population size expanded before the LGM period). Holocene expansion has been observed in human populations and could be attributed to the expansion of Neolithic agriculture (Gignoux et al., 2011;Li, Ye, et al., 2019;Zheng et al., 2011). Similar Holocene expansion patterns in insects have been shown for the soybean looper moth Chrysodeixis includens (Silva et al., 2020), seed beetle Callosobruchus maculatus (Kebe et al., 2017), squash bee Peponapis pruinosa (Lopez-Uribe et al., 2016) and Homalictus fijiensis (Dorey et al., 2021),  Table S6; Faust & Timon, 1995;Wang & Zhuang, 2001).
Here, the BSP analysis indicated that the most pronounced Ne increase for clade AL6 was from the Late Neolithic to ca. 1.5 kya.
Two main factors could contribute to the Holocene expansion of HAR population size: (i) the domestication and trading of the primary host; and (ii) anthropogenic niche construction processes.
Domestication could weaken plant resistance against insect herbivores and enhance the performance of herbivores, such as their survivorship . Artificial selection in peach led to rapid shifts in morphology and genetics that caused domesticated peaches to differ significantly from their wild progenitors (Cao et al., 2014). For example, domesticated peach leaves were larger and thinner than those of ancient wild species (Wang & Zhuang, 2001). As a sap-sucking insect, HAR feeds on the undersides of the leaves of peach (Blackman & Eastop, 2020 for HAR ( Figure S6c) and that the realized climatic niche underwent significant shifts along with the climate seasonality variables ( Figure 4). Overall, the expansion in the Ne of HAR might have been closely associated with the domestication of its host plant and anthropogenic niche construction processes.

| Implications for pest management
Many invasion events are challenging to recognize due to morphological indistinguishability between invasive and native species, known as cryptic invasions (Morais & Reichard, 2018). HAR is highly morphologically similar to the other two species of Hyalopterus broadly distributed in western Eurasia  and had long been considered a special type of Hyalopterus amygdali (Lozier et al., 2007;Mosco et al., 1997). The taxonomic status and distribution of HAR were not identified until recently Lozier et al., 2008 invasion events tends to be unequal across space, which is tightly associated with economic activity, population density or the human footprint (Essl et al., 2011;Gippet et al., 2019). Given the ongoing globalization of trade and travel, a project involving rapid response to and risk assessment of alien species is required to avoid more humanmediated translocation events. Additionally, Bernardo-Madrid et al.
(2019) reported that introductions are homogenizing the Eurasian mammal zoogeographical regions. However, insects, at least in terrestrial ecosystems, are generally the most pervasive and damaging group of animal invaders (Bradshaw et al., 2016). Biotic homogenization of insect communities may be more common, although this topic has received little research attention. Thus, we also call for more studies on the impacts of human dispersal and trans-Eurasian trade during the Holocene on agricultural pests, which will provide a new perspective on pest management and biological invasions.

| CON CLUS IONS
Based on extensive sampling and a multidisciplinary approach, the eco-evolutionary dynamics of the herbivorous insect HAR, a spe- attention needs to be paid to the effects of human-induced evolution on insect pests, which will not only provide new insights into the dispersal and expansion of pest populations but also help prevent species invasions due to the ongoing globalization of trade and travel.

ACK N OWLED G EM ENTS
We thank all the sample collectors for their assistance and Ms. Fen-Di Yang for making aphid slides of the voucher specimens.

CO N FLI C T O F I NTE R E S T S
The authors declare that there are no conflicts of interest.

PE E R R E V I E W
The peer review history for this article is available at https://publo ns.com/publo n/10.1111/ddi.13541.

DATA AVA I L A B I L I T Y S TAT E M E N T
Newly generated sequence data for COI, tRNA/COII, Cytb, 12S and EF-1α sequences in this study were deposited in GenBank of NCBI at https://www.ncbi.nlm.nih.gov, and the accession numbers are provided in supporting information.