Population Genetic Structure Reveals Two Lineages of Amynthas triastriatus (Oligochaeta: Megascolecidae) in China, with Notes on a New Subspecies of Amynthas triastriatus

Amynthas triastriatus (Oligochaete: Megascolecidae) is a widely distributed endemic species in Southern China. To shed light on the population genetic diversity and to elucidate the population differentiation and dispersal of A. triastriatus, a population genetic structure study was undertaken based on samples from 35 locations collected from 2010 to 2016. Two exclusive lineages within A. triastriatus—lineage A and lineage B—were revealed. Lineage A was mainly distributed at high altitudes while lineage B was mainly distributed at low altitudes in Southeast China. The genetic diversity indices indicated that the populations of A. triastriatus had a strong genetic structure and distinct dispersal histories underlying the haplogroups observed in this study. Combined with morphological differences, these results indicated a new cryptic subspecies of A. triastriatus. Lineage A was almost degenerated to parthenogenesis and lineage B had a trend to parthenogenesis, which suggested that parthenogenesis could be an internal factor that influenced the differentiation and dispersal of A. triastriatus. The divergence time estimates showed that A. triastriatus originated around Guangxi and Guangdong provinces and generated into two main lineages 2.97 Ma (95%: 2.17–3.15 Ma) at the time of Quaternary glaciation (2.58 Ma), which suggested that the Quaternary glaciation may have been one of main factors that promoted the colonization of A. triastriatus.


Introduction
Understanding the evolutionary history of a species is a fundamental goal of evolutionary biology [1] and many approaches are based on population genetic structure analysis, both in fauna and flora [2,3]. Population genetic structure studies combined with geographical distribution, time of differentiation and ancient geographical events are commonly used to determine the differentiation and dispersal of populations and to analyze their influencing factors.
Earthworms are very important macro-invertebrate detritivores that tend to remain in the same areas during long periods of time and show low ability to cross mountains; for these reasons they are often an ideal model for phylogeographical study [4][5][6][7][8]. Earthworms play important roles in the agroecosystem as key organisms that influence soil structure formation, soil carbon dynamics and biogeochemical cycles [9][10][11], and hence failure to recognize accurate species boundaries within  Table 1.

DNA Extraction, Amplification and Sequencing
Total genomic DNA was extracted using the E.Z.N.A. ® Mollusc DNA kit (Omega Bio-Tek, Norcross, GA, USA) according to the manufacturer's instructions. The quality of DNA samples was checked by electrophoresis using a 1% agarose gel and stored at −20 °C.
The mitochondrial genes were amplified using 50 μL of reaction mixture with 1 μL template DNA, 0.6 μL Trans Taq ® HIFI DNA Polymerase (Trans, Beijing, China), 4 μL dNTP, 5 μL HiFi Buffer, 2 μL primer and 35.4 μL ddH2O. The PCR of genes COI, COII, ATP8, 12S rRNA, 16S rRNA and NDI were carried out for 5 min at 94 °C followed by 32 cycles of 94 °C for 30 s, 50 °C for 30 s and 72 °C for 1 min, with an extension of 10 min at 72 °C. The PCR of gene ND6 was carried out for 5 min at 94 °C followed by 40 cycles of 94 °C for 30 s, 51.9 °C for 30 s and 72 °C for 1 min, with an extension of 10 min at 72 °C. The primers of mitochondrial COI, COII, 12S rRNA, 16S rRNA and NDI are listed in Table 2. The primers of genes ATP8 and ND6 were designed in this paper.
The PCR products were visualized on a 1% agarose gel and subsequently purified then sequenced by the Beijing Genomics Institute (Shanghai, China) using an ABI 3730 DNA analyzer. The sequences were compared with known earthworm sequences in GenBank using the BLAST search algorithm [35].

Sequence Alignment, Population Genetic Structure and Genetic Diversity, Time of Divergence and Colonization History
Population genetic structure and genetic diversity studies were conducted using single-gene datasets of COI with one Aporrectodea species as an outgroup. Time of divergence and colonization history were conducted using a combination of COI, COII, ATP8, ND6, 12S rRNA, 16S rRNA and NDI datasets with one Aporrectodea species as an outgroup. Sequences of COI and a combined dataset were aligned using the Clustal X 1.8 program [39].  Table 1.

DNA Extraction, Amplification and Sequencing
Total genomic DNA was extracted using the E.Z.N.A. ® Mollusc DNA kit (Omega Bio-Tek, Norcross, GA, USA) according to the manufacturer's instructions. The quality of DNA samples was checked by electrophoresis using a 1% agarose gel and stored at −20 • C.
The mitochondrial genes were amplified using 50 µL of reaction mixture with 1 µL template DNA, 0.6 µL Trans Taq ® HIFI DNA Polymerase (Trans, Beijing, China), 4 µL dNTP, 5 µL HiFi Buffer, 2 µL primer and 35.4 µL ddH 2 O. The PCR of genes COI, COII, ATP8, 12S rRNA, 16S rRNA and NDI were carried out for 5 min at 94 • C followed by 32 cycles of 94 • C for 30 s, 50 • C for 30 s and 72 • C for 1 min, with an extension of 10 min at 72 • C. The PCR of gene ND6 was carried out for 5 min at 94 • C followed by 40 cycles of 94 • C for 30 s, 51.9 • C for 30 s and 72 • C for 1 min, with an extension of 10 min at 72 • C. The primers of mitochondrial COI, COII, 12S rRNA, 16S rRNA and NDI are listed in Table 2. The primers of genes ATP8 and ND6 were designed in this paper.
The PCR products were visualized on a 1% agarose gel and subsequently purified then sequenced by the Beijing Genomics Institute (Shanghai, China) using an ABI 3730 DNA analyzer. The sequences were compared with known earthworm sequences in GenBank using the BLAST search algorithm [35][36][37][38].

Sequence Alignment, Population Genetic Structure and Genetic Diversity, Time of Divergence and Colonization History
Population genetic structure and genetic diversity studies were conducted using single-gene datasets of COI with one Aporrectodea species as an outgroup. Time of divergence and colonization history were conducted using a combination of COI, COII, ATP8, ND6, 12S rRNA, 16S rRNA and NDI datasets with one Aporrectodea species as an outgroup. Sequences of COI and a combined dataset were aligned using the Clustal X 1.8 program [39].
Models for sequence evolution and corresponding parameters were estimated using jModeltest v.2.1.7 [40]. The best-fitting models were used in the phylogenetic estimation for the maximum likelihood (ML) and Bayesian phylogeny estimation (BI) analyses. jModeltest revealed TPM3uf+G to be the best-fitting model for the COI gene (679 bp) with a gamma shape parameter of 0.2360, and GTR+G to be the best-fitting model for 4691 bp segments comprising mitochondrial genes COI, COII, ATP8, ND6, 12S rRNA, 16S rRNA and NDI with a gamma shape parameter of 0.3950.
The phylogenetic trees were reconstructed using different methods based on single-gene COI datasets of 65 A. triastriatus individuals and one individual of Aporrectodea trapezoides as outgroup. ML analyses were performed using PhyML v.3.0 [41]; clade support was assessed using bootstrapping with 1000 pseudo replicates. BI was performed using MRBAYES v.3.2 [42]; the parameters were set to 10,000,000 with trees sampled every 1000 generations, discarding the first 10% as burn-in. Posterior probabilities and bootstrap support for each branch were calculated from the sampled trees.
Population structure was evaluated in Arlequin v.3.5.2 [43] via AMOVA with 10,000 permutations for statistical confidence. A haplotype network was constructed using the median joining (MJ) method Network v.4.6.1.6 [44] to infer the relationships among haplotypes and their geographical distribution. Nucleotide mismatch distribution analyses of lineages A and B were calculated with DnaSP 5.0, and neutrality tests were used to test population equilibrium.
Basic statistics of mtDNA diversity, including nucleotide and haplotype diversity, Tajima's D [45] and Fu's Fs [46] neutrality tests and nucleotide mismatch distribution analysis with and within lineages were calculated with DnaSP 5.0 [47]. Bayesian Skyline Plot (BSP) was constructed with BEAST 1.8.2 [48] and Tracer 1.7 [49] to infer the timing of the population events.
The Bayesian tree and the time of divergence of A. triastriatus were estimated using BEAST 1.8.2 [48] based on 4691 bp segments comprising mitochondrial genes COI, COII, ATP8, ND6, 12S rRNA, 16S rRNA and NDI. These were from 18 individuals from 18 locations of lineage A, 17 individuals from 17 locations of lineage B and one individual of Aporrectodea trapezoides as an outgroup. A previous study has estimated 2.1 My −1 as the average substitution rate for COI, COII, ATP8, ND6, NDI, 12S and 16S genes [50], and this value was applied in the present study. Since earthworm lack a fossil record, the age of the calibration clade between A. triastriatus and Aporrectodea trapezoides was set at a mean age of 212 Ma (95%: 200-230 Ma); this represented the minimum age of Lumbricoidea and Megascolecoidea in our tree and is considered to be when the breakup and drift of Gondwana occurred [50][51][52][53]. Five replicate runs were performed with 10,000,000 MCMC steps, and the first 10% served as burn-in. The output was visualized in Tracer v1.6 to determine whether convergence and suitable effective sample sizes were achieved for all parameters. The maximum clade credibility tree was accessed using Tree Annotator v1.8.2 [48]. The tree was visualized and edited using Fig Tree v2.1.4.

Population Genetic Structure and Genetic Diversity
The posterior probabilities from the Bayesian analysis and the bootstrap values from the ML analysis were plotted on the ML tree ( Figure 2). A. triastriatus fell into two exclusive lineages, A and B, and this result was strongly supported in both the BI tree (1.00) and ML tree (100%). Lineage A consisted of four clades (clades 1-4) and lineage B consisted of two clades (clades 5, 6). Clade 1 was mainly distributed in Guizhou, Hubei, Chongqing and Sichuan provinces (10 of 15 individuals), clades 2 and 3 were mainly distributed in Anhui and Fujian provinces (six of 11 individuals) and clade 4 was mainly distributed in Guangxi, Guangdong and Hunan provinces (six of eight individuals). Clade 5 was mainly distributed in East China, such as Zhejiang, Anhui province, east of Jiangxi and Fujian province (18 of 22 individuals) and clade 6 was mainly distributed in the west of Southern China, such as Sichuan, Hunan, Guangdong, Guangxi and Guizhou provinces (six of nine individuals).
A1 of the Appendix. The mean base composition of the fragment showed a strong bias of A+T (A: 35.03%; C: 18.96%; G: 16.02%; T: 29.99%), as commonly found in invertebrate mitochondrial genomes [54]. The network showed that there were 31 nucleotide changes between lineages A and B, one within lineage B and 1-2 within lineage A. The sequenced region contained 16 polymorphic sites defining 14 haplotypes of lineage A, two polymorphic sites defining three haplotypes of lineage B and 49 polymorphic sites defining 17 haplotypes of A. triastriatus (Table 3). Haplotype and nucleotide diversities were used as measures of genetic diversity between and within lineages. General haplotype diversity was 0.871 for lineage A but only 0.123 for lineage B, and was 0.758 for A. triastriatus. The nucleotide diversity was low in lineage A and extremely low in lineage B. The genetic distances were 0.5% for lineage A, 0.1% for lineage B and 6.3% for A. triastriatus. These results showed that lineages A and B had exclusive haplotypes and suggested a great genetic differentiation exists between lineages.

Time of Divergence and Colonization History
As shown in Figure Table  4).
Hierarchical population structure in A. triastriatus was studied using AMOVA to check structure at the lineage level based on a fragment of the COI gene ( Table 5). The percentage of variation among lineages was 93.89%, and within lineages was just 6.11%. The FST value was 0.9389 (Nm = 0.02), indicating significant population genetic structure. The value of pairwise comparisons of ST among lineages was significantly high, suggesting a great genetic differentiation exists between lineages A and B.
Tajima's D values were negative within lineages, but positive and significant between lineages, indicating that the COI gene was possibly not a neutral evolution. Fu's Fs test is more sensitive to detect population expansion [46]; Fu's Fs test values were negative and nonsignificant within lineages, indicating that lineages A and B had experienced demographic expansion. Fu's Fs test values were positive and nonsignificant between lineages, indicating that A. triastriatus had experienced demographic expansion. The results of nucleotide mismatch distribution analysis based on fragments of COI, COII, ATP8, ND6, 12S rRNA, 16S rRNA and NDI genes showed the following: a transition between unimodal and bimodal distributions for COI and 12S, the same as expected for COII, 16S, ND1, ATP8 and ND6 of lineage A; unimodal distribution for COI and ND1, the same as expected for COII, 16S, ATP8 and ND6; and bimodal distribution for the 12S gene of lineage B. The results revealed that lineage A ( Figure 5a) and lineage B (Figure 5b) were not steady-state lineages, suggesting an older expansion. BSP revealed an explicit demographic history for the populations of A. triastriatus ( Figure 6): a flat population history during previous climate fluctuations of last glacial period (MIS 2-4), a sudden population decline during last glacial period, a rapid population expansion after last glacial period.  Table S1 of the Supplementary Materials. The mean base composition of the fragment showed a strong bias of A+T (A: 35.03%; C: 18.96%; G: 16.02%; T: 29.99%), as commonly found in invertebrate mitochondrial genomes [54]. The network showed that there were 31 nucleotide changes between lineages A and B, one within lineage B and 1-2 within lineage A. The sequenced region contained 16 polymorphic sites defining 14 haplotypes of lineage A, two polymorphic sites defining three haplotypes of lineage B and 49 polymorphic sites defining 17 haplotypes of A. triastriatus (Table 3). Haplotype and nucleotide diversities were used as measures of genetic diversity between and within lineages. General haplotype diversity was 0.871 for lineage A but only 0.123 for lineage B, and was 0.758 for A. triastriatus. The nucleotide diversity was low in lineage A and extremely low in lineage B. The genetic distances were 0.5% for lineage A, 0.1% for lineage B and 6.3% for A. triastriatus. These results showed that lineages A and B had exclusive haplotypes and suggested a great genetic differentiation exists between lineages. Table 3. Summary of genetic diversity measures and neutrality tests between and within lineages based on the COI gene, and the mean number of pairwise distances based on a fragment of the COI gene and on a combination of seven genes.

Time of Divergence and Colonization History
As shown in Figure 4, the root of the tree, corresponding to the diversification of A. triastriatus and Aporrectodea trapezoides, was estimated to occur about 215. 28 Table 4.   Table 4. Hierarchical population structure in A. triastriatus was studied using AMOVA to check structure at the lineage level based on a fragment of the COI gene ( Table 5). The percentage of variation among lineages was 93.89%, and within lineages was just 6.11%. The F ST value was 0.9389 (N m = 0.02), indicating significant population genetic structure. The value of pairwise comparisons of Φ ST among lineages was significantly high, suggesting a great genetic differentiation exists between lineages A and B. The results of nucleotide mismatch distribution analysis based on fragments of COI, COII, ATP8, ND6, 12S rRNA, 16S rRNA and NDI genes showed the following: a transition between unimodal and bimodal distributions for COI and 12S, the same as expected for COII, 16S, ND1, ATP8 and ND6 of lineage A; unimodal distribution for COI and ND1, the same as expected for COII, 16S, ATP8 and ND6; and bimodal distribution for the 12S gene of lineage B. The results revealed that lineage A (Figure 5a) and lineage B (Figure 5b) were not steady-state lineages, suggesting an older expansion. BSP revealed an explicit demographic history for the populations of A. triastriatus ( Figure 6

Application of Molecular Techniques in Taxonomy and Differentiation of Earthworms
Earthworms are important soil animals, but their simple structure and burrowing nature makes their taxonomy challenging and has led to a gross underestimation of the true level of earthworm biodiversity [55]. DNA taxonomy and associated molecular tools might be the only way to reveal the true level of biodiversity [56]. In recent years there has been increasing research use of sequence data along with morphological characters to distinguish earthworm species and to erect cryptic species [6,57,58].
Through analysising the population differentiation and diffusion process, population genetic structure studies and phylogenetic studies could provide an important basis to illustrate the classification statue of species. In previous study, phylogenetic trees and high haplotypes diversity indicated two distinct clades of Drawida japonica Michaelsen, 1892 collected in Shandong and Liaodong peninsulas [17]. A high genetic diversity suggests the presence of five cryptic allopatric species collected from the central Iberian Peninsula [15]. Lineages of A. triastriatus had a strong genetic structure and distinct demographic histories underlying the haplogroups observed in this study. The phylogenetic trees showed that A. triastriatus had two exclusive lineages, a differentiation was also shown by AMOVA, high ST index and the high percentage of variation between lineages. Taken together, these results indicated a high degree of isolation and a great genetic differentiation exists between two lineages.
Taxonomically speaking, COI gene divergence of Aporrectodea caliginosa, Lumbricus castaneus, Lumbricus terrestris and Satchellius mammalis was below 4% and ranged between 0.0% and 2.0% within Aporrectodea longa [12,59]. The interspecific COI gene divergence values of invasive Asian earthworms in the northeast United States ranged from 15.84% to 24.03%, and intraspecific COI gene divergence values ranged from 0.01% to 0.4% [29]. In this study, the COI gene divergence between the two lineages of A. triastriatus was 6.3%, which is neither in the intraspecific or interspecific range. However, there were also certain morphological differences between the two lineages. Thus, it is proposed to erect a new subspecies of A. triastriatus (Chen 1946), named Amynthas triastriatus usualis subsp. nov.; the description for this proposed subspecies is included as the final part of this paper.

Application of Molecular Techniques in Taxonomy and Differentiation of Earthworms
Earthworms are important soil animals, but their simple structure and burrowing nature makes their taxonomy challenging and has led to a gross underestimation of the true level of earthworm biodiversity [55]. DNA taxonomy and associated molecular tools might be the only way to reveal the true level of biodiversity [56]. In recent years there has been increasing research use of sequence data along with morphological characters to distinguish earthworm species and to erect cryptic species [6,57,58].
Through analysising the population differentiation and diffusion process, population genetic structure studies and phylogenetic studies could provide an important basis to illustrate the classification statue of species. In previous study, phylogenetic trees and high haplotypes diversity indicated two distinct clades of Drawida japonica Michaelsen, 1892 collected in Shandong and Liaodong peninsulas [17]. A high genetic diversity suggests the presence of five cryptic allopatric species collected from the central Iberian Peninsula [15]. Lineages of A. triastriatus had a strong genetic structure and distinct demographic histories underlying the haplogroups observed in this study. The phylogenetic trees showed that A. triastriatus had two exclusive lineages, a differentiation was also shown by AMOVA, high Φ ST index and the high percentage of variation between lineages. Taken together, these results indicated a high degree of isolation and a great genetic differentiation exists between two lineages.
Taxonomically speaking, COI gene divergence of Aporrectodea caliginosa, Lumbricus castaneus, Lumbricus terrestris and Satchellius mammalis was below 4% and ranged between 0.0% and 2.0% within Aporrectodea longa [12,59]. The interspecific COI gene divergence values of invasive Asian earthworms in the northeast United States ranged from 15.84% to 24.03%, and intraspecific COI gene divergence values ranged from 0.01% to 0.4% [29]. In this study, the COI gene divergence between the two lineages of A. triastriatus was 6.3%, which is neither in the intraspecific or interspecific range. However, there were also certain morphological differences between the two lineages. Thus, it is proposed to erect a new subspecies of A. triastriatus (Chen 1946), named Amynthas triastriatus usualis subsp. nov.; the description for this proposed subspecies is included as the final part of this paper.

Parthenogenesis and the Demographic and Dispersal of A. triastriatus
Parthenogenesis is a unisexual reproduction pattern corresponding to bisexual reproduction that is present in a few animal and plant species. Generally, the loss of sexual reproduction has been thought a dead end in evolution, leading to early extinction [60,61]. But in contrast to this perception, parthenogenesis has been observed in many species of earthworm (Lumbricidae and Megascolecidae) and parthenogenesis in earthworms is often related to polyploidy or aneuploidy [62][63][64][65][66][67][68][69][70]. Muldal (1952) points out that parthenogenesis is important as it makes the retention of polyploidy possible, and also favors the spread of polyploid forms into new areas, since even a single parthenogenetic individual may establish a population. The association between parthenogenesis and high polyploidy in earthworms produces an unexpected level of heterozygosity, an advantageous condition that provides resistance to environmental stress [65,69].
Parthenogenesis often results in polymorphism in earthworms [67], with morphological variability mainly relating to the reduction of reproductive structures such as seminal vesicles, spermathecae, prostates, and an empty seminal chamber. [69]. In the case of A. triastriatus, lineage A with a thin and lustreless seminal chamber and no prostate gland observed was almost degenerated to parthenogenesis, while lineage B with a plump and glossy seminal chamber and small prostate glands had a tendency to parthenogenetic reproduction. Compared to their sexual relatives, parthenogenetics often occur at high latitudes, high altitudes, on islands or island-like habitats, in xeric environments or in disturbed habitats [71][72][73][74][75]. Lineage A was mainly distributed at high altitudes in Southwest China, whereas lineage B was mainly distributed at low altitudes in Southeast China. Unisexual lineages can be highly ecologically successful, with broad geographical ranges that overlap and/or are more extensive than those of the sexual ancestor [76]. For example, the locations of A4 and B5, A8 and B7, A18 and B14 overlapped, and the AMOVA F ST values indicated that lineages A and B had a high degree of isolation and a significant genetic differentiation exists with lineage A having greater genetic diversity than lineage B. These findings supported those of Simon et al. (2003), that unisexual lineages with high genetic diversity can live in wider geographical ranges than their sexual ancestor.
Based on phylogenetic studies, parthenogenetic lineages usually occupy terminal nodes of phylogenetic trees; therefore, unisexual forms are usually assumed to arise from sexual congeners [76]. The phylogenetic trees produced from 321 species of Megascolecidae earthworms from China showed that all of these species fell into 14 groups, with A. triastriatus located at the terminal node of the terminal clade of the fourteenth group [52]. In the present study, lineage A was at the terminal node too. Additionally, the time of divergence for lineage A (0.48 Ma; 95%: 0.22-0.53 Ma) was later than for lineage B (0.72 Ma; 95%: 0.29-1.47 Ma), indicating that these unisexual forms did indeed arise from sexual congeners, which was in accordance with previous findings [28,53,77].

Differentiation and Colonization of A. triastriatus
The divergence time estimate of 215.28 Ma for outgroup Aporrectodea trapezoides (95%: 200.41-229.89 Ma) was not significantly different from previous estimates of the diversification between Lumbricidae and Megascolecidae at about 200-220 Ma [50][51][52][53]. Combining the results of population genetic structure and phylogenetic analysis, as well as paleogeographic events and paleoclimate, it was attempted to determine the differentiation and colonization routes of lineages A and B. In a previous study, the Megascolecidae family of China was separated into 14 groups after pheretima spread into China from the Indo-China Peninsula, then spread to the east or north [52]. Additionally, the split between A. triastriatus and GX201201-11A (an unknown species in the same cluster as A. triastriatus collected in Guangxi province, China) was estimated to have taken place 6.49 Ma (95%: 5.23-7.88 Ma) in Guangxi province [34].
It is proposed that A. triastriatus then formed into two main lineages (A and B) at around 2.97 Ma (95%: 2.17-3.15 Ma), which happened together with the Quaternary glaciation (2.58 Ma). The diversification of lineage A was divided into three stages.
The first stage (node 18 in Figure 4) occurred at the end of Zhonglianggan glaciation (MIS 12) and resulted from the Kunlun-Yellow River tectonic movement (1.2-0.6 Ma BP) that occurred at the transition from the Early to Middle Pleistocene [78,79]. Lineage A (Hap_14) dispersed to Guangdong province (Figure 7). The Heishiding Mountains in Guangdong province was a possible glacial refugium for Liquidambar species during the Quaternary glaciation [80], and the network analysis showed that hap_14 (GD56, collected in Heishiding mountains) was the hap that connected to lineage B. And BSP analysis shown that the Quaternary glaciation have influence on the diffusion of A. triastriatus, significantly several refuge should exist during this period. Therefore, the Heishiding Mountains is proposed as a potential refugium for lineage A.
(～0.15 Ma BP) occurred between the Middle and Late Pleistocene to raise the Tibetan Plateau to its present height, High-Asia experienced several temperature plunges, and glaciation existed in most of the mountains in western China [87][88][89][90]. Lineage A dispersed to the east to Fujian province, followed by dispersal to Guangxi province and north to Hunan and Jiangxi provinces. Humid and suitable climate conditions in South China promoted its colonization.
Unlike lineage A, lineage B mainly dispersed to the east. The diversification of lineage B was also divided into three stages. The first stage (nodes 2, 3, 4 and 6 in Figure 4) occurred during the Wangkun Glaciation (MIS [16][17][18][19][20] in which Kunlun-Yellow River tectonic movement not only strongly uplifted the Tibetan Plateau but also the surrounding mountains. Coverage of the Tibetan Plateau by ice sheet peaked during MIS 16 at 18 times the present glacial coverage [91,92]; the mid-Pleistocene Revolution (~0.9 Ma), the most significant transition boundary, led to a significant drop in global temperature and extension of the ice sheet [92]. In this period, lineage B dispersed to Guangdong, Fujian and Zhejiang provinces, far from the Tibetan Plateau. Guangdong (Heishiding), Fujian (Wuyi Mountain) and Zhejiang (Tianmu Mountain) were proposed to be refugia of flora and fauna [86,93]. A lineage B was mainly distributed in Zhejiang province and the surrounding area, it is suggested that this glaciation period blocked the colonization of additional locations by lineage B, resulting in refugia in Zhejiang province.   [78,[81][82][83][84]. Lineage A dispersed north to the mountainous terrain in the Sichuan, Chongqing and Guizhou provinces. Guizhou province, with its comparably stable ecological conditions during environmental fluctuations, is a refugium for Ginkgo biloba [85,86], and lineage A populations from Guizhou province contained 4 haplotypes with high levels of genetic diversity. Therefore, the Sichuan, Chongqing and Guizhou provinces are proposed as potential refugia for lineage A.
The third stage occurred during the last glacial period (MIS 2-4), in which the Gonghe Movement (~0.15 Ma BP) occurred between the Middle and Late Pleistocene to raise the Tibetan Plateau to its present height, High-Asia experienced several temperature plunges, and glaciation existed in most of the mountains in western China [87][88][89][90]. Lineage A dispersed to the east to Fujian province, followed by dispersal to Guangxi province and north to Hunan and Jiangxi provinces. Humid and suitable climate conditions in South China promoted its colonization.
Unlike lineage A, lineage B mainly dispersed to the east. The diversification of lineage B was also divided into three stages. The first stage (nodes 2, 3, 4 and 6 in Figure 4) occurred during the Wangkun Glaciation (MIS [16][17][18][19][20] in which Kunlun-Yellow River tectonic movement not only strongly uplifted the Tibetan Plateau but also the surrounding mountains. Coverage of the Tibetan Plateau by ice sheet peaked during MIS 16 at 18 times the present glacial coverage [91,92]; the mid-Pleistocene Revolution (~0.9 Ma), the most significant transition boundary, led to a significant drop in global temperature and extension of the ice sheet [92]. In this period, lineage B dispersed to Guangdong, Fujian and Zhejiang provinces, far from the Tibetan Plateau. Guangdong (Heishiding), Fujian (Wuyi Mountain) and Zhejiang (Tianmu Mountain) were proposed to be refugia of flora and fauna [86,93]. A lineage B was mainly distributed in Zhejiang province and the surrounding area, it is suggested that this glaciation period blocked the colonization of additional locations by lineage B, resulting in refugia in Zhejiang province.
The second stage (nodes 11 and 13) occurred during MIS 6 and 8, two glacial periods of the Penultimate Glaciation, with lineage B still dispersed in Zhejiang province. The third stage, as for lineage A, occurred during the last glacial period (MIS 2-4). One route of lineage B dispersal was to Anhui and Jiangxi provinces and another route was north to Guizhou and Guangxi provinces. Geological and climatic changes during this period had a profound effect on the colonization of lineage B, and it is also possible that the distant dispersal of B12 was affected by other factors such as birds or human beings.
The population structure of Drawida japonica Michaelsen, 1892 were analysed using 79 samples obtained from the Shandong and Liaodong peninsulas. Results show that the intraspecific genetic diversity of D. japonica has been influenced by geographic isolation [17]. The results of population genetic structure analysis of hormogastrid earthworms (Hormogaster elisae Álvarez 1977) collected from the central Iberian Peninsula and some invasive species (Amynthas corticis (Kinberg, 1867), Amynthas gracilis (Kinberg, 1867)) collected from a volcanic island show that environmental factors may have some influence on earthworms' genetic evolution [15,94]. The quaternary glaciation resulted in repeated drastic environmental changes that profoundly shaped the current distributions and genetic structure of many plants and animal species in temperate zones of the Northern Hemisphere [95][96][97]. As the youngest major glacial period, Quaternary glaciation ended in a recessional event extending from 0.015 to 0.01 Ma at mid-latitudes [98,99]. The present-day geographical distribution pattern of A. triastriatus formed about 0.01 Ma, with its dispersal process totally covered by the Quaternary glaciation period. This period had a major impact on climate and vegetation in Southeast Asia through the interaction of temperature, rainfall and topography [100,101]. We suggest that the Quaternary glaciation period promoted the differentiation and colonization of A. triastriatus. Drawida japonica, H. elisae, A. corticis, A. gracilis and A. triastriatus are wildspread species, ancient geography event, climate changes, geographic isolation (island, sea, river) and environment factors should play important role in the differentiation and diffusion of wildspread earthworm species.
The phylogenetic evaluation of the genus Amynthas was carried out using 77 Amynthas species from South China. It has been determined that at least one branch of Amynthas spread to the southeast of China, and another spread to the southwest of China [28]. In the present study, the spread of A. triastriatus was consistent with the research of Sun et al., with lineage A widely distributed in the southwest of China, and lineage B widely distributed in the southeast of China. Overall, the colonization routes of A. triastriatus were from south to north, which is in accord with the known dispersal direction of earthworms in China [52,102].

Conclusions
Amynthas triastriatus (Oligochaete: Megascolecidae) of the aeruginosus-group is a widely distributed endemic earthworm species in Southern China. Population genetic structure exhibited a very high degree of isolation and two exclusive lineages based on existing specimen and gene data. Lineage A was almost degenerated to parthenogenesis, lineage B had a tendency to parthenogenesis. There were distinguish genetic differentiations exist between lineages and lineage A had high haplotype diversity and strong genetic structure, additionally, lineage A mainly distributied at high alitude areas. These results indicated a new cryptic species of A. triastriatus; combined with morphological differences, a new subspecies of A. triastriatus is proposed. And the Quaternary glaciation may be one of main factors that promoted the colonization of A. triastriatus. In future, there may be more studies of differentiation and dispersal of A. triastriatus including more populations based on this study, and the population genetic structure analysis should play an important role in the earthworm taxonomy.

Description of a New Subspecies of Amynthas Triastriatus
Based on the above results, to unravel the classification puzzle of Amynthas triastriatus it is proposed that a new subspecies be established that encompasses the populations of lineage B. VII, two pairs presetal and two pairs postsetal on VIII, three pairs presetal on IX, 1/3 circumference ventrally apart from each other (paratypes: numbers of genital papillae variable from [14][15][16][17][18]. One pair of male pores in XVIII, 1/3 circumference apart ventrally, each on the top of a raised, round porophore surrounded by 3-5 circular ridges, two presetal and two postsetal genital papillae on median side. Two big collapse-topped genital papillae present on XVIII, 1/4 circumference ventrally apart from each other. (Figure 8A). Single female pore in XIV, milky.
Internal characters. Septa 5/6-7/8, thick and muscular, 10/11-12/13 slightly thickened, 8/9-9/10 absent. Gizzard bucket-shaped, in IX-X. Intestine enlarged distinctly from XV. One pair of intestinal caeca in XXVII, brown, simple, dorsal margin smooth, 17 short pointed saccules in ventral margin, extending anteriorly to XXIII. Four pairs of esophageal hearts in X-XIII. Ovaries in XIII. Two pairs of spermathecae in VIII-IX, about 1.7-2.0 mm long. Ampulla oval-shaped, stout duct as long as 1/5 ampulla. Diverticulum as long as 1/3 main pouch, terminal 1/2 dilated into an oval-shaped glossy seminal chamber. (Figure 8B). Holandric: two pairs of testis sacs, in X-XI undeveloped. Two pairs of seminal vesicles in XI and XII, the second pair more developed than the first pair. Holotype with one pair of thick, small, nubby lobate prostate glands and S-curved duct. Paratypes with one pair of prostate glands or only one side observed and other side S-curved duct ( Figure 8C). Ethymology. The species is named after the existence of prostate glands. Remarks. There are certain differences between the new subspecies and A. triastriatus. For instance, the new subspecies has small prostate glands that are absent in A. triastriatus, and the seminal chamber of the new subspecies is plump and glossy while thin and lustreless in A. triastriatus. Additionally, the first dorsal pore of the new subspecies is located at 11/12, but at 10/11 in A. triastriatus. Furthermore, the new subspecies has more papillae in both the spermathecal and male pore regions than A. triastriatus, and more finger sacs in the ventral intestinal caeca than A. triastriatus. Internal characters. Septa 5/6-7/8, thick and muscular, 10/11-12/13 slightly thickened, 8/9-9/10 absent. Gizzard bucket-shaped, in IX-X. Intestine enlarged distinctly from XV. One pair of intestinal caeca in XXVII, brown, simple, dorsal margin smooth, 17 short pointed saccules in ventral margin, extending anteriorly to XXIII. Four pairs of esophageal hearts in X-XIII. Ovaries in XIII. Two pairs of spermathecae in VIII-IX, about 1.7-2.0 mm long. Ampulla oval-shaped, stout duct as long as 1/5 ampulla. Diverticulum as long as 1/3 main pouch, terminal 1/2 dilated into an oval-shaped glossy seminal chamber. (Figure 8B). Holandric: two pairs of testis sacs, in X-XI undeveloped. Two pairs of seminal vesicles in XI and XII, the second pair more developed than the first pair. Holotype with one pair of thick, small, nubby lobate prostate glands and S-curved duct. Paratypes with one pair of prostate glands or only one side observed and other side S-curved duct ( Figure 8C).
Ethymology. The species is named after the existence of prostate glands.
Remarks. There are certain differences between the new subspecies and A. triastriatus. For instance, the new subspecies has small prostate glands that are absent in A. triastriatus, and the seminal chamber of the new subspecies is plump and glossy while thin and lustreless in A. triastriatus. Additionally, the first dorsal pore of the new subspecies is located at 11/12, but at 10/11 in A. triastriatus. Furthermore, the new subspecies has more papillae in both the spermathecal and male pore regions than A. triastriatus, and more finger sacs in the ventral intestinal caeca than A. triastriatus.