The genetic differentiation of a cricket (Velarifictorus micado) with two modes of life cycle in East Asia after the middle Pleistocene and the invasion origin of the United States of America

Abstract The cricket Velarifictorus micado is widely distributed in East Asia and colonized the United States of America (the USA) in 1959. It has two life cycles: egg and nymph diapause. We aimed to investigate the biogeographic boundary between them and determine when and why V. micado diverged. Mitochondrial fragments including COI and CytB were used for haplotype network, demographic analysis, and divergence time estimation in individuals of East Asia. We selected several samples from the USA to find out the colonization origin. The haplotype network indicated there were three lineages based on COI, NE lineage (the egg diapause and mainly distributed in the northern regions), SE lineage (the egg diapause and mainly distributed in the southern regions), and SN lineage (the nymph diapause and mainly distributed in the southern regions). The molecular chronograms indicated that the first divergence of V. micado into two main lineages, NE and southern lineages (SE and SN), was essentially bounded by the Yangtze River. It occurred around ~0.79 Ma (95% HPD: 1.13–0.46 Ma) in the Middle Pleistocene Transition. This was followed by the divergence of the southern lineage into two sublineages, SE and SN lineage, occurred around ~0.50 Ma (95% HPD: 0.71–0.25 Ma), corresponding to the time of development of glaciers in various parts of the Qinghai–Tibet Plateau (QTP) (0.73–0.46 Ma). SE lineage might originate from southwestern China based on the comparison between the haplotype network based on COI and CytB. Our study suggested that divergences of lineages have twice co‐occurred with tendency of cooling climatic in Asia after the Mid‐Pleistocene, and the life‐history strategy may play an important role in lineage diversification. Additionally, our results indicated that the USA populations were revealed at least twice separate Asian invasions. These both belonged to the egg diapause, which might provide a new perspective for invasion control.


| INTRODUC TI ON
The collision and continuous compression between the Indian plate and the Asian continent, a major event in the Phanerozoic eon, resulted in the emergence and uplift of the Himalayas and the Qinghai-Tibet Plateau (QTP) (59-50 Ma) (Hall, 2002;Hu et al., 2017;Murphy et al., 1997;Searle, 1987;Wang et al., 2014;Zahirovic et al., 2012).
Geographic changes resulted in the formation and development of Asian water system, such as the modern Yangtze River established during 23 Ma (Zheng et al., 2013), and driving the climate changes in Asia. The Yangtze River, a natural barrier, isolated the species distributed in northern and southern China. The climatic oscillations occurred frequently during the Pleistocene. The cold glacial and warm interglacial cycled regularly and repeatedly. Glaciers developed gradually in various parts of the QTP after the Kunhuang Movement (0.9 Ma) (Wu et al., 1999;Yao et al., 2000;Zheng, 2000;Zhou et al., 2002). Glacial cycles had an important impact on the distribution and genetic structure of a species (Hewitt, 2000(Hewitt, , 2004. Climate oscillation occurred in Asia during Pleistocene has affected the current distribution of species (Jia & Zhang, 2019;Ye et al., 2014;Qu et al., 2011), leading to genetic differentiation in species (Gu et al., 2013;Hazlitt et al., 2014;Rodrigues et al., 2014;Xu et al., 2009;Ye et al., 2018) and speciation in Asia (Bao et al., 2010;DeChaine et al., 2014;Lin et al., 2016;Shih et al., 2011).
Individuals colonized to new habitat due to the environmental changes in original habitat or population expansion. New behavioral, physiological, and genetic changes occurred in the adaption to new habitat (Avise et al., 1998). Studies suggested that life history is related to environmental factors (Marten et al., 2006;Papadopoulou et al., 2009). For example, the number of generations each year of Locusta migratoria (Orthoptera: Acrididae) increased with the decrease of latitudes in China (Chen, 1999;Guo et al., 1991;Tanaka & Zhu, 2008). The different life history within species played a more important role in speciation than geographic isolation, as in Teleogryllus commodus distributed in Southeastern and Northeastern Australia (Hogan, 1967), Pternemobius taprobanensis and Pternemobius fascipes distributed from the tropical to temperate region in East Asia (Masaki, 1983). Genetic data also indicated that individuals with different life cycles had been differentiated (Harrison & Bogdanowicz, 1995;Snyder et al., 2009).
Velarifictorus micado (Saussure, 1877) is widely distributed in East Asia, including China, Russia, Japan, Korea, Cambodia, Vietnam, Indonesia, and nearby islands, covering both the Palearctic and Oriental realms. The distribution range of this species has gradually expanded since its introduction to the United States in 1959 (Alexander & Walker, 1962;Bowles, 2018). It was also found that this cricket has two modes of life cycle. In a population that diapause as eggs, adults sing from August to October, and their eggs hatch after getting through the winter. In the other population that diapause as nymphs, adults sing in May-July hatch quickly after mating and get through the winter as nymphs. Different adaptation types of V. micado lead to differentiation of reproductive modes and might result in speciation among populations due to different life cycle patterns (He & Takeda, 2013). However, it is unknown that Accordingly, our research will use the cricket that widely distributed in temperate and tropical as the material. Aimed to find out suitable gene markers which can distinguish the two groups and the biogeographic boundary between them, investigate when and what drive V. micado diversification and the invasion origin of the USA, mitochondrial fragments COI and CytB were extracted from the large sampling specimens and used for phylogenetic analysis, time estimation, and demographic analysis. (ii) Evolutionary difference between two genes to speculate the origin of the ancient residents. A total of 346 individuals from 70 V. micado populations were collected, 55 populations from China, 8 from Japan, 1 from Korea, 1 from Vietnam, 1 from Cambodia, and 4 from the USA (Figure 1 and Figure 2). All materials were presented in 100% ethanol and stored in a freezer at −20°C, and genomic DNA was extracted from legs of the cricket using AxyPrep™ Multisource Genomic DNA Miniprep kit.

| PCR amplification and sequencing
Universal primers for cytochrome c oxidase unit 1 (COI) and cytochrome b (CytB) genes were designed in the previous study (Table 1). The PCR procedure for the two genes included an initial denaturation at 94°C for 4 min, followed by 35 cycles of 30 s at 94°C, 30 s at 45°C and 30 s at 72°C, ending with a final extension at 72°C for 5 min. Sequencing was performed on 3730xl DNA Analyzer, and sequencing was proofread in ATGC Ver 7.0.2 (Genetyx Corporation) and aligned in MEGA V. 7.0.14 (Kumar et al., 2016).

| Genetic polymorphism
Base substitution saturation test was performed in DAMBE (Xia & Xie, 2001) for phylogenetic analysis. Different haplotype, haplotype diversity (Hd), and nucleotide diversity (π) were calculated in DNASP 5.10.01 (Librado & Rozas, 2009). Defining groups of populations that were geographically homogeneous and maximally differentiated from each other was implemented in SAMOVA 1.0 (Dupanloup et al., 2002). Processes for K from 2 to 10 were running using 100 simulated annealing approach. A median-joining (MJ) haplotype network to test the relationships between groups was constructed in POPART (Leigh & Bryant, 2015), and each group has its own color.
Genetic differentiation (F ST ) among different lineages and the analysis of molecular variance (AMOVA) were calculated in ARLEQUIN 3.5.2.2 (Excoffier & Lischer, 2010). F ST was performed with pairwise difference genetic distance. AMOVA was used for test genetic structure and genetic variation analysis. Statistical significance of the estimates was assessed by 10,000 permutations.

| Historical demographic changes
Neutrality tests were used to speculate historical demographic change and implemented in ARLEQUIN 3.5.2 (Excoffier & Lischer, 2010), including Tajima' D (Tajima, 1989) and Fu' Fs (Fu, 1997) with 1,000 simulated samples. These statistics assume that the population has a constant size at mutation-drift equilibrium. A significantly negative value indicates there are excess of low frequency haplotypes in the population, suggesting that there might be expansion. Moreover, the value of Tajima' D is negatively correlated with the time since introduction. Namely, the higher (positive) the value of Tajima' D, the more recent introduction, indicating that these populations are recovering from the bottleneck effect.

| Divergence time estimation
Estimating divergence time was based on mitochondrial markers of COI in BEAST 2.5.0 (Bouckaert et al., 2014). Best substitution model was evaluated using Bayesian information criterion (BIC) selected in the program JModeltest v.2 (Darriba et al., 2012). A set of 15 specimens comprising 6 other Velarifictorus species and 9 ancient haplotypes of V. micado was employed for the estimation of divergence time (Table 2). A relaxed clock log normal was applied with 1.7% per site per lineage per million years for COI (Allegrucci et al., 2011(Allegrucci et al., , 2017Chobanov et al., 2016;Kiyoshi & Sota, 2006;Papadopoulou et al., 2010;Pons & Vogler, 2005;Shapiro et al., 2006). Divergence time was estimated using a Yule model.   & Walker, 1962). All populations in USA were integrated into the group USA. For the native populations, the value of F CT increased continuously from k = 2 to k = 10, and a distinct increase from 2 to 3 and from 5 to 6 using SAMOVA (Figure 3).

| Genetic polymorphism and haplotype network
The group comprising only one population when k = 6, thus, k = 3 was chose to defined the groups of populations. The three groups corresponded chiefly to three geographic regions. The group 1 contained populations mainly in northern China and Japan. The group 2 contained populations mainly in central China, Korea, and Japan. The group 3 was consisted of populations from southern China, Vietnam (VN), and Cambodia (KH) (Figure 2 and Table 3). The partitioning of total genetic variation in three groups using AMOVA indicated 71.45% diversity among groups, 21.76% within populations, and 6.79% among populations within groups (Table 4). Thirty-six unique haplotypes were derived from all 346 individuals. The distribution of different haplotypes based on COI was showed in Figure 4. The haplotype distributions were divided three starry shapes in the network, NE, SE, and SN, respectively. NE contained the individuals that produced with egg diapause and mainly distributed in the northern regions. SE included the populations with egg diapause and mainly distributed in the southern regions. SN was consisted with individuals that produced with nymph diapause and mainly distributed in the southern regions ( Figure 2 and Table 3). There was no sharing haplotype among the three lineages. Hap24, Hap 2, Hap 4, and Hap 11 haplotypes were the most frequent haplotypes, characterizing 25.72%, 25.14%, 17.92%, and 4.05% individuals, respectively. Hap 2 and Hap 11 were the ancestral haplotypes of NE, and Hap11 was disjoint from Hap 2, which suggested differentiation due to the long F I G U R E 2 Sampling localities of Velarifictorus micado. Black, red, and blue circles represent the NE (the individuals that produced with egg diapause and mainly distributed in the northern regions), SE (the individuals with egg diapause and mainly distributed in the southern regions), and SN (individuals that produced with nymph diapause and mainly distributed in the southern regions) lineages based on COI, respectively. Circle sizes are proportional to the number of individuals. The colorful regions represent the groups (1, 2, and 3) defined by SAMOVA   Table 3 and Table 9). The AMOVA results showed that the source of variation among groups categorized by three clades based on COI accounts for 81.25% (Table 7).

| Historical demographic changes
Based on COI, the negative and significant values of Tajima' D and Fu' Fs indicated past population expansion of SE and SN lineages.
The negative and significant value of D also implied the expansion in the past and the stable population size of NE lineage (Table 8).
However   Note: Genetic information included lineages, neutrality tests, and nucleotide polymorphism. Genetic information was based on COI. NE lineage was consisted of the individuals that produced with egg diapause and mainly distributed in the northern regions, SE lineage was consisted of the individuals with egg diapause and mainly distributed in the southern regions, and SN lineage was consisted of individuals that produced with nymph diapause and mainly distributed in the southern regions. Hap, haplotype distribution; Hd, haplotype diversity; π, nucleotide diversity; The bold black font represents the groups (1, 2, 3, and USA) divided by SAMOVA. *p < 0.05. **p < 0.01.

TA B L E 3 (Continued)
expansions recently causing the sympatric coexistence of lineages.
Besides the values of Fu' Fs were positive in ZJDP, ZJTM, SDDY, and ZJSX, there were at least two lineages in these populations.

| Divergence time estimation
The substitution models were selected under JModeltest v.2

| Life cycle and distribution
Velarifictorus micado has two ways to get through the winter, nymph diapause, and egg diapause (He & Takeda, 2013). Populations with egg diapause developed later than that with nymph diapause.
Adults of egg diapause population were found from July to October (Shandong, Heilongjiang, Henan, Guizhou and Yunnan provinces), and nymphs were from July to August (Shandong and Yunnan provinces). Adults with nymph diapause appear earlier, and nymphs could also be found later, because the next generation hatched before winter. The individuals with nymph diapause, as we found, were adults in March (Hainan province), April (Guangxi province), June (Hubei and Zhejiang provinces), and July (Guizhou province). The nymph was found in the spring or the later autumn, such as March (Hainan province), April (Guangxi province), August (Guizhou province), September, and October (Zhejiang province) ( Table 3). winter (Alexander, 1968). V. micado has two modes of (2) and (3).
The different life history reduced hybridization between diapausing and nondiapausing populations (Hogan, 1967), and one species with two life cycles, rather than geographic isolation, might result in speciation at certain latitudes (Hogan, 1967;Masaki, 1983 (Hogan, 1966). Walker proposed two hypotheses. The first one was consistent with Hogan, and the second was that the residents transferred to south and turned to produce multiple generations each year with the warmer climate (Walker, 1980). Our results suggested the northern egg-diapausing, the southern egg-diapausing, and the southern nymph-diapausing comprised a strongly supported monophyletic group, respectively. This was similar with previous studies that different life cycles of crickets comprised independent groups (Harrison & Bogdanowicz, 1995;Snyder et al., 2009). According to the biogeographic feature, the egg-diapausing distributed in most regions of East Asia ranging from 24° to 46°N. The nymph-diapausing distributed essentially in the south of Yangtze River in China, the northern limit reached to 32°N ( Figure 2). As other researches showed, with the decrease of latitudes, the number of generations each year of L. migratoria (Orthoptera: Acrididae) increased in China (Chen, 1999;Guo et al., 1991;Tanaka & Zhu, 2008). With the increase of latitude, the diapause mode of V. micado changed from nymph to egg diapause.
Our results suggested that there were both ancient and undifferentiated egg and nymph diapause haplotype in Yunnan, Guizhou, and Zhejiang provinces. However, the haplotype in Zhejiang was single.
Thus, the egg diapause in southern regions might originated from the southwestern China (Yunnan and Guizhou) ( Figure 5, Table 3 and  (Snyder et al., 2009). This phenomenon with two life histories bounded by the Yangtze River could be seen between Teleogryllus emma and Teleogryllus occipitalis (He et al., 2017;Lu et al., 2018). Wing dimorphism has been found in V. micado. The macropterous are capable of flying and the micropterous are flightless . However, there were only four macropterous individuals in our samples. It is possible that the migration ability of this cricket is weak. Obviously, there were at least two lineages in Jiangsu, Anhui, Shanghai, Zhejiang, and Yunnan provinces, so these are contact zones, which probably exists in other regions where has not be discovered due to the limited sampling time and size.

The wide sampling in East
According to the distribution of species, there is a broad faunal transition zone in the Quaternary between the Palearctic and Oriental realm in China (Zhang, 2002) and no strict biogeographic division (Norton et al., 2011). Although it is not strictly defined by the Yangtze River, our results were basically consistent with an suggestion of Wallace that the Yangtze River is the biogeographic boundary between the Oriental and Palearctic (Wallace, 1876). The sympatric coexistence of V. micado  Table 2 and Table 5). The divergence between them is maintained in the absence of obvious environmental difference and barriers to gene flow, such as the same photoperiod, temperature, humidity,, and habitats, while the genomic underpinning of ecological speciation often appear to have been found to be the result of a long period of allopatry (Bernatchez & Dodson, 1990;Feder et al., 2003;Foote & Morin, 2015;Gray et al., 2006;Kuehne et al., 2007;Lucek et al., 2018).

| The invasive origin of the USA
Our results confirmed that V. micado was transferred to America ( Figure 2 and Figure 4) (Bowles, 2018). V. micado was firstly found in the USA in 1959 (Alexander & Walker, 1962), and then, they were widely distributed around eastern and southeastern America (Peck et al., 1992;Walker, 1977). The recent research indicated the cricket had dispersed both northwards and westwards (Bowles, 2018). Our results revealed that there were at least twice invasions from Asia, one from eastern China (Hap 7) and another from Japan (Hap 5). Hap 8 might be introduced from Japan. Although it have not found in Japan due to the small sampling size, Hap 8 was related to Hap 5. They both got through winter as egg (Figure 4 and Table 3). Their distribution range in the United States may continue to expand due to their adaptation to dry and cold areas, for example, NE lineage was widely distributed in the northeastern China. The similar latitude and climate made the species origin from the USA easy to settle down to China, such as Corythucha ciliata, Homalodisca coagulata, and Anopheles quadrimaculatus (Li et al., 2007;Wang et al., 2008;Zhang et al., 2008).
Similarly, V. micado was easy to adapt to the environment of the USA.
Limited to sampling size, it is unknown if there are nymph populations in the USA (these may disperse southern). Our results suggested that the mitochondrial fragment COI could be easily used to check the origin of V. micado for prevention of invasion. Although it has not been found that this cricket has destructive damage to human or nature, the impact to ecosystem area of the new colony is unknown.

| Divergence time and major events
Haplotype network based on COI well supported that V. micado has di-  (Clark et al., 2006;Pisias & Moore, 1981;Raymo et al., 1997;Ruddiman et al., 1989), and the Kunhuang (Kunlun Huanghe) movement of 0.9 Ma BP (Li et al., 2001(Li et al., , 2004. During these global glaciers after the Kunhuang movement, glaciers developed gradually in various parts of the QTP (Wu et al., 1999;Yao et al., 2000;Zheng, 2000;Zhou et al., 2002). During MPT, the intensified aridification and the cooling climate have major impact on restructuring of vegetation and faunas in temperate northern East Asia. The steppe replaced the forest, some mammalian fauna extincted and new species produced due to the climate and related vegetation changes since MPT (Zhou et al., 2018). Subtropical zone reached at 42°N before the middle Pleistocene and retreated southwards after the mid-Pleistocene (Liu & Ding, 1983). Many species immigrated southwards, including boreal, thermophilic, and humid-preferring fauna. Species adapted to the new environment, while some remained and scattered in the locality during the Quaternary glaciation (Zhang, 2004). During the MPT, mid-Pleistocene Homo (MPH) (Bae, 2010;Wu & Poirier, 1995) (Wu et al., 1999).
The glaciation time in Yunshanping in southeastern QTP was 0.47-0.71 Ma (Yao et al., 2000;Zheng, 2000). The earliest glaciation time was dated to 0.46 Ma in Qilian Mountain in northeastern QTP (Zhou et al., 2002). Thus, glaciation and climatic shifts might trigger individuals in southern region diverging into two sublineages (SE and SN).
Our results suggested that there were both undifferentiated and ancient haplotype with egg and nymph diapause in Yunnan, Guizhou, and Zhejiang provinces. However, the haplotype in Zhejiang was single. Thus, the SE lineage might originate from the southwestern China (Yunnan and Guizhou) ( Figure 5, Table 3 and Table 9).
Our results indicated that lineage diversification of this cricket was driven by the climatic occasions after Mid-Pleistocene, and the different life cycles played an important role in that. But it is unknown whether  (Marten et al., 2006;Papadopoulou et al., 2009), and influence the gene flow, and thus lineage diversification (Fluker et al., 2014;Tibbets & Dowling, 1996;Whiteley et al., 2004). We had two hypotheses about the evolution of life cycle. If the evolution time of life cycle was much earlier than that of lineages, the egg-diapausing might be the ancient type, which was also supported by the Bayes tree ( Figure 6). However, if the evolution time of life cycle was similar to that of lineages, it is possible that the nymph diapause turned to egg diapause twice due to cooler climate after Mid-Pleistocene and diverged to three lineages. Namely, the nymph-diapausing might be the ancient type. Individuals in northern regions turned to egg diapause resulted from the climatic variability, and the individuals with nymph diapause were filtered out under natural selection in northern regions.
It should be noted that there are two questions remain to be

CO N FLI C T O F I NTE R E S T
The authors declare no conflict of interest. Zhu-Qing He https://orcid.org/0000-0003-4304-767X