Phylogeography and Genetic Structure of the Bush Cricket Decma fissa (Orthoptera, Tettigoniidae) in Southern China

© 2023 Academia Sinica, Taiwan Open Access Phylogeography and Genetic Structure of the Bush Cricket Decma fissa (Orthoptera, Tettigoniidae) in Southern China Qi Guo1, Qi-Di Zhu1, Zhi-Jun Zhou1, and Fu-Ming Shi1,* 1Key Laboratory of Zoological Systematics and Application of Hebei Province, College of Life Sciences, Hebei University, Baoding, Hebei 071002, China. *Correspondence: E-mail: [email protected] (Shi) E-mail: [email protected] (Guo); [email protected] (Zhu); [email protected] or [email protected] (Zhou) Received 29 November 2022 / Accepted 28 April 2023 / Published 12 July 2023 Communicated by Jen-Pan Huang Decma fissa is the most widely distributed species of the genus Decma occuring in southern China. This study presents the first phylogeographic work of D. fissa based on COI, Cytb and ITS sequence. We examined genetic diversity with ITS and mitochondrial sequence respectively, and phylogenetic work was based on the mitochondrial data. A high-level genetic diversity was revealed based on mitochondrial data but a low-level diversity was shown with ITS sequence. For the mitochondrial data, divergence time analysis displayed five lineages. Based on the Mantel test, geographic and genetic distances among D. fissa populations revealed a significant positive correlation. Bayesian skyline plot (BSP) analyses implied that none of three major lineages of D. fissa was seemingly affected by the last glacial maximum (LGM, 0.015–0.025 Mya). Ecological niche modeling was used to predict the distribution of D. fissa in four periods (LGM, Mid-Holocene, current and 2070) in China. Analysis of the ancestral area reconstruction indicated that D. fissa occurred in the South China area. Key words: Decma fissa, Phylogeography, Ecological niche models, Population structure, Demographic history Citation: Guo Q, Zhu QD, Zhou ZJ, Shi FM. 2023. Phylogeography and genetic structure of the bush cricket Decma fissa (Orthoptera, Tettigoniidae) in southern China. Zool Stud 62:32. doi:10.6620/ZS.2023.62-32. BACKGROUND The distribution and dispersing of insects are affected by many factors such as mutation, geographic isolation, natural selection, climate change and gene flow (Hewitt 2004; Garrido et al. 2012; Saeb and Al-Naqeb 2016; Tang et al. 2022). Understanding the historical processes of species can help us reveal their ability to confront environmental changes (Porretta et al. 2007; Lyons et al. 2012; Wei et al. 2013 2015). In particular, the repeated changes of the Quaternary climate (Clark et al. 2009) had a profound impact on the distribution pattern of existing species. Species had been reduced to refugia during the ice period, and then expanded again during the warm climate period (Hewitt 2000; Song et al. 2016). In addition, different species have different historical distributions and responses to environmental changes, so they would be affected differently (Liu et al. 2018). Unlike the climate of Europe and North America, due to the uplift of the Qinghai Tibet Plateau, the climate was warmer in eastern Asia in the LGM (Dai et al. 2011; Wei et al. 2013). Such a climate led to the absence of ice cover in low altitude areas. Although phylogeographic studies were carried out in many vertebrates (Dufresnes et al. 2016; Qiu et al. 2016; Liu et al. 2019), the distribution and genetic diversity of insects have been only partially studied in China (Ye et al. 2018; Zhou et al. 2021; Li et al. 2022; Wang et al. 2022; Tang et al. 2022). These studies show that insect species in eastern Asia might often survive in more than one refugia. Decma fissa (Orthoptera: Tettigoniidae: Meconematinae) is a long-wing species of Zoological Studies 62:32 (2023) doi:10.6620/ZS.2023.62-32 1 © 2023 Academia Sinica, Taiwan Meconematinae, which is a subfamily containing a high number of species (Cigliano et al. 2022). It was described as Xiphidiopsis fissa beloning to the genus Xiphidiopsis (Xia and Liu 1993). The genus Decma was established by Gorochov (1993) for Decma (Decma) stshelkanovtzevi. Gorochov et al (2005) moved Xiphidiopsis fissa into the genus Decma, and the new combination, Decma fissa, continues to be used today. Compared to the narrow distributions of other species of this genus (Liu and Yin 2004; Gorochov et al. 2005; Liu and Zhou 2007; Shi et al. 2013), the distribution of D. fissa ranges across several provinces in southern China (Xia and Liu 1993; Wei et al. 2016). In the present study, we predicted the potential distribution based on climate data and ecological modelling methods, and explored the dispersal pathway by ancestral area reconstruction. Our aims are to (a) determine the genetic structure and diversity; (b) study the demographic history among populations; and (c) identify the potential locations. MATERIALS AND METHODS Specimen Collection A total of 232 specimens were collected from 15 locations (Fig. 1). The location details are shown in table S1. All the individuals have been identified by morphological method and preserved in absolute ethanol and stored at -20℃ until DNA extraction. The key that we used for species identification is from Wang (2020). DNA extraction, amplification and sequencing Genomic DNA was extracted using TIANamp Genomic DNA Kit (Tiangen, Beijing) following the manufacturer’s protocol from muscle tissue of each adult individual. PCR was performed with Premix Taq™ (Takara, Beijing). A fragment of COI was amplified and sequenced using the primers COBU/ COBL (Huang et al. 2013), the primers Cytb-N/Cytb-J (Simon et al. 2006) were used for a fragment of Cytb, N Fig. 1. The geographic distribution collection points. page 2 of 15Zoological Studies 62:32 (2023) © 2023 Academia Sinica, Taiwan and the primers ITS-F/ITS-R (Weeker et al. 2001) were used for a fragment of ITS. PCR cycle profiles are: for COI, an initial denaturing of 3 min at 94℃ followed by 35 cycles of: 30s at 94℃, 30s at 49℃ and 90s at 72℃, then a step of 8 min at 72℃, finally forever at 4℃, for CYTB, an initial denaturing of 5 min at 95℃ followed by 35 cycles of: 30s at 95℃, 30s at 50℃ and 1 min at 72℃, then a step of 8 min at 72℃, finally forever at 4℃. For ITS, an initial denaturing of 3 min at 94℃ followed by 35 cycles of: 30s at 94℃, 30s at 46℃ and 90s min at 72℃, then a step of 8 min at 72℃, finally forever at 4℃. The information about primers is given in table 1. The PCR products were sent to Azenta (Tianjin, China) for sequencing in both directions after being analyzed with 1% agarose gels with ethidium bromide following electrophoresis. All sequences were edited with Seqman (DNASTAR, Lasergene 7.1) and aligned with BioEdit (Hall 1999). The combined data was concatenated with PhyloSuite v1.2.2 (Zhang et al. 2020). At last, a total of 232 COI, 232 Cytb and 229 ITS sequences were used in this work. Ecological niche models MaxEnt 3.3.3 (Phillips et al. 2006; Phillips and Dudík 2008; Warren et al. 2013) was used to predict the potential distribution with 23 distribution records including 15 collection points and 8 points from literature (Xia and Liu 1993; Liu and Yin 2004; Gorochov et al. 2005; Liu and Zhou 2007; Shi et al. 2013; Wei et al. 2016). All the distribution data from literature are provided in table S2. The ecological factors of the Last Glacial Maximum (LGM, about 22000 years ago), Mid-Holocene (about 6000 years ago), present (1960) and the future (2070) were downloaded from WorldClim (http://www.world clim. org) at 30s arc-min resolution and RCP4.5 (Thomson et al. 2011; Zhou et al. 2014). Firstly, we extracted the data of ecological factors with Arcgis 10.7 (Esri, Redlands, CA, USA) and counted the correlation between factors with IBM SPSS Statistics 20.0 (IBM Corp, Armonk, NY, USA). Then, we used ecological niche modeling (ENM) for current climate conditions using all factors. The model was performed with 10 replicate runs, 10000 maximum iterations and 25% for model training. The importance of factors was determined by percentage contribution in the result. The factors with high-relation (r > 0.8) and low percentage contribution were removed (Yang et al. 2013). Finally, nine factors were used to predict the distribution area of D. fissa in China: Bio2 (Mean Diurnal Range), Bio3 (Isothermality), Bio6 (Min Temperature of Coldest Month), Bio7 (Temperature Annual Range), Bio10 (Mean Temperature of Warmest Quarter), Bio12 (Annual Precipitation), Bio16 (Precipitation of Wettest Quarter), Bio17 (Precipitation of Driest Quarter), Bio18 (Precipitation of Warmest Quarter). We used the ENMeval package to optimize the regularization multiplier and feature class parameters in the R version 4.2.2 software (Muscarella et al. 2014; Kass et al. 2021). Based on the AICc calculated in the ENMeval, the best model setting for MaxEnt was the feature class (FC): linear (L), quadratic (Q), hinge (H), and regularization multiplier (RM) equal to 2. The accuracy of each model prediction was determined based on AUC scores. The score ranges from 0.5 to 1, and indicates excellent power of predicting when it is above 0.9 (Ye et al. 2014). Genetic diversity The number of haplotypes (Hap), haplotype diversity (Hd), nucleotide diversity (π), and variable sites were analyzed using DnaSP 5.10 (Librado and Rozas 2009). Population differentiation (FST) was implemented in Arlequin 3.5 (Excoffier and Lischer 2010) conducted with 10,000 permutations. FST value ranges 0 to 1, and the 0 implies complete panmixia (Hudson et al. 1992; Hudson and Richard 2000; Chabot et al. 2015). Mitochondrial haplotype network and divergence times Haplotype data of combined sequences were edited with DnaSP v5.10. The haplotype network was Table 1. Information about primers used in this study Gene Primers PCR annealing temperature COI COBU F: TYTCAACAAAYCAYAARGATATTGG 49℃ COBL R: TAAACTTCWGGRTGWCCAAARAATCA CYTB CYTB-N F: TTCTACTGGTCGRGCTCCAATTCA 50℃ CYTB-J R: GTTTTACCATGAGGTCAAATATC ITS ITS-F ITS-R TAGAGGAAGTAAAAGTCG GCTTAAATTCAGCGG 49℃ page 3 of 15Zoological Studies 62:32 (2023) © 2023 Academia Sinica, Taiwan constructed in PopArt 1.7 with Median Joining (Leigh and Bryant 2015). Divergence time was estimated with Beast v1.10.4 (Suchard et al. 2018). Due to the lack of fossil record, a uniform calibration prior, which is based on prior estimated divergence time, was set for the species Xizicus fascipes (GenBank code: NC018765.1) and Decma fissa at 42.84 Ma (Zhou et al. 2017). The analysis was run 100 million generations with TN + F + G4 model, uncorrelated relaxed clock, lognormal distribution, random starting tree, the Yule process, and sampling every 10,000 steps. The first 10% of the trees were discarded as burn-in. Tracer v1.7.1 (Rambaut et al. 2018) was used to assess convergence (all ESS parameters were > 200). The maximum clade credibility tree was produced in TreeAnnotator v1.10.4 and visualized with Figtree v4.1.3 (Rambaut 2012). Demographic history and population structure of mitochondrial data To trace the demographic changes of each population and three major lineages of D. fissa, neutrality tests and mismatch distribution analyses were conducted in Arlequin 3.5 (Excoffier and Lischer 2010). Two neutrality tests, Fu’s FS and Tajima’s D were used to evaluate historical demographic expansion with 1,000 simulated samples. Negative values of these two tests indicate that the population had undergone expansion, while positive values indicate the population experienced a bottleneck (Tajima 1989; Fu 1997). Mismatch distribution for three major lineages were conducted for the sum of squared deviations (SSD) and Harpending’s raggedness index (HRI) (Rogers and Harpending 1992; Harpending 1994; Excoffier 2004) by 1000 bootstrap replications. Genetic distances between populations were calculated with Mega X (Kumar et al. 2018) using Kimura’s 2-parameter and Tamura 3-parameter models. AMOVA (analysis of molecular variance) was implemented in Arlequin 3.5 conducted with 10,000 permutations. The spatial analysis of molecular variance (SAMOVA) was implemented by SAMOVA 2.0 (Dupanloup et al. 2002) with K = 2 to 15, and 10,000 iterations. The final K was determined by FCT (Heuertz et al. 2004). The population size of three major lineages over time was assessed using Bayesian skyline plots (BSP) analysis with concatenated mitochondrial genes (Drummond et al. 2005; Minin et al. 2008). This analysis was performed using BEAST v1.10.4 (Suchard et al. 2018). We focused on three major lineages (I, IV, and V) because other lineages (II and III) contained too few haplotypes. TN + F + G4 for lineage I and lineage V, HKY + F + G4 for lineage IV were determined to be the best substitution model with PhyloSuite (Zhang et al. 2020). Analyses were run using a strict molecular clock, assuming a divergence rate of 3.54% per million years for COI (Clock.rate parameter was set as a normal prior with an initial value = 0.0177, mean = 0.0177, and SD = 0.004) (Papadopoulou et al. 2010), and 4.22% per million years for Cytb (Clock.rate parameter was set as a normal prior with an initial value = 0.0211, mean = 0.0211, and SD = 0.004) (Pons and Vogler 2005). Analyses were run for 100 million generations and sampling every 10000 steps for lineages IV and V, and 200 million generations and sampling every 20000 steps for lineage I. Tracer v1.7.1 (Rambaut et al. 2018) was used to assess convergence (all ESS parameters were > 200) and visualization of median and 95% highest posterior probability density intervals (95% HPD). The population genetic structure of D. fissa was analyzed by Structure v. 2.3.4 (Pritchard et al. 2000). The number of clusters (K) was set as 2 to 15. Ten runs for K = 2 to 15 were analyzed, and a burn-in of 200000 followed by 1000000 Markov chain Monte Carlo (MCMC) iterations. The optimal K-value was determined with Structure Harvester 0.6.93 (Earl and Vonholdt 2012) based on the Delta K method. The results of 10 replicates for the chosen K-value were uploaded to Clumpak server (http://clumpak.tau.ac.il/ index) (Kopelman et al. 2015) to generate the final plots. Mantel test A Mantel test was conducted by Arlequin 3.5 with 1000 permutations to ensure that there was a relationship between the population differentiation (FST) based on mitochondrial combined sequence and geographic distance. Ancestral area reconstruction To trace the origin and diffusion process of D. fissa, Bayesian binary MCMC (BBM) analysis was implemented in RASP v4.2 (Yu et al. 2015). Five areas were considered based on the environmental conditions in the distribution range of D. fissa: (A) Central China, (B) South China, (C) Yunnan-Guizhou Plateau, (D) Southeast coastal area, and (E) Sichuan Basin. As input data, 20000 trees were produced by Beast 1.10.4 based on the mitochondrial dataset with HKY + F + G4 model (chosen with PhyloSuite). Clock rate and MCMC parameters were the same with BSP. 10000 random trees were selected to use for analysis after 4000 trees (20%) were discarded as burn-in. BBM analysis was implemented with default settings, except the number page 4 of 15Zoological Studies 62:32 (2023) © 2023 Academia Sinica, Taiwan of maximum areas, which was set to two because of the limited dispersal ability of the Meconematinae (Wang et al. 2019). RESULTS Distribution modeling Ecological niche models showed a high mean model fit (AUC > 0.9) indicating excellent performance of the models. Figure 2 shows the predicted areas of different periods. As for the present, the predicted areas were mainly the mountain regions including existing records. During the LGM, one high-suitability area is most obvious in the Sichuan basin. During the MidHolocene, a large low-suitability area existed in Yunnan province for the first time and the higher suitability areas were similar to the present-day predicted areas. Under RCP 4.5 in the year 2070, the high-suitability areas were slightly expanded, and other predicted areas were quite similar to the present-day predicted areas. Population genetic diversity and structure ITS sequences were 846 bp long, and 27 haplotypes were obtained based on 70 polymorphic sites. The haplotype diversity (Hd) and nucleotide diversity (π) are shown in table 2. For mitochondrial data, the combined sequences were 1215 bp long (COI 624 bp, and Cytb 591 bp). A total of 134 haplotypes were identified based on 165 polymorphic sites (91 sites from COI, and 74 sites from Cytb). The numbers of haplotypes, haplotype diversity (Hd), and nucleotide diversity (π) are shown in table 3. The genetic distances between populations based on Kimura’s 2-parameter and Tamura 3-parameter models are shown as a heatmap in figure 3 (details information Fig. 2. Potential distribution areas of D. fissa in different periods. Potential areas for (A) Current day; (B), Last Glacial Maximum; (C), MidHolocene; (D), year 2070 (RCP 4.5). N N N N page 5 of 15Zoological Studies 62:32 (2023) © 2023 Academia Sinica, Taiwan in Tables S3, S4). Pairwise FST between populations ranged from 0.02125 to 0.92767. The highest FST occurred between GXFCG and GZXN. The lowest FST occurred between GXFCG and YNHH (detailed information in Table S5, S6). The SAMOVA showed that the FCT value clearly increased from K = 2 to 3. Three groups (FCT = 0.427, p < .001) were delimited (S1: GDQY + GDZQ + GXGL + GXLB + GZTR + GZZY + HNCZ + HNCD + HNHH + HNZJJ + SCCD; S2: GZXN; S3: GXCZ + GXFCG + YNHH). AMOVA revealed that 59.68% genetic variation was found within populations, whereas 40.32% genetic variation was explained by differences among populations. Based on SAMOVA results, the most genetic variation was found within populations (45.71%), followed by the genetic variance among 3 groups (42.96%) and among populations within groups (11.60%). Results of the SAMOVA and AMOVA test are shown in table 4. We determined the population genetic structure of D. fissa. Plots of ΔK showed multiple peaks at K = 9 (Fig. 4). We found genetic similarities among populations GDQY, GDZQ, HNCZ and SCCD, and the individuals collected from populations GXCZ, GXFCG and YNHH demonstrate similarities. The geographic distance among populations ranged from 25.60 km (GDZQ vs. GDQY) to 1111.99 km (GDZQ vs. SCCD). Mantel tests indicated Fig. 3. The genetic distances among populations. A: based on Kimura’s 2-parameter; B: based on Tamura 3-parameter. Table 2. Diversity parameter based on ITS sequences pop code nHap Hd π Tajima’s DFu’s FS GZTR 39 4 0.358 0.00051 -1.30557 -1.34396 GZXN 7 3 0.667 0.00259 0.36328 1.50832 GZZY 6 2 0.333 0.00039 -0.93302 -0.00275 GXGL 15 4 0.619 0.00412 -1.34654 3.02493 GXCZ 6 1 0.000 0.00000 - - GXFCG 6 2 0.333 0.00039 -0.93302 -0.00275 GXLB 6 4 0.867 0.00426 -1.06907 0.56678 GDQY 17 5 0.691 0.00160 -0.79719 -0.57682 GDZQ 18 5 0.712 0.00107 0.11536 -1.50321 HNCZ 7 2 0.286 0.00068 -1.23716 0.85642 HNCD 14 3 0.604 0.00081 0.22615 0.08624 HNZZJ 15 2 0.248 0.00059 -0.51271 1.18757 HNHH 44 5 0.458 0.00536 -0.19705 6.24639 SCCD 10 2 0.356 0.00084 0.01889 1.52347 YNHH 19 4 0.380 0.00557 -0.32300 5.08703 Mean -0.52871 1.11051 Abbreviations: n, Number of samples; Hap, Number of haplotypes; Hd, Haplotype diversity; π, Nucleotide diversity. page 6 of 15Zoological Studies 62:32 (2023) © 2023 Academia Sinica, Taiwan a significant correlation between genetic differentiation (FST) and geographic distance (r = 0.3924, p < 0.001). Haplotype network and demographic history We obtained 134 haplotypes of combined sequences (designated as H1-H134). Among the 134 haplotypes, 14 haplotypes (H1, H3, H12, H17, H22, H24, H26, H42, H48, H49, H50, H76, H103 and H105) were shared by different population. H3 and H44 were both the most frequent haplotype shared by 10 samples, which suggests they are likely to be the ancestral haplotypes (Castelloe and Templeton 1994). H4 was only seen in population GZTR and H44 was shared by GZTR (n = 1), GDQY (n = 4) and GDZQ (n = 5). Divergence time analysis (MCC tree) revealed five mitochondrial lineages. The primary divergence within D. fissa commenced 1.31 Mya when two lineages (lineage IV and lineage V) split from the rest. Lineage IV and lineage V split around 1.01 Mya. The time of lineage III split from lineage I and two haplotypes (H73 and H74) was around 0.99 Mya (Fig. 5). The haplotype network (Fig. 6) showed phylogeographic relationships. H3 was the core of the “star-like” topology in lineage V. Tajima’s D test for lineage V (p < 0.01) was negative, while lineage I and lineage IV were negative but not significant. Fu’s FS test had significantly negative values (p < 0.001) for all three major lineages. Significantly negative values of Fu’s FS and Tajima’s D tests indicate a recent population expansion. The mismatch distribution analysis suggests that three lineages remained relatively stable. According to the BSP results, lineage I and lineage IV showed a population expansion in recent times (from around 0.025 and 0.05 Mya). Lineage V showed a constant population size for a long time, which suddenly decreased beginning around 0.01 Mya. Moreover, the BSP results showed that no lineage was affected by LGM (Fig. 7). Table 4. Results of SAMOVA and AMOVA analysis based on mitochondrial data SAMOVA group (K = 3) V% F-statistic Among groups 42.96 FCT = 0.427* Among populations within groups 11.60 FSC = 0.202* Within populations 45.71 FST = 0.543* AMOVA Among populations 40.32 FST = 0.40*** Within populations 59.68 Table 3. Diversity parameter based on combined mitochondrial sequences pop code nHap Hd π Tajima’s DFu’s FSSSD HRI GZTR 39 17 0.899 0.00737 0.07971 -0.63636 0.04929 0.04020 GZXN 7 6 0.952 0.00188 -1.03541 -2.91457 0.01626 0.10204 GZZY 6 6 - 0.00620 1.40404 -1.41226 0.02151 0.07111 GXGL 15 14 0.990 0.00756 -0.40346 -5.40267* 0.01730 0.04834 GXCZ 6 6 - 0.00455 -1.30088 -1.96438 0.05995 0.08889 GXFCG 6 3 0.733 0.00121 0.60031 0.46205 0.29990* 1.14667* GXLB 6 5 0.933 0.00916 0.10611 0.99560 0.07264 0.13778 GDQY 17 8 0.897 0.00577 -0.49237 1.30422 0.04859 0.06877 GDZQ 18 9 0.863 0.00353 -0.69798 -0.89728 0.02964 0.04891 HNCZ 7 4 0.714 0.00486 0.18346 2.10694 0.12366 0.32200 HNCD 14 10 0.945 0.00751 -0.27652 -0.65346 0.01574 0.02524 HNZZJ 15 15 - 0.00727 0.62230 -8.15016*** 0.00921 0.02068 HNHH 44 34 0.987 0.00862 -0.31133 -16.39408*** 0.00210 0.00427 SCCD 13 8 0.923 0.00684 -0.34641 0.65060 0.02633 0.04471 YNHH 19 10 0.912 0.00222 0.17169 -3.40730* 0.03212 0.09791 Lineage I 88 51 0.972 0.00673 -1.32809 -24.72986*** 0.00820 0.00695 Lineage IV 42 27 0.957 0.00767 -0.65331 -19.02936*** 0.01401 0.01154 Lineage V 79 38 0.957 0.00540 -1.96000** -24.89551*** 0.00396 0.01238 * p < 0.05. ** p < 0.01. *** p < 0.001. Abbreviations: n, Number of samples; Hap, Number of haplotypes; Hd, Haplotype diversity; π, Nucleotide diversity; SSD, Sum of squared deviations; HRI, Harpending’s raggedness index. page 7 of 15Zoological Studies 62:32 (2023) © 2023 Academia Sinica, Taiwan Ancestral area reconstruction The ancestral area reconstructions of D. fissa by BBM are shown in figure 8. BBM analyses indicated that D. fissa originated from the South China (area B) and dispersed to Central China (area A). Then, D. fissa dispersed to Yunnan-Guizhou Plateau (area C), Southeast coastal area (area D). The population of Sichuan Basin (area E) was from Central China and Yunnan-Guizhou Plateau (Fig. 9). Fig. 5. Divergence time estimation based on combined mitochondrial data of D. fissa. Fig. 4. Structure clustering results. (A) the posterior probability of each K; (B) the distribution of Delta K values; (C) Bayesian clustering results at K = 9; S1–3 was the groups defined by SOMOVA. page 8 of 15Zoological Studies 62:32 (2023) © 2023 Academia Sinica, Taiwan DISCUSSION The genetic diversity of the nuclear gene (ITS) is significantly lower than that of mitochondrial data. Because there is no significant sex-bias in migration capacity or adaptability in D. fissa, the reason for this condition may be the highly elevated mutation rate of mitochondrial DNA in natural populations or subsequent secondary contact after geographic isolation (Toews and Brelsford 2012). The low diversity also makes the ITS sequence unable to provide sufficient information in phylogeny. The 134 mtDNA haplotypes formed 5 haplotype lineages (lineage II only consisted of two haplotypes). Most haplotypes (120 of 134) were unique to one population, which indicated high level genetic differentiation among D. fissa populations. SAMOVA analysis showed that S2 (GZXN) and S3 (YNHH, GXFCG and GXCZ) formed a separate group respectively. The STRUCTURE analysis based on concatenated data also implied that S2 and S3 groups have unique genetic structures. The GZXN population had six unique haplotypes (H58 to H63), which connected with adjacent haplotype by more than 5 mutational steps. We also checked the difference between the results of SAMOVA and MCC tree. The S2 group defined by SOMAVA corresponds to a part of lineage I, and is connected with adjacent haplotype by more than 10 mutational steps. All haplotypes of the S2 group belong to lineage I. Lineages II, III, V and the rest of lineages I, IV correspond to the S1 group. High level haplotype diversity was observed in most populations, which may be a result of rapid population growth of ancestral populations (Yu et al. 2014). High genetic diversity often relates to the species with long evolutionary history and wide distribution (Li et al. 2019; Zhou et al. 2021). D. fissa is a longwinged Meconematinae species and its distribution is wider than most other species of Meconematinae in China. The Mantel test revealed a significant correlation between genetic differentiation (FST) and geographic distance. Dispersal ability is an important factor that affects the correlation between geographic distances and genetic diversity of the species (Wu et al. 2019). Significant FST values indicated that the D. fissa populations had undergone obvious genetic Fig. 6. The mitochondrial haplotype network of concatenated sequences. The dotted box represents the five lineages based on divergence time analysis. page 9 of 15Zoological Studies 62:32 (2023)