Patterns of Genetic Variation in the Eisenia nordenskioldi Complex (Oligochaeta: Lumbricidae) along an Elevation Gradient in Northern China

: Eisenia nordenskioldi is the dominant earthworm species in many tundra and boreal habitats. Nothing is known about the genetic diversity of this species along the elevation gradient in China. This study sampled 28 individuals in the E . nordenskioldi complex from Wuling Mountain, northern China, to examine their external morphology and genetic diversity. Mt. Wuling is the southern limit of the distribution of the E . nordenskioldi complex. The specimens from Mt. Wuling were classiﬁed into three groups along an elevation gradient. Mismatch distribution analysis suggested that the Pleistocene glaciations possibly did not signiﬁcantly affect the distribution of earthworm species in this region. We also found that elevation affected the genetic diversity, but not the external morphology of E . nordenskioldi . Given the altitudinal genetic diversity within the E . nordenskioldi complex, the phylogeography of this species provides important information for the zoogeographic reconstruction of the mountains in northern China. With the relatively limited sample size, the result is not conclusive, and further studies need to be conducted in the future to verify the results.


Introduction
The vertical zonation of mountain climates results in the formation of diverse, unique habitats for animals, plants, and microorganisms, giving rise to vertical differences in the biome along altitude gradients [1]. Xu et al. [2] showed that altitude did not affect the overall abundance of epigeic soil animals on Mt. Dongling but did affect the distribution of various functional feeding groups of animals at different altitudes. Few studies have examined earthworm diversity along an altitude gradient [3]. One suggested that the difference in earthworm species richness along an elevation gradient in the mountains of northeastern Puerto Rico was due to a combination of biotic and soil physical and chemical factors [4]. No study has examined the intra or inter-specific genetic diversity of Eisenia along an elevation gradient in northern China.
In northern China, Mt. Wuling is the main peak in the Yanshan Mountains, which have a typical warm-temperate, semi-humid, continental monsoon climate. It is a nature reserve with a total area of 143 km 2 and forest coverage of 76.2%. The vegetation and soil types vary along the elevation gradient [5].
Eisenia nordenskioldi, an earthworm species that is widespread in Northern Asia and adjacent regions, is known for its high morphological, karyotypic, and genetic variation [6,7]. The diagnostic features of this species include having purple dorsally and faint yellow ventrally; the body has no stripes; the intersegmental furrows in adult individuals are faint yellow; the clitellum is faint yellow, saddle-shaped in xxvii-xxxii; the spermathecal pores are paired in 9/10 and 10/11 ventrally; setae lumbricine, ab > bc, aa > bc; the spermathecae are tiny and ball-shaped; the gizzard is located in xvii-xviii. The Eisenia nordenskioldi complex contains two subspecies, the pigmented E. n. nordenskioldi and the unpigmented E. n. pallida [6,7]. Several genetic diversity and phylogeographic studies on E. nordenskioldi have been conducted previously [7][8][9][10][11][12][13][14]. In this study, we analyzed the genetic diversity of the E. nordenskioldi complex from Mt. Wuling along its elevation gradient and discussed its phylogeography.

Sample Collection and DNA Extraction
E. nordenskioldi is a common species widely distributed in tundra and boreal habitats in northeast Asia [6]. Nature reserves were selected with high, medium, and low latitudes in North China in order to discuss the genetic differentiation of this species at different altitudes. A total of 28 specimens of the E. nordenskioldi complex were collected by hand sorting in litter from four elevations in natural forest habitats on Mt. Wuling: 1000-1200, 1200-1400, 1400-1700, and 1700+ m (Table 1, Figure 1). Total DNA was extracted from the tail muscle of individuals fixed in 95% ethanol using a DNA extraction kit (Sangon, Shanghai, China) according to the manufacturer's instructions. Total genomic DNA was extracted using the Invitrogen Genomic DNA extraction kit according to the manufacturer's protocol. A COI gene fragment was amplified by PCR using primers LCO1490 (5 GGT CAA CAA ATC ATA AAG ATA TT 3 ) and HCO2198 (5 TAA ACT TCA GGG TGA CCA AAA AA 3 ), as described by [15]. The PCR products were sequenced directly using the Big Dye Terminator v3.1 Cycle Sequencing kit using these primers. Pairwise (p) distances were calculated using MEGA X [16]. Other samples of E. nordenskioldi were collected from nature reserves in northeast China (DQH, BS, SH, and LTDZ). The specimens from Russia were not collected, COI sequences from Russian earthworms were retrieved from GenBank (accession no. KU708313-708411 [9]), and those sequences of earthworms from Mongolia were taken from Blakemore [17]. The specimens collected from forests in northern China are stored in 95% ethanol in the Hebei Key Laboratory of Animal Diversity, Langfang Normal University, China. The net weight of the earthworm individuals was measured after ethanol fixing for two weeks.   Table 1 for coordinates).

Distribution of the E. nordenskioldi Complex along an Elevation Gradient
Group 1 (G1) was widely distributed not only in Mt. Wuling but also in northeast China. Group 2 (G2) was found only at low and mid-elevations on Mt Wuling. Groups 3 and 4 existed only in northeast China. Group 5 (G5) was distributed only at low elevation on Mt. Wuling. Some low-elevation locations from Mt. Wuling had their own set of haplotypes; these are described as cold intolerant with a narrow distribution (Figures 1 and 2). The topological relationships of the Mt. Wuling clusters were classified into three groups (G1, G2, and G5) along an elevation gradient: G1 could be described as cold tolerant as it is widely distributed in many zones (Mt. Wuling and northeast China) with very low temperature; G2 was adapted to either cold or warm habitats on Mt. Wuling, and gene flow was more frequent within this group; G5 existed only in a narrow lowaltitude zone (Figures 1 and 2). Gene flow estimated from the sequence data had values of DeltaSt = 0.05159, GammaSt = 0.08542, and Nm = 0.68. A haplotype network (Figure 3) was constructed using the 635-bp COI fragment of 128 individuals from different locations not only in Mt. Wuling but also in other adjacent regions (northeast China, Mongolia, and far east Russia). This indicated that the E. nordenskioldi complex in northern China comprises at least five exclusive groups, in accordance with the BI tree.

Distribution of the E. nordenskioldi Complex along an Elevation Gradient
Group 1 (G1) was widely distributed not only in Mt. Wuling but also in northeast China. Group 2 (G2) was found only at low and mid-elevations on Mt Wuling. Groups 3 and 4 existed only in northeast China. Group 5 (G5) was distributed only at low elevation on Mt. Wuling. Some low-elevation locations from Mt. Wuling had their own set of haplotypes; these are described as cold intolerant with a narrow distribution (Figures 1  and 2). The topological relationships of the Mt. Wuling clusters were classified into three groups (G1, G2, and G5) along an elevation gradient: G1 could be described as cold tolerant as it is widely distributed in many zones (Mt. Wuling and northeast China) with very low temperature; G2 was adapted to either cold or warm habitats on Mt. Wuling, and gene flow was more frequent within this group; G5 existed only in a narrow low-altitude zone (Figures 1 and 2). Gene flow estimated from the sequence data had values of DeltaSt = 0.05159, GammaSt = 0.08542, and Nm = 0.68. A haplotype network (Figure 3) was constructed using the 635-bp COI fragment of 128 individuals from different locations not only in Mt. Wuling but also in other adjacent regions (northeast China, Mongolia, and far east Russia). This indicated that the E. nordenskioldi complex in northern China comprises at least five exclusive groups, in accordance with the BI tree.    Table 1; Circle sizes are proportional to the number of individuals having this haplotype; COI sequences of Russian and Mongolian Groups were retrieved from Genbank.

Genetic Diversity and Differentiation of the Mt. Wuling E. nordenskioldi Complex
COI sequences were determined for 28 E. nordenskioldi specimens from Mt Wuling. The complete alignment included 635 sites (i.e., no length polymorphism was detected). There were 432 conserved sites and 203 variable sites, of which 182 were parsimony informative, and 21 were singleton sites. The overall average p distance was 0.174. The p distance was highest between L1 and H4 (0.208) and lowest between M3 and H4 (0.072). The p distance was highest within L1 (0.174). The genetic p distance was high both within and among the four elevation groups on Mt. Wuling. Figures 5 and 6 compare basic genetic diversity parameters among the four elevation groups, including the numbers of segregation sites, haplotypes and polymorphic sites, nucleotide and haplotype diversities, the total number of mutations, and average number of nucleotide differences. Mismatch distribution analysis suggested that the earthworm species had not experienced expansion (Figure 7). Table 2 shows the genetic diversity of E. nordenskioldi from different regions. Our E. nordenskioldi sequences fell into five groups in the Bayesian phylogeny, strongly supported by Bayesian posterior probabilities (Figures 1 and 2). Figures 2 and 3 present the data for G1 and G5 of the E. nordenskioldi complex showing both high haplotype and nucleotide diversities because of the elevation gradient. The probability obtained from a permutation test with 1000 replicates (PM test) was significant (p < 0.05) ( Table 3), indicating genetic differentiation among G1, G2, and G5 (Figures 1 and 2).
Length mm Width mm

Genetic Diversity and Differentiation of the Mt. Wuling E. nordenskioldi Complex
COI sequences were determined for 28 E. nordenskioldi specimens from Mt Wuling. The complete alignment included 635 sites (i.e., no length polymorphism was detected). There were 432 conserved sites and 203 variable sites, of which 182 were parsimony informative, and 21 were singleton sites. The overall average p distance was 0.174. The p distance was highest between L1 and H4 (0.208) and lowest between M3 and H4 (0.072). The p distance was highest within L1 (0.174). The genetic p distance was high both within and among the four elevation groups on Mt. Wuling. Figures 5 and 6 compare basic genetic diversity parameters among the four elevation groups, including the numbers of segregation sites, haplotypes and polymorphic sites, nucleotide and haplotype diversities, the total number of mutations, and average number of nucleotide differences. Mismatch distribution analysis suggested that the earthworm species had not experienced expansion (Figure 7). Table 2 shows the genetic diversity of E. nordenskioldi from different regions. Our E. nordenskioldi sequences fell into five groups in the Bayesian phylogeny, strongly supported by Bayesian posterior probabilities (Figures 1 and 2). Figures 2 and 3 present the data for G1 and G5 of the E. nordenskioldi complex showing both high haplotype and nucleotide diversities because of the elevation gradient. The probability obtained from a permutation test with 1000 replicates (PM test) was significant (p < 0.05) ( Table 3), indicating genetic differentiation among G1, G2, and G5 (Figures 1 and 2).          Abbreviations: S = Number of segregating sites, ps = S/n, Θ = ps/a1, π = nucleotide diversity, and D is the Tajima test statistic). NOTE: The analysis involved 174 nucleotide sequences. All positions containing gaps and missing data were eliminated. There was a total of 545 positions in the final dataset. Evolutionary analyses were conducted in MEGA X. Not significant p > 0.10.

Morphological and Genetic Variations within the E. nordenskioldi Complex According to Elevation
Eisenia nordenskioldi is known to have very high genetic diversity and contains several cryptic genetic lineages [7][8][9][10][11][12][13][14]. The E. nordenskioldi complex is separated into two species: E. nordenskioldi from northern and western Russia and Eisenia nordenskioldi from southern and southeastern Russia [13]. Both species need to be described and will be future work. The sequences of pigmented E. n. mongol and E. n. onon [17] were included within E. n. nordenskioldi from Russia, while unpigmented E. n. pallida from Korea [21] formed a separate branch within Eisenia nordenskioldi from southern and southeastern Russia. Our groups were morphologically similar to Eisenia from southern and southeastern Russia, characterized by body-color (purple dorsally and faint yellow ventrally), spermathecae shape (ball-shaped), and clitellum position (saddle-shaped in xxvii-xxxii).
It was indicated that altitudinal stratification has an effect on the genetic diversity of megascolecid worms from mountains in Taiwan, China [22]. The studies of Shekhotsov et al. [7][8][9][10][11][12][13][14] have not discussed the genetic diversity of E. nordenskioldi complex (Lumbricidae) in relation to its altitudinal distribution. In this work, we tackled the pattern of genetic variation in E. nordenskioldi along elevational gradients in northern China. Our results showed that elevation affected the genetic diversity, but not external morphological diversity, within E. nordenskioldi. The E. nordenskioldi complex from Mt. Wuling is not monophyletic and consists of three morphologically similar altitude groups (G1, G2, and G5) that have high genetic divergence (Figures 1 and 2).

Distribution of the E. nordenskioldi Complex from Mt. Wuling
Dong et al. [23] found two exclusive lineages for Amynthas triastriatus of the Megascolecidae in southeast China: Lineage A is distributed mainly at high altitudes and Lineage B mainly at low altitudes. Similarly, regarding the distribution of the E. nordenskioldi complex from Mt. Wuling (Figure 1), we hypothesized that G1 tolerates low temperatures in the tundra and high boreal habitats, while G5, which may be close to E. n. pallida, is intolerant to low temperatures and has a narrow regional distribution. We found that Mt. Wuling harbors high earthworm diversity. In the haplotype network (Figure 3), Mt. Wuling haplotypes were relatively close to Mongolian ones but far from Russian or northeast Chinese ones. This suggests a dispersal route from Siberia to Mongolia to Inner Mongolia to Mt. Wuling.

Phylogeography of the E. nordenskioldi Complex from Mt. Wuling
The Pleistocene glaciation resulted in dramatic shifts in animal habitats and a dramatic reduction in intra-specific diversity [24]. The cold Pleistocene glacial periods caused alpine insects to disperse into lowland regions [25]. Insects, such as bees and beetles, have significantly stronger dispersal ability compared with soil invertebrates, such as earthworms. Genetic differences were observed among different elevations on Mt. Wuling, which had their own haplotype sets (Figure 3). This indicated that the Pleistocene glaciations did not significantly affect the distribution of earthworm species. The analysis of allelic gene frequency and nucleotide mismatch implied that the events shaping the phylogeography of Mt. Wuling E. nordenskioldi occurred long before the last glacial maximum. Given this altitudinal genetic diversity within E. nordenskioldi, the phylogeography of this species provides important information for the zoogeographic reconstruction of the mountains in northern China.

Conclusions
Mt. Wuling is the southern limit of the distribution of the E. nordenskioldi complex. Our results showed that elevation has an effect on the genetic diversity of the E. nordenskioldi complex but not on the morphological diversity. The E. nordenskioldi complex in northern China has at least five exclusive genetic groups, and those of Mt. Wuling were classified into three groups (G1, G2, and G5) along an elevation gradient. Mismatch distribution analysis indicated that Pleistocene glaciations did not significantly affect the distribution of earthworm species in this region. Future work should examine the distribution and dispersal of the E. nordenskioldi complex in the Palearctic.