Population Analysis Based on Mito-nuclear Sequences: Implication of Dugesia japonica Decline in Taihang Mountains, China

Freshwater planarian Dugesia japonica is widespread in East Asia and has been used as a model animal in many research fields. Taihang Mountains is the demarcation line and indispensable treasure house of natural resources and ecological barrier in North China. In order to explore the genetic diversity, genetic structure and demographic history of D. japonica in Taihang Mountains, the mitochondrial Cytb and nuclear ITS-1 sequences of 120 individuals from 20 populations were analyzed. The results showed that 46 haplotypes were detected in total, including 10 shared haplotypes and 36 private ones. The overall haplotype diversity (Hd) and nucleotide diversity (π) were 0.961 and 0.00157, respectively. The phylogenetic relationship and haplotype network analysis indicated that there was no obvious correspondence between the haplotypes and the geographical locations, revealing no obvious phylogeographic pattern. Analysis of molecular variance (AMOVA) suggested that the genetic differentiation of D. japonica in Taihang Mountains was very great ( F ST =0.405, P <0.01). This significant genetic variation mainly occurred within populations (59.52%), followed by 36.54% deriving from among populations. As for demographic analysis, positive Tajima’s D and Fu’s Fs values (0.250 and 1.657, respectively) were found in concatenated Cytb/ITS-1 sequences as a whole for neutrality test. Moreover, the mismatch distribution was multimodal. Therefore, D. japonica in Taihang Mountains might have been undergoing a population decline, and conservation measures such as reducing human interference should be taken. We hope these findings will arouse conservation and management strategy regarding freshwater planarians and contribute to the biodiversity in the long run.


INTRODUCTION
F reshwater planarians which belong to phylum Platyhelminthes, class Turbellaria, order Tricladida are a group of bilaterian, dorsoventrally flattened, free-living invertebrate animals (Cebrià et al., 2018;Zhang et al., 2019a). Their striking regenerative ability has attracted the interest of scientists for centuries, and therefore led to their use as a classical model system of regeneration (Reddien and Alvarado, 2004). Meanwhile, due to their high susceptibility to xenobiotics, ease of acquisition and low maintenance cost, freshwater planarian has also been ideal test animal in toxicological research (Wu and Li, 2018). In addition, with the development of molecular marker and next generation sequencing technique, planarians with wide distribution and low dispersal capability has gradually become ideal model for exploring population genetics and phylogeography (Zhang et al., 2019a). Investigations in this field will not only be helpful to clarify the genetic variation within a species and its interaction with the environment, thereby deepening the understanding of biological evolution process, but also provide theoretical basis for the conservation and development of species and contribute to formulate reasonable and efficient conservation measures (Avise, 2000). As a common species distributed in East Asia, Dugesia japonica is ubiquitous in most provinces of China and has large population size. At present, D. japonica as a model species has been mainly used to study regenerative biology and toxicology (Pagán, 2017;Reddien, 2018;Zhang et al., 2018Zhang et al., , 2019bIreland et al., 2020). However, to our knowledge no work on the population genetics and phylogeography has been reported so far.
Taihang Mountains is located between the Loess Plateau and the North China Plain, which includes Beijing Municipality, Hebei, Shanxi and Henan provinces for an overall area of 136,500 km 2 . It extends from the northeast to southwest with the elevation ranging from sea level to nearly 3000 m (Geng et al., 2019). Due to the wide range of elevation and the typical earth-rocky features, Taihang Mountains possess highly heterogeneous features in O n l i n e F i r s t A r t i c l e topography, soil, climate and vegetation (Zhao et al., 2017;Fu et al., 2018;Geng et al., 2018). As the demarcation line of the second and third step of Chinese topography in north China, Taihang Mountains plays a vital role in ecological service, water conservation and climate regulation and is indispensable natural resources treasure house and ecological environment barrier (Hu et al., 2019;Zhang et al., 2019b). Molecular markers based on the variation of DNA molecules can directly reflect genetic polymorphism at the DNA level (Sunnucks, 2000). Due to possessing the features such as strictly maternal inheritance, faster evolution rate and lack of genetic recombination, mitochondrial DNA has been widely used for studying phylogeography and population genetics (Pakendorf and Stoneking, 2005). The commonly used mitochondrial cytochrome b gene (Cytb) is an appropriate marker for exploring the genetic diversity and genetic differentiation of populations (Cantatore et al., 1994). Nuclear DNA encompasses a large amount of genes and can theoretically meet the needs of phylogeographic research. Among them, the ribosomal internal transcribed spacer region (ITS) usually shows higher polymorphism and percentage of informative sites in studying intraspecies and closer interspecies relationships. At present, studies on population genetics and phylogeography of planarians using molecular markers Cytb and/or ITS have already been reported (Álvarez-presas et al., 2011;Lazaro et al., 2011;Solà et al., 2013;Riesgo et al., 2017).
In view of the unique geographical position and special topography of Taihang Mountains, we hypothesize that D. japonica in Taihang Mountains should have great  genetic differentiation and show a certain geographical pattern. In addition, with the development of tourism in Taihang mountain area, we wonder to know whether the human activity would affect the population structure and population dynamics of D. japonica. To this end, the genetic diversity, genetic structure, and demographic history of D. japonica in Taihang Mountains based on mitochondrial Cytb and nuclear ITS-1 were preliminarily analyzed in this paper.

Sample collection
Planarians were collected in rivers, streams or springs across Taihang mountainous area from 2018 to 2019. For each population some samples were fixed and preserved in absolute ethanol for molecular analysis, while some sexually mature individuals were killed with 1% nitric acid and then fixed in Bouin's fluid for morphological diagnosis (because species identification for planarians is mostly based on the anatomical characters of genital organs). If there were no sexually mature adults in a certain population, juveniles were brought indoors and reared to adulthood to guarantee reliable species identification. Only the samples identified as D. japonica were used in this study. Figure 1 and Table I show the geographical distribution and the detailed collection information of  D. japonica samples used in the current study, respectively. Voucher specimens were preserved in absolute ethanol and deposited in the Biological Museum of Henan Normal University (Xinxiang, China).

DNA extraction, PCR and sequencing
Total DNA was extracted from single planarian preserved in absolute ethanol by traditional phenolchloroform method (Sambrook and Russel, 2016) or the DNeasy Blood and Tissue Kit (Qiagen, Germany) following the manufacturer's protocol. After quality checking and concentration measuring, the obtained DNA was stored at 4 o C for immediate use or -2 o C for later use.
A 317 bp fragment of Cytb was amplified by polymerase chain reaction (PCR) using the primer pair CytbS (5'-ATGTCTCTTTGAGGAGCTACTGT-3') and CytbA (5'-CTAAAAAATACCACTCAGGCTTTAT-3') designed by ourselves. As for the nuclear gene ITS-1, a fragment of 669 bp was amplified using the published primers ITS9F (5'-GTAGGTGAACCTGCGGAAGG-3') and ITSR (5'-TGCGTTCAAATTGTCAATGATC-3') (Álvarez-presas et al., 2011). The PCR amplification conditions for Cytb were as follows: A total volume of 30 μL containing 15 μL TaqMix (dNTPs and Taq DNA polymerase plus, Vazyme Biotech, Nanjing, China), 12 μL ddH 2 O, 1μL 10 μM forward and reverse primers (synthesized by Sangon Biotech, Shanghai, China) and 1μL DNA template reacted in a GeneAmp PCR System 9700 (Applied Biosystems, Foster City, USA). The thermocycling procedure consisted of pre-denaturation at 94 o C for 5 min, 34 cycles of denaturation at 94 o C for 40 s, annealing at 45 o C for 30 s and extension at 72 o C for 1 min, followed by a final extension of 10 min at 72 o C. As for ITS-1, the annealing temperature was 45 o C, and other parameters were the same as Cytb. PCR products were confirmed as the target band by electrophoresis and then sent to Sangon Biotech (Shanghai, China) for sequencing with the corresponding amplifying primers.

DNA sequence analysis
Cytb and ITS-1 chromatograms including sense and antisense were checked manually and assembled using Seqman in DNASTAR Lasergene v7.1 to obtain single consensus sequence. For the protein-encoding gene Cytb, the sequence accuracy was confirmed by translating into the amino acids via the online software EMBOSS Transeq (https://www.ebi.ac.uk/Tools/st/emboss_transeq/). MEGA 7 (Tamura et al., 2011) was used to perform multiple sequence alignment. In order to avoid missing data, all sequences of Cytb and ITS-1 were truncated to the same length of 317 bp and 669 bp, respectively, and the trimmed sequences of two genes were then manually concatenated into one data matrix (Cytb/ITS-1) for further analyses. All sequences of Cytb and ITS-1 were submitted to GenBank and got the accession numbers.
In order to analyze the genetic diversity, number of haplotypes (N), haplotype diversity (Hd) and nucleotide diversity (π) were computed with the program DnaSP 5.0 (Librado and Rozas, 2009). Demographic history was explored by two different approaches in DnaSP 5.0. Firstly, Tajima's D and Fu's Fs statistics were used to test for neutrality. Then, the mismatch distribution to understand the demographic history was plotted. Using Arlequin 3.0 (Excoffier et al., 2005), an analysis of molecular variation (AMOVA) was conducted for 20 geographical populations which were divided into two groups (Sourthern and Northern), and pairwise F ST values to assess the genetic differentiation among populations were also quantified.
Phylogenetic topologies were reconstructed using Bayesian inference (BI) and maximum likelihood (ML) methods. The best-fit nucleotide substitution model TVM+I+G for the concatenated dataset was selected under Akaike Information Criterion (AIC) in the program jModeltest v2.1.7 (Darriba et al., 2012). Bayesian analysis was performed via MrBayes v.3.2.1. (Ronquist et al., 2012) based on the optimal model. Four Markov Chains (three heated and one cold) were run for ten million Markov chain Monte Carlo generations, sampling the chains every 500 generations. To verify the results, two independent runs were conducted. The initial 10% trees were discarded as burn-in samples. Remaining trees were used to generate a majority-rule consensus tree and Bayesian posterior probabilities (PP) (Huelsenbeck et al., 2001). PhyML 3.0 (Guindon and Gascuel, 2003) was used to perform ML analyses. We used the pruning and O n l i n e

F i r s t A r t i c l e
regrafting option to estimate better ML tree topology, and the non-parametric bootstrap with 1000 replicates to estimate the support for each internal branch of the ML phylogeny. Meanwhile, a haplotype network diagram was constructed using median-joining algorithm implemented in the Network 10.0 (Posada and Crandall, 2001).

Population genetic structure
The combined Cytb/ITS-1 sequences of 20 D. japonica populations sampling across Taihang Mountains were subjected to AMOVA analysis. The results suggested that the genetic differentiation index (F ST ) of D. japonica in Taihang Mountains was 0.405 with the statistical difference being extremely significant (P<0.01). The major molecular variation (59.52%) occurred within populations, followed by 36.54% deriving from among populations; whereas variation among groups accounted for only 3.94% (Table III). The pairwise F ST values computed from Cytb/ ITS-1 data were within the scope of -0.293 to 0.925, and the statistical differentiation reached significance (P<0.05) in 121 population pairs out of the total 190. The highest genetic variation was found between population HG and PShan (F ST =0.925, P<0.05), followed by that between HG and Yu (F ST =0.859, P<0.05) (Table IV).

O n l i n e F i r s t A r t i c l e
Population Analysis Based on Mito-nuclear Sequences

Phylogeny and haplotype network
Using D. ryukyuensis (GenBank Accession number: AB618488 for Cytb and FJ646909 for ITS-1), the sister species of D. japonica as the outgroup, phylogenetic trees based on Bayesian inference (BI) and maximum-likelihood (ML) methods were reconstructed (Fig. 2). Both trees presented the similar topology, and all the haplotypes were clustered into three haplogroups. Most of the haplotypes clustered into Clade I. Clade II included four haplotypes (H6, H37, H38 and H45) from LaiS, Qi, LQ, LC and ZQ, and clade III consisted of only two haplotypes (H36 and H46) from WT and LC, respectively. Meanwhile, the median-joining network of all haplotypes revealed the genetic structure congruent with the phylogenetic analysis ( Fig. 3). Three haplogroups including the same haplotypes as those in the phylogenetic tree were obviously clustered.

Demographic history
Neutrality tests were conducted using Tajima's D and Fu's Fs statistics, and similar results were obtained for 20 geographical populations as a whole or separately. When the 20 populations were considered as one group, the Tajima's D and Fu's Fs test values were 0.250 and 1.657, respectively. As for the single population, both Tajima's D and Fu's Fs results were found positive in population QY, PShan, XT1, XT2, XW, JY1, LQ, HG and YC (Table  II). Furthermore, the mismatch distribution over all populations was not distinctly unimodal but multimodal (Fig. 4).

DISCUSSION
Genetic diversity of species is the product of longterm evolution of organisms, and the level of genetic diversity is closely related to the viability and evolution potential of species (Humphries et al., 1995). Research on the level, formation mechanism and distribution pattern of species genetic diversity not only can be helpful to understand the evolutionary history of species, but also can provide an important basis for analyzing the evolutionary potential of species and predicting their development trends. Haplotype diversity (Hd) and nucleotide diversity (π) are two important indicators to measure genetic diversity (Féral, 2002). As far as D. japonica in Taihang Mountains are concerned, the overall haplotype diversity was 0.961 and the nucleotide diversity was 0.00157, which shows the high Hd and low π pattern typical for migratory species. Low nucleotide diversity together with high haplotype diversity might be indicative of genetic O n l i n e

F i r s t A r t i c l e
Population Analysis Based on Mito-nuclear Sequences bottleneck events (Wei et al., 2013;Sun et al., 2015). Both the phylogenetic tree and the haplotype network diagram showed the similar results. The haplotype topology did not display the distance relationship corresponding to the geographic location of the population, which shows there is no obvious pedigree geographic pattern. Studies on the population genetic structure not only can evaluate the variation level of populations and the relationship between different geographical groups, but also can determine the evolutionary significant units and management units so as to formulate the conservation and management strategies for resources (O'brien, 1994). Population genetic differentiation index (F ST ) is an important indicator to reflect the degree of genetic variation among populations (Allendorf, 1983). The degree of genetic differentiation is positively related to the F ST value, and it has been suggested that a value between 0.15 and 0.25 means great differentiation and values above 0.25 very great genetic differentiation (Wright, 1978;Balloux and Lugon-Moulin, 2002). In the current study, the F ST value drawn from AMOVA analysis was 0.405 (P<0.01) indicating very great genetic differentiation of overall D. japonica in Taihang Mountains. Consistently, most of the pariwise F ST values showed above 0.15 (122 out of the total 190), which suggests great variation between colonies. This result also sheds light on the reason behind the haplotype composition pattern of D. japonica in Taihang Mountains (less shared haplotypes and more private ones). Geographical isolation would separate populations, impede gene flow, hence generate genetic differentiation (Pyron and Burbrink, 2010;Liu et al., 2018). Since Cenozoic, Taihang Mountains have experienced multiple stages of geological movements and formed highly heterogeneous characteristics with a large variety of topography, soil, climate and vegetation (Geng et al., 2019). The unique and complex terrain may hinder genetic communication between populations, thereby promote the high differentiation of population genes (Liu et al., 2014). In addition, species with strong dispersal potential usually exhibits lower genetic variance between its different geographical populations. The limnetic D. japonica as a part of the benthos are slow-moving worms, which renders their long distance dispersal unlikely (Knakievicz, 2014). Therefore, the weak dispersal ability would be another reason responsible for the great genetic differentiation of D. japonica in Taihang Mountains.
Positive Tajima's D values indicate that populations have experienced a contraction while negative Tajima's D values occur in colonies that are undergoing demographic expansion or mutational selection (Tajima, 1989). When the 20 populations in Taihang Mountains were considered as one group, neutrality test results including Tajima's D and Fu's Fs were both positive for the combined Cytb/ITS-1sequences. Furthermore, the mismatch distribution over all populations was deemed multimodal and rejected the hypothesis of a sudden expansion (Rogers and Harpending, 1992). As a result, we can speculate that D. japonica in Taihang Mountains should have been experiencing population depression, even though the decline might be not very notable. This has also been verified by the analysis based on mitochondrial COI (unpublished data). Factors such as regional climate types, topography, natural environment and human activities will affect species diversity and population richness (Meng et al., 2017). D. japonica as a freshwater invertebrate mainly inhabits in the limnetic environment such as rivers, streams and springs. Taihang Mountains is the demarcation line in North China and has become one of the important biodiversity centers in China (Geng et al., 2019). But recently, with the rapid development of mountain tourism, the aquatic environment in which freshwater planarians inhabit has been polluted or at least interfered by human activity. Additionally, water drying up in the context of global warming might be another direct reason for the population degradation of D. japonica in Taihang Mountains.

CONCLUSIONS
In this study, the population genetic diversity, genetic structure, phylogeography and population dynamics of freshwater D. japonica in Taihang Mountains were investigated. On the whole, the genetic diversity analysis presents a high Hd and low π pattern typical for migratory species. Isolation effect exerted by the special topography in Taihang Mountains and weak dispersal ability of benthonic freshwater planarians are probably main reasons for the very great genetic differentiation of D. japonica in Taihang Mountains. The positive neutrality test results as well as the multimodal arrangement of mismatch distribution indicate that D. japonica in Taihang Mountains might have been undergoing population decline. Therefore, practical measures such as reducing human interference should be taken to strengthen the conservation of D. japonica in Taihang Mountains. We hope these findings will arouse conservation and management strategy regarding freshwater planarians and contribute to the biodiversity in the long run.

ACKNOWLEDGMENTS
This work was supported by the National Natural Science Foundation of China (No. 31702010).

Statement of conflict of interest
The authors have declared no conflict of interests.