Genetic structure and population dynamics of the silver pheasant (Lophura nycthemera) in southern China

Chaoying ZHU, Wenjing LU, Guanglong SUN, Jiang CHEN, Zhen ZHANG, Liangjie YU, Peng CHEN, Yuqing HAN, Zhifeng XU, Junpeng BAI, Dongqin ZHAO, Luzhang RUAN* School of Life Sciences, State Ministry of Education Key Laboratory of Poyang Lake Environment and Resource Utilization, Nanchang University, Nanchang, P.R. China Poyang Lake Construction Office, Water Resources Department of Jiangxi Province, Nanchang, Jiangxi Province, P.R. China Nanjing Institute of Environmental Sciences, Ministry of Ecology and Environment, Nanjing, P.R. China Key Laboratory of Animal Resistance of Shandong Province, College of Life Sciences, Shandong Normal University, Jinan, Shandong Province, P.R. China


Introduction
Pleistocene geologic and climatic events seriously affected the genetic patterns and phylogeographic structure of organisms (Cracraft, 1973;Hewitt CD, 2000;Wu et al., 2012;Wang et al., 2017), resulting in high diversity in tropical and subtropical areas of southern China, which was an important biodiversity site due to sufficient suitable habitats and abundant food resources during the glacial periods (Alain et al., 1998;Crowe et al., 2006;Huang et al., 2010;Cao et al., 2012). In these regions, birds could survive in different refugia and then may have differentiated into new lineages (Hewitt G, 2000) where gene flow could be prevented by some types of geographic barriers (Aleixo, 2004;Anthony et al., 2007;Nicolas et al., 2011).
The Pleistocene seriously affected the structure of most birds (Hewitt CD, 2000;Greer, 2013). Nevertheless, quite a few studies have discovered that the populations during the interglacial period were stable (Marín et al., 2013;Zhu et al., 2018), while other studies have shown that gene flow and genetic diversity are vital for population existence and evolution (Clement and Crandall, 2000;Trask et al., 2017).
High gene flow and population expansion can promote genetic diversity (Song and Lei, 2014;Ruan et al., 2018;Xu et al., 2019). However, low genetic diversity can lead to a population bottlenecks (Excoffier et al., 2009). Sexbiased dispersal can improve breeding opportunities for birds, thereby increasing gene flow, genetic diversity, and inclusive fitness (Clinton et al., 2007;Hamilton and May, 1977;Taylor, 1988;Pernetta et al., 2011), which are significantly beneficial to the evolutionary process of a population (Paris et al., 2016).
The silver pheasant (Lophura nycthemera) is a common species that is widely distributed in southern China and Southeast Asia, which prefers to inhabit shrub patches (Zhang et al., 2003). Based on mitochondrial markers, previous researchers studied the divergence time and evolutionary history of Lophura (Alain et al., 1998;Huang et al., 2010;Jiang et al., 2014). However, there are no previous studies of the sex-biased dispersal patterns or population structures of the silver pheasant in southern China. Considering previous results, we adopted mitochondrial and microsatellite approaches to explore the pattern of sex-biased dispersal of the silver pheasant and its influencing factors and to determine whether and how the interglacial period affected the population structure, gene flow, genetic diversity, and dynamics of the silver pheasant in southern China.

Sample collection
From 2010 to 2011, we collected 72 adult blood samples and 43 fresh feather samples from silver pheasants in 18 localities in southern China (101.72°E to 119.65°E and 25.86°N to 32.92°N). The samples were classified into seven populations, namely SC, AH, JX, FJ, ZJ, HB, and HN, based on the Euclidean distance and obvious barriers such as rivers and mountain regions near the sampling sites ( Figure  1; Table 1). We collected blood samples (approximately 0.1 mL) into heparinized vessels with a 26-gauge needle. All of the samples were stored in a refrigerator at -20 °C. For microsatellite data, the population SC was removed from gene flow, bottleneck, and pairwise difference analyses due to small sample size (n = 1; Tables 2 and 3).

DNA extraction and PCR amplification
We used the methods described by Hermet et al. (2013) and Luo and Xu (2004) to extract genomic DNA. PCR with primer pairs (Excoffier and Lischer, 2010) PHDH (5′-CATCTTGGCATCTTCAGTGCC-3′) and PHDL (5′-AGGACTACGGCTTGAAAAGC-3′) was conducted to amplify the partial sequence of the mitochondrial control region (D-loop) in the samples. The amplified volume was 15 µL. PCR amplification was performed as follows: denaturation at 95 °C for 4 min and then denaturation at 94 °C for 35 cycles, annealing of mitochondrial DNA at 53 °C for 30 s, elongation at 72 °C for 1 min, and elongation at 72 °C for 10 min as the final step. ABI Prism Large Dye Terminator Cycle sequencing chemistry was applied with the ABI 3730 automated DNA sequencer using an international sequencing technique for mitochondrial PCR products. Then we used SeqMan 8.0.2 (DNASTAR Inc., Madison, WI, USA) to calibrate, edit, and trim the chromatogram to a final length of 1053 bp. Nuclear microsatellite genotypes of 82 individuals were amplified by singleplex fluorescent PCR with primer pairs SR01, SR08, SR13 (Wang et al., 2009), SE02 (Jiang et al., 2006), PC125 (Hale et al., 2004), MCW98 (Crooijmans et al., 1996), and MCW146 (Crooijmans et al., 1997). The thermal cycling parameters were denaturation for 2 min at 94 °C and then 30 cycles for 30 s at 94 °C, annealing for 30 s at 52 °C, and elongation for 1 min at 72 °C. Each reaction had a 10-min elongation as the final step. In accordance with the Internal Lane Standard GeneScan 500 using ABI GeneScan Analysis [Rox], we also analyzed fragment sizes. Then the allele sizes were detected through ABI GeneMapper analysis. Sixteen samples (20%) were amplified twice to check the accuracy of results. We compared the complete mitochondrial sequences of this species from GenBank (Accession no. NC_012895.1) with our mtDNA sequence data to identify nucleotide differences in the control region.

Genetic differentiation and diversity
We used CLUSTALX 1.83 to align mitochondrial DNA sequences (Thompson et al., 1997) and DNASP 4.1 (Rozas et al., 2003) to calculate haplotype (H) and nucleotide diversities (π), while allelic richness (A) and expected (H E ) and observed heterozygosities (H O ) were calculated using ARLEQUIN 3.5 (Excoffier and Lischer, 2010). HP-Rare 1.1 (Kalinowski, 2005) was used to compute the rarefied allelic richness (A P ). Statistical deviations were examined using ARLEQUIN 3.5 from Hardy-Weinberg and genotypic linkage equilibrium. We also used ARLEQUIN 3.5 to analyze the genetic structure of the microsatellite and mitochondrial DNA data through analysis of molecular variance (AMOVA). In addition, we used Fand Φ-statistics to evaluate the population differences (Excoffier and Lischer, 2010). Genetic bottleneck tests via BOTTLENECK 1.2.02 (Cornuet and Luikart, 1996) were applied to examine the signatures of recent population bottlenecks. Each population was analyzed using twophase models (TPMs). Genetic bottleneck analysis with 10,000 repetitions was conducted to obtain P-values based on the Wilcoxon symbolic rank test, as recommended by Cornuet and Luikart (1996). It was expected that the rare alleles would rapidly be lost under the bottleneck conditions, consequently distorting the mode-shift (Luikart et al., 1998). Thus, we examined all population allele frequencies using the allele frequency distribution test. We also adopted the Garza-Williamson index (dividing the allelic range by the number of alleles; Garza and Williamson, 2001) as a measurement of the bottleneck experiment for all populations by using ARLEQUIN 3.5 (Excoffier and Lischer, 2010). A historical decline in sizes was defined when the M value was less than 0.68 (Alasaad et al., 2011).

Population expansion test
Using BEAST 1.6.1 (Drummond et al., 2003), we identified the nodes of the silver pheasant clades where the most recent common ancestor (TMRCA) lies. Assuming a strict molecular clock, we constructed a best substitution model (HKY+I+G) from MODELTEST 3.06 and among-site rate heterogeneity across all branches. We selected an uncorrelated lognormal model to explain the rate variation among all lineages. The fluctuations in population size were examined using the Bayesian skyline model. Independent Markov chain Monte Carlo (MCMC) analyses were run for 10 8 steps. The first 10% of data was ignored as burn-in, and then sampling was performed every 1000 generations. Demographic history time was  (Excoffier and Lischer, 2010), and mismatch distribution was determined by DNASP 4.1 (Rozas et al., 2003). The results of the neutrality test showed significantly negative values (P < 0.05) for Tajima's D and Fu's Fs tests, which were regarded as indicators of the recent population expansion. The raggedness index (RAG) was used to test the sum of square deviations (SSDs) for mismatch distribution analysis. The expansion time was estimated from the value of tau by using tau = 2 µkt (Rogers and Harpending, 1992), where tau, µ, and k are the mismatch distribution mode, the nucleotide mutation rate (1.19% per Myr) (Dong et al., 2013), and the number of nucleotides (k = 1053 bp), respectively.

Phylogeographic and population structure
Three mitochondrial D-loop sequences from Lophura hatinhensis (AJ300150), Lophura leucomelanos (AJ300153), and Lophura swinhoii (AJ300155) were treated as outgroups because they have a close phylogenetic relationship to the silver pheasant (Randi et al., 2001;Jiang et al., 2006). MODELTEST 3.06 was used to test the molecular evolution of nucleotide substitution for the mtDNA sequence data (Posada, 1998) with PAUP 4.0b10 (Swofford, 1998). We identified the best-fitted model with the Akaike information criterion (AIC; Posada, 1998). MrBayes 3.1.1 (Drummond et al., 2003) was used to construct the Bayesian trees in order to generate Bayesian posterior probabilities for phylogenetic deductions. ML trees were built with IQ-TREE (Trifinopoulos et al., 2016). We adopted two independent MCMC runs with three chains per run, totaling 3,000,000 generations after sampling every 100 steps. The convergence to a static distribution was tested by the splitting frequency between two independent runs (less than 0.01). TCS 1.21 (Clement and Crandall, 2000) was used to calculate genealogical relationships among haplotypes. A medianjoining network method (MJN; Bandelt et al., 1999) was also used to describe the relationship among the silver pheasant haplotypes. NETWORK 4.5 (http://www. fluxus-engineering.com) was used to compute the MJN. STRUCTURE version 2.3.3 was used to examine the patterns of population structure (Pritchard et al., 2000). We used twenty replicate runs to analyze K values from 1 to 6 with an initial 1 × 10 5 burn-in step with 2 × 10 6 analysis replicates following. STRUCTURE HARVESTER was used to evaluate the K value (Earl and Vonholdt, 2012).

Recent and historical gene flow test
The first-generation migration among populations was tested by GENECLASS 2.0 (Cornuet et al., 1999), and the long-term interisland gene flow rate was estimated by MIGRATE-N 3.2.6 (Beerli and Palczewski, 2010). MIGRATE-N 3.2.6 can estimate the mutation rate-scaled effective population size as θ = χN e µ, where χ, N e , and µ represent separate inheritance scaling factors (1 for mtDNA and 4 for microsatellite markers), mutation rate per generation per locus, and effective population size. The equation M = m/µ was used to test mutationscaled migration rates, where m represents the ratio of immigration per generation. We used a Brownian motion mutational model for microsatellite analyses. The θ value was bounded from 0 to 10 (delta = 0.1), while the M value was bounded from 0 to 30 (delta = 0.1). The static heating scheme was set including four chains (mathematical temperatures: 1.0, 1.5, 1.5, and 10 4 ). We repeated the entire analysis twice. The Mantel correlation coefficient (Mantel, 1967) was calculated by MANTEL 1.18. The estimates of genetic distance F ST /(1−F ST ) (Slatkin, 1995) were computed by ARLEQUIN 3.5 for both microsatellite and mitochondrial datasets (Excoffier and Lischer, 2010). The genetic distance measures of the pairwise populations were regarded as a linear correlation to pairwise logarithmic Euclidean distances (Fenster and Hardy, 2003), where 10,000 random permutations were carried out to calculate the significance of the correlation.

Genetic differentiation and diversity
After the sequential Bonferroni corrections, the Hardy-Weinberg test did not show any significant differences among microsatellite loci, while linkage disequilibrium tests indicated only one significant separate locus pair at the 0.05 level. This observation in population HN could be treated as occasional (Miller et al., 2012). Of all the values of expected heterozygosity, observed heterozygosity, and allelic richness, the highest values were detected in ZJ (H E = 0.717, H O = 0.917, A P = 1.900; Table 4). The lowest values of expected heterozygosity, allelic richness, and rarefied allelic richness were all detected in population HN (H E = 0.406, A = 2.000, A P = 1.450; Table 4). The lowest value of the observed heterozygosity was detected in population JX (H O = 0.549). There was no evidently consistent trend of genetic diversity among populations in the observations according to microsatellite data. Both the highest haplotype diversity and nucleotide diversities (H d = 0.914, π = 10.180 × 10 -3 ; Table 4) of mtDNA data were obtained in population SC. AMOVA showed that the main genetic variances were among the populations for mitochondrial data (56.63%, F ST = 0.566, P = 0.000; Table 5), whereas the main genetic variances were within populations for microsatellite data (95.65%, F ST = -0.293, P = 0.000). Pairwise F ST analysis indicated that there was relatively low genetic difference among populations. The lowest F ST estimates separately belonged to population pair HN/ZJ for mitochondrial DNA and population pair HB/ZJ for microsatellites (Table 2), whereas the highest F ST estimates belonged to population pair ZJ/FJ for mitochondrial DNA and population pair JX/HN for microsatellites. In analysis of genetic differences in both mtDNA and microsatellite datasets, populations JX, FJ, and HN were significantly different from the other populations.

Evidence for genetic bottlenecks
In TPM for bottleneck analysis, populations AH and ZJ revealed statistically significant results (P < 0.05). Even so, allelic shift tests suggested that only four populations of the seven, namely SC, JX, FJ, and HN, had an L-shaped distribution. A recent population bottleneck could explain the observed mode shift in HB, AH, and ZJ. The Garza-Williamson index for the HB, AH, FA, and ZJ populations was below 0.68 (Garza and Williamson, 2001), indicating that these populations underwent a historical population bottleneck (Williamson-Natesan, 2005).

Demographic history test
Tajima's D (Table 6) test was negative for all populations and significant for the whole population (P = 0.04). Fu's Fs test was negative for most populations except SC and AH and was extremely significantly negative for the whole population (P = 0.007). The tests of both RAG and SSD for mismatch distribution were all insignificant at the 0.05 level except for population SC, which implied that the curves fit the model test of population expansion well. In summary, both mismatch distribution analysis and the neutral test revealed that significant historical population expansion existed in the whole population. Calculated by t = tau/2µk (tau = 1.906, 95% confidence interval: 0.363-10.465), the whole population expansion could be evaluated to 76,053 years ago. The model (HKY+G+I) picked with the AIC was identified by MODELTEST to be the best for molecular evolution. The BSP analysis for the  whole population of China implied that the TMRCA was estimated to be approximately 970,300 ybp ( Figure 2).

Evidence for genetic structure and phylogeography
In a 1053-bp sequence for 115 individuals, we observed 33 variable sites, which defined 36 haplotypes. Most of the populations had only one private haplotype. Six shared haplotypes were observed: H3, H4, H19, H20, H24, and H29. Most haplotypes showed a weak geographic population structure, except those for population SC (Figure 3). The phylogenetic tree analysis showed no obvious geographical structure for most populations except population SC, which was consistent with the results from the MJN analysis ( Figure 4). Furthermore, STRUCTURE HARVESTER analysis showed that K = 2 or ΔK = 3 (log marginal likelihood = -1086.76, Figure 5) was the best clustering. Unfortunately, STRUCTURE analysis indicated that the whole population had no obvious geographical structure according to K = 2 or ΔK = 3.

Historical gene flow test
GENECLASS was applied to test the dispersal among populations. No individual was identified as a firstgeneration migrant (P < 0.01). Computed by MIGRATE-N 3.2.6, the gene flow analysis implied that the highest θ was in population FJ (0.891), while the highest M (2.930) was in neighbor populations JX to FJ. The lowest θ and M were observed in AH (0.524) and ZJ to JX (1.090), respectively (Table 3). In general, population JX was a main source of genetic migrants to the others, and the adjacent populations had apparently higher gene flows than the distant populations (Table 3). The Mantel test also indicated a positively insignificant correlation between the pairwise genetic distance and the Euclidean distance (r = 0.285, P = 0.8239 and r = 0.151, P = 0.2902 for mitochondrial and microsatellite, respectively).

Sex-biased dispersal test
Sex-biased dispersal was tested over the whole population, indicating that m AI (male m AI = -0.083, female m AI = 0.187, P = 0.170), F ST (male F ST = 0.043, female F ST = 0.066, P = 0.220), r (male r = 0.112, female r = 0.181, P = 0.080), and ν AI (male ν AI = 1.503, female ν AI = 0.851, P = 0.090) were not significant between females and males. In summary, all indexes in the sex-biased dispersal analyses showed no apparent sex-biased dispersal pattern in the silver pheasant.

High genetic diversity and gene flow lead to a weak population structure
Both the genetic diversity and gene flow of silver pheasant were high, whereas the population structure was comparatively weak. Both the mtDNA and microsatellite analyses in our investigation showed that genetic diversity was high in all the geographical populations (Table 4), possibly because of the wide distribution of the silver pheasant (Dong, 2011) and its habitat mainly distributed in warmer regions of southern China, a principal high-  level biodiversity area with abundant habitats (Crowe et al., 2006). The high genetic diversity in this area was also attributed to the long history (TMRCA = 970,300 ybp) of colonizing and evolving, and this observation was supported by the statistical analysis of historical demographic data (Table 6; Figure 2). Moreover, high gene flow and population expansion would promote the genetic diversity in Rallidae . Our results revealed high genetic exchange and low genetic differentiation among populations. AMOVA showed that the main genetic variances were among populations for mitochondrial data, whereas the main genetic variances were within populations for microsatellite data. Different evolutionary rates among markers (Brown and Swank, 1983) or male-biased gene flow (Gibbs et al., 2000;Miller et al., 2015) are both reasonable hypotheses, which could explain the differences in datasets. Pairwise F ST analysis indicated a slight genetic difference among populations. Birds usually have a higher gene flow than other organisms because of their ability to fly and move over longer distances (Miller et al., 2015). Nuclear genetic differentiation ascending with gene flow among populations explained that in adjacent populations, those with high gene flow always have low genetic differentiation . This study obtained consistent results with previous studies, which could be explained as high gene flow resulting in low genetic differentiation. Geographic features, such as potential river and basin barriers, could prevent gene flow in birds. Phylogenetic tree analyses suggested that no obvious geographical structure existed in the populations, except for population SC, which was in accordance with the results of Zhang et al. (2003). This phenomenon might be attributed to the unique topography of Sichuan, which is located north of the Yangtze River and is mainly dominated by basins, chiefly affected by rivers and basin barriers. STRUCTURE also supported a weak geographic structure. High gene flow and genetic diversity may correspond to a fragile population structure (Li et al., 2015;Ruan et al., 2018). However, there was no significant sex-biased dispersal pattern in the silver pheasant, indicating that sex-biased dispersal must not be an explanation for the weak population structure. Consequently, our results revealed that a weak population structure was mostly caused by high gene flow and genetic diversity.

Population dynamics during the interglacial period
In both Tajima's D and Fu's Fs tests, significantly negative findings were found in the whole population. Both the mismatch distribution analysis and the neutral test showed that the whole population had significant population expansion. The time of the whole population expansion could reach to 76,053 years ago, which indicated that the population in China expanded during the interglacial period, before the last glacial period. Several studies have stated that Galliformes may be susceptible to climate change with their low dispersal capability (e.g., Li et al., 2010;Wang et al., 2013), which likewise revealed that no sex-biased dispersal pattern occurred in the silver pheasant. Although a few local populations underwent a historical bottleneck, the whole population expanded in the last interglacial period, because the latest time that the global temperature was higher than the present is considered as the last interglacial period (Muhs et al., 2001). Therefore, the expansion population of the silver pheasant could be attributed to the abundant habitats and food resources in the last interglacial period. Fossils showed that the divergence time of Lophura was 38.9 Mya (Tuinen and Dyke, 2004), and the TMRCA of the entire population expansion occurred approximately 970,300 years ago (Figure 2), which was earlier than that of Lophura, indicating a long history for the current lineage for this species. Moreover, southern China has subtropical and tropical climates with complex topography and abundant habitats and food resources, suggesting that, typically, multiple refugia existed during the glacial periods (Song and Lei, 2014). Consequently, our study suggested that the long history of the current lineage of the silver pheasant might be related to the complex topography and enough refugia during the glacial periods in southern China. Briefly, the suitable climate in the last interglacial period provided rich habitats and food resources for the silver pheasant to survive, and the sophisticated topography with the refugia of southern China contributed to their reproduction and evolution.