Development and characterization of novel microsatellite markers for the Common Pheasant (Phasianus colchicus) using RAD-seq

The Common Pheasant (Phasianus colchicus) Linnaeus, 1758 is the most widespread pheasant in the world with a natural, geographic range spanning in temperate to subtropical regions of the Palearctic realm (Johnsgard 1999). This species exhibits a high-level of intra-specific differentiation in plumage coloration and patterns in males. Thirty subspecies forming five subspecies groups were defined mainly based on geographically distributed affinities and morphological characters (Cramp and Simmons 1980; Johnsgard 1999; Madge and McGowan 2002). The five subspecies groups are as follows: (1) the colchicus group (Black-necked Pheasants, west and south of the Caspian Sea, including the subspecies persicus, talischensis, colchicus and septentrionalis); (2) the principalis-chrysomelas group (White-winged Pheasants in Central Asia, including the subspecies principalis, zarudnyi, chrysomelas, bianchii, zerafschanicus and shawii); (3) the tarimensis group (Tarim Pheasant, tarimensis in Tarim Basin in southeastern Xinjiang, China); (4) the mongolicus group (Kirghiz Pheasants in northern Xinjiang, China and eastern Kazakhstan, comprising mongolicus and turcestanicus) and (5) the most subspecies-rich group, the torquatus group (Grey-rumped Pheasants, mostly found in China, containing 17 subspecies: decollatus, satscheuensis, pallasi, suehschanensis, torquatus, kiangsuensis, rothschildi, karpowi, strauchi, elegans, vlangalii, hagenbecki, edzinensis, alaschanicus, sohokhotensis, takatsukasae and formosanus) (Madge and McGowan 2002).

The Common Pheasant has a long history of captivity and being introduced as a common game species in western Europe, North America and Australia (Hill and Robertson 1988; Johnsgard 1999). This species deserves conservation management and sustainable use for several reasons. First of all, because natural populations of the Common Pheasant have been dramatically declining due to the loss of its natural habitats, hunting and other anthropogenic disturbances (Sotherton 1998), restocking of this bird is increasingly needed. For example, native habitat loss due to reclamation for agriculture caused a population decline in the subspecies principalis and persiscus in Iran (Solokha 1994). The subspecies turcestanicus is probably extinct now as a result of the aridification of the Aral Sea (Lepage 2007). As well, hybridized decendents between local subspecies and ex situ subspecies, due to artificial introduction of captive birds for hunting purposes, are evident in the wild (Braasch et al. 2011). Even worse is likely to occur, in so far as Common Pheasants in the wild may interbreed with a commercial, captive breed, the so-called “seven-color wild pheasant” which is a hybridized race between the Common Pheasant and its sister species, endemic to the Japan archipelagos, the Green Pheasant (Phasianus versicolor). For all these reasons genetic pollution in the wild Common Pheasant gene pool may occur. Last but not least, some range-restricted subspecies of the Common Pheasants inhabit isolated range and extreme environments such as arid regions, islands and mountains which preserve unique phenotypes and genotypes for future conservation and stocking (Braasch et al. 2011; Kayvanfar et al. 2017). For example, the formosanus subspecies of the Common Pheasant is endemic to Taiwan; other subspecies hagenbecki, alaschanicus and tarimensis are isolated and have adapted to semi-desert conditions (Johnsgard 1999). These conservation and management issues require evaluation using conservation approaches in genetics. Developing permanent genetic resources, such as autosomal microsatellites are of critical importance.

Microsatellites, also known as simple sequence repeats (SSRs), are a preferred type of markers in conservation genetics (Sunnucks 2000). Because of their heritable mode, SSRs usually have a higher mutation rate than that of mitochondrial and nuclear intronic markers and represent a very useful tool to genotype individuals and thus allow the quantifications of intraspecific genetic diversity, population structure and gene flow (Selkoe and Toonen 2006). Applications of SSRs are also reliable because of their relatively great abundance in genomes, high level of genetic polymorphism, co-dominant inheritance mode, analytical simplicity and repeatability of results across laboratories. So far, no species-specific microsatellites are available for the Common Pheasant although a previous study showed that cross-amplification of a very limited number of SSRs from other closely related Phasianinae species should be applicable for Common Pheasants (Baratti et al. 2001).

Recent advances in next generation sequencing (NGS) technologies enable the generation of large number of sequences efficiently and cost-effectively (reviewed in Ekblom and Galindo 2011). In addition, the so-called “Restriction-site Associated DNA” (RAD) method was consequently developed as a reliable means for genome complexity reduction (Baird et al. 2008). The concept is based on acquiring the sequence adjacent to a set of particular restriction enzyme recognition sites and then obtain sequences (RAD-seq) by NGS technology. Application of the RAD method, using the Illumina platform, has the advantage that it generates relatively long paired-end sequencing reads (100–150 bp), is cost-effective and sufficient to develop SSRs (Castoe et al. 2012). Another advantage is that a RAD-seq does not require a reference genome to be available and allows de novo assembly (Willing et al. 2011).

In this study we developed a set of autosomal microsatellites for the Common Pheasant using RAD-seq. We further designed multiplex PCR sets and tested for genetic polymorphism in a selected panel of 10 selected SSRs, which provide a tool for conservation genetic and studies of the evolution in the Common Pheasant.

To conduct RAD sequencing, we collected fresh tissues and blood samples from 60 individuals comprising seven subspecies in China (strauchi, vlangalii, kiangsuensis, karpowi, torquatus, elegans and tarimensis). Sample size for each subspecies varied between four to ten individual birds. Total genomic DNA was extracted using a QIAquick DNeasy kit (Qiagen, Hilden, Germany) following the manufacturer’s instructions.

Genomic DNA was digested with a restriction enzyme, ApeKI. Adapters P1 and P2 were ligated to the fragments. The P1 adapter contains a forward amplification primer site, an Illumina sequencing primer site and a barcode. Selected fragments were then subjected to end-repair and 3ʹ adenylated. The fragments are PCR amplified with P1- and P2-specific primers. Library was validating on the Agilent Technologies 2100 Bio-analyzer and the ABI StepOnePlus Real-Time PCR System. After adapter ligation and DNA cluster preparation, the samples were subjected to sequencing by a Hiseq 2000 sequencer (BGI, Shenzhen, China).

The raw data of 60 individuals have been processed by deleting adapter sequences and subsequently removing the reads, of which the rate of low quality (quality value ≤5 E) is more than or equal to 50%. All reads were then assigned to the individuals by the ambiguous barcodes and the specific recognition site (GWCC). Reads without a unique barcode and specific sequence were discarded. Final read length was trimmed to 82 nucleotides (minimum length). Then, the first four samples with high sequencing quality were selected to assemble the reference scaffolds. The assembly was performed using SOAPdenovo (Li et al. 2010), with scaffolds larger than 150 bp retained.

The collected polymorphic information of raw data from the 60 individuals was used to identify microsatellites by screening the sequence data for di-, tri-, tetra- and penta-nucleotide motifs with a minimum of ten repeats each. We applied MSATCOMMANDER 1.0.8 (Faircloth 2008) interfaces with the PRIMER3 software (http://bioinfo.ut.ee/primer3/), to allow the design of primers while minimizing potential structural or functional defects. The MSATCOMMANDER program was modified to ensure that the flanking region between the microsatellite and primer sequence would generate an amplicon size in the range of 100–250 bp, inclusive of the lengths of both primers (Brandt et al. 2014). After these procedures we randomly selected a panel of 30 novel di-nucleotide markers, together with five microsatellite markers isolated from the chicken genome (Baratti et al. 2001). We designed SSR primers using the online program Primer 3 (http://sourceforge.net/projects/primer3).

We continued to verify genetic polymorphism of the developed candidate microsatellites, by using another sample set, comprising 90 individuals from nine subspecies in China (vlangalii, satscheuensis, strauchi, elegans, decollatus, torquatus, kiangsuensis, karpowi and pallasi). We also included captive individuals of the Green Pheasant, a sister species of the Common Pheasant distributed in the Japanese archipelago. Sample size for each subspecies/species was restricted to ten individuals, allowing for unbiased comparisons of genetic polymorphism. Overall, we included 100 individuals for further analyses. Total genomic DNA was extracted using the same protocol with the samples for RAD-seq.

For each microsatellite locus, we calculated the number of alleles (N _A), observed (H _O) and the expected heterozygosities (H _E), as well as the polymorphism information content (PIC) with CERVUS v.3 (Kalinowski et al. 2007). Deviations from the Hardy–Weinberg equilibrium (HWE) and genotypic equilibrium between loci were tested with the same program. Significance levels were adjusted for multiple testing using the Bonferroni procedure (Rice 1989) if necessary.

In order to explore the detectability of population structures by the novel microsatellite set, we further identified the number of genetic clusters (K) among the 90 common and 10 Green Pheasants, using the Bayesian admixture model with the correlated allele frequencies option implemented in STRUCTURE v.2.3.4 (Pritchard et al. 2000; Falush et al. 2003). We performed one million Markov chain Monte Carlo (MCMC) repetitions and a burn-in of 200,000 repetitions with ten independent runs each for K = 1–7. The most likely number of genetic clusters was determined on the basis of the ad hoc statistics described in Evanno et al. (2005) using the STRUCTURE Harvester v.0.6.8 (Earl 2011).

About 127.98 G bases of raw data were generated for all the pooling lanes. After the raw data had been processed, about 123.27 G bases of clean data were retained. Then all reads with ambiguous barcodes were trimmed and about 114.91 G bases of clean data were kept for downstream analysis. The assembly generated 744,632 scaffolds larger than 150 bp, with a total length of 157 Mb. The average length was 211 bp with the largest of 1225 bp. The N50 statistic of the assembly was 254 bp. The microsatellite detection generated 7598 markers with 6419 di-nucleotide repeats, 766 tri-nucleotide repeats, 352 tetra-nucleotide repeats and 61 penta-nucleotide repeats (Additional file 1: Table S1). These data were further used for analysis of the population genome.

We detected significant departures from HWE at all loci, which is expected due to population structuring among individuals (Table 1). In addition, there was no evidence of genotypic disequilibrium after the Bonferroni correction. The number of alleles ranged from 2 to 8. The H _O values ranged from 0.03 to 0.98 and those of H _E from 0.049 to 0.905. The polymorphic information content (PIC) ranged between 0.048–0.893, with seven out of eleven loci having a PIC value around or above 0.50. In addition, we successfully amplified all 13 microsatellites with our 10 Green Pheasant samples. However, only six (PC6, PC8, PC9, MCW 97, MCW127 and MCW151) out of 13 loci showed polymorphism in the Green Pheasants.

The Bayesian clustering approach implemented in STRUCTURE suggested two genetic clusters (Fig. 1) as the most likely scenario based on the Evanno’s method. The means of the posterior probability, Ln P (D) (±SD) for different number of genetic groups (K), increased between K = 1–7 (Fig. 1b). The ∆K statistic reached a peak when K = 2 (Fig. 1c), suggesting that the result of two genetic clusters was the most likely scenario. This corresponds to the subdivision between the common and green pheasants. However, we also obtained a smaller peak when K = 5 (Fig. 1c) which most likely reflects a further population subdivision within Common Pheasants. We plotted the assignment of genetic clusters when K = 2–7 and found that the subspecies elegans, vlangalii and satscheuensis represented distinctive genetic clusters in contrast with the remaining subspecies (Fig. 1a).

Microsatellites are commonly used genetic markers in evolutionary, behavioral and conservation genetics due to their relatively high level of polymorphism and repeatability in genotyping. In a given species, microsatellite markers can be isolated and developed from scratch using methods such as magnetic beads (e.g. Wang et al. 2009) or by applying cross-species amplification using existing markers from closely related species (e.g. Dawson et al. 2010; Gu et al. 2012). The recent fast development of NGS has drastically promoted discoveries of novel microsatellites by using more cost-effective and less time-consuming strategies. Unlike laboratory-biased traditional methods, these approaches provide a massive amount of DNA sequencing reads as resources used for microsatellite detection by bioinformatic pipelines. Our results identified ten novel microsatellites in the Common Pheasant using an Illumina paired-end RAD sequencing strategy. Apart from RAD-seq, other sequencing strategies such as RNA-seq (e.g. Wang et al. 2012) and whole genome re-sequencing (e.g. Yang et al. 2015) have been consistently applied depending on the purpose of the research and on a budget.

We designed a panel of 13 multiplexed microsatellite markers that provided considerable polymorphic information to study population genetic structuring in common and green pheasants. In the Common Pheasant, the number of alleles ranged from 2 to 8, which is considered to be medium to high polymorphism (Botstein et al. 1980). However we found that seven loci in the green pheasants were monomorphic, probably owing to low level genetic diversity. The Green Pheasant has a smaller population and restricted range in its geographical distribution, compared to the wide ranging Common Pheasant. Another possibility is that, since only ten captive individuals of the Green Pheasant were included in the analysis, more tests with an enlarged sample set are likely required.

Using these novel microsatellites, we found a major genetic differentiation between the common and green pheasants. These two sister species of the genus Phasianus diverged about 2.8 million years ago based on a mitogenomic study by Li et al. (2015). We found further genetic population subdivisions within nine subspecies of the Common Pheasant in China. While no substantial genetic differentiation within the six subspecies in eastern China (decollatus, torquatus, kiangsuensis, karpowi, pallasi and strauchi), the three subspecies vlangalii, satscheuensis and elegans formed distinctive genetic clusters. These results corroborate the phylogeographic relationships revealed by mtDNA and nuclear intron data in an earlier study (Kayvanfar et al. 2017). In eastern China, subspecies are parapatric and their boundaries genetically vigorous, which is probably due to low genetic divergence and frequent gene flow. When it comes to western China, most subspecies are allopatric and their boundaries are distinct, probably indicating long-term isolation. However these first results are based only on subspecies within the torquatus group (Grey-rumped Pheasants). Samples of subspecies from the western Palearctic are needed to obtain a complete picture of population structuring within the Common Pheasant.

We can conclude that this novel set of microsatellites has proved to be useful for population genetics and conservation implications in the Common Pheasant as well as its sister species the Green Pheasant. As well, this study provides a checklist of 7598 candidate microsatellites (Additional file 1: Table S1) that can be used for genetic linkage mapping in the Common Pheasant and are probably useful to design markers for other closely related and endangered pheasant species (e.g. Chrysolophus, Lophura, Crossoptilon, Syrmaticus) (Gu et al. 2012).