Genetic Diversity and Population Structure of Spirobolus bungii as Revealed by Mitochondrial DNA Sequences

Simple Summary Spirobolus bungii plays a key role in soil ecolsystems, but there are few studies on this species, especially regarding genetics. In this study, 166 S. bungii individuals were collected from two cities in China (Nanjing and Tianjin), and COX2 and Cytb gene sequences of mitochondrial DNA and 18S rRNA gene fragments were sequenced. We conducted a population genetic analysis based on mtDNA sequences, revealing the population genetic diversity and genetic structure of S. bungii in the two cities, and confirming the hindering effect of geographical barriers on the gene flow of the populations. Abstract Soil macrofauna, such as Spirobolus bungii, are an important component of ecosystems. However, systematic studies of the genetic diversity, population genetic structure, and the potential factors affecting the genetic differentiation of S. bungii are lacking. We performed a population genetic study of 166 individuals from the mountains to the south of the Yangtze River, north of the Yangtze River in Nanjing city, and near Tianjin city, in order to investigate the correlations between geographical distance and genetic diversity. A total of 1182 bp of COX2 and Cytb gene sequences of mitochondrial DNA, and 700 bp of the 18S rRNA gene sequence were analyzed. There were two haplotypes and one variable site in the 18S rRNA gene, and 28 haplotypes and 78 variable sites in the COX2 and Cytb genes. In this study, the 18S rRNA gene was used for species identification, and mtDNA (concatenated sequences with Cytb and COX2) was used for population genetic analysis. Structure cluster analysis indicated that the genetic structures of the different populations of S. bungii tended to be consistent at small geographical scales. Phylogenetic trees revealed that the haplotypes were clearly divided into three branches: the area south of the Yangtze River, the area to the north of the Yangtze River in Nanjing, and the area in Tianjin. Large geographical barriers and long geographical distance significantly blocked gene flow between populations of S. bungii. Our results provide a basic theoretical basis for subsequent studies of millipede taxonomy and population genetic evolution.


Introduction
Soil macrofauna is an important part of biodiversity and include diverse macroinvertebrates living in the soil. They play a key role in regulating the physical, chemical, and microbiological properties of soils [1]. Soil macrofauna contribute to forest ecosystem sustainability by increasing organic matter dynamics and altering soil physical properties [2]. Due to these substantial impacts, soil macrofauna are described as soil engineers. Most studies on soil macrofauna have mainly focused on ecological functions, with only a few focusing on biological characteristics [3,4], such as taxonomy and genetics. In addition to annelids, platyhelminths, mollusks, and arthropods are important soil macrofauna. However, most genetic studies on soil macrofauna have been restricted to earthworms of the annelids. Furthermore, mitochondrial DNA is the most common molecular marker used in these studies [5,6]. In a previous study of soil macrofauna, mitochondrial DNA, nuclear DNA, and morphological techniques were considered as the methods of phylogenetic analysis [7][8][9][10].
As common soil macro-invertebrates, millipedes of the arthropods are one of the oldest terrestrial organisms and are widely distributed, with the earliest millipede fossils dating back to the late Silurian [11]. Millipedes contribute to soil processes, maintaining the soil fertility of mountain forests [12][13][14]. Millipedes are soil detritivores with a poor dispersal ability and are among the most important plant litter handlers in temperate forests [15]. The fragmented litter processed by millipedes is more easily degraded by other decomposers. Millipedes are used to monitor soil pollution owing to their important ecological functions [16,17]. In previous research, Jackson et al. used six genes (four mitochondrial and two nuclear) to present a phylogeny of the millipede family Xystodesmidae [9]. Piyatida et al. used two partial mitochondrial gene fragments (COI and 16S rRNA) and morphology to define four new Apeuthes species from Southeast Asia [10]. Spirobolus bungii is a common millipede species in the mountain forests of China. In this study, mitochondrial DNA was used as a marker to study the phylogeny and genetic diversity of S. bungii, to reveal the characteristics of the genetic diversity of millipedes. The 18S rRNA gene was used for species identification, to ensure that all the individuals studied were S. bungii.
This experimental study aimed to determine the population differentiation of S. bungii in eastern and northern China and to reveal the genetic structure and diversity of millipedes and their influencing factors. To study the genetic diversity of S. bungii, we needed to sample from two cities. First, considering the possibility of mass sampling, we selected Nanjing, East China. Second, S. bungii was found in the Yan Mountains of China [18]. Pan Mountain (Tianjin, China) is an offset of the Yan Mountains [19]; therefore, we chose Tianjin in North China as the second city for our study. Nanjing and Tianjin are developed cities with distinct differences [20,21]. Nanjing is located in a hilly area with green spaces of different sizes within the city [22,23]. The Yangtze River has run through the urban area of Nanjing since ancient times. Tianjin's urban area is mostly on a plain, with only Pan Mountain at 864.4 m above sea level, located in the north of the city, with no main river. As a large river, the Yellow River runs between the two cities. Therefore we hypothesized that the seven populations of S. bungii distributed in different regions could be divided into three groups: the southern bank of the Yangtze River (Group A), the northern bank of the Yangtze River (Group B), and in Tianjin (Group C).

Sample Collection
A total of 166 S. bungii individuals were collected from seven sites: Zijin Mountain (32 • Figure 1, Table 1). At each site, three soil samples were collected, along with individual organisms, between April and June 2021. After collection, all the biological samples were identified as S. bungii through identification of their physiological characteristics [24]. A portion of leg tissue was extracted from each animal, after which all individuals were released, and the collected tissue samples were stored at −80 • C until follow-up.

Soil Detection
The soil pH and electrical conductivity (EC) were measured in a 1:5 soil:water tion using an electronic pH meter (pHs-25, Shanghai Leici) and a conductivity met spectively. The soil moisture content was measured after oven-drying at 105 °C for The soil organic matter (SOM) was determined using the Walkley-Black potassiu chromate oxidation method [25]. The total nitrogen (TN) and available phosphorus (extracted with 0.5 mol L −1 NaHCO3) were analyzed using the standard Kjeldahl an sen methods, respectively. The available potassium (AK) extracted with 1 m NH4OAc was determined using an atomic absorption spectrophotometer (TAS-990 Beijing Persee, China). Nitrate was extracted using 1 mol L −1 KCl via ultraviolet sp photometry (TU-1901, Beijing Persee, China) [26].

DNA Extraction and MtDNA Amplification and Sequencing
The genome DNA were extracted from the tissue samples using DNAiso Reagent (TAKARA, Beijing, China). According to the nuclear and mitochondrial sequences of S.

Soil Detection
The soil pH and electrical conductivity (EC) were measured in a 1:5 soil:water solution using an electronic pH meter (pHs-25, Shanghai Leici) and a conductivity meter, respectively. The soil moisture content was measured after oven-drying at 105 • C for 24 h. The soil organic matter (SOM) was determined using the Walkley-Black potassium dichromate oxidation method [25]. The total nitrogen (TN) and available phosphorus (AP) (extracted with 0.5 mol L −1 NaHCO3) were analyzed using the standard Kjeldahl and Olsen methods, respectively. The available potassium (AK) extracted with 1 mol L −1 NH4OAc was determined using an atomic absorption spectrophotometer (TAS-990 AFG, Beijing Persee, China). Nitrate was extracted using 1 mol L −1 KCl via ultraviolet spectrophotometry (TU-1901, Beijing Persee, China) [26].

Sequence Analysis
The COX2 and Cytb genes of each individual millipede were concatenated using Editseq and Seqman. DnaSP5.0 [27] was utilized to determine the indices of genetic diversity, such as haplotype diversity (Hd), nucleotide diversity (pi), and the mean number of pairwise differences (K). The pairwise fixation index (F ST ), gene flow (haploid number of migrant, N M ) values, and analyses of molecular variance (AMOVA) using Arlequin version 3.5 [28] were employed to quantify the genetic structure of the S. bungii population, which was then analyzed using Structure 2.3.4 software [29]. The set population was K 2-8, each K value was repeated three times, and the lengths of the burn-in periods and MCMC repeats were 200,000 and 200,000, respectively.
The results were uploaded to Structure Harvester [30], to obtain the best K value. To investigate the genetic evolution of different populations, based on the haplotypes and individuals of the population, maximum likelihood (ML) and Bayesian inference (BI) trees were constructed using PhyloSuite v1.2.2 [31]. ML analyses were performed using the TPM2u + F + G4 model in the IQ-tree. The BI analyses were performed using the HKY + F + G4 model in MrBayes 3.2.6, and run for 10 million generations, with tree sampling every 1000 generations and a burn-in of 25% trees. The phylogenetic relationships between the mitochondrial DNA (mtDNA) haplotypes of S. bungii were estimated from an unrooted statistical parsimony network using Popart v1.7 [32]. TBtools software was used to create a soil heat map [33].

Genetic Diversity of S. bungii
A total of 1182 bp of COX2 (534 bp) and Cytb (648 bp) gene sequences of mtDNA and 700 bp of the 18S rRNA gene sequence from 166 individuals were analyzed (GenBank accessions: OL449413-OL449447). Gene sequence analyses revealed that there were two haplotypes and one variable site in the 18S rRNA gene, and 28 haplotypes and 78 variable sites in the COX2 and Cytb genes of the mtDNA ( Table 1). As molecular identification may be more accurate than morphological identification, we conducted molecular identification of the 18S rRNA gene, and the results verified that all collected samples were S. bungii.
Although the 18S rRNA gene has often been used for phylogenetic analysis among different species [34][35][36], it has not been commonly used for individual intraspecies analysis [37]. The 18S rRNA gene sequence of all collected biological samples was relatively conservative; thus, it was not used as the basis for phylogenetic analysis in our experiment. Except for the PS population (COX2: pi = 0.00094, S = 5, Hd = 0.368; Cytb: Pi = 0.00077, S = 6, Hd = 0.447), all population data indicated that Cytb had a higher genetic diversity than COX2. In the concatenated sequence, the number of haplotypes ranged from two to nine per sampling site, with ZJ having the fewest (n = 2) and PS having the most (n =9); 12 haplotypes were found in the southern bank of the Yangtze River, seven haplotypes were found in the northern bank of the Yangtze River, and nine haplotypes were found in Tianjin, and the haplotypes in each group were unique haplotypes. Haplotype diversity (Hd) ranged from 0.153 (ZJ) to 0.653 (PS), and the nucleotide diversity (pi) ranged from 0.00013 (ZJ) to 0.0085 (PS).

Population Genetic Structure of S. bungii
Gene flow occurs because of contact between different populations of a species. We studied the gene flow between each population, to explore the genetic diversity and population structure of S. bungii and its influencing factors ( Figure 2, Table 2). The range of pairwise F ST values between populations varied from 0.02778 (TS and ZJ populations) to 0.99368 (ZJ and LS populations) for the mtDNA, Cytb, and COX2 genes. The results showed that on the south bank of the Yangtze River, the level of gene flow among populations was relatively high (1.04729-8.74928), and genetic differentiation was not obvious (0.02778-0.19271). showed that on the south bank of the Yangtze River, the level of gene flow among populations was relatively high (1.04729-8.74928), and genetic differentiation was not obvious (0.02778-0.19271).  The structure analyses revealed the population genetic structure of S. bungii. The results showed that K = 4 (the maximum value of delta K) was the most likely number for the S. bungii genetic clusters (Figure 3a). However, according to the actual sampling, K = 3 may have been the most likely number of S. bungii genetic clusters ( Figure 3b). As shown in Figure 3a, the red gene cluster was almost entirely composed of samples from the ZJ, QX, FS, and TS populations; the mazarine blue gene cluster was primarily from individuals in the LS, LW, and PS populations, and the wathet blue gene cluster was mainly from individuals in the PS populations, with very few individuals in the other populations. Structure analyses demonstrated the consistency of the population genetic structure of S. bungii within a small geographical scale.  The structure analyses revealed the population genetic structure of S. bungii. The results showed that K = 4 (the maximum value of delta K) was the most likely number for the S. bungii genetic clusters (Figure 3a). However, according to the actual sampling, K = 3 may have been the most likely number of S. bungii genetic clusters ( Figure 3b). As shown in Figure 3a Hierarchical AMOVA from the geographical regions (northern and southern banks of the Yangtze River) indicated that most of the genetic variation was between groups (98.74%), whereas little genetic variation was within populations (1.02%) or among populations within groups (0.24%) ( Table 3). These results imply that the Yangtze River is a significant geographical barrier, leading to differences in the genetic diversity and population structure of S. bungii. However, as a control with populations within a small geographical range, AMOVA from the geographical regions (Nanjing and Tianjin) identified only 32.14% of the variation as being present between groups, whereas 66.40% was among populations within groups and 1.46% was within populations (Table 3).  Hierarchical AMOVA from the geographical regions (northern and southern banks of the Yangtze River) indicated that most of the genetic variation was between groups (98.74%), whereas little genetic variation was within populations (1.02%) or among populations within groups (0.24%) ( Table 3). These results imply that the Yangtze River is a significant geographical barrier, leading to differences in the genetic diversity and population structure of S. bungii. However, as a control with populations within a small geographical range, AMOVA from the geographical regions (Nanjing and Tianjin) identified only 32.14% of the variation as being present between groups, whereas 66.40% was among populations within groups and 1.46% was within populations (Table 3).

Phylogenetic Analyses and the TCS Haplotype Network
The ML tree (Figure 4a) showed that the haplotypes were clearly divided into three branches. Based on these results, three groups were considered to exist among the experimental individuals. The uppermost branch was the group of populations distributed on the southern bank of the Yangtze River (Group A), the middle branch was the group of populations located on the bank of the Yangtze River (Group B), and the bottom branch was the group of populations spread across Tianjin (Group C). The BI tree (Figure 4b) showed that Group B was more closely related to Group C than to Group A. Thus, the division of S. bungii populations into three groups according to their geographical distribution was supported by the results of the phylogenetic tree. division of S. bungii populations into three groups according to their geographical distribution was supported by the results of the phylogenetic tree.
A dominant role of intrapopulation variations in intragroup variations was found in the analyses of these two geographical regions. This result indicated that two populations within the same geographical region (S. bungii populations distributed along the Yangtze River in the Nanjing areas) were more distant than those distributed in the Tianjin area. The influence of the Yangtze River on the genetics of S. bungii was more profound than that of the Yellow River.

Relationship between Environmental Factors and Haplotype Distribution
The physical and chemical properties of the soils in the different mountains of Nanjing differed to some extent ( Figure 5). The relationships between different haplotypes and the corresponding soil physical and chemical properties were reflected through color gradients and similarities in the soil heat map (Figure 6). Hap1 was the most common haplotype in Group A and was distributed across every population in this group (ZJ, QX, FS, and TS).
The corresponding soil physical and chemical properties represented the average values for this group. Hap28 was a unique haplotype of Ziijin Mountain (ZJ), and the corresponding soil samples had high levels of soil moisture, whereas the soil pH was below the average for Group A. Hap21, Hap22, Hap23, and Hap24 were unique haplotypes of Qixia Mountain (QX), where the available soil potassium was high. Hap2 and Hap3 were distributed across Fang Mountain (FS), and most of the physical and chemical soil properties, except for pH, were lower than the average for this group. Hap25, Hap26, and Hap27 were the normal haplotypes of Tang Mountain (TS), where the available soil phosphorus, potassium, and soil moisture showed extremely low levels.

Differences in Genetic Structure between Soil Macrofauna
Dispersal is an important ecological process, which is widely believed to play a role in determining the community structure of species, together with local factors [38,39]. There are two forms of dispersal in animals: namely, active dispersal and passive dispersal. Active dispersal has strict requirements for dispersers [40]. Individual size and the motor organs are generally thought to be related to dispersal ability [41]. Within a certain geographical range, the dispersal ability of a species affects its genetic structure and the gene flow of the population [42]. The effect of dispersal ability on genetic structure has been demonstrated in some soil animals. The limited dispersal ability of soil annelids

Differences in Genetic Structure between Soil Macrofauna
Dispersal is an important ecological process, which is widely believed to play a role in determining the community structure of species, together with local factors [38,39]. There are two forms of dispersal in animals: namely, active dispersal and passive dispersal. Active dispersal has strict requirements for dispersers [40]. Individual size and the motor organs are generally thought to be related to dispersal ability [41]. Within a certain geographical range, the dispersal ability of a species affects its genetic structure and the gene flow of the population [42]. The effect of dispersal ability on genetic structure has been demonstrated in some soil animals. The limited dispersal ability of soil annelids (such as earthworms) hinders the gene exchange between their populations. Their population genetic diversity is positively correlated with habitat size [7,40,43]. Some soil insects (such as beetles and ants) are excellent at dispersal. Their wings allow them to travel further for genetic exchange with other populations [44,45].
Millipedes, as soil macrofauna, are large and have numerous pairs of legs. In the present study, we observed higher levels of gene flow (1.04729-8.74928; Table 2) on the southern bank of the Yangtze River (Group A). According to the results, S. bungii shows lower levels of population genetic differentiation than some soil animals proven to have limited dispersal ability (such as earthworms) [46]. Therefore, we speculate that the dispersal ability of S. bungii in Nanjing is relatively strong and that the habitat patches connecting habitats serve as bridges between populations, making gene exchange possible ( Figure 2). On an evolutionary time scale, exchanges between populations reduce the level of population genetic differentiation.
We believe that S. bungii in Nanjing has a certain dispersal ability, which is, however, only one of the factors influencing the genetic diversity and structure of the population in the area. Local factors (natural and artificial habitat spatial heterogeneity) may also play a role [41]. Our current understanding of the relationship between the dispersal ability of soil animals, population differentiation, and genetic diversity is rather limited. Therefore, future research should focus on the relationship between the population genetic structure and diversification factors.

Gene Flow between S. bungii Populations Blocked by Geographical Barriers
S. bungii population analyses revealed a significant relationship between phylogeny and geography, demonstrating the presence of a phylogeographic structure. S. bungii exhibited a high level of genetic differentiation, with many unique haplotypes and without shared haplotypes, which showed that the S. bungii populations were differentiated from each other. The AMOVA results (Table 2) also corresponded to the results of the phylogenetic tree (ML and BI trees) (Figure 4a,b); namely, groups A and B were geographically close but genetically distant, and groups B and C were geographically far apart but genetically close. We speculate that the reason for this might be related to the geological timeframes of the formation of the Yangtze and Yellow rivers.
Studies on the basalt age in the Nanjing area have shown that the Yangtze river formed in the Miocene (~23 million years ago) [47][48][49]. It has been generally acknowledged that the Yellow river formed during the Pleistocene epoch (~1.15 million years ago) [50,51]. The results of the millipede phylogenetic analyses appear to be explained by the geological timeframes of the formation of the Yangtze and Yellow rivers. Recent studies have indicated that large geographical barriers have an obvious blocking effect on genetic diversity and gene exchange within populations [52][53][54]. At present, the isolation of soil fauna in the Yangtze and Yellow rivers has not been extensively described, but some studies have shown that the Yangtze and Yellow rivers block the gene exchange of some terrestrial animals [55,56].
Therefore, we considered that the obvious differences in the haplotypes of S. bungii might be intimately and directly related to the formation of the Yangtze and Yellow rivers. Their formation could have carved large populations into smaller ones, leading to neutral evolutionary processes such as genetic drift. Prior to the Miocene, there may have been no obvious genetic differentiation in the eastern coastal areas of China. Following the formation of the Yangtze River, the populations to the north (the predecessor species of Groups B and C) could not contact the populations to the south (the predecessor species of Group A). Over thousands of years of evolution, the genetic distance between the two populations increased, and genetic differentiation may have become evident. The formation of the Yellow River in the early Pleistocene might have hindered the gene exchange of S. bungii once again, and genetic differentiation may have appeared between the populations distributed to the north of the Yangtze River and in the Tianjin region. Due to the later divergence time, the genetic distance between Groups B and C was smaller than that between Groups A and B.

Environmental Factors on the Genetic Diversity of S. bungii
In the present study, we investigated the relationships between the environmental factors represented by the basic physical and chemical properties of the soil and the genetic diversity of S. bungii. Soil is the main living space that supports the activities of soil-dwelling animals [57]. Biological soil invertebrate communities influence the litter decomposition rate [58,59], and environmental factors lead to the differentiation of soil animal communities [60,61]. Thus, it is important to study the relationships between soil-residing animals, represented by S. bungii and the soil environment.
We speculated that several physical and chemical soil properties may influence the S. bungii genotypes. The results ( Figure 6) revealed some heterogeneity in the soils of the different mountain locations. Taking the samples from Group A as an example, there were clear differences in the soil properties of the four mountains. The soil available phosphorus content of Tang Mountain was significantly lower than the average level, whereas the total nitrogen content of Qixia Mountain was lower, and the soil moisture content of Zijin Mountain was significantly higher. However, the distribution of haplotypes was relatively concentrated, and a universally common haplotype was present in each population.
Although there was no evident correlation between haplotype distribution and soil properties, habitat conditions had an undeniable influence on the genetic diversity of S. bungii. Mitochondria have the advantages of a simple structure, low molecular weight, and maternal inheritance [62]; therefore, we selected gene fragments in the mitochondrial DNA of S. bungii to study their genetic diversity. Our results suggested that the relationship between soil properties and millipede mitogenomes is not clear; however, the effects of other habitat conditions on genetic diversity cannot be neglected. Therefore, in subsequent experiments, we aim to use further methods, such as simple sequence repeats, and more genes, such as nuclear genes, to assess the potential impacts of habitat conditions on the genetic diversity of S. bungii.

Conclusions
This study focused on the genetic diversity and population structure of S. bungii in Nanjing and Tianjin. The results revealed that the gene flow between populations of S. bungii increased the consistency of the genetic structure. However, large geographical barriers, such as between the Yangtze and Yellow rivers, significantly blocked the gene flow between S. bungii populations, leading to high levels of genetic differentiation. Moreover, the correlation between the soil physical and chemical properties and haplotype distribution of S. bungii is not obvious. Our study fills research gaps in the genetic diversity of S. bungii at different geographical scales and provides a theoretical basis for subsequent studies of millipede phylogenetic analyses, community composition, and genetic diversity.