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

Abstract

Background

The Common Pheasant (Phasianus colchicus) Linnaeus, 1758 is the most widespread pheasant in the world and widely introduced as a game bird. Increasing needs for conservation genetics and management of both wild and captive populations require permanent genetic resources, such as polymorphic microsatellites in order to genotype individuals and populations.

Methods

In this study, 7598 novel polymorphic microsatellites for the Common Pheasant were isolated using a RAD-seq approach at an Illumina high-throughput sequencing platform. A panel of ten novel microsatellites and three existing ones from the chicken genome were multiplexed and genotyped on a set of 90 individuals of Common Pheasants (representing nine subspecies and ten individuals each) and 10 individuals of the Green Pheasant (P. versicolor).

Results

These 13 microsatellites exhibited moderate to high levels of polymorphism, with the number of alleles per locus ranging from 2 to 8 and expected heterozygosities from 0.049 to 0.905. The first analysis of the genetic structure of subspecies/populations using a Bayesian clustering approach, implemented in STRUCTURE, showed two genetic clusters, corresponding to both the Green and the Common Pheasant, with further evidence of subpopulation structuring within the Common Pheasants.

Conclusion

These markers are useful genetic tools for sustainable uses and evolutionary studies in these two Phasianus pheasants and probably other closely related game birds.

Background

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.

Methods

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.

We checked their polymorphism on 2.5% agarose gels. In the end, ten novel and three extant markers with polymorphic signals were retained (Table 1). The selected 13 loci were further arranged into three PCR multiplex sets (Table 2) and each forward primer was labeled with a fluorescent dye. Each amplification was carried out in a 10-µL reaction volume containing 5 µL of PCR mix (QIAGEN Multiplex Kit), 1 µL of a primer mix and 1 µL of template DNA. The PCR conditions were as follows: initial denaturation at 95 °C for 5 min, followed by 35 cycles of denaturation at 94 °C for 30 s, annealing at 58 °C for 45 s and at 72 °C for 90 s and a final extension at 72 °C for 10 min. Products were isolated and detected on an ABI Prism 3730XL Genetic Analyzer (Applied Biosystems, service provided by Invitrogen, Shanghai, China). Fragment lengths were checked in comparison to an internal standard size (GeneScanTM-600LIZ, Applied Biosystems), using GeneMarker software v.2.2 (Soft Genetic).

Table 1 Characteristics of 13 microsatellites in a sample set of 90 individuals of the Common Pheasant (Phasianus colchicus)

Table 2 Characteristics of three multiplex sets used in genotyping of the Common Pheasant (Phasianus colchicus)

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).

Results

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).

Fig. 1

a Population structuring of Common and Green Pheasants inferred from STRUCTURE that used admixture model with correlated allele frequencies when K = 1–7. Each line represents an individual and the length of each color represents the possibility of the individual assigned into a given genetic cluster. The abbreviation stands for subspecies of Common Pheasant (VLA: vlangalii, SAT: satscheuensis, STR: strauchi, ELE: elegans, DEC: decollatus, TOR: torquatus, KIA: kiangsuensis, KAR: karpowi, PAL: pallasi, PV: Green Pheasant). b Although the means of the posterior probability Ln P (D) (±SD) increased per K, c two genetic clusters were suggested by the maximum value of ΔK and the order rate of change in posterior likelihood Ln P (X/K) per K

Discussion

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).