Phylogeography and demographic history of Gyrodactylus konovalovi (Monogenoidea: Gyrodactylidae), an ectoparasite on the East Asia Amur minnow (Cyprinidae) in Central China

Abstract Gyrodactylus konovalovi is an ectoparasite on the Amur minnow (Rhynchocypris lagowskii) that is widely distributed in the cold fresh waters of East Asia. In the present study, the phylogeography and demographic history of G. konovalovi and the distribution of its host in the Qinling Mountains are examined. A total of 79 individual parasites was sequenced for a 528 bp region of the mitochondrial NADH dehydrogenase subunit 5 (ND5) gene, and 25 haplotypes were obtained. The substitution rate (dN/dS) was 0.068 and indicated purifying selection. Haplotype diversity (h) and nucleotide diversity (π) varied widely in the Qinling Mountains. Phylogenetic trees based on Bayesian inference (BI), maximum likelihood (ML), and maximum parsimony (MP) methods and network analysis revealed that all haplotypes were consistently well‐supported in three different lineages, indicating a significant geographic distribution pattern. There was a significant positive correlation between genetic differentiation (F st) and geographic distance. The results of mismatch distribution, neutrality test and Bayesian skyline plot analyses showed that whole populations underwent population contraction during the Pleistocene. Based on the molecular clock calibration, the most common ancestor was estimated to have emerged in the middle Pleistocene. Our study suggests for the first time that a clearly phylogeography of G. konovalovi was shaped by geological events and climate fluctuations, such as orogenesis, drainage capture changes, and vicariance, during the Pleistocene in the Qinling Mountains.


Gyrodactylus konovalovi
In this study, we first examined the phylogeography and demographic history of the G. konovalovi based on the mitochondrial

| Ethics statement
This study was approved by the Animal Care and Use Committee of Shaanxi Normal University. The host is not evaluated in IUCN red list status (https ://www.iucnr edlist.org). None of the species (fish or parasite) sampled is endangered or protected in China (Yue & Chen, 1998). Host sampling was permitted by the local level authority in scientific research.

| Sample collection
Specimens of R. lagowskii were sampled from 48 localities, which covered the most Qinling Mountains from May to October in 2016 and 2017. Fish were rapidly euthanized by a blow to the head and stored in 96% ethanol. Gyrodactylus species were collected from the skin and fins of unique host under a stereomicroscope in the laboratory, and one parasite specimen per host individual was used to avoid pseudoreplication. Subsequently, the Gyrodactylus specimens were examined microscopically, and species identification was performed based on the morphology of the opisthaptor (Ergens, 1976). A total of 79 individuals of G. konovalovi from eleven localities were identified (Figure 2 and Table 1). Finally, each parasite specimen and host specimen were individually stored in 96% ethanol at 4°C. Voucher specimens of the parasites and hosts were deposited in the Fish Disease Laboratory, College of Life Sciences, Shaanxi Normal University, Xi'an, China, 710062.

| DNA extraction, PCR amplification, and direct sequencing
Total genomic DNA of G. konovalovi was extracted from a single individual using the TIANamp Micro DNA kit (Tiangen Biotech, Beijing, China) following the manufacturer's protocol. The ND5 gene sequences of the 79 individuals of G. konovalovi from the eleven populations were amplified by polymerase chain reaction (PCR) amplification using a newly specific primers designed based on the congeneric species Gyrodactylus brachymystacis and Gyrodactylus parvae (Ye, Easy, King, Cone, & You, 2017), the forward primer gkND5-F (5'-ATTAGAAAGAGAGCAGTGT-3') and the reverse primer gkND5-R (5'-ATTTGTAGATGATAGCAAG-3') were used to amplify a segment F I G U R E 1 Gyrodactylus konovalovi Ergens, 1976 (a) Whole mount specimen. Milli-Q water. The following cycling conditions were applied: initial denaturation for 1 min at 93°C followed by 35 cycles of denaturation for 10 s at 92°C, annealing for 1.5 min at 50°C and extension for 2 min at 60°C with a final extension for 6 min at 72°C. All PCR fragments were initially purified with a PCR purification kit (BGI Biotech, Shenzhen, China), subsequently subjected to electrophoresis in a 1% agarose gel and finally sequenced with PCR forward primer with an ABI Prism ® 3,730 automated sequencer (Applied Biosystems, Foster City, USA).

| Population genetic diversity
A total of 79 ND5 gene sequences were visually inspected and manually edited using BioEdit v7.0.9.0 (Hall, 1999) and then aligned with MUSCLE (Edgar, 2004)  The molecular diversity indices for the number of haplotypes (H), haplotype diversity (h), and nucleotide diversity (π) were calculated in DnaSP v5.10.1 (Librado & Rozas, 2009). A substitution model for the haplotype dataset was determined using the Akaike Information Criterion (AIC) in jModeltest v2.2.10 (Posada, 2008). The TrN model of evolution with the gamma shape parameter and proportion of invariable sites (TrN + I + G) was selected for the analysis of molecular variances (AMOVA) and phylogenetic analysis.

| Phylogenetic and network analyses
Phylogenetic relationships among the mitochondrial ND5 haplotypes were reconstructed using Bayesian inference (BI), maximum likelihood (ML), and maximum parsimony (MP) methods. For the BI and ML analyses, the TrN + I + G model was the most appropriate model of nucleotide substitution. The congeneric species Gyrodactylus parvae You, Easy & Cone, 2008 was selected as the outgroup (Ye et al., 2017). MP analysis was implemented in PAUP 4.0b10a (Swofford, 2002). Heuristic searches with tree-bisectionreconnection were executed for 1,000 random addition replicates F I G U R E 2 Map of sampling localities for Gyrodactylus konovalovi populations. The locality codes are given in Table 1. The lineages labeled red, blue, and green represent for lineages a, b, and c, respectively TA B L E 1 Sampling information and haplotype diversity based on the ND5 gene for eleven populations of Gyrodactylus konovalovi  (2), Hap20 (7) 0.389 ± 0.164 0.00074 ± 0.00031 HDT Ningshan Co.
Hap7 (2) with all characters treated as unordered and equally weighted. ML analysis was conducted using RAxML7.2.8 (Stamatakis, Ludwig, & Meier, 2005), with bootstrap analysis performed with 1,000 replicates. BI analysis was performed using MrBayes3.1.2 (Ronquist & Huelsenbeck, 2003), and one set of four chains was allowed to run simultaneously for 8 million generations. The trees were sampled every 1,000 generations, with the first 25% being discarded as burn-in. Stationarity was considered to be reached when the average standard deviation of split frequencies was below 0.01. A haplotype network was then constructed using the median-joining (MJ) network approach with MP Calculation in Network 5.0.1.1 (Bandelt, Forster, & Rohl, 1999).

| Population genetic structure
Genetic differentiation among populations was assessed using were analyzed to test for isolation by distance (IBD) (Slatkin, 1993).
The strength and significance of the relationships between genetic differentiation and geographic distance were assessed using linear regression in Graph Pad Prism v 5.0 (www. graphpad. com).

| Divergence time estimation
Divergence time estimation based on the strict molecular clock model was performed in BEAST v1.6.1 (Drummond & Rambaut, 2007). The divergence times of lineages were estimated using the TN93 + I + G model and a Yule speciation tree prior because of the absence of fossil or geological data for calibration. Given that the generation time for gyrodactylids may be as short as 2 days , which is considerably shorter than that of other taxa of monogeneans and may result in more rapid divergence at mitochondrial loci. The intermediate mutation rate of 13% per million years ago (Mya) was adopted, as it has been employed for mitochondrial gene analysis in Gyrodactylidae species (Lumme et al., 2016;Pettersen et al., 2015). Analyses were performed for 8 million generations while sampling every 1,000th tree, and the first 10% of trees sampled were treated as burn-in.

| Demographic history
Three methods were used with the ND5 haplotype dataset to trace the demographic history of three lineages. The neutrality test between Tajima's D (Tajima, 1989) and Fu's Fs (Fu & Li, 1993) was used to test for neutral evolution in Arlequin v3.5.1.2 (Excoffier & Lischer, 2010), with significantly negative values indicating population expansion.
The mismatch distribution was calculated in Arlequin v3.5.1.2 to test for signals of demographic expansion with a smooth unimodal curve (Harpending, 1994), and the significance of the sum of squared deviations (SSD) and Harpending's raggedness index (Hri) evaluated out by bootstrap resampling (10,000 replicates). The beginning time of expansion (t) was calculated in accordance with a previous study (Lumme et al., 2016). Bayesian skyline plot (BSP) analysis was applied in BEAST v1.6.1 (Drummond & Rambaut, 2007) to describe demographic history by assessing the time variation of effective population size. This analysis was performed using the TN93 + I + G substitution model with a mutation rate of 13% per Mya and 8 million generations.

| RE SULTS
A total of 79 sequences from eleven geographic populations of G.
konovalovi was obtained ( Figure 2). The sequence alignment pro-  Network analysis revealed that all haplotypes were consistently well-supported by three different lineages (Figure 4)

| Genetic diversity
Total haplotype diversity (h) and genetic diversity (π) across the samples were 0.955 ± 0.007 and 0.04984 ± 0.00458, respectively, indicating high level of genetic diversity. The h was highest in the YP population, whereas π was highest in the TTG population. The h and π were lowest in the XYB population, except for three populations (HDT, LCH and LYG), and the relationship between genetic diversity   (Hewitt, 1996). The genetic diversity of the marine parasite Mazocraeoides gonialosae also displayed a decreasing tendency from south to north along the coasts of the South and East China Seas (Li et al., 2011). However, populations with high haplotype diversity and low nucleotide diversity of M. gonialosae and G. thymalli and a star-shaped haplotype network are indicative of a classic postglacial expansion after a period of low effective population size, with rapid population growth enhancing the retention of new mutations, accumulating haplotype diversity, but lacking enough time to accumulate nucleotide diversity (Grant & Bowen, 1998;Hewitt, 1996).
Gene flow was considerably higher among oviparous and viviparous monogenean parasites that were studied along the coasts of the South and East China Seas and the Glomma River in Norway (Li et al., 2011;Pettersen et al., 2015). Such restricted gene flow and most unique haplotypes might cause mutations to persist in isolated populations that have been separated for a long time. A small population might lead to a strong genetic drift, and the reduced genetic diversity due to bottleneck increased genetic drift might increase inbreeding and genetic homogeneity to decrease population fitness (Huyse et al., 2017;Mills, 2007).
The substitution rate (dN/dS) was 0.068, indicating strong purifying selection with genes of relatively low nucleotide diversity (π) against nonsynonymous substitutions in accordance with results from other gyrodactylids based on COI and ND5 genes, such as G. thymalli, G. salaris, G. teuchis, and G. arcuatus (Lumme et al., 2016;Meinilä, Kuusela, Ziętara, & Lumme, 2004;Pettersen et al., 2015). The low values of the substitution rate (dN/dS) noted among all the protein-coding genes below 0.1 demonstrate that each gene is subject to strong purifying selection, in which the biological functions of proteins are safeguarded against deleterious mutations (Huyse, Buchmann, & Littlewood, 2008); this feature is a typical of mitochondrial protein-coding genes (Rand & Kann, 1998). However, the dN/dS ratio of COI was greater than one (dN/dS > 1), which may indicate that this gene has been under positive selection from G. corydori in Brazil (Bueno-Silva et al., 2011).

| Phylogenetic and network analyses
The topologies of the phylogenetic trees and network analyses yielded similar results, both indicating well-supported different lineages, revealing obvious phylogeographical pattern in the Qinling Mountains.
The pattern clearly conformed to the geographic structure of the Qinling Mountains. The Qinling Mountains have played important roles in influencing the phylogeography of amphibian and fish species and forming obvious geographic structures (Bowen et al., 2016;Liu et al., 2015;Lumme et al., 2016;Meng, Li, & Qiao, 2014;Wang et al., 2012Wang et al., , 2013Yu et al., 2014), suggesting a historical vicariance pattern.
The orogeny of the Qinling Mountains may have initiated the divergence into independent lineages (Yan, Wang, Chang, Ji, & Zhou, 2010). Thus, barriers to gene flow produced by complex geological history appear to be responsible for driving the high level of species diversity in these mountains in China. These mountains represent a natural boundary between the north and the south, dividing the Chinese temperate and subtropical climatic zones (Ding et al., 2013) and resulting in differentiated terrestrial and freshwater fauna (Li, 1981;Zhang, 2011). The host R. lagowskii and other gyrodactylids also showed high genetic differentiation and strong geographical structure in East Asia, Europe and Brazil (Bueno-Silva et al., 2011;Hahn et al., 2015;Hassan et al., 2015;Huyse et al., 2017;Lumme et al., 2016;Min & Yang, 1986;Pettersen et al., 2015;Xue et al., 2017). The M. gonialosae results showed a high rate of gene flow between different populations and a lack of genetic structure because of the continuous distribution in marine environments (Li et al., 2011). Furthermore, there was a significant positive correlation between genetic differentiation (F st ) and geographic distance fit the IBD model in accordance with other studies (Meng et al., 2014;Pettersen et al., 2015;Yu et al., 2014), suggesting that the distribution of genetic variation is due to geographical separation, rather than natural selection. The population genetic differentiation within F. taihangnica and Odorrana schmackeri rejected the IBD model (Li et al., 2011Wang et al., 2013), indicating high gene flow and panmixia in these species. The remarkable spatial genetic differentiations among F. taihangnica lineages likely reflect the effects of historical isolations induced by tectonic changes and Pleistocene climatic changes . More distant sites within the hydrological network were more genetically differentiated in the mountains. The host of this ectoparasite primarily occurs in the clear cold freshwater from middle steam to upstream (Kang, Min, & Yang, 2000;Nishida et al., 2015). This specific ecological habitat results in small population sizes and restricts gene flow among populations (with each population possessing unique haplotypes).
Therefore, these populations most likely survive in isolated habitats and differentiate through genetic drift and selection.

| Divergence time estimation
A strict molecular clock was estimated for the haplotypes (Figure 3). The time since the most recent common ancestor of the F I G U R E 5 Plots of genetic differentiation estimates of F st against geographic distance (km) between populations within the ND5 dataset of Gyrodactylus konovalovi. The linear regression overlays the scatter plots (r = .7222, p < .001) whole ingroup was dated to 0.87 Mya. Lineage A diverged 0.59 Mya, and lineages B and C diverged 0.52 Mya as estimated using a mutation rate of 13% per Mya. All divergence times occurred during the middle Pleistocene. The divergence time of the populations during the middle Pleistocene are in accordance with those for the gyrodactylids G. corydori, G. teuchis, G. arcuatus, and G. gondae (Bueno-Silva et al., 2011;Hahn et al., 2015;Huyse et al., 2017;Lumme et al., 2016).
In contrast, the divergence time of G. thymall from the Glomma River in Norway was estimated as very recent several thousand years during the Holocene after the LGM (Pettersen et al., 2015). The orogeny of the Qinling Mountains may have initiated divergence into the independent lineages (Yan et al., 2010), while the estimated divergence times in the early to middle Pleistocene might correspond to the rapid uplift (by more than 500-1,500 m) of these mountains that was influenced by the Qinghai-Tibet Plateau movement (Zhang & Fang, 2012). The phylogeographic distribution of freshwater fishes was affected by climate fluctuation during the Pleistocene (Gao et al., 2012), and the phylogeographic distribution of R. oxycephalus was affected by geological events and Pliocene climate fluctuations (Yu et al., 2014).

| Demographic history
The neutrality test between Tajima's D and Fu's Fs values was positive and not significant (Table 4). However, the values of the SSD and Hri index of all lineages except lineage C and of the total population did not reject the hypothesis that sudden expansion and mismatch distribution were multimodal ( Figure 6). The neutrality test and mismatch distribution, rejected the hypothesis of population expansion, indicating that the populations were stable and did not undergo expansion. The values from the neutrality test for the lineages were both positive and nonsignificant, and the mismatch distribution was multimodal, suggesting a bottleneck effect. The BSP suggested that the effective population size increased slowly for lineages A and C and the total population after approximately 0.01 Mya during the Holocene; however, the total population underwent a steep decline after 0.25 Mya, corresponding to the Pleistocene glacial epoch ( Figure 7). Therefore, all lineages and total population were discordant in the demographic scenarios based on three methods. The results of the mismatch distribution, neutrality test and BSP showed that the total population underwent contraction during the end of  (Huang, 1982;Wang et al., 2013). The values from the neutrality test for the lineages were positive and nonsignificant, and the mismatch distribution was multimodal, suggesting a bottleneck effect which is a sharp reduction in the size due to climate change after the glacial epoch (Lande, 1988).
A slightly different form of bottlenecks can occur if a small group becomes geographically separated from the main population, such as through a founder events; such populations are often small and show increased sensitivity to genetic drift, an increased inbreeding, and relatively low genetic variation (Lynch, Conery, & Burger, 1995). Ancient population bottlenecks or founder effects may have reduced genetic diversity. Therefore, the climate fluctuation and tectonic changes that led to vicariance and drainage capture during the Pleistocene affected by the phylogeographic distribution of G. konovalovi in the Qinling Mountains.

| CON CLUS IONS
The results of this study showed that the climate fluctuation and tectonic changes led to formed vicariance and drainage capture during the Pleistocene, which promoted lineages divergence and influ- konovalovi and its host R. lagowskii at the population level and the magnifying glass hypothesis.

ACK N OWLED G M ENTS
The authors would like to thank Fangluan Gao (College of Plant Protection, Fujian Agriculture and Forestry University) for assistance with the analysis software. The authors would also like to thank American Journal Experts for editing the manuscript for English language. This work was supported by the National Natural

CO N FLI C T O F I NTE R E S T
None Declared.

AUTH O R CO NTR I B UTI O N S
You P designed the research; Chen T, Chen J, Tang L, Chen X, Yan J, and You P performed the research; Chen T analyzed the data; Chen T and You P wrote the paper.