Multilocus evidence provides insight into the demographic history and asymmetrical gene flow between Ostrinia furnacalis and Ostrinia nubilalis (Lepidoptera: Crambidae) in the Yili area, Xinjiang, China

Abstract Tianshan Mountains provide a model for studying biological evolution and speciation. Here we assess the evolutionary history of Ostrinia furnacalis (ACB) and Ostrinia nubilalis (ECB), which are sympatric in the Yili River Valley in Xinjiang, China. Our study is based on the historical gene flow analyses of two species by using three mitochondrial DNA (mtDNA, COI, COII, Cytb) and four nuclear DNA (nuDNA, EF‐1α, Wingless, RPS5, CAD) markers obtained from representatives of HC (Huocheng), YN (Yining), XY (Xinyuan), and MNS (Manasi). Our results reveal that there is an asymmetrical gene flow pattern between the four populations. The population migratory pathways between these different populations show inflow into HC and YN, outflow from XY, and that MNS maintained a flow balance. Bayesian divergence time dating based on the COI gene suggests that the genetic divergence between the two species in this area may have occurred in Holocene at 0.008 Mya. Neutrality tests (Tajima's D, Fu's F s), and mismatch distribution test results suggest that population expansion events may not have occurred in the recent past. The demographic history and gene flow pattern between ACB and ECB may follow the “mountain isolation” hypothesis. The ML and BI trees of the mtDNA haplotype dataset show that ECB haplotypes are grouped together in a distinct clade and are clearly separate from ACB haplotypes. However, the geographical pattern of haplotype distribution is less clear for both ACB and ECB, supporting that there has been frequent gene flow among the geographic populations in the Tianshan Mountains. These findings indicate that the Tianshan Mountains are less likely a barrier to gene flow of the two species.


| INTRODUC TI ON
Deciphering genetic variations and demographic history of the existing distribution patterns of different populations within and between species have been a major focus in biogeography for decades (Avise, 1992;Yao et al., 2021). Several biogeographic hypotheses asserting that diversification and speciation were driven by isolation of different populations, e.g., "glacial refugia" hypothesis (Avise, 1992;Smith, 1965) and "mountain isolation" hypothesis (Avise, 1998).
These hypotheses assume that the different populations were partially or completely isolated by glaciers or mountains during the glacial period, and the genetic differences of isolated populations gradually accumulated with genetic drift and/or natural selection (Willis et al., 2004). Therefore, the complex historical trajectories linked to species divergence have posited to understand current population genetic processes and to predict the potential for populations to respond to selection and divergent environmental conditions. However, the study that focuses on testing biogeographic hypotheses is scarce in biodiversity hotspots in Central Asia, particularly in the northwestern of China (Yao et al., 2021).
The Tianshan Mountains of Xinjiang experienced complex orogenic and climatic events during the Pleistocene (0.0117-2.588 Mya) and the Holocene (present to 0.0117 Mya) are located in the northwestern China and provide an excellent model for inferring biogeographical processes associated with species divergence according to proposed biogeographic hypotheses (Sun et al., 2004).
The topography and landforms of the Tianshan Mountains became very complicated, ranging from 300 to 7500 m above sea level (Hewitt, 1996) and were heavily influenced by the uplift of the Qinghai-Tibet Plateau. A series of large and small mountain ranges are running north-south and form the North, Middle, and South Tianshan (e.g., Alatau Mountain, Boluokenu Mountain, Narathi Mountain) with a number of intermountain basins (e.g., Junggar Basin, Yili Basin, Tarim Basin) and the Yili River Valley. Many glacial refuges are found within these complex landforms (Médail & Diadema, 2009), which are considered as important glacial refuge hotspots. Some species may have experienced frequent distribution expansion and contraction during glacial-interglacial cycles, leading to potential gene flow between isolated populations in refugia on the Tianshan Mountains (Brower & Desalle, 1998;Qu et al., 2011Qu et al., , 2014Song et al., 2009). However, the Tianshan Mountains may act as a barrier to gene flow between geographically isolated populations and close relatives, resulting in spatial genetic differentiation and the current complex population structure of a sibling species complex. For example, the Qinghai-Tibet Plateau, the south of the Tianshan Mountains, and its surrounding regions contributed to the divergence between Primula nutans and Primula fasciculata (Ren et al., 2018).
The ACB (Asian corn borer), O. furnacalis (Guenée), and the ECB (European corn borer), O. nubilalis (Hübner) (Crambidae: Pyraustinae), are worldwide maize pests that cause substantial yield losses in corn production Frolov et al., 2007;Mutuura & Munroe, 1970). These two corn borers both belong to the O. nubilalis species group (trilobed uncus of male genitalia, which is a structure derived from the 10th abdominal tergite to grasp the female during copulation; see Yang et al., 2021: 830, Clade III in Figure 1), one of the most evolutionarily and ecologically interesting but taxonomically difficult groups in Lepidoptera. The species group includes 10 species and 23 subspecies worldwide (Frolov et al., 2007;Mutuura & Munroe, 1970;Yang et al., 2021). Incongruence between molecular phylogenetic relationships and the traditional classification of Ostrinia has been puzzling for a long time, leading to a number of members including ACB and ECB being morphologically indistinguishable and making accurate species identification extremely difficult (Hoshizaki et al., 2008;Kim et al., 1999;Mutuura & Munroe, 1970;Wang et al., 2017;Yang et al., 2021). On the other hand, this species complex is a good example to understand biogeographic patterns and evolutionary histories owing to their distributional complexity of the world. In these species, the ECB has spread from Europe to the central regions of Asia, as far east as Uzbekistan and to the western edge of Xinjiang (Xinjiang Uyghur Autonomous Region) in China, while the ACB has been reported to damage maize in eastern Asia, Northern Australia (Nafus & Schreiner, 1991), and in the major maize-growing regions of Eastern China (Hu & Sun, 1979;Tang et al., 1988).
Early studies indicated that only ECB occurred in the Yili area of the Western Xinjiang (Li et al., 1982;Huang et al., 2017;Tang et al., 1988).
Recent studies now indicate that ACB also occurs in this area (Li et al., 2010(Li et al., , 2013Yang et al., 2008Yang et al., , 2011. There are now some regions in the Yili Kazak Autonomous Prefecture Area of Xinjiang where ACB and ECB co-occur (Wang et al., 2017). As a result, it is the only known area in the world where these two corn borer species, ACB and ECB, have made contact. This provides the rare opportunity to study the evolution of sympatric sibling species in situ (Wang et al., 2017).
Previous studies on the genus Ostrinia showed that assortative mating is evident among sympatric ECB and Ostrinia scapulalis (the Adzuki bean borer) in Europe based on pheromone race-specific gene marker genetic analysis (Malausa et al., 2005). Coates et al. (2018) utilized a single-nucleotide polymorphism (SNP) at the pgfar locus to study the genomic mechanisms of sympatric ecological and sexual divergence of ECB, illuminating the genetic basis of strain-specific adaptive traits and gene flow between ECB strains. Wang et al. (2017) proposed that there was introgression of genes and the presence of hybrid individuals between ACB and ECB within the Yili area, suggesting that incomplete lineage sorting and historical gene flow may shape the evolutionary trajectory of these two species in this area. Specimens were collected using both sweeping net and light trap methods and have been stored as pinned dry collection at the Entomological Museum, Northwest A&F University, Yangling, Shaanxi Province, China (NWAFU). Some of the adult specimens were preserved in 95% ethanol at −20°C at the NWAFU. We selected 15-19 samples from each population for DNA extraction (Appendix S1).
Voucher specimens and corresponding mounted slides were deposited at the NWAFU. All newly generated sequences have been submitted to GenBank with the following accession numbers:

| Genetic polymorphism and population structure analyses
To evaluate the variation in genetic diversity among the four sampled geographical populations, we calculated the number of haplotypes (N h ), haplotype diversity (H d ), and nucleotide diversity (π) for each population of different species using DnaSP v 6.12.03 (Rozas & Librado, 2017) based on the datasets of the combined three mitochondrial genes (MTD) and the combined four nuclear genes (NUD), respectively. The mtDNA haplotype dataset (MTHD) and nuDNA haplotype dataset (NUHD) were defined separately using DnaSP v 6.12.03 (Rozas & Librado, 2017) based on the MTD and NUD. The defined MTHD and NUHD were used in the subsequent phylogenetic analyses. The haplotype networks of the two species were constructed using TCS (Templeton, Crandall, and Sing's parsimonious network, Templeton et al., 1992) network strategies in PopART (Bandelt et al., 1999;Clement et al., 2002) based on MTHD.
An analysis of molecular variance (AMOVA) was used to investigate the genetic variation between and within populations in the Arlequin 3.5 software (Excoffier & Lischer, 2010) with 1000 permutations based on MTD and NUD, respectively.
Genetic differentiation (F st ) of different geographical populations for each species was calculated in the DnaSP v 6.12.03 software (Rozas & Librado, 2017) based MTHD and NUHD, respectively.

| Phylogenetic analyses
A sequence saturation test was carried out for MTHD and NUHD, respectively, by using the DAMBE software (Xia & Lemey, 2009).
The phylogenetic relationships of different haplotypes were reconstructed with maximum likelihood (ML) and Bayesian inference (BI) based on MTHD and NUHD. Both ML and BI were used for each dataset. The species Ostrinia latipennis (Warren) (H27) was selected as the outgroup in the phylogenetic analyses. ML and BI analyses were executed in the PhyloSuite v 1.2.2 software (Zhang et al., 2020). ML analyses were conducted with 1000 bootstraps and were run 10 times starting from random seeds under the GTRGAMMA model. The best partition schemes and substitution models for each gene of BI analyses were estimated for the MTHD and NUHD using PartitionFinder v 2.1.1 (Lanfear et al., 2017) under

| Historical demography and divergence dating analyses
The historical demographic dynamics of different geographical populations for each species was estimated via the effective population size using multiple approaches based on MTD and NUD, respectively. First, Tajima's D (Tajima, 1989) and Fu's F s (Fu, 1997) statistics were used to assess whether the nucleotide polymorphisms deviated from expectations under the neutral theory with 10,000 coalescent simulations in the Arlequin 3.5 software (Excoffier & Lischer, 2010). in Tracer v 1.7.2 . The nucleotide substitution model of COI marker was selected through MrModeltest v 2.3 (Nylander, 2004) as the "HKY + I" model. Due to the lack of dated fossils or geological events useful for molecular clock dating, we used the widely accepted mutation rates for the insect mitochondrial COI marker (0.0177 site −1 Ma −1 ; Papadopoulou et al., 2010).
Samples were taken every 40,000 generations, and 400 million generations were performed, with the first 10% of samples discarded as burn-in. Tracer v 1.7.2  was used to monitor the stability of the chain and the convergence of the two parallel runs to determine whether the effective sample size (ESS) of each parameter reached the recommended value of >200.
TreeAnnotator    (Govindajuru, 1989). In order to investigate the gene flow migratory pathways, we used Bayesian stochastic search variable selection (BSSVS) analyses (Su et al., 2015) based on the MTD.

| Population structure and phylogenetic analyses
We identified 27 haplotypes based on MTHD (Figure 2c DAMBE analysis showed that haplotype sequences were far from saturated and could be used for subsequent analyses (p < .5, Iss (0.0301) ≪ (0.7851)). For the MTHD dataset, phylogenetic analyses inferred with BI and ML methods were slightly different (Figure 2a,b).
The topology inferred by the BI method showed that ECB haplotypes were grouped together in a distinct clade and clearly separated from the ACB clade, while the ML analysis showed that most ECB haplotypes (Hap10, Hap9, Hap5, Hap11) were grouped together except for the Hap23, which was nested into the ACB clade. In addition, haplotypes from the same locality were not lumped together supporting no geographical structuring (Figure 2). In contrast, phylogenetic analyses inferred with ML and BI methods based on NUHD showed that ACB and ECB were not recovered as monophyletic lineages (Appendix S7).
The discordant results between MTHD and NUHD suggested that incomplete lineage sorting might be present between these two species.

| Divergence time analyses
The BEAST analyses suggested that the first divergence time between the two species of ACB and ECB was estimated at 0.008 Mya

| Gene flow analyses
The N m was applied to identify the gene flow frequency within and between the two species across their different geographical populations (YN, XY, HC, MNS) based on MTD and NUD (Table 1)  of inflow migratory pathways. In contrast, XY had a large number of outflow migratory pathways, whereas MNS had almost equal inflows and outflows. In general, we found that population migratory route direction between different populations was outflow in XY and inflow in YN, HC ( Figure 5). In addition, the statistics of BF and Indicator based on MTD were summarized in

| Genetic differentiation and demographic history of ACB and ECB in sympatric regions within the Yili area
Our findings indicate that there is high genetic differentiation between the two species (ACB and ECB) in Xinjiang. However, the degree of genetic differentiation among different geographical populations within each species is low. These results are consistent with previous studies reported that a low level of genetic differentiation is present among Chinese populations of ACB Li et al., 2014). Furthermore, our molecular variance analysis (AMOVA) indicates that a large proportion of the total genetic variance is attributed to variations within populations (Appendix S5). Moreover, genetic distance and Mantel test results show that there is no significant correlation between geographical distance and genetic distance, suggesting that genetic variation among populations of corn borers may not be associated with geographical distance but most likely related to biological characteristics and host plants differentiations among populations proposed by Wang et al., 2018). By contrast, Coates et al. (2011) and Kim et al. (2011) suggested that significant genetic differentiation may present among intermountain regions due to topographic barriers in the Eastern United States based on microsatellite and SNP markers.
Our neutrality tests results (Tajima's D and Fu's F s , shown in Appendix S4 and Figure 3) indicate that recent population expansion has not occurred within each of four populations (Fu, 1997;Harpending et al., 1998;Tajima, 1989). In addition, the mismatch distribution among haplotypes is multimodal, further supporting there being no population expansion within these populations. However, limitation of

| Gene flow and migratory pathways between the different populations of ACB and ECB
According to our gene flow analysis (Table 1), our results show that a high level of gene flow exists between populations for both ACB and ECB, which corresponds to the conclusion that a low degree of genetic differentiation is present among populations within ACB and ECB due to their strong ability to spread (Li et al., 2010Yang et al., 2015). In addition, a low degree of gene flow was found between ACB and ECB within the sympatric Yili area in Xinjiang, which is consistent with the gene flow estimates based on concatenated COI-COII mitochondrial haplotypes between locations where ACB and ECB co-occur (Wang et al., 2017). However, this is inconsistent with previous studies on other species pairs of the genus Ostrinia Coates et al., 2013;Malausa et al., 2005).
According to high-throughput SNP and microsatellite markers and based on concatenated COI & COII, Wang et al. (2017) confirmed that the Yili area is a hybrid zone between ACB and ECB and found that there is gene flow from the invading ACB into ECB. Therefore, we hypothesize that reinforcement might be the major driver of reproduction isolation of ACB and ECB in sympatry. However, the molecular mechanisms of hybrids and the possible causes of the reinforcement and introgressive hybridization of ACB and ECB in sympatry are unclear yet, and further investigation based on genomic data, which provides greater inferential resolution for resolving finescale phylogeographic patterns, is needed.
According to migratory pathways between the different popula-  (Avise, 1992;Smith, 1965  this may lead to fragmentation of habitats and genetic differentiation among isolated geographical populations (Bansal et al., 2013;Schluter, 2009;Zheng et al., 2002). In general, our results preliminarily reveal that the evolution pattern of ACB and ECB in this sympatric region is likely consistent with the Holocene climatic oscillations and "mountain isolation" hypothesis (Hewitt, 1996;Pyron & Burbrink, 2010). However, our proposed hypotheses should be tested through multiple lines of methods such as deep phylogeographic and population genetic structure investigation on ACB and ECB in this area.

| CONCLUSION
Our results highlight that there is an asymmetrical gene flow pattern among four different geographical populations of ACB and ECB in the Tianshan Mountains where they are less likely a barrier to gene flow of the two species. The geological factors and climatic oscillations during the Holocene may have driven genetic patterns of two species. However, to fully understand the genetic variation and demographic history of corn borers in Xinjiang, and the gene flow between ACB and ECB in the sympatric Yili area, it will be necessary to expand the range of specimen collection and use additional genes in further investigation.

AUTH O R CO NTR I B UTI O N S
Bing Li: Software (lead); writing -original draft (lead); writing -review and editing (lead). Zhaofu Yang: Project administration (lead).

ACK N OWLED G M ENTS
We thank Akito Y. Kawahara (Florida Museum of Natural History,  -15).

CO N FLI C T O F I NTE R E S T
The authors declare no competing interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are openly available in